[gdal-dev] [PATCH] Support Mercator_2SP in GeoTIFF
Even Rouault
even.rouault at mines-paris.org
Mon Nov 4 11:12:56 PST 2013
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