[postgis-users] Interesting Datum Conversion Issue

[Frank Warmerdam]

James G Wilkinson wrote:
> A situation has occurred that I am having difficulty getting my hands 
> around.
> I have a single point with the following attributes:
> 
> lon = -115.10101
> lat = 36.15526
> PROJ4 entry in SPATIAL_REF_SYS
> +proj=longlat +ellps=GRS80 +datum=NAD83 +no_defs
> 
> In ARC/Info, when I project this point to the following:
> 
> projection: LAMBERT
> datum: NONE
> units: meters
> spheroid: sphere
> 1st parallel: 33.0
> 2nd parallel: 45.0
> central meridian: -118.0
> latitude of origin: 37.0
> false easting: 0.0
> false northing: 0.0
> 
> I get back a point at  y-coord = -89372 m.
> 
> The above projection is entered into my SPATIAL_REF_SYS table
> with the following PROJ4 parameters (SRID=888889):
> +proj=lcc +lat_1=33.0 +lat_2=45.0 +lat_0=37.0 +lon_0=-118.0 +x_0=0 
> +y_0=0 +a=6370997.0 +b=6370997.0 +units=m
> 
> Now with PostgreSQL v.7.4.8, PostGIS v.0.9.1, and PROJ4 v.4.4.7,
> this same point using SELECT Y ( TRANSFORM (the_geom, 888889) ) FROM 
> spatial.pt_test
> returns -89074 (which is close enough for my work.....for now).
> 
> When I migrated to PostgreSQL v.8.1.3, PostGIS v.1.1.2, and PROJ4 v.4.4.9,
> my SELECT returns -109641 (!! and this is what baffles me !!).
> 
> Any help that folks can lend will be greatly appreciated.

Jim,

I believe the problem is that PROJ.4 is attempting to correct for the
difference in ellipsoid when it should really just treat your NAD83
lat/long location as a value on the defined sphere (without conversion).

What it is doing now is taking the location in NAD83, converting it to
a geocentric x/y/z relative to the center of the earth.  Then converting
that location onto the sphere.  Because of the substantial difference
between GRS80 ellipsoid and a sphere this produces a large apparent
shift in lat/long between the two earth models.

My suggestion is that if you want to treat your lat/long coordinate as
being on the sphere, then use the input coordinate system:
   +proj=latlong +a=6370997.0 +b=6370997.0

I will add that the PROJ.4 behavior, while seemingly "proper", has been
no end of problem and I will likely alter it in some future version.  It
seems the common practice in the GIS industry is to do what ESRI does.
If not provided with datum shift information to just ignore the difference
in ellipsoids.

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    | President OSGF, http://osgeo.org