[OSRS-PROJ] Datum Shifting
Frank Warmerdam
warmerda at home.com
Thu Jul 6 21:04:41 EDT 2000
Folks,
Well, I have been beavering away on datum shifting support for a couple
days, and have made substantial progress.
At the coordinate system definition level I have added:
+datum=name
The name is looked up in a list of datums, and used to imply an
ellipsoid, and a datum definition (+towgs84= or +nadgrids=)
+towgs84=xshift,yshift,zshift[,xrot,yrot,zrot,scale]
Defines a shift to WGS84 for a datum.
+nadgrids=file[,file...]
Gives a list of files to be used to perform grid shift to WGS84
for this datum (only in the output format nad2bin for now).
+proj=latlong
Defines a coordinate system as being geodetic, rather than projected.
The PJ data structure has been extended with the following fields:
int datum_type; /* PJD_UNKNOWN/3PARAM/7PARAM/GRIDSHIFT/WGS84 */
double datum_params[7];
int is_latlong; /* proj=latlong ... not really a projection at all */
Note that for the latlong coordinate system a forward, and inverse function
have been implemented that do nothing. Just mapping a location in angular
radians back onto itself.
The following prototypes have been added to projects.h. The point_offset
argument indicates the stepsize in doubles between entries in the x, y, and
z arrays. For data held in xyzxyz format the point_offset would be 3. For
data in separate x, y, z arrays the point_offset would be 1.
int pj_transform( PJ *src, PJ *dst, long point_count, int point_offset,
double *x, double *y, double *z );
int pj_datum_transform( PJ *src, PJ *dst, long point_count, int point_offset,
double *x, double *y, double *z );
int pj_geocentric_to_geodetic( double a, double ra,
long point_count, int point_offset,
double *x, double *y, double *z );
int pj_geodetic_to_geocentric( double a, double ra,
long point_count, int point_offset,
double *x, double *y, double *z );
int pj_compare_datums( PJ *srcdefn, PJ *dstdefn );
int pj_apply_gridshift( const char *, int,
long point_count, int point_offset,
double *x, double *y, double *z );
void pj_deallocate_grids();
The only aspect that I plan to document for public use is pj_transform()
which operates as described before. However, the geocentric/geodetic
conversion functions and some others might be useful to some applications
in isolation. Each of the functions that returns an int is returning an
error code (0 is success). I have also learn about the pj_errno flag which
is also set properly by these functions.
The pj_init() function works as it always did, but now calls pj_datum_set()
to parse datum parameters (+datum=, +towgs84= and +nadgrids=). The datum
parameters are all optional. The pj_inv() and pj_fwd() functions are unchanged.
The pj_transform() is the entry point to take advantage of datum shifting.
Note that pj_transform() does 3/7 parameter geocentric shifts as well as
supporting NAD27/83 grid shifting using the existing nad_cvt() API.
I have changed proj to support the -ld options to list datums, but have
otherwise left it unchanged. I have implemented a new cs2cs program with
works someone similarly to the proj program, but with less options (no binary,
Chebyshev polynomials, error estimates, or verbose out). For now I am using
the +to parameter as a separator between the from and to coordinate systems.
For example:
warmerda at cs46980-c[384]% cs2cs -ld
__datum_id__ __ellipse___ __definition/comments______________________________
WGS84 WGS84 towgs84=0,0,0
GGRS87 GRS80 towgs84=-199.87,74.79,246.62
Greek_Geodetic_Reference_System_1987
NAD83 GRS80 towgs84=0,0,0
North_American_Datum_1983
NAD27 clrk66 nadgrids=conus
North_American_Datum_1927
warmerda at cs46980-c[386]% cs2cs -v +proj=latlong +datum=NAD27 +to +proj=latlong
+datum=GGRS87
# ---- From Coordinate System ----
#Lat/long (Geodetic)
#
# +proj=latlong +datum=NAD27 +ellps=clrk66 +nadgrids=conus
# ---- To Coordinate System ----
#Lat/long (Geodetic)
#
# +proj=latlong +datum=GGRS87 +ellps=GRS80 +towgs84=-199.87,74.79,246.62
-117 33
116d59'54.943"W 32d59'53.873"N -154.525
warmerda at cs46980-c[390]% cs2cs -v +units=us-ft +init=nad27:406 +to +units=m
+proj=utm +zone=11 +datum=WGS84
# ---- From Coordinate System ----
#Lambert Conformal Conic
# Conic, Sph&Ell
# lat_1= and lat_2= or lat_0
# +units=us-ft +init=../nad/nad27:406 +proj=lcc +datum=NAD27 +lon_0=-116d15
# +lat_1=33d53 +lat_2=32d47 +lat_0=32d10 +x_0=609601.2192024384 +y_0=0
# +no_defs +ellps=clrk66 +nadgrids=conus
# ---- To Coordinate System ----
#Universal Transverse Mercator (UTM)
# Cyl, Sph
# zone= south
# +units=m +proj=utm +zone=11 +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
400000 100000
82688.56 3585908.61 0.000
As you can see from the ``cs2cs -ld'' command, I have a rather limited set of
datums hardcoded so far. I am interested in how I should name them, and what
datums I should use. Possible sources are EPSG, Mapinfo's datum list, PCI's
datum list, Geotrans's list (several of these are based on DMA TR 8350.2-B)
but this still leaves a question of what short names to use.
Anyways, I am being called to bath the kids, so I have to go now. The code
is checked into CVS in case anyone wants to look at it.
Best regards,
---------------------------------------+--------------------------------------
I set the clouds in motion - turn up | Frank Warmerdam, warmerda at home.com
light and sound - activate the windows | http://members.home.com/warmerda
and watch the world go round - Rush | Geospatial Programmer for Rent
----------------------------------------
OSRS PRoject PROJ Discussion List
To Subscribe: send a message to majordomo at remotesensing.org with 'subscribe osrs-proj' in the body
To Unsubscribe: send a message to majordomo at remotesensing.org with 'unsubscribe osrs-proj' in the body
To Report Problems: send a message to owner-osrs-proj at remotesensing.org
More information about the Proj
mailing list