Opened 14 years ago

Closed 9 years ago

#3731 closed defect (invalid)

Equirectangular ProjCenterLatGeoKey handling wrong for Mars imagery?

Reported by: melanopsis Owned by: warmerdam
Priority: normal Milestone:
Component: OGR_SRS Version: 1.7.2
Severity: normal Keywords: GTiff Equirectangular
Cc:

Description (last modified by melanopsis)

We have some imagery from Mars that we try to process using a custom software being developed in our department. Our pre-processor program relies heavily on gdal and it seems that there is a problem with coordinate calculations in DMS. I've attached both listgeo and gdalinfo outputs to this message and the sample data can be found at anatolia.geology.ucdavis.edu/misc/mars_clip.tif

Thank you

Attachments (2)

listgeo.txt (2.0 KB ) - added by melanopsis 14 years ago.
gdal_info.txt (1.2 KB ) - added by melanopsis 14 years ago.

Download all attachments as: .zip

Change History (17)

by melanopsis, 14 years ago

Attachment: listgeo.txt added

by melanopsis, 14 years ago

Attachment: gdal_info.txt added

comment:1 by melanopsis, 14 years ago

Description: modified (diff)

comment:2 by melanopsis, 14 years ago

Version: unspecified1.7.2

comment:3 by Even Rouault, 14 years ago

Resolution: invalid
Status: newclosed

You didn't specify which result you believed was right, but I'd presume that GDAL 1.7.2 is right. You presumably use libgeotiff < 1.3.0, since libgeotiff 1.3.0 had a fix to produce the same result as GDAL. See http://trac.osgeo.org/geotiff/changeset/1652/trunk/libgeotiff/geotiff_proj4.c. With libgeotiff 1.3.0, listgeo outputs the same coordinates in DMS as gdalinfo.

See also http://trac.osgeo.org/gdal/ticket/2706

comment:4 by melanopsis, 14 years ago

Thank you for the quick response and the links. Frankly, I thought listgeo's coordinates were correct since they were in agreement with both ENVI and Google Mars. Then I tried ArcGIS and ArcGIS gave different coordinates. Using the coordinate values in meters, I made some quick calculations and my results were in agreement with ArcGIS. So in summary;

UL coordinates:

gdalinfo 1.7.2 --> 28d20'45.59"W, 44d13'0.13"N listgeo with libgeotiff 1.3.0 --> 28d20'45.59"W, 44d13'0.13"N listgeo with libgeotiff 1.2.5 --> 18d36'47.46"W, 24d13'0.13"N ENVI 4.7 --> 18d36'47.46"W, 24d13'0.13"N Google Mars --> 18d36'47.46"W, 24d13'0.13"N ArcGIS 9.3 --> 28d20'45.59"W 24d13'0.13"N

I would appreciate if you could investigate this further.

Thank you

comment:5 by melanopsis, 14 years ago

Resolution: invalid
Status: closedreopened

comment:6 by Even Rouault, 14 years ago

Hum, I don't understand why you reopened the ticket if your results are in agreement with arcgis, and if arcgis and gdal results are identical...

in reply to:  6 comment:7 by melanopsis, 14 years ago

Replying to rouault:

Hum, I don't understand why you reopened the ticket if your results are in agreement with arcgis, and if arcgis and gdal results are identical...

ArcGIS and gdal are not in agreement. latitude value is 20 degrees off.

comment:8 by Even Rouault, 14 years ago

Ah right I looked too quickly. The 20 degree difference in latitude is likely due to the +lat_0=20 that is found in the translation of the SRS to proj.4 strings by GDAL : +proj=eqc +lat_ts=0 +lat_0=20 +lon_0=180 +x_0=0 +y_0=0 +a=3393833.2607584 +b=3393833.2607584 +units=m +no_defs

Well I'm afraid I can't help more on those projection issues. Perhaps someone else will be interested in investigating into this.

in reply to:  8 comment:9 by melanopsis, 14 years ago

Replying to rouault:

Ah right I looked too quickly. The 20 degree difference in latitude is likely due to the +lat_0=20 that is found in the translation of the SRS to proj.4 strings by GDAL : +proj=eqc +lat_ts=0 +lat_0=20 +lon_0=180 +x_0=0 +y_0=0 +a=3393833.2607584 +b=3393833.2607584 +units=m +no_defs

Well I'm afraid I can't help more on those projection issues. Perhaps someone else will be interested in investigating into this.

Well, yes, that's the puzzling part. For some reason the value of lat_0 is added to lat value of the point.

