[gdal-dev] NetCDF rotated grid support

[Rousseau Lambert2, Louis-Philippe (EC)]

Great thanks a lot !

I also think that something like GRIB2 would be nice for NetCDF. I will take note of this for future improvement to the NetCDF driver :)

LP

On 2020-04-23 14:46, Even Rouault wrote:

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/19dd3bd9/attachment.html>