Just an informational question.
The grdinfo -L2 command uses a area weighting scheme for global geographic grids, according to a note in the grdinfo man page:
Note: If the grid is geographic then each node represents a physical area that decreases with increasing latitude. We therefore report spherically weighted statistical estimates for such grids.
Peering into the gmt_stat.c source code, it can be seen that this is accomplished in the gmtstat_get_geo_cellarea function. What I’m trying to get a handle on, is the value used when the parameter R2 is evaluated: R2 = pow (0.001 * GMT->current.proj.mean_radius, 2.0).
Specifically, how can I ask GMT to tell me what the value of GMT->current.proj.mean_radius is?
Ask how?
That value is stored in the API structure and is set from the …default settings or otherwise from a gmtset ... setting.
The GMT defaults parameter decides what kind of radius:
PROJ_MEAN_RADIUS
with default authalic. Read the gmt.conf manual on other options. It then computes that from the selected ellipsoid (PROJ_ELLIPSOID). As for the actual computation you would have to peak inside or google what these various man radii mean.
Appreciate the replies, but I think you missed what I was asking. Or, that GMT is unable to output an internally used value.
Is there a command that can be issued to have GMT report the value of the parameter, GMT->current.proj.mean_radius, given the selected PROJ_ELLIPSOID?
Somehow, GMT, internally, is using a value - what is that value? I understand it is dependent on the selection of PROJ_ELLIPSOID and that it is authalic. I’m using WGS-84.
I don’t need to see how it is computed, just want the value that GMT is using, internally, for the grdinfo -L2 commands. I need this to check that another software package is computing these things properly.
A little more digging led me to this equation in gmt_map.c (line 2286):
r = sqrt (0.5 * a * a + 0.5 * b * b * atanh (GMT->current.proj.ECC) / GMT->current.proj.ECC)
Given the WGS-84 ellipsoid parameters, I can probably compute this value, outside of GMT.
Thanks Paul and Joaquim! John
1 Like
I didn’t and that’s why I asked Ask how?. You’d need an FFI to do that, which is particular troublematic for the API (and its doughter, GMT) structures that change for each GMT version.
An FFI - hmmm. Fatal Familial Insomnia? 
I’m surmising that GMT is unable to offer values of mnay of its internal variables. No worries.
Using the above equation, from gmt_map.c, I get 6367477.33489 meters, for WGS-84. Presuming that current.proj.ECC is equivalent to the square of the first eccentricity (= 2f - f^2, where f = 1/298.257223563).
I am unable to locate that particular equation in the literature (yet).
Foreign Function Interface … or fprintf(stderr, ...).
No shared lib gives up its secret numbers if we don’t pull them out of it.
I could say Julia but the problem in this case is what I said. the API and GMT struct contents change too often to make it worth wrapping them.
Just need to make a correction to my above reply.
The parameter, current.proj.ECC, is not the square of the first eccentricity (e^2 is commonly used in various geometric geodesy equations); it is simply the eccentricity, e. Thus the corrected WGS-84 radius, using GMT’s internal formulation is: 6371007.2342 meters.
Using a formulation found in Jordan’s Handbook of Geodesy, I get a value of: 6371007.0743 meters.
But I’m not going to worry about 16cm difference when it comes to relative cell weighting! Radius probably doesn’t matter that much, when weighting cells, anyway.