Hi Even,Ok, that clears things up a little. I guess the real question I have at this point is - which is correct in this circumstance?
The original GeoTIFF? ArcGIS? GDAL's conversion? None of the above? I get the feeling it's all a bit subjective. >From your comments, I'm surprised this doesn't happen more often with various >packages\files. I thought projections were nice and reasonably well defined >things. Thanks again, Jonathan ---- On Tue, 17 Nov 2015 10:33:37 -0800 Even Rouault<even.roua...@spatialys.com> wrote ---- Le mardi 17 novembre 2015 19:13:52, Jonathan Moules a écrit : > Hi Even, > Thanks for the reply. > > > I guess GDAL chose Mercator_1SP because that's the one that's explicitly > defined in the projection text. No, it is not explicitly defined, and GDAL ignores PCSCitationGeoKey to determine the projection method. It will use ProjCoordTransGeoKey = CT_Mercator for that, and in GDAL 1.11, it can only map it to Mercator_1SP > > > I've had a quick look at the two files using listgeo, and highlighted the > different lines with **** prefix below (it seems GDAL removes 6 tags, and > adds one). > > > The before: > Keyed_Information: > GTModelTypeGeoKey (Short,1): ModelTypeProjected > GTRasterTypeGeoKey (Short,1): RasterPixelIsArea > GTCitationGeoKey (Ascii,9): "Mercator" > GeographicTypeGeoKey (Short,1): GCS_WGS_84 > GeogCitationGeoKey (Ascii,7): "WGS 84" > ****GeogGeodeticDatumGeoKey (Short,1): Datum_WGS84 > ****GeogLinearUnitsGeoKey (Short,1): Linear_Meter > GeogAngularUnitsGeoKey (Short,1): Angular_Degree > ****GeogEllipsoidGeoKey (Short,1): Ellipse_WGS_84 > GeogSemiMajorAxisGeoKey (Double,1): 6378137 > ****GeogSemiMinorAxisGeoKey (Double,1): 6356752.31424518 > GeogInvFlatteningGeoKey (Double,1): 298.257223563 > ProjectedCSTypeGeoKey (Short,1): User-Defined > ****PCSCitationGeoKey (Ascii,9): "Mercator" > ProjectionGeoKey (Short,1): User-Defined > ProjCoordTransGeoKey (Short,1): CT_Mercator > ProjLinearUnitsGeoKey (Short,1): Linear_Meter > ****ProjStdParallel1GeoKey (Double,1): 60 > ProjNatOriginLongGeoKey (Double,1): 0 > ProjNatOriginLatGeoKey (Double,1): 0 > ProjFalseEastingGeoKey (Double,1): 0 > ProjFalseNorthingGeoKey (Double,1): 0 > End_Of_Keys. > > ------------ > The after: > Keyed_Information: > GTModelTypeGeoKey (Short,1): ModelTypeProjected > GTRasterTypeGeoKey (Short,1): RasterPixelIsArea > GTCitationGeoKey (Ascii,9): "Mercator" > GeographicTypeGeoKey (Short,1): GCS_WGS_84 > GeogCitationGeoKey (Ascii,7): "WGS 84" > GeogAngularUnitsGeoKey (Short,1): Angular_Degree > GeogSemiMajorAxisGeoKey (Double,1): 6378137 > GeogInvFlatteningGeoKey (Double,1): 298.257223563 > ProjectedCSTypeGeoKey (Short,1): User-Defined > ProjectionGeoKey (Short,1): User-Defined > ProjCoordTransGeoKey (Short,1): CT_Mercator > ProjLinearUnitsGeoKey (Short,1): Linear_Meter > ProjNatOriginLongGeoKey (Double,1): 0 > ProjNatOriginLatGeoKey (Double,1): 0 > ProjFalseEastingGeoKey (Double,1): 0 > ProjFalseNorthingGeoKey (Double,1): 0 > ****ProjScaleAtNatOriginGeoKey (Double,1): 1 > End_Of_Keys. > End_Of_Geotiff. > > Actually when trying that with GDAL 2.0, you would now get : PROJCS["Mercator", 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["Mercator_2SP"], PARAMETER["standard_parallel_1",60], PARAMETER["central_meridian",0], PARAMETER["false_easting",0], PARAMETER["false_northing",0], UNIT["metre",1, AUTHORITY["EPSG","9001"]]] due to the changes done in https://trac.osgeo.org/gdal/ticket/5791 The new logic to determine Mercator_1SP vs _2SP is best explainted by the code itself ;-) {{{ /* If a lat_ts was specified use 2SP, otherwise use 1SP */ if (psDefn->ProjParmId[2] == ProjStdParallel1GeoKey) { if (psDefn->ProjParmId[4] == ProjScaleAtNatOriginGeoKey) CPLError( CE_Warning, CPLE_AppDefined, "Mercator projection should not define both StdParallel1 and ScaleAtNatOrigin.\n" "Using StdParallel1 and ignoring ScaleAtNatOrigin.\n" ); oSRS.SetMercator2SP( adfParm[2], adfParm[0], adfParm[1], adfParm[5], adfParm[6]); } else oSRS.SetMercator( adfParm[0], adfParm[1], adfParm[4], adfParm[5], adfParm[6] ); }}} That new behaviour is more consistant with the geotiff keys (and should be the one you get with ArcGIS initialy), although it perhaps reveals a problem with the encoding of the geotiff file itself if you say that the Mercator_1SP interpretation with scale=1 leads to correct geopositioning. > > I'm not clear on why the things that changed were changed, but I can see > that GDAL removed the ProjStdParallel1GeoKey and value. Should not that > have been displayed by GDALSRSInfo for the Before file, even if it was > wrongly set for this given projection? GDAL is all about abstraction. It interprets at its best the information from the geotiff keys to a well-known projection method, and thus can ignore silently some keys that it doesn't need. > > > Thanks, > Jonathan > > > ---- On Tue, 17 Nov 2015 09:18:49 -0800 Even > Rouault&lt;even.roua...@spatialys.com&gt; wrote ---- > > Le mardi 17 novembre 2015 17:55:15, Jonathan Moules a écrit : > &gt; Hi List, > &gt; I have a Geotiff which includes this projection: > &gt; > &gt; PROJ.4 : '+proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +datum=WGS84 +units=m > &gt; +no_defs ' > &gt; > &gt; > &gt; OGC WKT : > &gt; PROJCS["Mercator", > &gt; GEOGCS["WGS 84", > &gt; DATUM["WGS_1984", > &gt; SPHEROID["WGS 84",6378137,298.257223563, > &gt; AUTHORITY["EPSG","7030"]], > &gt; AUTHORITY["EPSG","6326"]], > &gt; PRIMEM["Greenwich",0], > &gt; UNIT["degree",0.0174532925199433], > &gt; AUTHORITY["EPSG","4326"]], > &gt; PROJECTION["Mercator_1SP"], > &gt; PARAMETER["central_meridian",0], > &gt; PARAMETER["scale_factor",1], > &gt; PARAMETER["false_easting",0], > &gt; PARAMETER["false_northing",0], > &gt; UNIT["metre",1, > &gt; AUTHORITY["EPSG","9001"]]] > &gt; > &gt; > &gt; > &gt; > &gt; If I load this raster into ArcGIS, it displays in the wrong place (a > few &gt; thousand kilometres North). > &gt; > &gt; > &gt; I then run it through gdal_translate (GDAL 1.11.1), with no flags: > &gt; * gdal_trainslate input.tif output.tif > &gt; > &gt; For output.tif, GDALSRSInfo shows that the projection is identical, > but now &gt; the file loads correctly in ArcGIS. The same file works fine > in QGIS both &gt; before and after the "translation". > &gt; > &gt; > &gt; Looking at the projection info in ArcGIS, it displays one difference: > &gt; Before (not working): > &gt; Standard_parallel_1 = 60 > &gt; > &gt; > &gt; After (working): > &gt; Standard_parallel_1 = 0 > &gt; > &gt; > &gt; But I don't see anything about those in either of the GDALSRSInfo > outputs. &gt; > &gt; > &gt; So my questions: > &gt; - What is gdal_translate doing to the file to "fix" it? > &gt; - If it is something to do with Standard Parallel 1 - why isn't this > &gt; component of the projection exposed by GDAL? > > Yes, in Mercator_1SP, there's no Standard Parallel 1, this is for > Mercator_2SP. > > See > http://www.remotesensing.org/geotiff/proj_list/mercator_1sp.html > http://www.remotesensing.org/geotiff/proj_list/mercator_2sp.html > > I guess your original geotiff file has some unusual formulation which is > detected as Mercator_1SP by GDAL, and probably Mercator_2SP by ArCGIS. > > You could try with the listgeo utility that comes with libgeotiff to > display the geotiff keys. > > &gt; > &gt; Thoughts welcome. > &gt; Thanks, > &gt; Jonathan -- Spatialys - Geospatial professional services http://www.spatialys.com
_______________________________________________ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev