[OSRS-PROJ] Datum Shifting

Frank Warmerdam warmerda at home.com
Thu Jul 6 21:04:41 EDT 2000


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:

	The name is looked up in a list of datums, and used to imply an
        ellipsoid, and a datum definition (+towgs84= or +nadgrids=)
        Defines a shift to WGS84 for a datum. 

        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). 

        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  
       NAD83 GRS80        towgs84=0,0,0                 
       NAD27 clrk66       nadgrids=conus                

warmerda at cs46980-c[386]% cs2cs -v +proj=latlong +datum=NAD27 +to +proj=latlong
# ---- 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