[FWTools] projecting with ogr2ogr v. ArcInfo

[Frank Warmerdam]

On 12/15/05, Bryan Keith <bryan at geomega.com> wrote:
> Hello,
>
> I'm attempting to project data from geographic coordinates to a
> Minnesota county coordinate system (Hennepin County) using ogr2ogr and
> ArcInfo.  The results are slightly different, and I'm not sure why.  X
> coordinates are about the same (within 0.1 feet), but Y coordinates
> differ by ~3.1 feet.  This difference seems consistent though I haven't
> verified that thoroughly.  I include a small subset of the points in
> this e-mail.
>
> Input coordinates are long-lat:
>
> -93.17356183         45.04758353
> -93.36960993         45.08857070
> -93.27394885         44.97432735
>
> Output from ogr2ogr:
>
> 554221.66861940 193587.15836056
> 503544.69316568 208462.19508218
> 528309.75142366 166825.34464421
>
> Output from ArcInfo:
>
> 554221.6767         193584.0575
> 503544.6937         208459.0942
> 528309.7556         166822.2438

Bryan,

I believe that Arc/Info is just doing the reprojection
step on the custom ellipsoid, but not taking any action
to try and translate that to the NAD83 datum.

If I just do this with PROJ.4 directly, I get the same
results as Arc/Info.

% proj +proj=lcc +lat_1=44.88333333333333 +lat_2=45.13333333333333
+lat_0=44.79111111111111 +lon_0=-93.38333333333334 +x_0=152400.3048
+y_0=30480.06096 +a=6378418.941 +b=6357033.31   +units=us-ft +no_defs
-93.17356183         45.04758353
554221.68       193584.06

But if I do a translation from NAD83 to that coordinate system I get
the values you saw.

% cs2cs +proj=latlong +datum=NAD83 +to +proj=lcc
+lat_1=44.88333333333333 +lat_2=45.13333333333333
+lat_0=44.79111111111111 +lon_0=-93.38333333333334 +x_0=152400.3048
+y_0=30480.06096 +a=6378418.941 +b=6357033.31 +units=us-ft +no_defs
-93.17356183         45.04758353
554221.67       193587.16 -281.47

I believe this is because PROJ.4 is doing an ellipsoid to ellipsoid
conversion.

So, if you just want to do the reprojection with no datum shifting
you should specify "+proj=latlong +a=6378418.941 +b=6357033.3"
as your source coordinate system.

Note, this automatic ellipsoid to ellipsoid conversion (when there is
no explicit datum shift info) is something I am likely to remove from
PROJ.4 in the next revision since it is causing so much problem.  It
definately seems to *not* be what other software packages do.

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