gdalwarp to Stereographic projection using proj4 string produces Oblique_Stereographic projection

[pds]

When attempting to reproject a GeoTIFF raster to Stereographic SRS using a Proj4 string, the actual resulting Raster has an Oblique_Stereographic WKT.

Command used:

gdalwarp -t_srs '+proj=stere +lat_ts=-37 +lat_0=-90 +lon_0=145 +k_0=1.0 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs' melb-small.tif savedResults/Melb_USGS/refactored_gdal/melb-small_PSt.tif

Relevant projection/SRS gdalinfo on resulting raster:

PROJCS["unnamed",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Oblique_Stereographic"],
    PARAMETER["latitude_of_origin",-37],
    PARAMETER["central_meridian",145],
    PARAMETER["scale_factor",1],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]]]

[pds]

Update: the same thing happens even if I specify a WKT asking for Stereographic:

gdalwarp -t_srs 'PROJCS["unnamed", GEOGCS["WGS 84", DATUM["WGS_1984", SPHEROID["WGS 84",6378137,298.257223563, AUTHORITY["EPSG","7030"]], AUTHORITY["EPSG","6326"]], PRIMEM["Greenwich",0], UNIT["degree",0.0174532925199433], AUTHORITY["EPSG","4326"]], PROJECTION["Stereographic"], PARAMETER["latitude_of_origin",-37.5], PARAMETER["central_meridian",145], PARAMETER["scale_factor",1], PARAMETER["false_easting",0], PARAMETER["false_northing",0], UNIT["metre",1, AUTHORITY["EPSG","9001"]]]' melb-small.tif savedResults/Melb_USGS/refactored_gdal/melb-small_St.tif

Output WKT:

PROJCS["unnamed",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Oblique_Stereographic"],
    PARAMETER["latitude_of_origin",-37.5],
    PARAMETER["central_meridian",145],
    PARAMETER["scale_factor",1],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]]]

[pds]

Sample input file used

[pds]

Note: cause seems to be in GeoTiff driver, line 704:

         case CT_Stereographic:
            oSRS.SetOS( adfParm[0], adfParm[1],
                        adfParm[4],
                        adfParm[5], adfParm[6] );
            break;

[pds]

Possible patch to overcome this issue writing Stereographic proj to GeoTIFF

[etourigny]

The patch works for me.

I'm not commiting the change, as it has to be commited upstream also.

[Even Rouault]

Note : GDAL is its own "upstream" for frmts/gtiff/gt_wkt_srs.cpp.

Reading ​http://www.remotesensing.org/geotiff/proj_list/stereographic.html and ​http://www.remotesensing.org/geotiff/proj_list/oblique_stereographic.html might help understanding why CT_Stereographic had been historically mapped to SetOS(). Although the pages would suggest indeed to map CT_Stereographic to WKT Stereographic and CT_ObliqueStereographic to WKT ObliqueStereographic. The difficulty seems to be rather when translating this to/from proj4 string.

[etourigny]

It is a bit confusing indeed...

If the WKT string specifies Stereographic then it should indeed be recognized as such, especially since it is valid WKT.

This causes a problem in the netcdf driver with the CF "Stereographic" projection, as "ObliqueStereographic" is not recognized.

The PROJ.4 representation is not as well defined.

Here is how it should be (I think):

Polar stereographic: '+proj=stere +lat_ts=-37 +lat_0=-90 +lon_0=145 +k_0=1.0 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs'

stereographic: '+proj=stere +lat_0=-37 +lon_0=145 +k_0=1.0 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs'

oblique stereographic: '+proj=sterea +lat_0=-37 +lon_0=145 +k_0=1.0 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs'

There a 3 problems

1) in frmts/gtiff/gt_wkt_srs.cpp

Stereographic is set to ObliqueStereographic The proposed patch solves that problem

2) in ogr/ogr_srs_proj4.cpp (lines 549 - 567)

The distinction between Stereo and ObliqueStereo is in the 'k' parameter, whereas it should be in the proj argument (stere/sterea):

Furthermore, the "k" parameter is always set to 1 for Stereographic, although one should be able to change it ("Defaults to 1.0 when not available").

Attaching a small fix, in addition to Patrick's proposed patch.

'+proj=stere +lat_0=-37 +lon_0=145 +k=1 +x_0=...' gets mapped to OpliqueStereographic, (k parameter) although it should be Stereographic

'+proj=stere +lat_0=-37 +lon_0=145 +x_0=...' gets mapped to Stereographic (OK but no k parameter)

'+proj=sterea +lat_0=-37 +lon_0=145 +k=1 +x_0=...' gets mapped to ObliqueStereographic (OK but have to specify the k parameter)

To illustrate my point, an ObliqueStereo WKT is represented by a PROJ.4 string ( as output by OSRExportToProj4() ) as +proj=sterea .

Example from a modified gdalinfo output

Coordinate System is:
PROJCS["unnamed",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Oblique_Stereographic"],
    PARAMETER["latitude_of_origin",-37],
    PARAMETER["central_meridian",145],
    PARAMETER["scale_factor",0.99],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]]]
PROJ.4 string is:
'+proj=sterea +lat_0=-37 +lon_0=145 +k=0.99 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs '

3) in ogr/ogr_srs_proj4.cpp

The old name of Scaling factor 'k' is used instead of the new 'k_0' standard (c.f. the remotesensing.org pages mentioned above).

Should I file a bug about this?

Ideally both would be accepted as input, but k_0 would be set on output of OSRExportToProj4().

