[gdal-dev] NetCDF rotated grid support

[Even Rouault]

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>
To:	"crs.swg at lists.opengeospatial.org" <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>
To:	"crs.swg at lists.opengeospatial.org" <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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20200415/1e55b456/attachment-0001.html>