[gdal-dev] NetCDF rotated grid support

Thu Apr 23 09:14:47 PDT 2020

Hi Even,

I can confirm I also find that working with those projections is a nightmare ;)

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.

I found that I also had the same issue with rotated_pole.nc when trying to display the file in QGIS using newer version of gdal (3.0.4) and Proj (6.3.1). The NetCDF which should be on Europe, is misplaced and centered on Africa (using 0,0 as origin I guess).

But, I tried to open the same file (rotated_pole.nc) in QGIS using gdal (2.4.1) and proj (5.2.0) and the file was correctly placed, thus using those versions the rotated pole grid is correctly supported/read. The issue I'm having with rotated grid seems to be only in recent version of gdal/proj.

And if I look at the projection reported int the gdalinfo I get:

  *   gdal 3.0.4 and proj 6.3.0:
     *   Driver: netCDF/Network Common Data Format Files: rotated_pole.nc Size is 137, 108 Coordinate System is: PROJCRS["unnamed", BASEGEOGCRS["unknown", DATUM["unnamed", ELLIPSOID["Spheroid",6367470,594.313048347956, LENGTHUNIT["metre",1, ID["EPSG",9001]]]], PRIMEM["Greenwich",0, ANGLEUNIT["degree",0.0174532925199433, ID["EPSG",9122]]]], CONVERSION["unnamed", METHOD["Rotated_pole"]], CS[Cartesian,2], AXIS["easting",east, ORDER[1], LENGTHUNIT["metre",1, ID["EPSG",9001]]], AXIS["northing",north, ORDER[2], LENGTHUNIT["metre",1, ID["EPSG",9001]]]] Data axis to CRS axis mapping: 1,2 PROJ.4 string is: '+proj=ob_tran +o_proj=longlat +lon_0=18 +o_lon_p=0 +o_lat_p=39.25 +a=6367470 +b=6367470 +to_meter=0.0174532925199 +wktext' Origin = (-35.470000000560987,23.650000304819269) Pixel Size = (0.440000001121970,-0.439999999286972) Metadata: NC_GLOBAL#Conventions=CF-1.4 projection#false_easting=0 projection#false_northing=0 projection#grid_north_pole_latitude=39.25 projection#grid_north_pole_longitude=-162 projection#north_pole_grid_longitude=0 projection#semi_major_axis=6367470 projection#semi_minor_axis=6356756 tg#grid_mapping=projection ... Corner Coordinates: Upper Left ( -35.4700000, 23.6500003) Lower Left ( -35.4700000, -23.8699996) Upper Right ( 24.8100002, 23.6500003) Lower Right ( 24.8100002, -23.8699996) Center ( -5.3299999, -0.1099997) Band 1 Block=137x1 Type=Int16, ColorInterp=Undefined Computed Min/Max=-3778.000,1476.000
  *   gdal 2.4.1 and proj 5.2.0:
     *   Driver: netCDF/Network Common Data Format Files: rotated_pole.nc rotated_pole.nc.aux.xml Size is 137, 108 Coordinate System is: PROJCS["unnamed", GEOGCS["unknown", DATUM["unknown", SPHEROID["Spheroid",6367470,594.3130483479559]], PRIMEM["Greenwich",0], UNIT["degree",0.0174532925199433]], PROJECTION["Rotated_pole"], EXTENSION["PROJ4","+proj=ob_tran +o_proj=longlat +lon_0=18 +o_lon_p=0 +o_lat_p=39.25 +a=6367470 +b=6367470 +to_meter=0.0174532925199 +wktext"]] PROJ.4 string is: '+proj=ob_tran +o_proj=longlat +lon_0=18 +o_lon_p=0 +o_lat_p=39.25 +a=6367470 +b=6367470 +to_meter=0.0174532925199 +wktext' Origin = (-35.470000000560987,23.650000304819269) Pixel Size = (0.440000001121970,-0.439999999286972) Metadata: NC_GLOBAL#Conventions=CF-1.4 projection#false_easting=0 projection#false_northing=0 projection#grid_north_pole_latitude=39.25 projection#grid_north_pole_longitude=-162 projection#north_pole_grid_longitude=0 projection#semi_major_axis=6367470 projection#semi_minor_axis=6356756 tg#grid_mapping=projection ... Corner Coordinates: 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) Band 1 Block=137x1 Type=Int16, ColorInterp=Undefined Min=-3778.000 Max=1476.000 Computed Min/Max=-3778.000,1476.000

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 ?

Thanks

LP

On 2020-04-15 10:13, Even Rouault wrote:

Louis-Philippe,

Support for rotated grid is a nightmare due to different conventions used for GRIB, netCDF or PROJ ob_tran method. The GDAL GRIB and netCDF drivers normally do what is needed map from the native product convention to the PROJ one.

Below 2 emails I wrote some time ago the OGC CRS SWG group which might bring some light (or deep obscurity!) to this topic:

1)

Rotated pole/grid

