Opened 10 years ago

Closed 10 years ago

#2478 closed defect (fixed)

lat error from unfamiliar keyword "pseudo_standard_parallel_1"

Reported by: thare Owned by: warmerdam
Priority: highest Milestone: 1.6.0
Component: OGR_SRS Version: svn-trunk
Severity: normal Keywords: equirectangular pseudo_standard_parallel_1
Cc: shane@…, dgrichard, Even Rouault

Description (last modified by warmerdam)

There is what appears to be a new keywrod listed as "pseudo_standard_parallel_1" for the Equirectangular projection. I think this has caused a problem in GDAL? There is the original "latitude_of_origin" which appears to be correct. Now the "pseudo..." keyword is the same value as the "latitude_of_origin" and this keyword seems to be the cause for the lat range to off by exactly the same amount (in the example off by 10 degrees).

-Trent


Correct center lat,lon range of image:
Latitude (centered): 10.6 Longitude (East): 156.4

Gdal 1.6.0dev reports the center as:
Center      (-1374994.500,  630301.250) (156d26'27.39"E, 20d38'7.68"N)

File: http://hirise-pds.lpl.arizona.edu/download/PDS/RDR/PSP/ORB_008100_008199/PSP_008199_1910/PSP_008199_1910_COLOR.JP2 (355MB)

Using GDAL 1.6.0dev built inside FWTools 2.2.1.
-------------

> gdalinfo PSP_008199_1910_COLOR.jp2

...

PROJCS["Equirectangular MARS",
    GEOGCS["GCS_MARS",
        DATUM["unknown",
            SPHEROID["unnamed",3395582.0270805,0]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433]],
    PROJECTION["Equirectangular"],
    PARAMETER["latitude_of_origin",10],
    PARAMETER["central_meridian",180],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    PARAMETER["pseudo_standard_parallel_1",10],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]]]
Origin = (-1376386.875000000000000,637256.375000000000000)
Pixel Size = (0.250000000000000,-0.250000000000000)

