I've never used or built GeoTools, so I didn't actually run the software to test it. I did examine the code for geotools and am fairly confident they are handling 1SP/2SP in geotiff the same way I am. If anyone uses geotools, I can provide a geotiff in 2sp.
The ticket is the same issue and I think I've resolved it in the manner proposed by Frank Warmerdam in http://trac.osgeo.org/gdal/ticket/1797#comment:4 The issue with the file in the ticket is that it specifies a latitude of origin other than the equator and does not specify a standard latitude 1. Does the latitude in the file refer to the origin, i.e. where y = 0, or the latitude of true scale, e.g. where k=1? Given the keys available in geotiff, it seems to clearly be the best that ProjNatOriginLatGeoKey should be the origin and ProjStdParallel1GeoKey should be lat_ts. However, proj4 doesn't support a latitude of origin other than 0 for Mercator, or even have a way to specify it for that matter. Thus one might consider it to not be a parameter for the projection. This leaves lat_ts as the sole latitude to be specified. One might then place it in ProjNatOriginLatGeoKey, as that is usually the first latitude parameter used in a projection. IMHO, this is clearly wrong and the result of a careless programmer just copying a pattern from another projection without considering exactly what they were doing. Another possibility when reading GeoTIFFs that specify lat_ts, is to normalize them to 1SP, i.e. lat_ts=0, and a k other than 1. However, calculating k from lat_ts requires some non-trivial math that involves the eccentricity of the ellipsoid and it seems far more appropriate to do that in proj4/ogr than in libgeotiff. It also seem like a bad idea for maintaining accuracy. For instance, NOAA chart 18440 says on the chart, "Scale 1:150,000 at Lat 47o 40'" and I produce WKT projection with, PARAMETER["standard_parallel_1",47.667] for that chart. The same number as in the chart. If that's converted to a scale factor it's not going to have the same parameter anymore. And while theoretically exact, floating point math and ASCII decimal numbers are going to result in converting lat_ts -> k -> lat_ts not producing the same value as input. On Mon, Nov 4, 2013 at 11:12 AM, Even Rouault <[email protected]> wrote: > Trent / others, > > Although your patch comes with a detailed analysis, I find myself a bit light > on the topic to blindly apply it. > > A few observations : > > * I found an old ticket that might be related > http://trac.osgeo.org/gdal/ticket/1797. Would your patch solve the issue > raised in that ticket ? > > * You mention that you followed what Geotools did. Did you actually check the > interoperability of generated files ? > > * Do people have knowledge of other software (not based on GDAL) that would > read / produce GeoTIFFs with a Mercator_2SP projection, and could provide > samples ? > > * Now on the patch itself, there's a GDAL part (gt_wkt_srs.cpp) and a > libgeotiff one (libgeotiff/geo_normalize.c). The later should be submitted in > GeoTIFF trac http://trac.osgeo.org/geotiff/newticket and approved there, > before > being downstreamed into GDAL > > * Concerning the patch geo_normalize.c, the use of NAN / isnan() could cause > portability issues on some platforms. I'd suggest using boolean flags instead. > > I've added geotiff mailing list in CC, in the hope of broadening the target > audience. And Daniele too if he has some memories of his experience with > fixing > http://jira.codehaus.org/browse/GEOT-2163 (particularly if he had the > opportunity to check interoperability with other software) > > Best regards, > > Even > > Le lundi 04 novembre 2013 02:48:26, Trent Piepho a écrit : >> Mercator_[12]SP are both the same projection in GeoTIFF. Lacking a >> definition in the official specification, the intention appears to be that >> for Mercator_2SP the latitude of true scale (lat_ts) should be specified in >> ProjStdParallel1GeoKey and for Mercator_2SP the scale at origin (k) should >> be specified in ProjScaleAtNatOriginGeoKey. >> >> Current behavior when creating a GTiff with 2SP is to drop lat_ts and >> create a default k value of 1.0, which will product an incorrect >> projection. When reading, a default k of 1.0 is used and lat_ts is >> ignored, also changing the projection. >> >> This patches changes the behavior to supply ProjScaleAtNatOriginGeoKey for >> Mercator_1SP or ProjStdParallel1GeoKey for Mercator_2SP when writing. >> >> For reading either or both values are supplied when present by libgeotiff's >> normalization code. If neither is present, a default k of 1.0 is created. >> >> GDAL's GTiff driver is changed to create a 2SP projection is lat_ts is >> present, 1SP otherwise. It emits a warning of both last_ts and scale are >> present (and ignores scale). >> --- >> gdal/frmts/gtiff/gt_wkt_srs.cpp | 28 ++++++++++--- >> gdal/frmts/gtiff/libgeotiff/geo_normalize.c | 59 >> ++++++++++++++++++++++++++- 2 files changed, 80 insertions(+), 7 >> deletions(-) >> >> diff --git a/gdal/frmts/gtiff/gt_wkt_srs.cpp >> b/gdal/frmts/gtiff/gt_wkt_srs.cpp index 5934df5..0afb101 100644 >> --- a/gdal/frmts/gtiff/gt_wkt_srs.cpp >> +++ b/gdal/frmts/gtiff/gt_wkt_srs.cpp >> @@ -737,9 +737,21 @@ char *GTIFGetOGISDefn( GTIF *hGTIF, GTIFDefn * psDefn >> ) break; >> >> case CT_Mercator: >> - oSRS.SetMercator( adfParm[0], adfParm[1], >> - adfParm[4], >> - adfParm[5], adfParm[6] ); >> + /* 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] ); >> >> if (psDefn->Projection == 1024 || psDefn->Projection == 9841) >> // override hack for google mercator. { >> @@ -1510,14 +1522,18 @@ int GTIFSetFromOGISDefn( GTIF * psGTIF, const char >> *pszOGCWKT ) GTIFKeySet(psGTIF, ProjNatOriginLongGeoKey, TYPE_DOUBLE, 1, >> poSRS->GetNormProjParm( SRS_PP_CENTRAL_MERIDIAN, 0.0 ) ); >> >> - GTIFKeySet(psGTIF, ProjScaleAtNatOriginGeoKey, TYPE_DOUBLE, 1, >> - poSRS->GetNormProjParm( SRS_PP_SCALE_FACTOR, 1.0 ) ); >> - >> GTIFKeySet(psGTIF, ProjFalseEastingGeoKey, TYPE_DOUBLE, 1, >> poSRS->GetProjParm( SRS_PP_FALSE_EASTING, 0.0 ) ); >> >> GTIFKeySet(psGTIF, ProjFalseNorthingGeoKey, TYPE_DOUBLE, 1, >> poSRS->GetProjParm( SRS_PP_FALSE_NORTHING, 0.0 ) ); >> + >> + if( EQUAL(pszProjection,SRS_PT_MERCATOR_2SP) ) >> + GTIFKeySet(psGTIF, ProjStdParallel1GeoKey, TYPE_DOUBLE, 1, >> + poSRS->GetNormProjParm( SRS_PP_STANDARD_PARALLEL_1, >> 0.0 ) ); + else >> + GTIFKeySet(psGTIF, ProjScaleAtNatOriginGeoKey, TYPE_DOUBLE, 1, >> + poSRS->GetNormProjParm( SRS_PP_SCALE_FACTOR, 1.0 ) >> ); } >> >> else if( EQUAL(pszProjection,SRS_PT_OBLIQUE_STEREOGRAPHIC) ) >> diff --git a/gdal/frmts/gtiff/libgeotiff/geo_normalize.c >> b/gdal/frmts/gtiff/libgeotiff/geo_normalize.c index 4761ea0..75705ee >> 100644 >> --- a/gdal/frmts/gtiff/libgeotiff/geo_normalize.c >> +++ b/gdal/frmts/gtiff/libgeotiff/geo_normalize.c >> @@ -1495,8 +1495,65 @@ static void GTIFFetchProjParms( GTIF * psGTIF, >> GTIFDefn * psDefn ) break; >> >> /* -------------------------------------------------------------------- */ >> - case CT_LambertConfConic_1SP: >> case CT_Mercator: >> +/* -------------------------------------------------------------------- */ >> + if( GTIFKeyGet(psGTIF, ProjNatOriginLongGeoKey, >> + &dfNatOriginLong, 0, 1 ) == 0 >> + && GTIFKeyGet(psGTIF, ProjFalseOriginLongGeoKey, >> + &dfNatOriginLong, 0, 1 ) == 0 >> + && GTIFKeyGet(psGTIF, ProjCenterLongGeoKey, >> + &dfNatOriginLong, 0, 1 ) == 0 ) >> + dfNatOriginLong = 0.0; >> + >> + if( GTIFKeyGet(psGTIF, ProjNatOriginLatGeoKey, >> + &dfNatOriginLat, 0, 1 ) == 0 >> + && GTIFKeyGet(psGTIF, ProjFalseOriginLatGeoKey, >> + &dfNatOriginLat, 0, 1 ) == 0 >> + && GTIFKeyGet(psGTIF, ProjCenterLatGeoKey, >> + &dfNatOriginLat, 0, 1 ) == 0 ) >> + dfNatOriginLat = 0.0; >> + >> + >> + if( GTIFKeyGet(psGTIF, ProjStdParallel1GeoKey, >> + &dfStdParallel1, 0, 1 ) == 0 ) >> + dfStdParallel1 = NAN; >> + >> + if( GTIFKeyGet(psGTIF, ProjScaleAtNatOriginGeoKey, >> + &dfNatOriginScale, 0, 1 ) == 0 ) >> + { >> + /* Default only if dfStdParallel1 isn't defined */ >> + if( isnan(dfStdParallel1) ) >> + dfNatOriginScale = 1.0; >> + else >> + dfNatOriginScale = NAN; >> + } >> + >> + /* notdef: should transform to decimal degrees at this point */ >> + >> + psDefn->ProjParm[0] = dfNatOriginLat; >> + psDefn->ProjParmId[0] = ProjNatOriginLatGeoKey; >> + psDefn->ProjParm[1] = dfNatOriginLong; >> + psDefn->ProjParmId[1] = ProjNatOriginLongGeoKey; >> + if( !isnan(dfStdParallel1) ) >> + { >> + psDefn->ProjParm[2] = dfStdParallel1; >> + psDefn->ProjParmId[2] = ProjStdParallel1GeoKey; >> + } >> + if( !isnan(dfNatOriginScale) ) >> + { >> + psDefn->ProjParm[4] = dfNatOriginScale; >> + psDefn->ProjParmId[4] = ProjScaleAtNatOriginGeoKey; >> + } >> + psDefn->ProjParm[5] = dfFalseEasting; >> + psDefn->ProjParmId[5] = ProjFalseEastingGeoKey; >> + psDefn->ProjParm[6] = dfFalseNorthing; >> + psDefn->ProjParmId[6] = ProjFalseNorthingGeoKey; >> + >> + psDefn->nParms = 7; >> + break; >> + >> +/* -------------------------------------------------------------------- */ >> + case CT_LambertConfConic_1SP: >> case CT_ObliqueStereographic: >> case CT_TransverseMercator: >> case CT_TransvMercator_SouthOriented: > > -- > Geospatial professional services > http://even.rouault.free.fr/services.html _______________________________________________ gdal-dev mailing list [email protected] http://lists.osgeo.org/mailman/listinfo/gdal-dev
