[postgis-users] Interesting Datum Conversion Issue

[James G Wilkinson]

Mr. Warmerdam,

Many thanks for the follow up.  It addresses many of my concerns.

On another front, is it possible that anyone out there has developed the 
particular
files to do this datum translation (i.e., geographic NAD83 to LCC 
datumless--sphere;
LCC NAD83 to geographic datumless--sphere)?  I suspect not, but it never 
hurts to ask.

Can anyone point me to some good literature (paper or web-based with 
web-based
preferred) on how to develop three-point and seven-point datum 
translations that I
can then plug into PROJ4?  Anything that I may develop I will be happy 
to share
with the community.

Regards,

Jim Wilkinson

Frank Warmerdam wrote:

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