[etourigny]

PROJ.4 makes the clear distinction between 3 oblique projections (stere, sterea, rouss) with different formula:

1) "Stereographic" (+proj=stere), with it 4 variants (OBLIQ, EQUIT,S_POLE and N_POLE), with equations for the 4 variants (file PJ_stere.c) from USGS (Snyder)

2) "Oblique Stereographic Alternative" (+proj=sterea), also known as "Oblique Stereographic" (EPSG) and "Double Stereographic" (ESRI), with a different formula than the OBLIQ variant from USGS (file PJ_sterea.c)

3) "Roussilhe Stereographic" (+proj=rouss), from Roussilhe (1922), not supported by GDAL (file proj_rouss.c)

Therefore, it makes good sense to separate "Stereographic" (+stere) from "Oblique_Stereographic" ("Double Stereographic") (+sterea) , both in WKT and PROJ.4 forms, so that the appropriate transformation is used by libproj. The translation is rather straightforward, and the 2 patches resolve this issue.

Quoting Gerald Evenden in ​http://osgeo-org.1803224.n2.nabble.com/Implementation-s-for-Oblique-Stereographic-td2062321.html

Libproj4 provides for three (3) different Oblique Stereographics projections: "+proj=stere" for the USGS version (Snyder), "+proj=sterea" for the Double Projection method and "+proj=rouss" based upon Roussilhe's 1922 paper.

Unfortunately the bug report (PROJ.4 ???) mentioned in that thread has disappeared: ​http://bugzilla.remotesensing.org/show_bug.cgi?id=980

[etourigny]

related: bug #1428 (Oblique_ vs Double_Stereographic)

[etourigny]

I think we that we should apply the following map (no changes to PS)

ESRI WKT                   OGC WKT                      PROJ.4

Double_Stereographic       Oblique_Stereographic      +proj=sterea
Stereographic              Stereographic              +proj=stere
Stereographic_South_Pole   Polar_stereographic        +proj=stere +lat_ts=-37 +lat_0=-90 ...
Stereographic_North_Pole   Polar_stereographic        +proj=stere +lat_ts=37 +lat_0=90 ...

Attached patch fixes proj<->ogc_wkt and proj->ogc_wkt->esri_wkt , need esri files for Double_Stereographic and Stereographic to test.

[etourigny]

Tested latest patch with ESRI WKT file generated with the following:

gdal_translate -a_srs 'ESRI::PROJCS["Double_Stereographic",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],PROJECTIONDouble_Stereographic,PARAMETER["latitude_of_origin",-37],PARAMETER["central_meridian",145],PARAMETER["scale_factor",1],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["Meter",1]] ' melb-small.tif sta-esri.tif

the ESRI WKT was generated with morphtoESRI(), it would be best to have a real ESRI (one Stereographic and one Double_Stereographic) to test.

[etourigny]

patch with esri import

[pds]

Replying to etourigny:

Tested latest patch with ESRI WKT file generated with the following:

gdal_translate -a_srs 'ESRI::PROJCS["Double_Stereographic",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],PROJECTIONDouble_Stereographic,PARAMETER["latitude_of_origin",-37],PARAMETER["central_meridian",145],PARAMETER["scale_factor",1],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["Meter",1]] ' melb-small.tif sta-esri.tif

the ESRI WKT was generated with morphtoESRI(), it would be best to have a real ESRI (one Stereographic and one Double_Stereographic) to test.

There is ESRI documentation online on what parameters are needed for each proj system, see e.g. ​http://webhelp.esri.com/arcgiSDEsktop/9.3/index.cfm?TopicName=Double_Stereographic

When I was testing the Orthographic proj in ESRI, somewhere found a pointer to a more specific API reference to each of their PE strings, but haven't been able to locate this again.

[etourigny]

Fixed in trunk (r23246) and added new autotest osr_esri_20, which tests for remapping (ESR-OGC-PROJ.4) for Stereographic, Oblique_Stereographic and Polar_Stereographic.

Modified autotest/osr/osr_esri.py gdal/frmts/gtiff/gt_wkt_srs.cpp gdal/ogr/ogr_srs_esri.cpp gdal/ogr/ogr_srs_proj4.cpp

Hermann, thanks for the files, I used the projection definitions from 3 of your files, with minor modifications. Can you please test the fix with actual data?

[peifer]

Replying to etourigny:

Can you please test the fix with actual data?

I tested with 29 shapefiles, using (0). The shapefiles originate from the Netherlands and Romania, their .prj information is along the lines of "RD New.prj" and "Stereo 1970.prj" respectively.

The earlier error message is gone now (1) and the source SRS is properly detected, see (2). It is good to know that this long-standing issue has now been fixed. Thanks.

(0) ogr2ogr -t_srs epsg:4326 ...

(1) ERROR 6: No translation for Double_Stereographic to PROJ.4 format is known.

(2) OGRCT: Source: +proj=sterea +lat_0=...

[etourigny]

ok, will close this issue. Just to be sure 100% can you confirm that the transformation is adequate numerically (i.e. the points in your shapefiles end up where they should be), in addition to the source srs being properly detected?

regards, Etienne

[peifer]

Replying to etourigny:

can you confirm that the transformation is adequate numerically (i.e. the points in your shapefiles end up where they should be)

I have no appropriate data for doing any high precision checking (in terms of centimetres and millimetres), but the transformation results look OK within +/- 10 metres.

[etourigny]

great! many thanks for your input