About the ERROR from psbasemap/grdimage [ERROR]: "South" (min elevation) is outside 0-90 degree range

[Hanbing59]

Hi,

I just got errors when running PlotMap for a skyplot.

psbasemap [ERROR]: "South" (min elevation) is outside 0-90 degree range
psbasemap [ERROR]: General map projection error
grdimage [ERROR]: "South" (min elevation) is outside 0-90 degree range
grdimage [ERROR]: General map projection error
psbasemap [ERROR]: "South" (min elevation) is outside 0-90 degree range
psbasemap [ERROR]: General map projection error
psxy [ERROR]: "South" (min elevation) is outside 0-90 degree range
psxy [ERROR]: General map projection error

I checked the GMT plot script generated by PlotMap, found a negative value for the min elevation set by the region variable:

region=-180/180/-1.15146e-13/90

The region variable serves as a setting for the -R option which requires the min elevation is in the range 0-90 for Polar coordinates projections (-Jp|P).

Then I checked the file PlotMap.cpp, found it relates to the calculation of variable minB, which can be set by the user or inferred from the input data file if not set (default) according to the documentation. In my case, I leave it empty, which should mean "not set and inferred from the input data" according to my understanding on the documentation.

As I am not very familiar with C++ at this point and can not inspect the source code thoroughly to find the precise code that makes this rounding errors, by results from AI coding agent:

For gridded data, points are treated as cell centers by default, and the code expands boundary by half-grid spacing when gridlineRegistered is false (this is my case) and the boundary is set by lat - bufferLat. So if the minimum latitude in file is 1.5 and inferred increment is 3, the inferred map boundary is expected to be 1.5 - 1.5 = 0. The -1.15146e-13 is just floating-point roundoff near zero (from coordinate/increment computations).

And the practical options:

1. Set gridlineRegistered=1 if the file coordinates are node values (not cell centers).
2. Set minPhi explicitly.
3. Code-side improvement: clamp tiny values near zero before writing region (e.g., abs(val) < 1e-12 -> 0).

Now back to my thoughts:

1. set gridlineRegistered=yes works. But as I understand, for many data, it can be regarded as either point values or cell means. For example, the antenna PCV data, regarded as point values or cell means both makes sense. This should not be the root reason.
2. Set minPhi explicitly also works but not the best practice from programming perspective.
3. clamp tiny values near zero before writing region. This may work but only should be in the right place. I need more time to investigate the code.

So far, that's all I have. I will update my progress here. If I am wrong somewhere, or whoever has some comments, you are more than welcome to leave several words here.

Regards,
Hanbing

[Hanbing59]

Update:

as a reply to the answer from @tmayerguerr (#396 (comment)), in the step GnssAntennaDefinition2Skyplot, the generated grid data has already been limited with the latitude range [0, 90].

[Hanbing59]

Hi @tmayerguerr,

I think I found the root reason for this error. (Unfortunately, not yet the solution)

Let me explain what I've found.

When we use the program GnssAntennaDefinition2Skyplot to generate a GriddedData file. Due to some kind of rounding errors, the coordinates in the GriddedData file are not precise perfectly, even though the setting values for deltaLambda and deltaPhi of the grid should result in perfectly spaced grids without rounding of dividing remainders. For example, in my case, I set deltaLambda=1 and deltaPhi=0.2 (I tried larger values, may have same problem) for a geograph grid with a border set as

-180 180 0 90 0

The obtained GriddedData file looks like (only the ending several lines shown here):

As you can see, the min elevation is not exactly 0.1(which it should be) but very close, i.e. 0.09999.... This should be not a problem, generally.

However, when this GriddedData servers as input to the program PlotMap to plot a skyplot. If user don't set the variable minPhi for min.latitude (which is allowed and probably the case for most of the time), it will infer minPhi from the input grid data. It does this inferring by searching for the min value of the latitude vector subtracted by half of the latitude space, which is 0.1 in my case. In another word, it does the following subtraction to get the value for minPhi:

9.999999999988502808e-02 - 0.1

which gives the small negative result like -1.15146e-13. And this would be written into the GMT plotting scripts as the region variable, causing the error:

"South" (min elevation) is outside 0-90 degree range

because of using polar coordinates projections.

Hope I've made this clear. I think there was another user having similar concerns about rounding errors in GROOPS, like #211 ?

I did not check the grid sampling code yet to see if it does use the similar algorithm for calculating the coordinates as used in the time series.

Regards,
Hanbing