[Proj] Verification of cs2cs results for geographical coordinate system transformation

[Byeong Kim]

Hi all,
 
This is a rephrased posting of what I posted on Nov 2. I am not sure how to remove the old posting. Sorry for my duplication.
*---------------------------------------------------------------------------------------------------------------------------------------------
 
I am trying to figure out the location difference by using a Lon/Lat value mistakenly. 
Suppose you have a collaborator who does not know anything about datum or geocoordinate transformation.  You asked him a lat/lon pair for one of you study subject without specifying the ellipsoid and somehow s/he did not ask you back about it either.
The number s/he gave you is 
-84 33. 0 (Lon Lat Height)
It turned out s/he assumed that this number is based on a perfect sphere (likely R=6370997.0m).
Later, however, you found out that it is based on Clarke1866. In ESRI's term it can be GCS_Assumed_Geographic_1 and this ellipsoid is frequently used for many US Census data. 
Since you use the lat/lon value to do Lambert Conformal Projection, now you will be curious about what the possible error in (X,Y) will be.
This what I am trying to figure out.
 
To answer the question, I wrote a BAT file (I am using a PROJ.4 utilities on XP after compiling it with MinGW and datum shift tables) to compute it. 
 
Steps are following: 
 
0. I created a text file containing a pair of lat/lon value. The name of file is MonLonLat that contains
-84 33. 0
 
1. I transformed the given lat/lon into a Lat/Lon on a Perfect Sphere with cs2cs
cs2cs +proj=longlat +ellps=clrk66 +datum=NAD27 +towgs84=-3,142,183,0,0,0,0 +to +proj=latlong +ellps=sphere +R=6370997 +towgs84=0,0,0 +no_def -f "%%f2" -w9 MonLonLat > MonLonLat.Clarke1866ToSphere
With this line, I got a file MonLonLat.Clarke1866ToSphere what contains
-84.0000002     32.8244522 834.0192882
 
2. Then I performed LCCprojection for both pairs with proj.
 
2a. Using a lat/lon from Clarke1866 ellipsoid
proj +proj=lcc +lat_1=60.0 +lat_2=30.0 +lat_0=40.0 +lon_0=-100.0 +x_0=0 +y_0=0 +ellps=sphere +R=6370997 -f "%%f2" -w9  MonLonLat
This gives me
1464716.5725192 -614514.0442942 0
 
2b. Using a lat/lon from a Perfect Sphere
proj +proj=lcc +lat_1=60.0 +lat_2=30.0 +lat_0=40.0 +lon_0=-100.0 +x_0=0 +y_0=0 +ellps=sphere +R=6370997 -f "%%f2" -w9 MonLonLat.Clarke1866ToSphere
This gives me
1468546.7535622 -633426.1117132 834.0192882
 
2a case represents the event that you use a lat/lon value based on Clarke1866 to do Lambert Conformal Projection as if the number is based on a Perfect Sphere. That is, this is the wrong case.
2b case represents the event that you correct a lat/lon value based on Clarke1866 to be a lat/lon on a Perfect Sphere before doing a projection. That is, it is the correct way I believe.
 
What surprised me is that the error (or difference) is very large for my application.  As you see, the displacement is about 4km in east-west and about 20 km in north-south.
 
I am wondering if any of my calculation is wrong or is there any chance proj4 may have a bug in cs2cs utility.


-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://lists.maptools.org/pipermail/proj/attachments/20061106/899938a0/attachment.html