Thanks for your help

comment:10 by melanopsis, 14 years ago

to follow up on this, I revisited the calculations that I made earlier and it seems that I miscalculated the longitude value (a colleague also confirmed these results). The final values are;

18d36'47.46"W, 24d13'0.13"N.

so now with gdal, in addition to lat value being 20 degrees off, It seems that lon value is also about 10 degrees off (gdal --> 28d20'45.59"W, 44d13'0.13"N).

comment:11 by warmerdam, 14 years ago

Summary: Coordinates in DMS calculated wrongly for Mars imageryEquirectangular ProjCenterLatGeoKey handling wrong for Mars imagery?

comment:12 by warmerdam, 14 years ago

Component: defaultOGR_SRS
Keywords: GTiff Equirectangular added

GeoTIFF coordinate system essentials:

      ProjCoordTransGeoKey (Short,1): CT_Equirectangular
      ProjLinearUnitsGeoKey (Short,1): Linear_Meter
      ProjStdParallel1GeoKey (Double,1): 0                
      ProjFalseEastingGeoKey (Double,1): 0                
      ProjFalseNorthingGeoKey (Double,1): 0                
      ProjCenterLongGeoKey (Double,1): 180              
      ProjCenterLatGeoKey (Double,1): 20         

Current transformation to WKT and PROJ.4:

PROJCS["unnamed",
    GEOGCS[,
        DATUM["unknown",
            SPHEROID["unnamed",3393833.2607584,0]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433]],
    PROJECTION["Equirectangular"],
    PARAMETER["latitude_of_origin",20],
    PARAMETER["central_meridian",180],
    PARAMETER["standard_parallel_1",0],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]]]

 +proj=eqc +lat_ts=0 +lat_0=20 +lon_0=180 +x_0=0 +y_0=0 +a=3393833.2607584 +b=3393833.2607584 +units=m +no_defs

However, I see if instead of treating the 20degree ProjCenter as a lat_0 value, I treat it instead as a latitude of true scale (+lat_ts) I get the values you expect.

cs2cs +proj=eqc +lat_ts=20 +lat_0=0 +lon_0=180 +x_0=0 +y_0=0 \
      +a=3393833.2607584 +b=3393833.2607584 +units=m +no_defs +to +proj=latlong +datum=WGS84
8983007.125 1434441.625
18d36'47.456"W  24d13'0.13"N 0.000

So the issue revolves around how the ProjCenterGeoKey in a geotiff file should be interpreted. Is it at "latitude origin" or is it a "latitude of true scale".

Unfortunately the geotiff specification is grossly negligent in specifying the parameters for various projections so we are left operating based on best practice and common usage. The GDAL practice is to look for the latitude of true scale as a ProjStdParallel1GeoKey and for the latitude origin in ProjCenterLatGeoKey. This behavior has, I believe, been the result of adjustments made for past tickets though I haven't looked them up at this time.

So, I'd like to throw a question back - why do you think the geotiff formulation is appropriate?

I'd note that ArcGIS and, I think, Google Earth both use various versions of GDAL for interpreting geotiff files so your comparisons may be just comparisons of different versions of GDAL and how the GDAL coordinate system was mapped into the projection engines of Google Earth and ArcGIS rather than truely independent test points.

comment:13 by melanopsis, 14 years ago

Well, frankly, I had no idea where the problem might be and I assumed the dataset we have was compiled properly. It seems that latitude_of_origin and latitude_of_true scale are also interchangeably used which in my opinion is wrong. In fact, for this projection lat_0 should always set to be zero (since eqc doesn't require lat_0 defined) and only lat_ts/lat_1 and lon_0 should be specified. According to wikipedia;

x = lambda * cos(phi_1)
y = phi

where

lambda is the longitude from the central meridian of the projection
phi is the latitude
phi_1 are the standard parallels (north and south of the equator) where the scale of the projection is true. 

Only solution that I can think of is to check whether lat_0 == 0 and if not, set it to zero and move the value from lat_0 to lat_ts.

Thank you for your help, it was very informative...

comment:14 by Jukka Rahkonen, 9 years ago

Was the conclusion that GDAL had no bug but the reason for the troubles were due to confusion in which GeoTIFF tags should be used for defining the Martial projection? In that case the ticket can be closed again.

comment:15 by melanopsis, 9 years ago

Description: modified (diff)
Resolution: invalid
Status: reopenedclosed

Yes, this can be closed now.

Note: See TracTickets for help on using tickets.