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