[Gdal-dev] Re: a bug in OGRCoordinateTransformation

[Frank Warmerdam]

On 6/24/05, Yang, Zhiqiang <zhiqiang.yang at oregonstate.edu> wrote:
>  
> Hi, 
>   
> I am trying to do some transformation with OGRCoordinateTransformation, but
> the results do not seem correct when compared with result from Erdas Imagine
> 8.7. Here is the details. 
>   
> Source projection is 
> 'PROJCS["UTM Zone 10, Northern
> Hemisphere",GEOGCS["WGS_1984",DATUM["WGS_1984",SPHEROID["WGS
> 84",6378137,298.2572235630016],TOWGS84[0,0,0,0,0,0,0]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-123],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["meters",1]]'
>  
> Target projection is 
> `PROJCS["Lambert Azimuthal
> Equal-area",GEOGCS["Undefined",DATUM["Undefined",SPHEROID["Sphere
> of Radius
> 6370997m",6370997,0]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433]],PROJECTION["Lambert_Azimuthal_Equal_Area"],PARAMETER["latitude_of_center",45],PARAMETER["longitude_of_center",-100],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["meters",1]]'
>  
> Coordinate to test: (360055.5, 5132807.5) 
> Results from OGRCoordinateTransformation: (-1874278.666369, 420070.981388) 
> Results from Imagine: (-1870537.414966, 439891.156463) 
>   
> The difference between the two results is so big that
> OGRCoordinateTransformation more likely to have bug in it. 

Yang,

The src coordinate system is internally converted to the PROJ.4 
representation:

   +proj=utm +zone=10 +ellps=WGS84 +datum=WGS84 +units=m +no_defs

The target coordinate system is converted to:

   +proj=laea +lat_0=45 +lon_0=-100 +x_0=0 +y_0=0 +a=6370997
+b=6370997 +units=m +no_defs

The most likely complication here is that the target coordinate system has
a very different earth model.  Currently, by default, OGR (and PROJ.4)
will attempt to translate between earth models accurately.    That is,
the lat/long coordinate in on the WGS84 ellipsoid is converted into an
x/y/z geocentric location, and then this is converted into a lat/long/height
on the spherical earth model. 

This results in a big latitude shift (and a substantial elevation change).
However, in my experience, when people work with spherical earth
models, they generally do not want the lat/long location converted via
geocentric location.  They would instead prefer the lat/long be used without
adjustment on the sphere. 

Well I thought I would be able to demonstrate that this is what was going on
in this case, but I don't seem to be able to get a matching result!  The
intermediate steps are:

1) Convert from source UTM coordinates to WGS84.

warmerda at gdal2200[27]% testepsg -t src.wkt WGS84 360055.5 5132807.5
(360055.500000,5132807.500000) -> (-124.818312,46.334384)

2) Convert this lat/long to output LAEA projected coordinates, just
taking the lat/long location directly on the sphere.

warmerda at gdal2200[28]% testepsg -t '+proj=latlong +a=6370997
+b=6370997 +no_defs' dst.wkt -124.818312 46.334384
(-124.818312,46.334384) -> (-1867682.518646,440660.120598)

But these results are not particularly close to the results from Imagine
either!  

I think to better establish why there is a difference you will need to 
do your imagine transformation in steps so we can see what step
differs.  

1) convert utm to wgs84 and compare. 
2) convert wgs84 to sphere and compare. 
3) convert spherical lat/long to LAEA and compare. 

PS. Are you really using the WKT directly with Imagine or 
converting it via some mechanism?  Perhaps this translation is
not working well? 

Best regards,
-- 
---------------------------------------+--------------------------------------
I set the clouds in motion - turn up   | Frank Warmerdam, warmerdam at pobox.com
light and sound - activate the windows | http://pobox.com/~warmerdam
and watch the world go round - Rush    | Geospatial Programmer for Rent