grdvolume (some) wrong results
|Target version:||Candidate for next minor release|
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
#2 Updated by Paul over 2 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 2 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.