added datum support

git-svn-id: http://svn.osgeo.org/metacrs/proj/trunk@852 4e78687f-474d-0410-85f9-8d5e500ac6b2

1 files changed, 65 insertions, 4 deletions

diff --git a/src/pj_init.c b/src/pj_init.c
index 68b77c9c..944b85b2 100644
--- a/src/pj_init.c
+++ b/src/pj_init.c
@@ -10,6 +10,7 @@ static const char SCCSID[]="@(#)pj_init.c	4.13   95/09/05 GIE REL";
 	static paralist
 *start;
 extern FILE *pj_open_lib(char *, char *);
+
 	static paralist *
 get_opt(FILE *fid, char *name, paralist *next) {
 	char sword[52], *word = sword+1;
@@ -26,8 +27,14 @@ get_opt(FILE *fid, char *name, paralist *next) {
 				first = 0;
 			else if (!first && word[1] == '>')
 				break;
-		} else if (!first && !pj_param(start, sword).i)
-			next = next->next = pj_mkparam(word);
+		} else if (!first && !pj_param(start, sword).i) {
+                        /* don't default ellipse if datum is set */
+                        if( strncmp(word,"ellps=",6) != 0 
+                            || !pj_param(start, "tdatum").i )
+                        {
+                            next = next->next = pj_mkparam(word);
+                        }
+                }
 	if (errno == 25)
 		errno = 0;
 	return next;
@@ -64,7 +71,17 @@ get_init(paralist *next, char *name) {
 		errno = 0; /* unknown problem with some sys errno<-25 */
 	return next;
 }
-	PJ *
+
+/************************************************************************/
+/*                              pj_init()                               */
+/*                                                                      */
+/*      Main entry point for initialing a PJ projections                */
+/*      definition.  Note that the projection specific function is      */
+/*      called to do the initial allocation so it can be created        */
+/*      large enough to hold projection specific parameters.            */
+/************************************************************************/
+
+PJ *
 pj_init(int argc, char **argv) {
 	char *s, *name;
 	void *(*proj)(PJ *);
@@ -73,6 +90,7 @@ pj_init(int argc, char **argv) {
 	PJ *PIN = 0;
 
 	errno = pj_errno = 0;
+
 	/* put arguments into internal linked list */
 	if (argc <= 0) { pj_errno = -1; goto bum_call; }
 	for (i = 0; i < argc; ++i)
@@ -81,6 +99,7 @@ pj_init(int argc, char **argv) {
 		else
 			start = curr = pj_mkparam(argv[i]);
 	if (pj_errno) goto bum_call;
+
 	/* check if +init present */
 	if (pj_param(start, "tinit").i) {
 		paralist *last = curr;
@@ -89,36 +108,62 @@ pj_init(int argc, char **argv) {
 			goto bum_call;
 		if (curr == last) { pj_errno = -2; goto bum_call; }
 	}
+
 	/* find projection selection */
 	if (!(name = pj_param(start, "sproj").s))
 		{ pj_errno = -4; goto bum_call; }
 	for (i = 0; (s = pj_list[i].id) && strcmp(name, s) ; ++i) ;
 	if (!s) { pj_errno = -5; goto bum_call; }
+
 	/* set defaults, unless inhibited */
 	if (!pj_param(start, "bno_defs").i)
 		curr = get_defaults(curr, name);
 	proj = pj_list[i].proj;
+
 	/* allocate projection structure */
 	if (!(PIN = (*proj)(0))) goto bum_call;
 	PIN->params = start;
+        PIN->is_latlong = 0;
+
+        /* set datum parameters */
+        if (pj_datum_set(start, PIN)) goto bum_call;
+
 	/* set ellipsoid/sphere parameters */
 	if (pj_ell_set(start, &PIN->a, &PIN->es)) goto bum_call;
+
 	PIN->e = sqrt(PIN->es);
 	PIN->ra = 1. / PIN->a;
 	PIN->one_es = 1. - PIN->es;
 	if (PIN->one_es == 0.) { pj_errno = -6; goto bum_call; }
 	PIN->rone_es = 1./PIN->one_es;
+
+        /* Now that we have ellipse information check for WGS84 datum */
+        if( PIN->datum_type == PJD_3PARAM 
+            && PIN->datum_params[0] == 0.0
+            && PIN->datum_params[1] == 0.0
+            && PIN->datum_params[1] == 0.0
+            && PIN->a == 6378137.0
+            && ABS(PIN->es - 0.006694379990) < 0.000000000050 )/*WGS84/GRS80*/
+        {
+            PIN->datum_type = PJD_WGS84;
+        }
+        
 	/* set PIN->geoc coordinate system */
 	PIN->geoc = (PIN->es && pj_param(start, "bgeoc").i);
+
 	/* over-ranging flag */
 	PIN->over = pj_param(start, "bover").i;
+
 	/* central meridian */
 	PIN->lam0=pj_param(start, "rlon_0").f;
+
 	/* central latitude */
 	PIN->phi0 = pj_param(start, "rlat_0").f;
+
 	/* false easting and northing */
 	PIN->x0 = pj_param(start, "dx_0").f;
 	PIN->y0 = pj_param(start, "dy_0").f;
+
 	/* general scaling factor */
 	if (pj_param(start, "tk_0").i)
 		PIN->k0 = pj_param(start, "dk_0").f;
@@ -130,6 +175,7 @@ pj_init(int argc, char **argv) {
 		pj_errno = -31;
 		goto bum_call;
 	}
+
 	/* set units */
 	s = 0;
 	if (name = pj_param(start, "sunits").s) { 
@@ -144,6 +190,7 @@ pj_init(int argc, char **argv) {
 		PIN->fr_meter = 1. / PIN->to_meter;
 	} else
 		PIN->to_meter = PIN->fr_meter = 1.;
+
 	/* projection specific initialization */
 	if (!(PIN = (*proj)(PIN)) || errno || pj_errno) {
 bum_call: /* cleanup error return */
@@ -160,7 +207,18 @@ bum_call: /* cleanup error return */
 	}
 	return PIN;
 }
-	void
+
+/************************************************************************/
+/*                              pj_free()                               */
+/*                                                                      */
+/*      This is the application callable entry point for destroying     */
+/*      a projection definition.  It does work generic to all           */
+/*      projection types, and then calls the projection specific        */
+/*      free function (P->pfree()) to do local work.  This maps to      */
+/*      the FREEUP code in the individual projection source files.      */
+/************************************************************************/
+
+void
 pj_free(PJ *P) {
 	if (P) {
 		paralist *t = P->params, *n;
@@ -170,7 +228,10 @@ pj_free(PJ *P) {
 			n = t->next;
 			pj_dalloc(t);
 		}
+
 		/* free projection parameters */
 		P->pfree(P);
 	}
 }
+
+