OGR does not support PROJ.4 +proj=ob_tran , rotated_pole is not supported by the netcdf driver

[etourigny]

It seems OGR and geotiff (maybe proj.4 too?) need to be modified to support the ob_tran projection (General Oblique Transformation).

The problem was initially reported in bug #4251:

rotated pole grids ARE supported by PROJ.4 (through ob_tran transformation, in a somewhat counter-intuitive way, see my recent post at:  http://osgeo-org.1803224.n2.nabble.com/forward-inverse-projection-logic-for-ob-tran-td6874351.html). So the CRS of the data is, in PROJ.4 notation:

+proj=ob_tran +o_proj=latlon +o_lon_p=-162 +o_lat_p=39.25 +lon_0=180 

However, GDAL/OGR's limited projection support via gdal/ogr/ogr_srs_proj4.cpp doesn't include ob_tran transformation, which is indeed a pity. 

As rotated_pole is quite common in netcdf-cf files created by RCMs (mostly from European centers), it is important to support it within gdal. In fact, ob_tran may be the solution to support any other projections we might not support yet in the netcdf driver?

Perhaps Pat and 'peifer' would like to work on this? Adding Frank and Even in cc as they might have something to say.

It seems from this message that in fact ob_tran is not fully supported by proj.4? ​http://osgeo-org.1803224.n2.nabble.com/forward-inverse-projection-logic-for-ob-tran-td6874351.html

From this post it seems we need to use the inverse projection from rotated pole to lat/lon. ​http://osgeo-org.1803224.n2.nabble.com/Re-rotated-pole-ob-tran-help-needed-tt5149504.html#none

Perhaps this workaround suggested by Frank, using a WKT EXTENSION[] node works?: ​http://osgeo-org.1803224.n2.nabble.com/Does-GDAL-handle-rotated-spheres-td2030419.html

more info:

​http://www.osgeo.org/pipermail/gdal-dev/2005-June/005839.html ​http://osgeo-org.1803224.n2.nabble.com/quot-standards-quot-process-td2846302.html

Related to bugs #4251 and #4284 that deal with the netcdf driver specifically.

[peifer]

Replying to etourigny:

Perhaps Pat and 'peifer' would like to work on this?

I would if I could. But I'm not a programmer and already too old to become one :-(. I could perhaps help out with testing, or some such. Hermann (osgeo id: peifer)

[rhattersley]

Replying to etourigny:

From this post it seems we need to use the inverse projection from rotated pole to lat/lon. ​http://osgeo-org.1803224.n2.nabble.com/Re-rotated-pole-ob-tran-help-needed-tt5149504.html#none

It is possible to use the forward projection instead of the inverse projection. The key step is to note that for any rotated pole definition it is possible to derive a definition which encodes the inverse rotation.

Hence, it is possible to convert the CF-style rotated pole specification (i.e. the location of the rotated pole expressed in standard lat/lon) into the proj.4 style rotated pole specification (the location of the standard North pole in the rotated coords).

For example, when north_pole_grid_longitude is zero, the proj.4 parameters are given by:

• o_lon_p = 0
• o_lat_p = grid_north_pole_latitude
• lon_0 = 180 + grid_north_pole_longitude

And by rummaging in the proj.4 source code, I have determined a way to obtain the resulting rotated pole coords in degrees by using the to_meter parameter.

The resulting proj.4 definition strings look like:

+proj=ob_tran +o_proj=latlon +o_lon_p=... +o_lat_p=... +lon_0=... +a=1 +to_meter=0.0174532925199

NB. This definition ceases to work if you, or some intervening software layer, adds +units=....

[etourigny]

Thanks for you suggestion.

Do you know what to do if north_pole_grid_longitude is different from 0? In other workds, what would be the inverse transform (and forward transform following your logic)?

Your formula works like the example in the previously linked post:

$ proj -f "%f" +proj=ob_tran +o_proj=latlon +o_lon_p=0 +o_lat_p=42.5 +lon_0=263  +a=1 +to_meter=0.0174532925199 
-127.163970947 12.5568675995 
-33.880011	-28.379995 

$ invproj -f "%f" -m 57.295779506 +proj=ob_tran +o_proj=latlon +o_lon_p=83 +o_lat_p=42.5 +lon_0=180
-127.163970947 12.5568675995 
-33.880011	-28.379995 

I will also attach a small file for which both methods work. Can you provide numerical examples?

Finally, would you recommend using the workaround posted in ​http://osgeo-org.1803224.n2.nabble.com/Does-GDAL-handle-rotated-spheres-td2030419.html , but using a GEOGCS and your suggested proj.4 string in EXTENSION attribute?

Adding full support for ob_tran is perhaps too involved, so this workaround might be sufficient.

[etourigny]

subset of data, original is at ​http://www.narccap.ucar.edu/data/table4/orog_HRM3.nc

[rhattersley]

Replying to etourigny:

Do you know what to do if north_pole_grid_longitude is different from 0? In other workds, what would be the inverse transform (and forward transform following your logic)?

Simply:

• o_lon_p = north_pole_grid_longitude

I will also attach a small file for which both methods work. Can you provide numerical examples?

What kind of numerical examples are you after?

Finally, would you recommend using the workaround posted in ​http://osgeo-org.1803224.n2.nabble.com/Does-GDAL-handle-rotated-spheres-td2030419.html , but using a GEOGCS and your suggested proj.4 string in EXTENSION attribute?

Sorry - I've not tried using this with WKT/GDAL so I can't really comment. (I'm currently working via pyproj and cs2cs.) I have an interest in extending to WKT (and quite possibly GDAL) though, so I'll post any relevant findings.

[etourigny]

Replying to rhattersley:

Simply:

• o_lon_p = north_pole_grid_longitude

OK. Any idea how the "normal" (inverse) transformation would work in this case?

To be clear in your suggested configuration we have:

• o_lon_p = north_pole_grid_longitude
• o_lat_p = grid_north_pole_latitude
• lon_0 = 180 + grid_north_pole_longitude

(e.g. +proj=ob_tran +o_proj=latlon +o_lon_p=0 +o_lat_p=42.5 +lon_0=263 )

In the "normal" inverse configuration we have

• o_lon_p=grid_north_pole_longitude
• o_lat_p=grid_north_pole_latitude
• lon_0 = 180

(e.g. +proj=ob_tran +o_proj=latlon +o_lon_p=83 +o_lat_p=42.5 +lon_0=180 )

where does north_pole_grid_longitude fit in here?

Thanks

[peifer]

I might be overlooking something, but isn't the issue as simple as follows:

• either one uses ob_tran transformation according to the book [1], i.e. the position of the rotated pole is given, in normal lon/lat coordinates. Then proj/invproj behave in a somewhat counter-intuitive way;
• or one uses somewhat counter-intuitive transformation parameters (as suggested by rhattersley), i.e. the position of the standard pole, in rotated lon/lat coordinates. Then proj/invproj behave "normally";
• both options do not really help us any further, as gdal/ogr/ogr_srs_proj4.cpp doesn't support ob_tran
• neither will any WKT EXTENSION help: "+wktext cannot help when OSR doesn't recognize the projection. It only helps to keep the extra options that OSR doesn't handle by itself." (Quote from Evan Rouault at [2])

[1] ​http://home.comcast.net/~gevenden56/proj/manual.pdf
[2] ​http://osgeo-org.1803224.n2.nabble.com/hammer-parameters-td5911561.html

[Even Rouault]

Herman,

I don't remember in which context I said originally the sentence you quote, but it deserves further explanation because it might be true or wrong depending on what you really want to do.

What Frank suggests in ​http://osgeo-org.1803224.n2.nabble.com/Does-GDAL-handle-rotated-spheres-td2030419.html actually works : you can build any WKT recognized by OSR, and add a EXTENSION PROJ4 to it with a string. When OSR transform the spatial reference object to a proj4 string, it will then just extract the string you provided. You can for example have a WKT that says "I'm a UTM projection" and add a PROJ4 EXTENSION that says "I'm a mercator projection". The latter will win when the SRS is translated to proj4 world.

But, to go back to my quote, if you pass "+proj=ob_tran +o_proj=eqc +lon_0=-40 +o_lat_p=22 +x_0=639367 +y_0=1473324 +a=6371000 +b=6371000 +wktext" to OSRImportFromProj4(), it will fail "of course" because there's no code in OSR to map ob_tran to something OSR knows.

Hope things are clearer...

I've not followed all the story, but it looks like it is for some netCDF files that are expressed in the ob_tran SRS. If a somehow normalized form of expression ob_tran as WKT exists, then I'd suggest to implement it into OSR. That would be the cleaner approach. If there's none, and there's no approaching form of WKT that could be used as an approximation, then I'd suggest to build a fake WKT (perhaps with a LOCAL_CS top node) that cannot be used in any other context than reprojection with GDAL.

[peifer]

Replying to rouault:

... it will fail "of course" because there's no code in OSR to map ob_tran to something OSR knows.

Unfortunately, this is exactly our issue. Thanks for your suggestions. I'm not so sure if there is something like a "normalised ob_tran SRS". ob_tran is an oblique transformation of basically any projection.

[etourigny]

Replying to rouault:

What Frank suggests in ​http://osgeo-org.1803224.n2.nabble.com/Does-GDAL-handle-rotated-spheres-td2030419.html actually works : you can build any WKT recognized by OSR, and add a EXTENSION PROJ4 to it with a string. When OSR transform the spatial reference object to a proj4 string, it will then just extract the string you provided. You can for example have a WKT that says "I'm a UTM projection" and add a PROJ4 EXTENSION that says "I'm a mercator projection". The latter will win when the SRS is translated to proj4 world.

If I understand correctly, you can put anything in the WKT before the EXTENSION? (say WGS 84, or just a LOCAL_CS node?)

Or it it better/essential that you encode the "target" srs in the wkt? (i.e. in our case a lat/lon CRS, and in Frank's suggestion a Plate_Carree CRS).

I've not followed all the story, but it looks like it is for some netCDF files that are expressed in the ob_tran SRS.

Something along those lines. The actual files have a rotated pole crs which is supported in proj through +ob_tran

  char rotated_pole
    rotated_pole:grid_mapping_name = "rotated_latitude_longitude" ;
    rotated_pole:grid_north_pole_latitude = 32.5 ;
    rotated_pole:grid_north_pole_longitude = 170. ;

If a somehow normalized form of expression ob_tran as WKT exists, then I'd suggest to implement it into OSR. That would be the cleaner approach.

As far as I understand it and as Hermann pointed out, ob_tran is a transformation rather than a CRS, how can that fit into WKT? Are there any similar examples in OGR/WKT?

How about the WKT Parameterized Transforms, specifically the "Affine Transform"?

If there's none, and there's no approaching form of WKT that could be used as an approximation, then I'd suggest to build a fake WKT (perhaps with a LOCAL_CS top node) that cannot be used in any other context than reprojection with GDAL.

[Even Rouault]

If I understand correctly, you can put anything in the WKT before the EXTENSION? (say WGS 84, or just a LOCAL_CS node?)

Yes

Or it it better/essential that you encode the "target" srs in the wkt? (i.e. in our case a lat/lon CRS, and in Frank's suggestion a Plate_Carree CRS).

I have no strong mind on this. There are pros/cons. Pros : it gives an idea of the target SRS... Cons : If you encode the "target" SRS and the user doesn't do any warping, then it might be confusing, for example if he just translates to GTiff. The PROJ4 EXTENSION will be lost, and the SRS reported on the GTiff will then be wrong. Hence my suggestion to use a WKT that obviously shows that it is not a really WKT.

As far as I understand it and as Hermann pointed out, ob_tran is a transformation rather than a CRS, how can that fit into WKT? Are there any similar examples in OGR/WKT?

Not that I am aware of.

How about the WKT Parameterized Transforms, specifically the "Affine Transform"?

I had to actually Google to know what it is. AFAIK, there's no support for it in current OSR implementation. Could perhaps be a good enhancement ! But not sure if it is widely supported elsewhere.

[rbarnes]

​This gis.stackexchange question links here.

[Even Rouault]

In 40777:

netCDF: add support for reading files in rotated pole projection (fixes #4285)