[postgis-users] Interesting Datum Conversion Issue

Thu May 25 00:23:03 PDT 2006

we make datum shift with postgis with custom srids, like this... search in
the lis for this email and this subject: *Datum Shift Problems with PostGIS

Regards and contact me for any questions... the following problem was solved
and currently we make 3 an 7 parameters datum transformations

REGARDS*

Datum shift Problems with PostGIS

Dear List Users, In order to make datum change with proj from command line,
i have the following example that works fine.

Datum Change from Sad69/utm 19s to wgs84/utm 19s

# cs2cs +proj=utm +zone=19 +south +ellps=GRS67 +towgs84=-75,-1,-44 +unit=m
+to +proj=utm +zone=19 +south +datum=WGS84

In Order to make the same procces, but using Postgis we have added a record
en the spatial_ref_sys table, with a custom epsg code to get a personalized
transformation, and in en the proj4text field the following sentence (like
the parameters of cs2cs) has been added.

Record added to spatial_ref_sys Table

srid : 50000
auth_name : EPSG
auth_srid: 50000
srtext : Sad69 a wgs84
proj4text : +proj=utm +zone=19 +south +ellps=GRS67 +towgs84=75,1,44 +unit=m
+to +proj=utm +zone=19 +south +datum=WGS84 +no_defs

To run this custom datum transformation we use the transform() function
build in a query inserted as parameter in the pgsql2shp command.

pgsql2shp -h localhost -p 5432 -u pgsql -P pgsql -f shapefile_name.shp
db_name
"select transform(setsrid(the_geom,29179),50000) as the_geom from
nombre_tabla"

(we have our cartography with SRID=-1, reason why we use the setsrid
function)

we get a shapefile correctly transformed to wgs84/utm 19s Datum

Here is the first Question: Why we need to change the sign in the
transformation parameters, +towgs84=-75,-1,-44 in command line and
+towgs84=75,1,44 in the proj4text fiel of the spatial_ref_sys table ?????

we had many problems with this, we solve it, but we don't know why we need
to change it if postgis works with proj.

Now, we try to make the inverse process, in order to make a datum shift from
WGS84 utm 19s to sad69 utm 19s, and the proj documentation tells that the -I
parameter is the only that we need to make the inverse datum shift. With
cs2cs command the -I parameter works fine, but the inclusion of this
parameter in the proj4text field in the spatial_ref_sys table does not make
a correct transformation, giving us around 1000 meters of displacement.

Second Question: How can we make this inverse transformation with postgis
(wgs84 to psad56) ?

2006/5/15, James G Wilkinson <jgw at alpinegeophysics.com>:
>
> 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,
>
> _______________________________________________
> postgis-users mailing list
> postgis-users at postgis.refractions.net
> http://postgis.refractions.net/mailman/listinfo/postgis-users
>

-- 
René F. Viáncos S.
Director de Geomática y TIC
Unidad de Ambiente, Región Territorio y Espacio, UNARTE
Facultad de Ciencias
Universidad de Chile
Tel (56-2) 632 62 09
Cel (56 9) 933 72 66
rviancos at uchile.cl
rviancos at gmail.com
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/postgis-users/attachments/20060525/10eebf02/attachment.html>