From: Even Rouault <even.rouault at spatialys.com><mailto:even.rouault at spatialys.com>

To: "crs.swg at lists.opengeospatial.org"<mailto:crs.swg at lists.opengeospatial.org> <crs.swg at lists.opengeospatial.org><mailto:crs.swg at lists.opengeospatial.org>

Date: 03/10/2019 14:26

Hi again,

(was initially a part of my previous email about "Expressing a ProjectedCRS of

a DerivedGeographicCRS", but this deserves its own thread)

Standardization of rotated pole is another interesting/annoying topic I've

discussed with a few people privately. This tends to be tricky as I always

stumble on different conventions regarding the semantics of parameters and the

maths behind.

For posterity, and my own recollection, here are my latest findings:

1) PROJ's ob_tran (https://proj.org/operations/projections/ob_tran.html)

formulas are a straightforward implementation from Snyder's "Map projections -

a working manual", in particular equations (5-7) and (5-8b) for the forward

path, at page 31 of https://pubs.usgs.gov/pp/1395/report.pdf

(although I've some disagreement with the 'new pole' terminology currently

used in the above HTML help page:

https://github.com/OSGeo/PROJ/pull/1653 )

2) But they are different from the conventions used for PDS3 or ISIS3

planetary image formats, a bit like if the forward and inverse paths had been

switched, so you need to 'switch' parameters, do crazy things like 180 -

pole_latitude, etc when crossing different conventions... Interestingly one

doc for PDS3 mentions at the middle of page 10 of https://pds-imaging.jpl.nasa.gov/documentation/Cassini_BIDRSIS.PDF that "The derivation

of these equations is given in reference [12], p. 29—32, but a

critical pair of signs appear to be exchanged between the formulae

for the forward and inverse transformations. Reference [13], p. 6,

gives the correct forward transformation, and reference [14], pp. 77—

78, give both the forward and inverse transformations"

"reference [12] p29-32" is Snyder's "Map projections - a working manual" I

mentionned before. So this somewhat acknowledges the existence of different

conventions (not sure if there are really "wrong" and "correct" ones)

reference [13] is "Bugayevskiy, L. M., and Snyder, J. P., 1995, Map

projections: A reference manual, Taylor & Francis, London and Bristol,

Pennsylvania, 328 pp."

and [14] "Yang, Q. H., Snyder, J. P., and Tobler, W. R., 2000, Map projection

transformations: Principles and applications, Taylor & Francis, London and

Bristol, Pennsylvania, 367 pp"

PDS/ISIS3 implementation of [13] and [14] is in

https://github.com/USGS-Astrogeology/ISIS3/blob/3.8.0/isis/src/base/objs/

ObliqueCylindrical/ObliqueCylindrical.cpp#L311

3) When dealing with netCDF rotated pole, I've also seen issues with

conventions.

An implementation of that convention is in netcdf-java:

https://github.com/Unidata/netcdf-java/blob/

3ce72c0cd167609ed8c69152bb4a004d1daa9273/cdm/core/src/main/java/ucar/unidata/

geoloc/projection/RotatedPole.java

4) GRIB2 also uses a yet another different convention (with south pole

coordinates as parameters)

An implementation of that convention is also in netcdf-java in

https://github.com/Unidata/netcdf-java/blob/

3ce72c0cd167609ed8c69152bb4a004d1daa9273/cdm/core/src/main/java/ucar/unidata/

geoloc/projection/RotatedLatLon.java

I was also pointed by a data producer of GRIB2 rotated pole datasets that they

used https://www.mathworks.com/matlabcentral/fileexchange/43435-rotated-grid-transform?focused=3795951&tab=function&ue&nocookie=true as the inspiration for

their maths, and luckily it gives consistent result with the Java

RotatedLatLon code.

All in all, a lot of confusion which I think depends if you consider north

pole vs south pole, or the coordinates of the unrotated pole in the rotated

grid, or the coordinates of the new rotated pole in the unrotated grid...

As far as I experimented, there's always a way of using a given implementation

(in my case ob_tran) with the other conventions by tweaking the values of the

3 angular parameters. So, despite the diversity of conventions, considering

standardizing one could still be useful.

Even

2)

Re: Rotated pole/grid

From: Even Rouault <even.rouault at spatialys.com><mailto:even.rouault at spatialys.com>

To: "crs.swg at lists.opengeospatial.org"<mailto:crs.swg at lists.opengeospatial.org> <crs.swg at lists.opengeospatial.org><mailto:crs.swg at lists.opengeospatial.org>

Date: 22/11/2019 02:35

Hi,

reviving this not-so-exciting topic :-)...

For the purpose of dealing with GRIB rotated long/lat grids,

https://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_doc/grib2_temp3-1.shtml,

I'm about to use a "Pole rotation (GRIB convention)" conversion method with

the following 3 parameters, whose naming and semantics comes mostly from the above GRIB spec:

- "Latitude of the southern pole (GRIB convention)"

Latitude of the point from the unrotated

CRS, expressed in the unrotated CRS, that will become the south pole of the

