Making sense of grdvolume output

Q&A
1

I’ve used grdvolume and its output is unexpected. I probably don’t understand the documentation, but the documentation is difficult to follow. (For one, the output format is undocumented.)

I have a mask file in netCDF4 format. Dimensions are lon, lat, and z. z is either 0 or 1.

This command:

grdvolume file.grd -Sn -R$area -C0

produces:

0	14310	14085	1

I assume the first entry is the contour level, the second entry is the volume, the third entry is the area, and the last entry is volume/area. What I don’t understand is how volume is calculated.

From the documentation:

grdvolume reads a 2-D grid file and calculates the volume contained between the surface and the plane specified by the given contour (or zero if not given) and reports the area, volume, and maximum mean height (volume/area).

What is the “surface”? A naive definition of “surface” is a contour level of 0, in which case I expect the volume of the command with -C0 to be 0. (What’s also confusing about this sentence is that the order of area, volume, and volume/area don’t match the order of the output.)

More mystery: the documentation indicates that if I leave out the -C flag, a contour value of 0 is assumed. But this command:

grdvolume file.grd -Sn -R$area

returns:

0	176344	14088	0

The volumes are dramatically different.

The mask file has z values of either 0 or 1. I try this:

grdvolume file.grd -Sn -R$area -C1

and get this:

1	0	0	NaN

How is the area zero? 1s exist in the data.

Furthermore, this:

grdvolume file.grd -Sn -R$area -C-10

produces this:

-10	176368	1777766	10

despite the data lacking any value of -10. And this:

grdvolume file.grd -Sn -R$area -C-2

produces this:

-2	176368	366825	2

which has a different area from the “-10” contour.

2

z(x,y) is the surface. You are looking at the volume between the plane defined by a contour level and z(x,y). So, if you choose a contour of 1 then there is nothing above your plane since you said all z values are either 0 or 1. For a contour of -10 you are adding the huge slab between the plane z = -10 and z = 0.
If you are able to share your subset (gmt grdcut file.grd -R$area -Gsubset.grd) then I could have a look. I am not sure why -C0 and no -C would give different values. Either a bug or a documentation failure.

3

Thanks for the response. Part of my confusion is that the output is structured counter to my assumption—it’s (contour level, area, volume, volume/area) and not (contour level, volume, area, volume/area). I’d like to recommend that inputs and outputs be documented in a clear way, as they are for Python libraries (e.g., numpy docs separate parameters and returns).

So, if you choose a contour of 1 then there is nothing above your plane since you said all z values are either 0 or 1.

This makes sense. But why is there no area? The surface and the plane intersect.

For a contour of -10 you are adding the huge slab between the plane z = -10 and z = 0.

To be clear, shouldn’t the statement be “between the plane at contour = -10 and z = 0 where z = 0 and contour = -10 and z = 1 where z = 1”?

Here are the data:

https://drive.google.com/open?id=1O1sR6dIhzQx-hOBWzFhiQ_IbI1ausSey

4

Thanks for the comments. I have updated grdvolume documentation in master to be clear about the output format.
I will try to find some time this weekend to look at your grid before I comment on some of the results.

5

Your case -C1: yes true the plane and surface intersect, but the algorithm is looking for volume above the plane and it is zero; it is computing the area in that vein as well, finding nothing.
A couple of other points: You may be (or not) assuming that your 0/1 grid surface has vertical walls but of course we cannot represent both 0 and 1 at the same node. So each of your one-nodes surrounded by 0s will look more like pyramids:

map
map2048×1929 67.3 KB

The center node is one, and grdcontour shows the 0.1 contours. If this was the only one in a grid then we would compute slices of this pyramid. In verifying the results I found a bug affecting contributions to the volume for some cases near a contour - this has been fixed master. Also note that if you simply want your one-nodes to represent square plateaus of height one and evaluate the volume and area of those then it is probably simpler to use grdmath and add up the contributions from those nodes:

gmt grdmath your.grd AREA MUL SUM = result.grd

and run gmt grdinfo result.grd to see the z_min = zmax value for the volume (here == area since z = 1). You will need to convert this result (in km^2) to your nautical miles squared.

6

OK, thanks. One more recommendation. The documentation uses the word “between” when referring to the calculation of the volume above the contour value and below the surface. I’d replace “between” with an explicit “above” or “below”—“between” is ambiguous because it could refer to either. The discussion of inner and outer volumes is buried in the last 2 sentences for the -C flag.

7

OK, now using above and below in those places.