[Geotiff] [gdal-dev] [PATCH] Support Mercator_2SP in GeoTIFF

Even Rouault even.rouault at mines-paris.org
Mon Nov 4 14:12:56 EST 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 Geotiff mailing list