[gdal-dev] NetCDF rotated grid support

[Even Rouault]

Louis-Philippe,

> I did some more tests with the rotated_pole.nc in gdal GitHub repo
> (https://github.com/OSGeo/gdal/tree/master/autotest/gdrivers/data) just to
> make sure the issue was not with my data.
> 
> With the new versions, the EXTENSION is missing from the WKT, which I guess
> makes the reprojection fail when I try to view it in Qgis ?

Kind of... I confirm the issue. A complex interaction between GDAL and PROJ due to the 
netCDF driver using a non-standard way of expressing a rotated pole CRS.

Well, there isn't really a standardized way, but what I did for GRIB2 recently was more in the 
right direction (hopefully!)

So to address the issue, I've done 2 things:
- a fix for GDAL (in time for 3.1.0): https://github.com/OSGeo/gdal/pull/2438
- and one for PROJ: https://github.com/OSGeo/PROJ/pull/2185

The GDAL fix should probably be sufficient without the PROJ one, but I won't bet too much 
on that...

Ultimately the netCDF driver should be improved similarly to the GRIB one, that is creating a 
DERIVINGCONVERSION["Pole rotation (netCDF convention)"] instead of a hack.

So with both I now get with gdalinfo:

[...]
Coordinate System is:
GEOGCRS["unnamed",
    BASEGEOGCRS["unknown",
        DATUM["unknown",
            ELLIPSOID["unknown",6367470,0,
                LENGTHUNIT["metre",1,
                    ID["EPSG",9001]]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8901]]],
    DERIVINGCONVERSION["unknown",
        METHOD["PROJ ob_tran o_proj=longlat"],
        PARAMETER["lon_0",18,
            ANGLEUNIT["degree",0.0174532925199433,
                ID["EPSG",9122]]],
        PARAMETER["o_lon_p",0,
            ANGLEUNIT["degree",0.0174532925199433,
                ID["EPSG",9122]]],
        PARAMETER["o_lat_p",39.25,
            ANGLEUNIT["degree",0.0174532925199433,
                ID["EPSG",9122]]]],
    CS[ellipsoidal,2],
        AXIS["longitude",east,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433,
                ID["EPSG",9122]]],
        AXIS["latitude",north,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433,
                ID["EPSG",9122]]]]
[...]
Upper Left  ( -35.4700000,  23.6500003) ( 55d 6'46.50"W, 56d15'18.83"N)
Lower Left  ( -35.4700000, -23.8699996) ( 16d 4'19.39"W, 18d42'19.70"N)
Upper Right (  24.8100002,  23.6500003) ( 78d43'49.95"E, 63d51'22.90"N)
Lower Right (  24.8100002, -23.8699996) ( 42d35'19.21"E, 22d45'11.78"N)
Center      (  -5.3299999,  -0.1099997) (  9d37'52.99"E, 50d20'18.47"N)

The corner transformed in (non-rotated) geographic coordinates are now the same as GDAL 
2.4, whereas before the fix, they were centered close to (0,0) as you mentionned.

Even

-- 
Spatialys - Geospatial professional services
http://www.spatialys.com
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20200423/5ee11b02/attachment.html>