[gdal-dev] [PATCH] Support Mercator_2SP in GeoTIFF
Trent Piepho
tpiepho at gmail.com
Mon Nov 4 15:04:55 PST 2013
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
<even.rouault at mines-paris.org> 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
More information about the gdal-dev
mailing list