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

[Greg Troxel]

"Byeong Kim" <Byeong_Kim at dnr.state.ga.us> writes:

> 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.

>From your email address and -84 33, I'm going to assume that you are
talking about locations in or near the state of Georgia in the US.

It seems highly unlikely that anyone near Georgia (who is not
well-versed in datums and ellipsoids) would be able to come up with
numbers based on a perfect sphere for locations near Georgia.  For you
to do it, you had to use proj!  People that don't understand the
issues almost always either read numbers off a GPS receiver or look at
maps.  USGS topo maps are typically in NAD27 (with the Clarke 1866
ellipsoid).  Around Massachusetts the displacements from NAD83 are
typically a bit under 50m.  People without datum clue usually have not
changed receivers from WGS84, and if so it's probably NAD27 to match
paper maps.  At least in the age of GPS when you get lat/lon without a
datum it's usually in WGS84.

> [lots of description of transforms omitted]

I think I understand what you are trying to do, but I think you are
going about it wrong.  I'm also not perfectly clear on what happened.

There are basically two issues, datum and ellipsoid.  Datum includes
ellipsoid, of course, but you could have two datums that use the same
ellipsoid but have different coordinates for a given point because the
origin is different.

Then, when you transform from geodetic (lat/long) to projected, you
are not changing the datum, but just projecting into a new coordinate
system.  While the math is hairy, from a geodesy point of view this is
far simpler than datum shifts.

If you take a given coordinate pair, and project using different
ellipsoids, you will get very different values.  This is because the
length along the ellipsoid from the equator to your latitude is very
different.  However, if you look at the lat/long values for the same
oint in say NAD27 and NAD83/WGS84 (which are essentially the same if
you're not worrying about the last ~1m), and compute the different
between them *assuming they are in the same datum*, you'll send up
with ~50m.  This is because the datums are not aligned.  But if you
then transform the NAD27 values and the NAD83 values into UTM, where
each UTM conversion uses the datum's ellipsoids (which is what is
normally done), you'll get values that are much farther off, because
of the ellipsoid.


So basically you should be taking your lat/lon values, converting them
all to the same datum, and *only then* projecting them all using the
*same ellipsoid*.  I would expect you to be using your state plane
coordinate system, but that's transverse mercator in two zones for
both nad27 and nad83 according to the files that come with proj.

If you have in the past taken values and projected them using
different ellipsoids, and then tried to relate them to each other,
that's just plain broken and you really should reprocess your original
data.

If you've merely taken values in NAD27 and NAD83 without worrying
about which was which, and projected them into LCC with the Clarke
1866 ellipsoid (or any other *single* ellipsoid, and you want to know
how much error there is, you should

1) decide on the "right" datum in which you will compute errors.  Call
   this D_R.

2) Make a list of 'wrong" datums for which your input might be in.
   Call one D_W, and repeat for all.

2) take a lot of input points, and for each

  a) project to LCC with the D_R's ellipsoid.
  b) assume the point is in D_W, and convert to lat/lon in D_R.  Then
  project to LCC with D_R's ellipsoid.

Then, compute cartesian distance from a-b output, and look at those
over the input space.  You should see a fairly smooth change in error.

As to what you did

> 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

So this is hugely different, about 0.8 degreees of latitude which is
~50 miles.  It surprises me that the differences are so large, but no
one uses a spherical earth model -- because it is such a poor
approximation to reality.

In doing that step, you've conflated two issues.  One is the ellipsoid
used for the geodetic projection (from XYZ in ECEF), and the other is
datum shifts.  But the datum shift values are small, and they are even
smaller in North America.

> 2. Then I performed LCCprojection for both pairs with proj.
> 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 would expect more like 89 km in northing given 0.8 degrees, but I
suspect the datum/ellipsoid shift errors and the project erors from
using differnet ellipsoids may largely cancel, but this is
sufficiently far off from what one ought to be doing that there may be
something else going on.

> I am wondering if any of my calculation is wrong or is there any
> chance proj4 may have a bug in cs2cs utility.

None of this is evidence of problems with proj4.