Opened 6 years ago

Closed 6 years ago

Last modified 6 years ago

#4267 closed defect (fixed)

gdalwarp to Stereographic projection using proj4 string produces Oblique_Stereographic projection

Reported by: pds Owned by: etourigny
Priority: normal Milestone:
Component: default Version: svn-trunk
Severity: normal Keywords: stereographic, oblique_stereographic, proj4
Cc: peifer

Description

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"]]]

Attachments (6)

melb-small.tif (15.1 KB) - added by pds 6 years ago.
Sample input file used
geoTiffStereographic-patch.diff (489 bytes) - added by pds 6 years ago.
Possible patch to overcome this issue writing Stereographic proj to GeoTIFF
patch2.txt (1.2 KB) - added by etourigny 6 years ago.
patch-stereo.txt (3.0 KB) - added by etourigny 6 years ago.
patch-stereo2.txt (4.4 KB) - added by etourigny 6 years ago.
patch with esri import
stereo.zip (16.2 KB) - added by peifer 6 years ago.

Download all attachments as: .zip

Change History (21)

comment:1 Changed 6 years ago by 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"]]]

Changed 6 years ago by pds

Attachment: melb-small.tif added

Sample input file used

comment:2 Changed 6 years ago by 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;

Changed 6 years ago by pds

Possible patch to overcome this issue writing Stereographic proj to GeoTIFF

comment:3 Changed 6 years ago by etourigny

The patch works for me.

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

comment:4 Changed 6 years ago by 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.

comment:5 Changed 6 years ago by 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().

Changed 6 years ago by etourigny

Attachment: patch2.txt added

comment:6 Changed 6 years ago by 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

comment:7 Changed 6 years ago by etourigny

related: bug #1428 (Oblique_ vs Double_Stereographic)

comment:8 Changed 6 years ago by etourigny

Cc: peifer added

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.

Changed 6 years ago by etourigny

Attachment: patch-stereo.txt added

comment:9 Changed 6 years ago by etourigny

Owner: changed from warmerdam 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.

Changed 6 years ago by etourigny

Attachment: patch-stereo2.txt added

patch with esri import

comment:10 in reply to:  9 Changed 6 years ago by 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.

Changed 6 years ago by peifer

Attachment: stereo.zip added

comment:11 Changed 6 years ago by 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?

comment:12 in reply to:  11 Changed 6 years ago by 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=...

comment:13 Changed 6 years ago by etourigny

Resolution: fixed
Status: newclosed

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

comment:14 in reply to:  13 Changed 6 years ago by 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.

comment:15 Changed 6 years ago by etourigny

great! many thanks for your input

Note: See TracTickets for help on using tickets.