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.roua...@mines-paris.org] Sent: Thursday, July 10, 2014 2:35 PM To: gdal-dev@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 _______________________________________________ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev