Bug #895

grdvolume (some) wrong results

Added by Joaquim over 5 years ago. Updated over 5 years ago.

Status:ClosedStart date:2016-04-29
Priority:NormalDue date:
Assignee:Paul% Done:


Target version:Candidate for next minor release
Affected version:5.x-svn Platform:


Compute the lower half sphere (r = 10) with center at surface and top half as well

gmt grdmath -R-10.1/10.1/-10.1/10.1 -I0.1 X Y HYPOT 10 DIV ACOS SIN 10 MUL NEG = bot_half.grd
gmt grdmath bot_half.grd NEG = top_half.grd

compute top half's surface and volume

gmt grdvolume top_half.grd
0 408.04 2094.31958366 5.13263303515

The volume is rather correct (analytical val = 2094.3951...) but surface is way off.
Now compare with lower half 'outer' volume

gmt grdvolume bot_half.grd -Cr-10/0
0 310 2090.47541863 0

Now surface is close but volume is less precise (but it should be exactly the same)

Try to set a unit (meters)

gmt grdvolume top_half.grd -Se
0 5.0190580961e+012 2.58161011191e+013 5.14361472309

here the problem comes from the fact that it 'thinks' the grid is in geogs


#1 Updated by Paul over 5 years ago

  • Status changed from New to In Progress
  • Assignee set to Paul
  • % Done changed from 0 to 50

Partial fix for now: -S will now check if the grid is geographic and if not will refuse and tell users that if the grid is geographic that they need to supply -fg.

#2 Updated by Paul over 5 years ago

  • % Done changed from 50 to 80

Perhaps we can consider an under-the-hood for the reason for these problems: The level you are requesting (0) is a large flat area. If I select the contour 1e-6 for the top and -1e-6 for the lower then I get

0.00000100    316.79999133    2094.16297681    6.61036311
0.00000000    316.79999133    2094.16297681    0.00000000

which are acceptable, no? However, there are times when we want to use exactly the zero contour so perhaps educating the user via comments in the documentation is the way to go.
Also note that the -Cr option seem to not set the 1st and 4th output column.

#3 Updated by Paul over 5 years ago

  • Status changed from In Progress to Resolved
  • Target version set to Candidate for next minor release
  • % Done changed from 80 to 100

I think the safe solution (which I now have implemented) is to add small to each contour. Thus, if you want to zero contour you are getting the zero+small contour which is 1e-6 the data range (or possibly the contour increment). This avoids the pitfall of the -C0 above where we end up adding a large part of the area outside the sphere as part of the area inside the contour. With this change your commands give the answers I posted above (using -C0 and -Cr-10/0). I also fixed the missing height calculation for -Cr, and added a new test for grdvolume (sphere_volume.sh) based on the above In r16375.

#4 Updated by Paul over 5 years ago

  • Status changed from Resolved to Closed

Closed as fixed - reopen if there are other problems related to this.

Also available in: Atom PDF