[gdal-dev] NetCDF config option for west > east?

[Markus Metz]

On Sun, Nov 25, 2018 at 2:24 AM Even Rouault <even.rouault at spatialys.com>
wrote:
>
> On samedi 24 novembre 2018 23:09:34 CET Markus Metz wrote:
> > I have created a PR for this, but there are still issues:
> >
> > is
> > Pixel Size = (0.099999991285714,-0.100000004000000)
> > Corner Coordinates:
> > Upper Left  ( -25.0000122,  69.9999970)
> > Lower Left  ( -25.0000122,  29.9999954)
> > Upper Right (  44.9999817,  69.9999970)
> > Lower Right (  44.9999817,  29.9999954)
> > Center      (   9.9999847,  49.9999962)
> >
> > should be
> > Pixel Size = (0.100000000000000,-0.100000000000000)
> > Corner Coordinates:
> > Upper Left  ( -25.0000000,  70.0000000)
> > Lower Left  ( -25.0000000,  30.0000000)
> > Upper Right (  45.0000000,  70.0000000)
> > Lower Right (  45.0000000,  30.0000000)
> > Center      (  10.0000000,  50.0000000)
>
> Yes, the issue is that the data type for the x variable in the .nc file is
> float32, so it has not enough precision. You could try to round the
coodinates
> with some heuristics, like (untested !)
>
> for the resolution:
>
> double K = std::pow(10, 4 - std::round(std::log(res) / std::log(10)));
> if( std::fabs(K * res - std::round(K * res)) < 1 )
> {
>         res = std::round(K * res) / K;
> }
>
> for the coordinates:
>
> // Case where coordinates are aligned on a multiple of the resolution
> if ( std::fabs(x / res - std::round(x / res)) < 1e-3 )
> {
>     x = std::round(x / res) * res
> }
> // Case were they are aligned with a half-pixel shift
> else if (std::fabs((x - res/2)/res) - std::round((x - res/2)/res)) < 1e-3
)
> {
>     x = std::round((x - res/2)/res) * res + res/2;
> }

Many datasets with degrees as coordinate units have resolution and extents
that can be expressed as integer or 0.5 seconds, e.g. SRTM at 1, 3, and 30
arc seconds. Some of these deviations caused by fp precision limits go away
if you convert to arc seconds and then round to the nearest e.g. 1e-6 arc
second. GRASS GIS is using this approach successfully.

About the PR: if the tested nc files are formally legal with their
longitude jump at 359.95, 0.04998779, then this could be regarded as a bug
fix, right?

Markus M
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20181125/d15d6af3/attachment.html>