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

[Martin Chapman]

Even,

Thanks for the quick response.  Yeah, I goofed, I meant to say units, guess 
I am just getting old.  Anyway, thanks for pointing me in the right 
direction!

Best regards,
Martin

-----Original Message-----
From: Even Rouault [mailto:even.rouault at mines-paris.org]
Sent: Thursday, July 10, 2014 2:35 PM
To: gdal-dev at lists.osgeo.org
Cc: Martin Chapman
Subject: Re: [gdal-dev] State Plane to WGS84 Transformation problem with 
false easting and northing values

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],AUTH
> ORIT
> Y["EPSG","4269"]],PROJECTION["Lambert_Conformal_Conic_2SP"],PARAMETER[
> "sta
> ndard_parallel_1",35.46666666666667],PARAMETER["standard_parallel_2",3
> 4.03
> 333333333333],PARAMETER["latitude_of_origin",33.5],PARAMETER["central_
> meri
> dian",-118],PARAMETER["false_easting",21527734.72222222],PARAMETER["fa
> lse_
> 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],AUTH
> ORIT
> Y["EPSG","4269"]],PROJECTION["Lambert_Conformal_Conic_2SP"],PARAMETER[
> "sta
> ndard_parallel_1",35.46666666666667],PARAMETER["standard_parallel_2",3
> 4.03
> 333333333333],PARAMETER["latitude_of_origin",33.5],PARAMETER["central_
> meri
> dian",-118],PARAMETER["false_easting",6561666.666666666],PARAMETER["fa
> lse_ northing",1640416.666666667],UNIT["US survey
> foot",0.3048006096012192,AUTHORITY["EPSG","9003"]],AUTHORITY["EPSG","2
> 6945
> "]]
>
>
>
> 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