Modified PJ structure to hold a_orig, es_orig, ellipsoid definition before

adjustment for spherical projections.
Modified pj_datum_transform() to use the original ellipsoid parameters,
not the ones adjusted for spherical projections.
Modified pj_datum_transform() to not attempt any datum shift via
geocentric coordinates if the source *or* destination are raw ellipsoids
(ie. PJD_UNKNOWN).  All per PROJ bug #1602, GDAL bug #2025.


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

5 files changed, 78 insertions, 15 deletions

diff --git a/ChangeLog b/ChangeLog
index 73c70727..22cffc98 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,3 +1,20 @@
+2007-11-25  Frank Warmerdam  <warmerdam@pobox.com>
+
+	* pj_transform.c: Do ellipsoid comparisons using the _orig ellipse
+	values rather than the adjusted one.  Use these original values for
+	any conversion to/from geocentric coordinates. 
+
+	Also, only do pj_datum_transform if neither the source nor destination 
+	is PJD_UNKNOWN.  This means we will no longer attempt via-geocentric 
+	adjustments for coordinate systems lacking a datum definition (having 
+	only an ellipsoid.
+
+	* projects.h, pj_init.c: added a_orig and es_orig values in the PJ
+	structure so we can distinguish between the originally requested 
+	ellipsoid, and the ellipsoid after adjustment for spherical projections
+
+	Todays changes courtesy of bug 1602. 
+
 2007-09-28  Frank Warmerdam  <warmerdam@pobox.com>
 
 	* nad/esri.extra: Add "900913" code for google mercator.
diff --git a/nad/td_out.dist b/nad/td_out.dist
index ac6d82ac..d08b0d0d 100644
--- a/nad/td_out.dist
+++ b/nad/td_out.dist
@@ -14,19 +14,19 @@ Test MD used where available
 79d58'00.000"W 36d58'00.000"N 0.0	79d57'59.128"W	36d58'0.501"N 0.000
 ##############################################################
 Test raw ellipse to raw ellipse
-79d58'00.000"W 37d02'00.000"N 0.0	79d58'W	37d1'50.514"N 699.067
-79d58'00.000"W 36d58'00.000"N 0.0	79d58'W	36d57'50.52"N 699.407
+79d58'00.000"W 37d02'00.000"N 0.0	79d58'W	37d2'N 0.000
+79d58'00.000"W 36d58'00.000"N 0.0	79d58'W	36d58'N 0.000
 ##############################################################
 Test NAD27 to raw ellipse
-79d00'00.000"W 35d00'00.000"N 0.0	78d59'59.107"W	34d59'58.563"N 718.018
+79d00'00.000"W 35d00'00.000"N 0.0	79dW	35dN 0.000
 ##############################################################
 Between two 3parameter approximations on same ellipsoid
 0d00'00.000"W 0d00'00.000"N 0.0	0dE	0dN 4.000
 79d00'00.000"W 45d00'00.000"N 0.0	78d59'59.821"W	44d59'59.983"N 0.540
 ##############################################################
 3param to raw ellipsoid on same ellipsoid
-0d00'00.000"W 0d00'00.000"N 0.0	0dE	0dN 5.000
-79d00'00.000"W 45d00'00.000"N 0.0	78d59'59.776"W	44d59'59.978"N 0.675
+0d00'00.000"W 0d00'00.000"N 0.0	0dE	0dN 0.000
+79d00'00.000"W 45d00'00.000"N 0.0	79dW	45dN 0.000
 ##############################################################
 Test simple prime meridian handling.
 0d00'00.000"W 0d00'00.000"N 0.0	1dW	0dN 0.000
diff --git a/src/pj_init.c b/src/pj_init.c
index 4b120080..205b468a 100644
--- a/src/pj_init.c
+++ b/src/pj_init.c
@@ -30,6 +30,15 @@
  ******************************************************************************
  *
  * $Log$
+ * Revision 1.19  2007/11/26 00:21:59  fwarmerdam
+ * Modified PJ structure to hold a_orig, es_orig, ellipsoid definition before
+ * adjustment for spherical projections.
+ * Modified pj_datum_transform() to use the original ellipsoid parameters,
+ * not the ones adjusted for spherical projections.
+ * Modified pj_datum_transform() to not attempt any datum shift via
+ * geocentric coordinates if the source *or* destination are raw ellipsoids
+ * (ie. PJD_UNKNOWN).  All per PROJ bug #1602, GDAL bug #2025.
+ *
  * Revision 1.18  2006/10/12 21:04:39  fwarmerdam
  * Added experimental +lon_wrap argument to set a "center point" for
  * longitude wrapping of longitude values coming out of pj_transform().
@@ -282,6 +291,9 @@ pj_init(int argc, char **argv) {
 	/* set ellipsoid/sphere parameters */
 	if (pj_ell_set(start, &PIN->a, &PIN->es)) goto bum_call;
 
+        PIN->a_orig = PIN->a;
+        PIN->es_orig = PIN->es;
+
 	PIN->e = sqrt(PIN->es);
 	PIN->ra = 1. / PIN->a;
 	PIN->one_es = 1. - PIN->es;
diff --git a/src/pj_transform.c b/src/pj_transform.c
index 3c2bed1a..40a7bf8d 100644
--- a/src/pj_transform.c
+++ b/src/pj_transform.c
@@ -30,6 +30,15 @@
  ******************************************************************************
  *
  * $Log$
+ * Revision 1.23  2007/11/26 00:21:59  fwarmerdam
+ * Modified PJ structure to hold a_orig, es_orig, ellipsoid definition before
+ * adjustment for spherical projections.
+ * Modified pj_datum_transform() to use the original ellipsoid parameters,
+ * not the ones adjusted for spherical projections.
+ * Modified pj_datum_transform() to not attempt any datum shift via
+ * geocentric coordinates if the source *or* destination are raw ellipsoids
+ * (ie. PJD_UNKNOWN).  All per PROJ bug #1602, GDAL bug #2025.
+ *
  * Revision 1.22  2007/09/11 20:32:25  fwarmerdam
  * mark the transient error array const
  *
@@ -197,7 +206,7 @@ int pj_transform( PJ *srcdefn, PJ *dstdefn, long point_count, int point_offset,
             }
         }
 
-        if( pj_geocentric_to_geodetic( srcdefn->a, srcdefn->es,
+        if( pj_geocentric_to_geodetic( srcdefn->a_orig, srcdefn->es_orig,
                                        point_count, point_offset, 
                                        x, y, z ) != 0) 
             return pj_errno;
@@ -294,7 +303,7 @@ int pj_transform( PJ *srcdefn, PJ *dstdefn, long point_count, int point_offset,
             return PJD_ERR_GEOCENTRIC;
         }
 
-        pj_geodetic_to_geocentric( dstdefn->a, dstdefn->es,
+        pj_geodetic_to_geocentric( dstdefn->a_orig, dstdefn->es_orig,
                                    point_count, point_offset, x, y, z );
 
         if( dstdefn->fr_meter != 1.0 )
@@ -464,8 +473,8 @@ int pj_compare_datums( PJ *srcdefn, PJ *dstdefn )
     {
         return 0;
     }
-    else if( srcdefn->a != dstdefn->a 
-             || ABS(srcdefn->es - dstdefn->es) > 0.000000000050 )
+    else if( srcdefn->a_orig != dstdefn->a_orig 
+             || ABS(srcdefn->es_orig - dstdefn->es_orig) > 0.000000000050 )
     {
         /* the tolerence for es is to ensure that GRS80 and WGS84 are
            considered identical */
@@ -598,6 +607,10 @@ int pj_geocentric_from_wgs84( PJ *defn,
 
 /************************************************************************/
 /*                         pj_datum_transform()                         */
+/*                                                                      */
+/*      The input should be long/lat/z coordinates in radians in the    */
+/*      source datum, and the output should be long/lat/z               */
+/*      coordinates in radians in the destination datum.                */
 /************************************************************************/
 
 int pj_datum_transform( PJ *srcdefn, PJ *dstdefn, 
@@ -611,16 +624,26 @@ int pj_datum_transform( PJ *srcdefn, PJ *dstdefn,
     pj_errno = 0;
 
 /* -------------------------------------------------------------------- */
+/*      We cannot do any meaningful datum transformation if either      */
+/*      the source or destination are of an unknown datum type          */
+/*      (ie. only a +ellps declaration, no +datum).  This is new        */
+/*      behavior for PROJ 4.6.0.                                        */
+/* -------------------------------------------------------------------- */
+    if( srcdefn->datum_type == PJD_UNKNOWN
+        || dstdefn->datum_type == PJD_UNKNOWN )
+        return 0;
+
+/* -------------------------------------------------------------------- */
 /*      Short cut if the datums are identical.                          */
 /* -------------------------------------------------------------------- */
     if( pj_compare_datums( srcdefn, dstdefn ) )
         return 0;
 
-    src_a = srcdefn->a;
-    src_es = srcdefn->es;
+    src_a = srcdefn->a_orig;
+    src_es = srcdefn->es_orig;
 
-    dst_a = dstdefn->a;
-    dst_es = dstdefn->es;
+    dst_a = dstdefn->a_orig;
+    dst_es = dstdefn->es_orig;
 
 /* -------------------------------------------------------------------- */
 /*      Create a temporary Z array if one is not provided.              */
diff --git a/src/projects.h b/src/projects.h
index 660c4a3b..9270755c 100644
--- a/src/projects.h
+++ b/src/projects.h
@@ -28,6 +28,15 @@
  ******************************************************************************
  *
  * $Log$
+ * Revision 1.27  2007/11/26 00:21:59  fwarmerdam
+ * Modified PJ structure to hold a_orig, es_orig, ellipsoid definition before
+ * adjustment for spherical projections.
+ * Modified pj_datum_transform() to use the original ellipsoid parameters,
+ * not the ones adjusted for spherical projections.
+ * Modified pj_datum_transform() to not attempt any datum shift via
+ * geocentric coordinates if the source *or* destination are raw ellipsoids
+ * (ie. PJD_UNKNOWN).  All per PROJ bug #1602, GDAL bug #2025.
+ *
  * Revision 1.26  2007/03/11 17:03:18  fwarmerdam
  * support drive letter prefixes on win32 and related fixes (bug 1499)
  *
@@ -289,8 +298,10 @@ typedef struct PJconsts {
         int is_geocent; /* proj=geocent ... not really a projection at all */
 	double
 		a,  /* major axis or radius if es==0 */
-		e,  /* eccentricity */
+                a_orig, /* major axis before any +proj related adjustment */
 		es, /* e ^ 2 */
+                es_orig, /* es before any +proj related adjustment */
+		e,  /* eccentricity */
 		ra, /* 1/A */
 		one_es, /* 1 - e^2 */
 		rone_es, /* 1/one_es */
@@ -298,7 +309,7 @@ typedef struct PJconsts {
 		x0, y0, /* easting and northing */
 		k0,	/* general scaling factor */
 		to_meter, fr_meter; /* cartesian scaling */
-
+    
         int     datum_type; /* PJD_UNKNOWN/3PARAM/7PARAM/GRIDSHIFT/WGS84 */
         double  datum_params[7];
         double  from_greenwich; /* prime meridian offset (in radians) */