rotated CRS.

- "Longitude of the southern pole (GRIB convention)"

Longitude of the point from the unrotated

CRS, expressed in the unrotated CRS, that will become the south pole of the

rotated CRS.

- "Axis rotation (GRIB convention)"

The angle of rotation about the new polar

axis (measured clockwise when looking from the southern to the northern pole)

of the coordinate system, assuming the new axis to have been obtained by

first rotating the sphere through southPoleLongInUnrotatedCRS degrees about

the geographic polar axis and then rotating through

(90 + southPoleLatInUnrotatedCRS) degrees so that the southern pole moved

along the (previously rotated) Greenwich meridian.

The explanation of the axis rotation parameter is rather obscure to me. Would need

some scetch... Anyway, in practice, the few datasets I've accessed to

(like the GRIB2 ones at http://collaboration.cmc.ec.gc.ca/cmc/cmos/sample_grib/ or

the GRIB1 one at https://trac.osgeo.org/gdal/ticket/7104 ) use axis rotation = 0

I explicitly chose the " (GRIB convention)" suffix for parameter names as well,

as I'm pretty sure that the semantics of similarly named parameters

for other formulations of Pole rotation could be different.

So an example of a resulting DerivedGeographicCRS:

GEOGCRS["Coordinate System imported from GRIB file",

BASEGEOGCRS["Coordinate System imported from GRIB file",

DATUM["unnamed",

ELLIPSOID["Sphere",6371229,0,

LENGTHUNIT["metre",1,

ID["EPSG",9001]]]],

PRIMEM["Greenwich",0,

ANGLEUNIT["degree",0.0174532925199433,

ID["EPSG",9122]]]],

DERIVINGCONVERSION["Pole rotation (GRIB convention)",

METHOD["Pole rotation (GRIB convention)"],

PARAMETER["Latitude of the southern pole (GRIB convention)",-31.758312,

ANGLEUNIT["degree",0.0174532925199433,

ID["EPSG",9122]]],

PARAMETER["Longitude of the southern pole (GRIB convention)",-92.402969,

ANGLEUNIT["degree",0.0174532925199433,

ID["EPSG",9122]]],

PARAMETER["Axis rotation (GRIB convention)",0,

ANGLEUNIT["degree",0.0174532925199433,

ID["EPSG",9122]]]],

CS[ellipsoidal,2],

AXIS["latitude",north,

ORDER[1],

ANGLEUNIT["degree",0.0174532925199433,

ID["EPSG",9122]]],

AXIS["longitude",east,

ORDER[2],

ANGLEUNIT["degree",0.0174532925199433,

ID["EPSG",9122]]]]

And the maths are below:

{{{

#include <math.h>

#include <stdio.h>

#define PI 3.141592653589793

// Implementation of https://github.com/Unidata/netcdf-java/blob/3ce72c0cd167609ed8c69152bb4a004d1daa9273/cdm/core/src/main/java/ucar/unidata/geoloc/projection/RotatedLatLon.java

void rotate(double lon, double lat, double rot1, double rot2, double cosDlat, double s, double& rlon, double &rlat)

{

double e = (lon - rot1) / 180 * PI; // east

double n = lat / 180 * PI; // north

double cn = cos(n);

double x = cn * cos(e);

double y = cn * sin(e);

double z = sin(n);

double x2 = cosDlat * x + s * z;

double z2 = -s * x + cosDlat * z;

double R = sqrt(x2 * x2 + y * y);

double e2 = atan2(y, x2);

double n2 = atan2(z2, R);

rlon = e2 / PI * 180 - rot2;

rlat = n2 / PI * 180;

}

// Unrotated -> rotated

void forward_transform(double lon, double lat, // Input

double south_pole_lat,

double south_pole_lon,

double axis_rotation,

double& rlon, double& rlat // Output

)

{

double dlat_rad = (south_pole_lat - (-90)) / 180 * PI;

double sinDlat = sin(dlat_rad);

double cosDlat = cos(dlat_rad);

rotate(lon, lat, south_pole_lon, axis_rotation, cosDlat, sinDlat, rlon, rlat);

}

int main()

{

double south_pole_lat = -39.3;

double south_pole_lon = 18;

double axis_rotation = 45;

double lon = 12;

double lat = 55;

double rlon, rlat;

forward_transform(lon, lat,

south_pole_lat, south_pole_lon, axis_rotation,

rlon, rlat);

printf("%.8f %8f\n", rlon, rlat); // -48.44759099 4.439720

// Equivalent to +proj=ob_tran +o_proj=longlat +o_lon_p=-axis_rotation +o_lat_p=-south_pole_lat +lon_0=south_pole_lon

//

//$ echo "12 55" | cs2cs -d 8 +proj=longlat +to +proj=ob_tran +o_proj=longlat +o_lon_p=-45 +o_lat_p=39.3 +lon_0=18

// -48.44759099 4.43972046 0.00000000

}

}}}

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/25f45852/attachment-0001.html>