Meters to Feet for false east/northing being done twice

[larrybiehl]

I have a geotiff file which is in NAD83 / Indiana East (ft US), EPSG Code 2965, where the false eastings and northings are being return as 1076386.7350175 and 2690966.839184167. This is 3.280833 times larger than they should be. The information from GDALGetProjectionRef is: PROJCS["NAD83 / Indiana East(ftUS)",GEOGCS["NAD83",DATUM["North_American_Datum_1983",SPHEROID["GRS190",6378137,298.2572221010002,AUTHORITY["EPSG","7019"]],AUTHORITY["EPSG","6269"],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433],AUTHORITY["EPSG","4269"]],PROJECTIONTransverse_Mercator,PARAMETER["latitude_of_origin",37.5],PARAMETER["central_meridian",-85.66666666666667],PARAMETER["scale_factor",0.999966667],PARAMETER["false_easting",1076386.7350175],PARAMETER["false_northing",2690966.839184167],UNIT["us_survey_feet",0.3048006096012192],AUTHORITY["EPSG","2965"]] .

I have traced this to:

int unitCode = 0;

GTIFKeyGet(hGTIF, ProjLinearUnitsGeoKey, &unitCode, 0, 1 ); if(unitCode != KvUserDefined) {

adfParm[5] /= psDefn->UOMLengthInMeters; adfParm[6] /= psDefn->UOMLengthInMeters;

}

/* -------------------------------------------------------------------- */ /* Translation the fundamental projection. */ /* -------------------------------------------------------------------- */

switch( psDefn->CTProjection ) {

case CT_TransverseMercator:

oSRS.SetTM( adfParm[0], adfParm[1],

adfParm[4], adfParm[5], adfParm[6] );

break;

in GTIFGetOGISDefn in the gt_wkt_srs.cpp file. adfParm[5] and adfParm[6] is converted from meters to feet here. They are also converted again in oSRS.SetTM at:

SetProjection( SRS_PT_TRANSVERSE_MERCATOR );

SetNormProjParm( SRS_PP_LATITUDE_OF_ORIGIN, dfCenterLat ); SetNormProjParm( SRS_PP_CENTRAL_MERIDIAN, dfCenterLong ); SetNormProjParm( SRS_PP_SCALE_FACTOR, dfScale ); SetNormProjParm( SRS_PP_FALSE_EASTING, dfFalseEasting ); SetNormProjParm( SRS_PP_FALSE_NORTHING, dfFalseNorthing );

---

It appears to be related to the citation in the file. This file has the following citation:

IMAGINE GeoTIFF Support ERDAS Desktop 2011 Version 11.0.3 11.0.0.896 Projection Name = State Plane Units = us_survey_feet GeoTIFF Units = us_survey_feet|

If this citation is not included, the calculations are correct. There is only one conversion. There is one in the code given above. The following code which gets called by SetNormProjParm comes up with a conversion factor of 1 for this case:

void OGRSpatialReference::GetNormInfo(void) const

{

if( bNormInfoSet )

return;

/* -------------------------------------------------------------------- */ /* Initialize values. */ /* -------------------------------------------------------------------- */

OGRSpatialReference *poThis = (OGRSpatialReference *) this;

poThis->bNormInfoSet = TRUE;

poThis->dfFromGreenwich = GetPrimeMeridian(NULL); poThis->dfToMeter = GetLinearUnits(NULL); poThis->dfToDegrees = GetAngularUnits(NULL) / CPLAtof(SRS_UA_DEGREE_CONV); if( fabs(poThis->dfToDegrees-1.0) < 0.000000001 )

poThis->dfToDegrees = 1.0;

}

The info from GDALGetProjectionRef which is correct is:

PROJCS["NAD83 / Indiana East (ftUS)",GEOGCS["NAD83",DATUM["North_American_Datum_1983",SPHEROID["GRS 1980",6378137,298.2572221010002,AUTHORITY["EPSG","7019"]],AUTHORITY["EPSG"," 6269"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433],AUTHORITY["EPSG" ,"4269"]],PROJECTIONTransverse_Mercator,PARAMETER["latitude_of_origin",37.5],PARAMETER["central_meridian",-85.66666666666667],PARAMETER["scale_factor",0.999966667],PARAMETER["false_easting",328083.333],PARAMETER["false_northing",820208.333],UNIT["US survey foot",0.304 8006096012192,AUTHORITY["EPSG","9003"]],AUTHORITY["EPSG","2965"]]

I have not figured how to correct for this.

I have attached two files. One which illustrates the problem and one which does not.

[larrybiehl]

Contains samples of file which is okay one that illustrates the problem.

[larrybiehl]

I have found what may be a solution to this, but do not know if it is good for all situations. linearUnitIsSet seems to be set by info in citations. A fix which handles my situation is:

if(unitCode != KvUserDefined && !linearUnitIsSet)

{

adfParm[5] /= psDefn->UOMLengthInMeters; adfParm[6] /= psDefn->UOMLengthInMeters;

}

instead of

if(unitCode != KvUserDefined)

{

adfParm[5] /= psDefn->UOMLengthInMeters; adfParm[6] /= psDefn->UOMLengthInMeters;

}

[Jukka Rahkonen]

I made a test with GDAL 2.0 and it suffers still from the same issue. When I compared the tags of the GeoTIFFs with listgeo the only difference is that the one which fools GDAL to get the "false_easting" and "false_northing" values wrong has this tag:

      GTCitationGeoKey (Ascii,157): "IMAGINE GeoTIFF Support.\nERDAS Desktop 201
1 Version 11.0.3 11.0.0.896.\nProjection Name = State Plane.\nUnits = us_survey_
feet.\nGeoTIFF Units = us_survey_feet"

[Even Rouault]

No longer reproducible with GDAL 2.2