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 )
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)
Change History (17)
by , 14 years ago
Attachment: | listgeo.txt added |
---|
by , 14 years ago
Attachment: | gdal_info.txt added |
---|
comment:1 by , 14 years ago
Description: | modified (diff) |
---|
comment:2 by , 14 years ago
Version: | unspecified → 1.7.2 |
---|
comment:3 by , 14 years ago
Resolution: | → invalid |
---|---|
Status: | new → closed |
comment:4 by , 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 , 14 years ago
Resolution: | invalid |
---|---|
Status: | closed → reopened |
follow-up: 7 comment:6 by , 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...
comment:7 by , 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.
follow-up: 9 comment:8 by , 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.
comment:9 by , 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 , 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 , 14 years ago
Summary: | Coordinates in DMS calculated wrongly for Mars imagery → Equirectangular ProjCenterLatGeoKey handling wrong for Mars imagery? |
---|
comment:12 by , 14 years ago
Component: | default → OGR_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 , 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 , 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 , 9 years ago
Description: | modified (diff) |
---|---|
Resolution: | → invalid |
Status: | reopened → closed |
Yes, this can be closed now.
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