[OSRS-PROJ] Re: remotesensing.org PROJ reading datum gridshift files unnecessarily?

[Paul Kelly]

May I suggest the attached patch to pj_transform.c as a starting point to
fixing the issues described in bug 368. The two main issues it addresses
are:

1. 
A grid shift is now only applied when both the source and destination
co-ordinate systems have datum parameters specified. This solves the
probem whereby the following two command-lines gave different output:

echo '-103 44' | cs2cs +proj=longlat +ellps=clrk66 +nadgrids=conus +to \
+proj=longlat +ellps=clrk66

echo '-103 44' | cs2cs +proj=longlat +ellps=clrk66 +to +proj=longlat \ 
+ellps=clrk66

(the first one incorrectly applied a shift to the input values)

2.
Transformation unconditionally goes through geocentric co-ordinates and
performs an ellipsoid conversion, even if datum parameters are specified
for one system but not the other or for neither. This solves the problem
whereby the following two command-lines gave different output:

echo '10 10' | cs2cs  +proj=latlong +ellps=clrk66 +to +proj=latlong \
+ellps=WGS84

echo '10 10' | cs2cs +proj=latlong +ellps=clrk66 +towgs84=0,0,0 +to \
+proj=latlong +ellps=WGS84 +towgs84=0,0,0

(the first one did not perform an ellipsoid conversion while the second
one did; after the patch is applied they both do)

As discussed before, it is not clear whether a simple ellipsoid conversion
should be done without warning if different ellipsoids but incomplete or
missing datum parameters are specified. The alternative solution would be
for pj_transform() to fail to transform the points in this case. However
on balance I feel the former behaviour is better because
(a) It does not change the existing behaviour too much, and
(b) It may be hard to test if two ellipsoids are the same (comparing
floating point values), e.g. deep inside GRASS there may be a rounding
error converting between eccentricity squared and inverse flattening, so
the parameters presented to PROJ may appear to be for different
ellipsoids. But after the conversion to and from geocentric co-ordinates
on these 'different' ellipsoids, the resulting error would presumably be
negligible.

Ideally, if we want to use pj_transform() to do a simple projection to or
from latitude/longitude, we would not need to specify an ellipsoid for the
second co-ordinate system (it would be assumed to be the same as the
first), but I realise that this is not possible the way pj_init() is
currently used to define the co-ordinate systems (it must have an
ellipsoid specified). The problem I have in mind here is that the method
used in GRASS to determine the second ellipsoid *may* result in a
negligible difference in the parameters that are supplied to PROJ such
that a floating point equality comparison may make the two ellipsoids
appear different. However as described above, this does not cause a
problem with the current PROJ or with the attached patch.

I hope it will be of some use

Paul Kelly

-------------- next part --------------
--- pj_transform.c.orig	Mon Mar 31 16:10:45 2003
+++ pj_transform.c	Fri Jul 25 20:45:01 2003
@@ -485,62 +485,62 @@
 /*	If this datum requires grid shifts, then apply it to geodetic   */
 /*      coordinates.                                                    */
 /* -------------------------------------------------------------------- */
-    if( srcdefn->datum_type == PJD_GRIDSHIFT )
+    if( srcdefn->datum_type != PJD_UNKNOWN
+        && dstdefn->datum_type != PJD_UNKNOWN )
     {
-        pj_apply_gridshift( pj_param(srcdefn->params,"snadgrids").s, 0, 
-                            point_count, point_offset, x, y, z );
-        CHECK_RETURN;
+        if( srcdefn->datum_type == PJD_GRIDSHIFT )
+        {
+            pj_apply_gridshift( pj_param(srcdefn->params,"snadgrids").s, 0, 
+                                point_count, point_offset, x, y, z );
+            CHECK_RETURN;
 
-        src_a = SRS_WGS84_SEMIMAJOR;
-        src_es = 0.006694379990;
-    }
+            src_a = SRS_WGS84_SEMIMAJOR;
+            src_es = 0.006694379990;
+        }
 
-    if( dstdefn->datum_type == PJD_GRIDSHIFT )
-    {
-        dst_a = SRS_WGS84_SEMIMAJOR;
-        dst_es = 0.006694379990;
+        if( dstdefn->datum_type == PJD_GRIDSHIFT )
+        {
+            dst_a = SRS_WGS84_SEMIMAJOR;
+            dst_es = 0.006694379990;
+        }
     }
-        
-/* ==================================================================== */
-/*      Do we need to go through geocentric coordinates?                */
-/* ==================================================================== */
-    if( srcdefn->datum_type == PJD_3PARAM 
-        || srcdefn->datum_type == PJD_7PARAM
-        || dstdefn->datum_type == PJD_3PARAM 
-        || dstdefn->datum_type == PJD_7PARAM)
-    {
+       
 /* -------------------------------------------------------------------- */
 /*      Convert to geocentric coordinates.                              */
 /* -------------------------------------------------------------------- */
-        pj_geodetic_to_geocentric( src_a, src_es,
-                                   point_count, point_offset, x, y, z );
-        CHECK_RETURN;
+    pj_geodetic_to_geocentric( src_a, src_es,
+                               point_count, point_offset, x, y, z );
+    CHECK_RETURN;
 
 /* -------------------------------------------------------------------- */
 /*      Convert between datums.                                         */
 /* -------------------------------------------------------------------- */
-        if( srcdefn->datum_type != PJD_UNKNOWN
-            && dstdefn->datum_type != PJD_UNKNOWN )
-        {
-            pj_geocentric_to_wgs84( srcdefn, point_count, point_offset,x,y,z);
-            CHECK_RETURN;
+    if( srcdefn->datum_type != PJD_UNKNOWN
+        && dstdefn->datum_type != PJD_UNKNOWN 
+        && (srcdefn->datum_type == PJD_3PARAM 
+        || srcdefn->datum_type == PJD_7PARAM
+        || dstdefn->datum_type == PJD_3PARAM 
+        || dstdefn->datum_type == PJD_7PARAM) )
+    {
+        pj_geocentric_to_wgs84( srcdefn, point_count, point_offset,x,y,z);
+        CHECK_RETURN;
 
-            pj_geocentric_from_wgs84( dstdefn, point_count,point_offset,x,y,z);
-            CHECK_RETURN;
-        }
+        pj_geocentric_from_wgs84( dstdefn, point_count,point_offset,x,y,z);
+        CHECK_RETURN;
+    }
 
 /* -------------------------------------------------------------------- */
 /*      Convert back to geodetic coordinates.                           */
 /* -------------------------------------------------------------------- */
-        pj_geocentric_to_geodetic( dst_a, dst_es,
-                                   point_count, point_offset, x, y, z );
-        CHECK_RETURN;
-    }
+    pj_geocentric_to_geodetic( dst_a, dst_es,
+                               point_count, point_offset, x, y, z );
+    CHECK_RETURN;
 
 /* -------------------------------------------------------------------- */
 /*      Apply grid shift to destination if required.                    */
 /* -------------------------------------------------------------------- */
-    if( dstdefn->datum_type == PJD_GRIDSHIFT )
+    if( srcdefn->datum_type != PJD_UNKNOWN
+        && dstdefn->datum_type == PJD_GRIDSHIFT )
     {
         pj_apply_gridshift( pj_param(dstdefn->params,"snadgrids").s, 1,
                             point_count, point_offset, x, y, z );