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

[Martin Chapman]

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

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20140710/be7e2798/attachment.html>