[gdal-dev] State Plane to WGS84 Transformation problem with false easting and northing values

[Even Rouault]

Martin,

Yes this issue has been reported several times (last time was 
http://lists.osgeo.org/pipermail/gdal-dev/2014-May/039251.html) and is a GDAL 
bug that the user code shouldn't compensate for. I think it is 
https://trac.osgeo.org/gdal/ticket/3901 / 
https://trac.osgeo.org/gdal/ticket/4954

The 0.3048 ratio you see is the length of a foot in meters (not sure why you 
mention "unit radians" here). Basically the conversion from foot to meters is 
wrongly applied twice in the geotiff SRS decoding.

21527734.72222222 * 0.3048 * 0.3048 = 2000000

Something I've looked at a bit, but did not pursue because of the lack of 
incentive to do so.

Even

Le jeudi 10 juillet 2014 22:12:04, Martin Chapman a écrit :
> Frank, Even or whoever,
> 
> 
> 
> I have a source geotiff that has the following projection:
> 
> 
> 
> PROJCS["State
> Plane",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],AUTHORIT
> Y["EPSG","4269"]],PROJECTION["Lambert_Conformal_Conic_2SP"],PARAMETER["sta
> ndard_parallel_1",35.46666666666667],PARAMETER["standard_parallel_2",34.03
> 333333333333],PARAMETER["latitude_of_origin",33.5],PARAMETER["central_meri
> dian",-118],PARAMETER["false_easting",21527734.72222222],PARAMETER["false_
> northing",5381933.680555555],UNIT["feet",0.3048006096012192],AUTHORITY["EP
> SG","26945"]]
> 
> 
> 
> I want to convert the bounds to WGS84.  When I setup my coordinate
> transform object the OGRProj4CT::InitializeNoLock() function in ogrct.cpp
> calls poSRSSource->exportToProj4() and returns the following:
> 
> 
> 
> +proj=lcc +lat_1=35.46666666666667 +lat_2=34.03333333333333 +lat_0=33.5
> +lon_0=-118 +x_0=6561666.666666666 +y_0=1640416.666666667 +datum=NAD83
> +units=us-ft +no_defs
> 
> 
> 
> Then, when calling TransformEx() on the bounds the resulting coordinates
> are wrong.  It appears that somewhere in the transformation process the
> UNIT RADIANS value is not being factored correctly with the false eating
> and northing values.
> 
> 
> 
> If I translate the source image into a new geotiff with the exact same
> parameters as the source geotiff then the new projection is as follows:
> 
> 
> 
> PROJCS["NAD83 / California zone
> 5",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],AUTHORIT
> Y["EPSG","4269"]],PROJECTION["Lambert_Conformal_Conic_2SP"],PARAMETER["sta
> ndard_parallel_1",35.46666666666667],PARAMETER["standard_parallel_2",34.03
> 333333333333],PARAMETER["latitude_of_origin",33.5],PARAMETER["central_meri
> dian",-118],PARAMETER["false_easting",6561666.666666666],PARAMETER["false_
> northing",1640416.666666667],UNIT["US survey
> foot",0.3048006096012192,AUTHORITY["EPSG","9003"]],AUTHORITY["EPSG","26945
> "]]
> 
> 
> 
> And the exportToProj4() returns this:
> 
> 
> 
> +proj=lcc +lat_1=35.46666666666667 +lat_2=34.03333333333333 +lat_0=33.5
> +lon_0=-118 +x_0=2000000 +y_0=500000.0000000001 +datum=NAD83 +units=us-ft
> +no_defs
> 
> 
> 
> Which produces a correct result when using TransformEx().
> 
> 
> 
> When I get the SRS_PP_FALSE_EASTING and SRS_PP_FALSE_NORTHING from both
> projections using the OGRSpatialReference::GetProjParm() function I notice
> the difference between the two is a factor of the GEOGCS unit radians:
> 
> 
> 
> SOURCE IMAGE:
> 
> false_easting:       21527734.72222
> 
> false_northing:      5381933.680556
> 
> 
> 
> TRANSLATED IMAGE:
> 
> false_easting:       6561666.666667
> 
> false_northing:      1640416.666667
> 
> 
> 
> 
> 
> The translated offsets divided by the source offsets seem to be a factor
> of the unit radians: 0.3048
> 
> 
> 
> I think the reason that the translate to geotiff corrects the problem is
> that the geotiff keys for the projection are calculated in a different
> manner than when using the OGRCoordinateTransform object.
> 
> 
> 
> Any help on where to apply the fix in the GDAL code if you think one is
> required would be much appreciated.  If you think I need to accommodate
> this in my code any help would be greatly appreciated.
> 
> 
> 
> Best regards,
> 
> Martin

-- 
Geospatial professional services
http://even.rouault.free.fr/services.html