lat error from unfamiliar keyword "pseudo_standard_parallel_1"

[thare]

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)

...

[warmerdam]

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.

[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

[warmerdam]

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.

[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 ?

[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.

[dgrichard]

setEquirectangular() proposed fix

[dgrichard]

exportToProj4() proposed fix

[dgrichard]

GTiff driver proposed fix

[dgrichard]

GTiff driver proposed fix

[dgrichard]

ISIS2 driver proposed fix

[dgrichard]

ISI3 driver proposed fix

[dgrichard]

Geoconcept driver proposed fix

[dgrichard]

setEquirectangular2() proposed fix (intermediate)

[dgrichard]

setEquirectangular2() proposed fix (intermediate)

[dgrichard]

SWIG osr interface SetEquirectangular2() proposed fix

[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 :

• change all calls to SetEquirectangular() to SetEquirectangular2();
  
• test driver per driver with driver's maintainer;
  
• remove old SetEquirectangular() from interfaces and rename SetEquirectangular2() into SetEquirectangular().
  

What do maintainers/users think about this proposal ?

[dgrichard]

ISIS driver proposed fix

[dgrichard]

IDRISIS proposed fix

[dgrichard]

ILWIS proposed fix

[dgrichard]

HDF4 proposed fix

[dgrichard]

NITF proposed fix

[dgrichard]

HFA proposed fix

[dgrichard]

Panorama proposed fix

[dgrichard]

PCI proposed fix

[dgrichard]

USGS proposed fix

[dgrichard]

Fromepsg proposed fix

[dgrichard]

Pymod interface with GDAL proposed fix

[dgrichard]

MrSID proposed fix

[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.

[Even Rouault]

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

[warmerdam]

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.

[warmerdam]

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