Corner Coordinates:
Upper Left  (-1376386.875,  637256.375) (156d25'1.51"E, 20d45'10.17"N)  //Note the lat range is off by 10 degrees
Lower Left  (-1376386.875,  623346.125) (156d25'1.51"E, 20d31'5.19"N)
Upper Right (-1373602.125,  637256.375) (156d27'53.28"E, 20d45'10.17"N)
Lower Right (-1373602.125,  623346.125) (156d27'53.28"E, 20d31'5.19"N)
Center      (-1374994.500,  630301.250) (156d26'27.39"E, 20d38'7.68"N)

...

Attachments (23)

ogrspatialreference.diff (1.9 KB) - added by dgrichard 10 years ago.
setEquirectangular() proposed fix
ogr_srs_proj4.diff (1.5 KB) - added by dgrichard 10 years ago.
exportToProj4() proposed fix
gt_wkt_srs.diff (1.8 KB) - added by dgrichard 10 years ago.
GTiff driver proposed fix
geo_normalize.diff (1.2 KB) - added by dgrichard 10 years ago.
GTiff driver proposed fix
isis2dataset.diff (744 bytes) - added by dgrichard 10 years ago.
ISIS2 driver proposed fix
isis3dataset.diff (675 bytes) - added by dgrichard 10 years ago.
ISI3 driver proposed fix
geoconcept_syscoord.diff (1.8 KB) - added by dgrichard 10 years ago.
Geoconcept driver proposed fix
ogr_srs_api.diff (747 bytes) - added by dgrichard 10 years ago.
setEquirectangular2() proposed fix (intermediate)
ogr_spatialref.diff (800 bytes) - added by dgrichard 10 years ago.
setEquirectangular2() proposed fix (intermediate)
osr.diff (815 bytes) - added by dgrichard 10 years ago.
SWIG osr interface SetEquirectangular2() proposed fix
pdsdataset.diff (659 bytes) - added by dgrichard 10 years ago.
ISIS driver proposed fix
IdrisiDataset.diff (526 bytes) - added by dgrichard 10 years ago.
IDRISIS proposed fix
ilwiscoordinatesystem.diff (769 bytes) - added by dgrichard 10 years ago.
ILWIS proposed fix
hdf4imagedataset.diff (609 bytes) - added by dgrichard 10 years ago.
HDF4 proposed fix
nitfdataset.diff (543 bytes) - added by dgrichard 10 years ago.
NITF proposed fix
hfadataset.diff (1014 bytes) - added by dgrichard 10 years ago.
HFA proposed fix
ogr_srs_panorama.diff (648 bytes) - added by dgrichard 10 years ago.
Panorama proposed fix
ogr_srs_pci.diff (568 bytes) - added by dgrichard 10 years ago.
PCI proposed fix
ogr_srs_usgs.diff (646 bytes) - added by dgrichard 10 years ago.
USGS proposed fix
ogr_fromepsg.diff (726 bytes) - added by dgrichard 10 years ago.
Fromepsg proposed fix
gdal.diff (681 bytes) - added by dgrichard 10 years ago.
Pymod interface with GDAL proposed fix
mrsiddataset.diff (932 bytes) - added by dgrichard 10 years ago.
MrSID proposed fix
fix_old_equirectangular.patch (2.9 KB) - added by Even Rouault 10 years ago.
Patch to revert to behaviour of EquiRectangular? to the one of GDAL 1.5.0

Download all attachments as: .zip

Change History (33)

comment:1 Changed 10 years ago by warmerdam

Cc: dgrichard added
Description: modified (diff)
Keywords: GeoJP2 added; GDAL 1.6.0dev removed
Status: newassigned
Version: unspecifiedsvn-trunk

Trent,

Would there happen to be any smaller files demonstrating this problem?

It appears SetEquidistantRectangular2() was introduced in r13482 for ticket #2134. I'm adding dgrichard to the cc.

comment:2 Changed 10 years ago by thare

Here is an ISIS2 example. It appears to effect all image using Equirectangular. http://astrogeology.usgs.gov/Projects/MDIM21/data/ME70N180E.CUB.GZ 0.4MB

thanks, Trent

In this example the correct lat range is:

MINIMUM_LATITUDE = 65.0 MAXIMUM_LATITUDE = 75.0 WESTERNMOST_LONGITUDE = 0.0 EASTERNMOST_LONGITUDE = 360.0 MAP_RESOLUTION = 16.0 MAP_SCALE = 3.68827819912536 MAP_PROJECTION_TYPE = "EQUIRECTANGULAR_CYLINDRICAL" CENTER_LONGITUDE = 180.0 CENTER_LATITUDE = 60.0

since the projection issue I think adds it over -+90 the values get set to zero in gdalinfo: E:\GIS\Mars\MRO\HIRISE_test_jp2>gdalinfo me70N180E.cub Driver: ISIS2/USGS Astrogeology ISIS cube (Version 2) Files: me70N180E.cub Size is 2880, 160 Coordinate System is: PROJCS["EQUIRECTANGULAR_CYLINDRICAL MARS",

GEOGCS["GCS_MARS",

DATUM["D_MARS",

SPHEROID["MARS_localRadius",3381164.391679916,0]],

PRIMEM["Reference_Meridian",0], UNIT["degree",0.0174532925199433]],

PROJECTIONEquirectangular?, PARAMETER["latitude_of_origin",60], PARAMETER["central_meridian",180], PARAMETER["false_easting",0], PARAMETER["false_northing",0], PARAMETER["pseudo_standard_parallel_1",60]]

Origin = (-5311120.606740518500000,4425933.838950431900000) Pixel Size = (3688.278199125360100,-3688.278199125360100) Corner Coordinates: Upper Left (-5311120.607, 4425933.839) ( 0d 0'0.00"W,135d 0'0.00"N) Lower Left (-5311120.607, 3835809.327) ( 0d 0'0.00"W,125d 0'0.00"N) Upper Right ( 5311120.607, 4425933.839) ( 0d 0'0.00"E,135d 0'0.00"N) Lower Right ( 5311120.607, 3835809.327) ( 0d 0'0.00"E,125d 0'0.00"N) Center ( 0.000, 4130871.583) (180d 0'0.00"E,130d 0'0.00"N) Band 1 Block=2880x1 Type=Byte, ColorInterp?=Undefined

comment:3 Changed 10 years ago by warmerdam

Component: GDAL_RasterOGR_SRS
Keywords: equirectangular pseudo_standard_parallel_1 added; GeoJP2 removed
Priority: normalhigh

I have reproduced the problem with the file Trent references.

GDAL trunk:

Coordinate System is:
PROJCS["EQUIRECTANGULAR_CYLINDRICAL MARS",
    GEOGCS["GCS_MARS",
        DATUM["D_MARS",
            SPHEROID["MARS_localRadius",3381164.391679915,0]],
        PRIMEM["Reference_Meridian",0],
        UNIT["degree",0.0174532925199433]],
    PROJECTION["Equirectangular"],
    PARAMETER["latitude_of_origin",60],
    PARAMETER["central_meridian",180],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    PARAMETER["pseudo_standard_parallel_1",60]]
Origin = (-5311120.605468750000000,4425933.837890625000000)
Pixel Size = (3688.278198242187500,-3688.278198242187500)
OGRCT: Source: +proj=eqc +lat_ts=60 +lat_0=60 +lon_0=180 +x_0=0 +y_0=0 +a=338116
4.391679915 +b=3381164.391679915 +units=m +no_defs 
OGRCT: Target: +proj=longlat +a=3381164.391679915 +b=3381164.391679915 +no_defs 
Corner Coordinates:
Upper Left  (-5311120.605, 4425933.838) (  0d 0'0.00"E,135d 0'0.00"N)
Lower Left  (-5311120.605, 3835809.326) (  0d 0'0.00"E,125d 0'0.00"N)
Upper Right ( 5311120.605, 4425933.838) (  0d 0'0.00"W,135d 0'0.00"N)
Lower Right ( 5311120.605, 3835809.326) (  0d 0'0.00"W,125d 0'0.00"N)
Center      (       0.000, 4130871.582) (180d 0'0.00"E,130d 0'0.00"N)

GDAL 1.5:

Coordinate System is:
PROJCS["EQUIRECTANGULAR_CYLINDRICAL MARS",
    GEOGCS["GCS_MARS",
        DATUM["D_MARS",
            SPHEROID["MARS_localRadius",3381164.391679915,0]],
        PRIMEM["Reference_Meridian",0],
        UNIT["degree",0.0174532925199433]],
    PROJECTION["Equirectangular"],
    PARAMETER["latitude_of_origin",60],
    PARAMETER["central_meridian",180],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0]]
Origin = (-5311120.605468750000000,4425933.837890625000000)
Pixel Size = (3688.278198242187500,-3688.278198242187500)
OGRCT: Source: +proj=eqc +lat_ts=60 +lon_0=180 +x_0=0 +y_0=0 +a=3381164.39167991
5 +b=3381164.391679915 +units=m +no_defs 
OGRCT: Target: +proj=longlat +a=3381164.391679915 +b=3381164.391679915 +no_defs 
Corner Coordinates:
Upper Left  (-5311120.605, 4425933.838) (  0d 0'0.00"E, 75d 0'0.00"N)
Lower Left  (-5311120.605, 3835809.326) (  0d 0'0.00"E, 65d 0'0.00"N)
Upper Right ( 5311120.605, 4425933.838) (  0d 0'0.00"W, 75d 0'0.00"N)
Lower Right ( 5311120.605, 3835809.326) (  0d 0'0.00"W, 65d 0'0.00"N)
Center      (       0.000, 4130871.582) (180d 0'0.00"E, 70d 0'0.00"N)

Note the different latitudes in the corner coordinates, and the different PROJ.4 definition used in the two cases. I am thinking that we should *not* be introducing this new PSEUDO_STD_PARALLEL_1 parameter in "plain" equirectangular set with OGRSpatialReference::SetEquirectangular?(), and/or we should not be introducing +lat_ts in the exportToProj4() method.

Didier, I'll wait for some explanation from you before backing out some of your changes.

comment:4 Changed 10 years ago by dgrichard

I have add SetEquirectangular2() (and kept SetEquirectangular?()) sometime after the introduction of the generalized equidistant cylindrical projection in PROJ.4 (within GDAL, it was named Equirectangular). This projection adds an optional lat_0 parameter. I then added SetEquirectangular2() for taking this parameter into account. Using SetEquirectangular2() was a way for not modifying SetEquirectangular?(). The former gets called when lat_0!=lat_ts and both different from 0 (See ogr_srs_proj4.cpp). In all examples above, they seem to be equal and not 0 which is quite strange cause I have checked eqc projection in PROJ.4 and there was no lat_0 defined. Could it be a side effect ?

Looking at the WKT part, I see nothing wrong in the definition. Looking at the corner coordinates, it seems that the lat_0 value has been substracted from the latitude. This parameter is added by SetEquirectangular2().

If we remove PSEUDO_STD_PARALLEL_1 and +lat_0, we lost the generalization of equidistant cylindrical projection.

I am not an image expert, but I believe that there are some codes in image drivers not taking into account the new projection and then mixing things up. Shortest way : comment lines calling SetEquirectangular2() in ogr_srs_proj4.cpp to force using SetEquirectangular?(), it should fixe the error. Toughest way : find out where lat_0 is substracted and why ?

comment:5 Changed 10 years ago by dgrichard

Seems that setEquirectangular() method holds the bug in setting SRS_PP_PSEUDO_STD_PARALLEL_1 with SRS_PP_LATITUDE_OF_ORIGIN directly. As a side effect, it seems that exportToProj4() has the same problem in setting SRS_PP_PSEUDO_STD_PARALLEL_1 directly with SRS_PP_LATITUDE_OF_ORIGIN instead of setting it to 0 when not defined. Proposed patches are attached to the ticket.

Changed 10 years ago by dgrichard

Attachment: ogrspatialreference.diff added

setEquirectangular() proposed fix

Changed 10 years ago by dgrichard

Attachment: ogr_srs_proj4.diff added

exportToProj4() proposed fix

Changed 10 years ago by dgrichard

Attachment: gt_wkt_srs.diff added

GTiff driver proposed fix

Changed 10 years ago by dgrichard

Attachment: geo_normalize.diff added

GTiff driver proposed fix

Changed 10 years ago by dgrichard

Attachment: isis2dataset.diff added

ISIS2 driver proposed fix

Changed 10 years ago by dgrichard

Attachment: isis3dataset.diff added

ISI3 driver proposed fix

Changed 10 years ago by dgrichard

Attachment: geoconcept_syscoord.diff added

Geoconcept driver proposed fix

Changed 10 years ago by dgrichard

Attachment: ogr_srs_api.diff added

setEquirectangular2() proposed fix (intermediate)

Changed 10 years ago by dgrichard

Attachment: ogr_spatialref.diff added

setEquirectangular2() proposed fix (intermediate)

Changed 10 years ago by dgrichard

Attachment: osr.diff added

SWIG osr interface SetEquirectangular2() proposed fix

comment:6 Changed 10 years ago by dgrichard

I have made a deeper survey on the problem :

The use of SetEquirectangular?() is in line with PROJ.4.6.0 eqc.
The use of SetEquirectangular2() is in line with PROJ trunk eqc (generalized version with both lat_ts and lat_0).

SetEquirectangular?() switched lat_ts and lat_0 (the latter being always equal to 0). All drivers using this method have to be changed. I have modified GTiff and ISIS2/3 to demonstrate the work (See attached patches) :

$ CPL_DEBUG="ON" gdalinfo ME70N180E.CUB

ISIS2: using projection EQUIRECTANGULAR_CYLINDRICAL
ISIS2: local radius: 3381164.391680
GDALRaw: RawRasterBand(0x80584d0,1,0x8058f08,
              Off=27136,PixOff=1,LineOff=2880,Byte,1)

GDAL: GDALOpen(ME70N180E.CUB) succeeds as ISIS2.
Driver: ISIS2/USGS Astrogeology ISIS cube (Version 2)
Files: ME70N180E.CUB
Size is 2880, 160
Coordinate System is:
PROJCS["EQUIRECTANGULAR_CYLINDRICAL MARS",
    GEOGCS["GCS_MARS",
        DATUM["D_MARS",
            SPHEROID["MARS_localRadius",3381164.391679916,0]],
        PRIMEM["Reference_Meridian",0],
        UNIT["degree",0.0174532925199433]],
    PROJECTION["Equirectangular"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",180],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    PARAMETER["pseudo_standard_parallel_1",60]]
Origin = (-5311120.605468750000000,4425933.837890625000000)
Pixel Size = (3688.278198242187500,-3688.278198242187500)
OGRCT: Source: +proj=eqc +lat_ts=60 +lat_0=0 +lon_0=180 +x_0=0 +y_0=0 +a=3381164.391679916 +b=3381164.391679916 +units=m +no_defs
OGRCT: Target: +proj=longlat +a=3381164.391679916 +b=3381164.391679916 +no_defs
Corner Coordinates:
Upper Left  (-5311120.605, 4425933.838) (  0d 0'0.00"E, 75d 0'0.00"N)
Lower Left  (-5311120.605, 3835809.326) (  0d 0'0.00"E, 65d 0'0.00"N)
Upper Right ( 5311120.605, 4425933.838) (  0d 0'0.00"W, 75d 0'0.00"N)
Lower Right ( 5311120.605, 3835809.326) (  0d 0'0.00"W, 65d 0'0.00"N)
Center      (       0.000, 4130871.582) (180d 0'0.00"E, 70d 0'0.00"N)
Band 1 Block=2880x1 Type=Byte, ColorInterp=Undefined
  NoData Value=0
GDAL: GDALClose(ME70N180E.CUB)
GDAL: GDALDeregister_GTiff() called.

The plan would be to :

What do maintainers/users think about this proposal ?

Changed 10 years ago by dgrichard

Attachment: pdsdataset.diff added

ISIS driver proposed fix

Changed 10 years ago by dgrichard

Attachment: IdrisiDataset.diff added

IDRISIS proposed fix

Changed 10 years ago by dgrichard

Attachment: ilwiscoordinatesystem.diff added

ILWIS proposed fix

Changed 10 years ago by dgrichard

Attachment: hdf4imagedataset.diff added

HDF4 proposed fix

Changed 10 years ago by dgrichard

Attachment: nitfdataset.diff added

NITF proposed fix

Changed 10 years ago by dgrichard

Attachment: hfadataset.diff added

HFA proposed fix

Changed 10 years ago by dgrichard

Attachment: ogr_srs_panorama.diff added

Panorama proposed fix

Changed 10 years ago by dgrichard

Attachment: ogr_srs_pci.diff added

PCI proposed fix

Changed 10 years ago by dgrichard

Attachment: ogr_srs_usgs.diff added

USGS proposed fix

Changed 10 years ago by dgrichard

Attachment: ogr_fromepsg.diff added

Fromepsg proposed fix

Changed 10 years ago by dgrichard

Attachment: gdal.diff added

Pymod interface with GDAL proposed fix

Changed 10 years ago by dgrichard

Attachment: mrsiddataset.diff added

MrSID proposed fix

comment:7 Changed 10 years ago by Even Rouault

Frank, Didier,

I'm attaching a patch 'fix_old_equirectangular.patch' that basically reverts to the behaviour of plain EquiRectangular? to what is was before r13868. 3 things are done:

  • OGRSpatialReference::SetEquirectangular?() doesn't define SRS_PP_PSEUDO_STD_PARALLEL_1 anymore (goes back to what is was before r13868)
  • OGRSpatialReference::importFromProj4() only calls SetEquirectangular2() if lat_0 is found, otherwise it calls SetEquirectangular?()
  • OGRSpatialReference::exportToProj4() doesn't specify lat_0 if SRS_PP_PSEUDO_STD_PARALLEL_1 is not found in the projection parameter list

Disclaimer, I propose those changes from just a "logical" analysis of the changes done by r13868, but without really understanding how the projection works and the meaning of the lat_0 and lat_ts parameters, and if the changes done in proj are backward compatible or not.

I've tested it with ME70N180E.CUB dataset, proj 4.6.0 and with that patch, we get back to GDAL 1.5.0 behaviour.

Changed 10 years ago by Even Rouault

Patch to revert to behaviour of EquiRectangular? to the one of GDAL 1.5.0

comment:8 Changed 10 years ago by warmerdam

Cc: Even Rouault added
Priority: highhighest

I'm elevating the priority to Highest, and intend to treat this as an RC blocker.

I am reviewing the equirectangular / equidistant cylindrical docs from EPSG, ESRI, PROJ.4, and hopefully CSMap in an effort to find an appropriate WKT formulation for this projection method, and correct mappings to other formats. I believe the geotiff guidance on this has been in error and this has impacted GDAL.

comment:9 Changed 10 years ago by warmerdam

I have migrated the core issue of Equirectangular parameters to #2706 to break it from other topics in this ticket.

comment:10 Changed 10 years ago by warmerdam

Resolution: fixed
Status: assignedclosed
Note: See TracTickets for help on using tickets.