Skip to content

Commit

Permalink
GTIFGetDefn(): add missing normalization of angular units to degree
Browse files Browse the repository at this point in the history 
Contrary to what the documentation of GTIFDefn::ProjParm[] mentionned,
we failed to normalize angular measures to degrees when reading them
from projection parameters (ProjXXXXXGeoKey's) when ProjCoordTransGeoKey
was present (but we did normalize them when reading them from the PROJ
database when there were was only a EPSG PCS code)

Relates to OSGeo/gdal#10154 and OSGeo/gdal#10158
  • Loading branch information
@rouault
rouault committed Jun 7, 2024
1 parent 9053d20 commit 46244ce
Show file tree
Hide file tree
Showing 7 changed files with 157 additions and 8 deletions.
33 changes: 29 additions & 4 deletions libgeotiff/geo_normalize.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2241,15 +2241,16 @@ static void GTIFFetchProjParms( GTIF * psGTIF, GTIFDefn * psDefn )
break;
}

for( int iParam = 0; iParam < psDefn->nParms; iParam++ )
{
switch( psDefn->ProjParmId[iParam] )
{

/* -------------------------------------------------------------------- */
/* Normalize any linear parameters into meters. In GeoTIFF */
/* the linear projection parameter tags are normally in the */
/* units of the coordinate system described. */
/* -------------------------------------------------------------------- */
for( int iParam = 0; iParam < psDefn->nParms; iParam++ )
{
switch( psDefn->ProjParmId[iParam] )
{
case ProjFalseEastingGeoKey:
case ProjFalseNorthingGeoKey:
case ProjFalseOriginEastingGeoKey:
Expand All @@ -2263,6 +2264,30 @@ static void GTIFFetchProjParms( GTIF * psGTIF, GTIFDefn * psDefn )
}
break;

/* -------------------------------------------------------------------- */
/* Normalize any angular parameters into degrees. In GeoTIFF */
/* the angular projection parameter tags are normally in the */
/* units of GeogAngularUnit. Note: this conversion is only done */
/* since libgeotiff 1.7.4 */
/* -------------------------------------------------------------------- */

case ProjStdParallel1GeoKey:
case ProjStdParallel2GeoKey:
case ProjNatOriginLongGeoKey:
case ProjNatOriginLatGeoKey:
case ProjFalseOriginLongGeoKey:
case ProjFalseOriginLatGeoKey:
case ProjCenterLongGeoKey:
case ProjCenterLatGeoKey:
case ProjStraightVertPoleLongGeoKey:
case ProjRectifiedGridAngleGeoKey:
if( psDefn->UOMAngleInDegrees != 0
&& psDefn->UOMAngleInDegrees != 1.0 )
{
psDefn->ProjParm[iParam] *= psDefn->UOMAngleInDegrees;
}
break;

default:
break;
}
Expand Down
12 changes: 9 additions & 3 deletions libgeotiff/geo_normalize.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -121,9 +121,15 @@ typedef struct {
int nParms;

/** Projection parameter value. The identify of this parameter
is established from the corresponding entry in ProjParmId. The
value will be measured in meters, or decimal degrees if it is a
linear or angular measure. */
is established from the corresponding entry in ProjParmId.
In GeoTIFF keys, the values of the projection parameters are expressed
in the units of ProjLinearUnitsGeoKey (for linear measures) or
GeogAngularUnitsGeoKey (for angular measures).
However, the value returned in ProjParam[] will be normalized to meters
or decimal degrees.
Note: until libgeotiff 1.7.3, the conversion to degrees for angular
measures was *not* done when ProjCoordTransGeoKey is present.
*/
double ProjParm[MAX_GTIF_PROJPARMS];

/** Projection parameter identifier. For example ProjFalseEastingGeoKey.
Expand Down
4 changes: 3 additions & 1 deletion libgeotiff/test/Makefile.am
View file
Original file line number Diff line number Diff line change
Expand Up @@ -55,7 +55,9 @@ EXTRA_DIST = testlistgeo testlistgeo_out.dist \
data/lambert_cylindrical_equal_area.tif \
data/ProjectedCSTypeGeoKey_4087_equidistant_cylindrical.tif \
data/equidistant_cylindrical.tif \
data/pixel_is_point_wgs84.tif
data/pixel_is_point_wgs84.tif \
data/epsg_27563_only_pcs_code.tif \
data/epsg_27563_allgeokeys.tif

check-local:
$(TESTLISTGEO) $(LISTGEOEXE)
Expand Down
Binary file added libgeotiff/test/data/epsg_27563_allgeokeys.tif
View file
Binary file not shown.
Binary file added libgeotiff/test/data/epsg_27563_only_pcs_code.tif
View file
Binary file not shown.
8 changes: 8 additions & 0 deletions libgeotiff/test/testlistgeo
View file
Original file line number Diff line number Diff line change
Expand Up @@ -244,6 +244,14 @@ mv testlistgeo_out.dist.normalized.tmp testlistgeo_out.dist.normalized
sed "s/ProjCenterLongGeoKey: 46.437229 ( 46d26'14.02\"E)/ProjCenterLongGeoKey: 46.437229 ( 46d26'14.03\"E)/" < testlistgeo_out.dist.normalized > testlistgeo_out.dist.normalized.tmp
mv testlistgeo_out.dist.normalized.tmp testlistgeo_out.dist.normalized

echo "Testing listgeo epsg_27563_only_pcs_code.tif" >> ${OUT}
$EXE ${DATA_DIR}/epsg_27563_only_pcs_code.tif >>${OUT}
echo "" >>${OUT}

echo "Testing listgeo epsg_27563_allgeokeys.tif" >> ${OUT}
$EXE ${DATA_DIR}/epsg_27563_allgeokeys.tif >>${OUT}
echo "" >>${OUT}

# do 'diff' with distribution results
# after cleaning for avoid spurios result due
# to different build dir
Expand Down
108 changes: 108 additions & 0 deletions libgeotiff/test/testlistgeo_out.dist
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2005,3 +2005,111 @@ Upper Right ( 2d 0' 0.00"E, 54d 0' 0.00"N)
Lower Right ( 2d 0' 0.00"E, 50d 0' 0.00"N)
Center ( 0d 0' 0.00"E, 52d 0' 0.00"N)

Testing listgeo epsg_27563_only_pcs_code.tif
Geotiff_Information:
Version: 1
Key_Revision: 1.0
Tagged_Information:
ModelTiepointTag (2,3):
0 0 0
827294.141443773 523985.644431661 0
ModelPixelScaleTag (1,3):
0.500636999995913 0.500637000001734 0
End_Of_Tags.
Keyed_Information:
GTModelTypeGeoKey (Short,1): ModelTypeProjected
GTRasterTypeGeoKey (Short,1): RasterPixelIsArea
GTCitationGeoKey (Ascii,33): "NTF (Paris) / Lambert Sud France"
GeogCitationGeoKey (Ascii,12): "NTF (Paris)"
GeogAngularUnitsGeoKey (Short,1): Angular_Grad
ProjectedCSTypeGeoKey (Short,1): Code-27563 (NTF (Paris) / Lambert Sud France)
ProjLinearUnitsGeoKey (Short,1): Linear_Meter
End_Of_Keys.
End_Of_Geotiff.

PCS = 27563 (NTF (Paris) / Lambert Sud France)
Projection = 18093 (Lambert Sud France)
Projection Method: CT_LambertConfConic_1SP
ProjNatOriginLatGeoKey: 44.100000 ( 44d 6' 0.00"N)
ProjNatOriginLongGeoKey: 0.000000 ( 0d 0' 0.00"E)
ProjScaleAtNatOriginGeoKey: 0.999877
ProjFalseEastingGeoKey: 600000.000000 m
ProjFalseNorthingGeoKey: 200000.000000 m
GCS: 4807/NTF (Paris)
Datum: 6807/Nouvelle Triangulation Francaise (Paris)
Ellipsoid: 7011/Clarke 1880 (IGN) (6378249.20,6356515.00)
Prime Meridian: 8903/Paris (2.337229/ 2d20'14.03"E)
Projection Linear Units: 9001/metre (1.000000m)

Corner Coordinates:
Upper Left ( 827294.141, 523985.644) ( 2d59' 3.48"E, 46d58'37.71"N)
Lower Left ( 827294.141, 523980.638) ( 2d59' 3.47"E, 46d58'37.54"N)
Upper Right ( 827299.148, 523985.644) ( 2d59' 3.71"E, 46d58'37.70"N)
Lower Right ( 827299.148, 523980.638) ( 2d59' 3.71"E, 46d58'37.54"N)
Center ( 827296.645, 523983.141) ( 2d59' 3.59"E, 46d58'37.62"N)

Testing listgeo epsg_27563_allgeokeys.tif
Geotiff_Information:
Version: 1
Key_Revision: 1.0
Tagged_Information:
ModelTiepointTag (2,3):
0 0 0
827294.141443773 523985.644431661 0
ModelPixelScaleTag (1,3):
0.500636999995913 0.500637000001734 0
End_Of_Tags.
Keyed_Information:
GTModelTypeGeoKey (Short,1): ModelTypeProjected
GTRasterTypeGeoKey (Short,1): RasterPixelIsArea
GTCitationGeoKey (Ascii,33): "NTF (Paris) / Lambert Sud France"
GeographicTypeGeoKey (Short,1): GCS_NTF_Paris
GeogCitationGeoKey (Ascii,12): "NTF (Paris)"
GeogGeodeticDatumGeoKey (Short,1): Code-6807 (Nouvelle Triangulation Francaise (Paris))
GeogPrimeMeridianGeoKey (Short,1): PM_Paris
GeogAngularUnitsGeoKey (Short,1): Angular_Grad
GeogAngularUnitSizeGeoKey (Double,1): 0.0157079632679489
GeogEllipsoidGeoKey (Short,1): Ellipse_Clarke_1880_IGN
GeogSemiMajorAxisGeoKey (Double,1): 6378249.2
GeogInvFlatteningGeoKey (Double,1): 293.466021293627
GeogPrimeMeridianLongGeoKey (Double,1): 2.5969213
GeogTOWGS84GeoKey (Double,7): -168 -60 320
0 0 0
0
ProjectedCSTypeGeoKey (Short,1): Code-27563 (NTF (Paris) / Lambert Sud France)
PCSCitationGeoKey (Ascii,33): "NTF (Paris) / Lambert Sud France"
ProjectionGeoKey (Short,1): Code-18093 (Lambert Sud France)
ProjCoordTransGeoKey (Short,1): CT_LambertConfConic_1SP
ProjLinearUnitsGeoKey (Short,1): Linear_Meter
ProjLinearUnitSizeGeoKey (Double,1): 1
ProjNatOriginLongGeoKey (Double,1): 0
ProjNatOriginLatGeoKey (Double,1): 49
ProjFalseEastingGeoKey (Double,1): 600000
ProjFalseNorthingGeoKey (Double,1): 200000
ProjScaleAtNatOriginGeoKey (Double,1): 0.999877499
VerticalUnitsGeoKey (Short,1): Linear_Meter
End_Of_Keys.
End_Of_Geotiff.

PCS = 27563 (NTF (Paris) / Lambert Sud France)
Projection = 18093 (Lambert Sud France)
Projection Method: CT_LambertConfConic_1SP
ProjNatOriginLatGeoKey: 44.100000 ( 44d 6' 0.00"N)
ProjNatOriginLongGeoKey: 0.000000 ( 0d 0' 0.00"E)
ProjScaleAtNatOriginGeoKey: 0.999877
ProjFalseEastingGeoKey: 600000.000000 m
ProjFalseNorthingGeoKey: 200000.000000 m
GCS: 4807/NTF (Paris)
Datum: 6807/Nouvelle Triangulation Francaise (Paris)
Ellipsoid: 7011/Clarke 1880 (IGN) (6378249.20,6356515.00)
Prime Meridian: 8903/Paris (2.337229/ 2d20'14.03"E)
TOWGS84: -168,-60,320,0,0,0,0
Projection Linear Units: 9001/metre (1.000000m)

Corner Coordinates:
Upper Left ( 827294.141, 523985.644) ( 2d59' 3.48"E, 46d58'37.71"N)
Lower Left ( 827294.141, 523980.638) ( 2d59' 3.47"E, 46d58'37.54"N)
Upper Right ( 827299.148, 523985.644) ( 2d59' 3.71"E, 46d58'37.70"N)
Lower Right ( 827299.148, 523980.638) ( 2d59' 3.71"E, 46d58'37.54"N)
Center ( 827296.645, 523983.141) ( 2d59' 3.59"E, 46d58'37.62"N)

0 comments on commit 46244ce

Please sign in to comment.