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

[Paul Kelly]

Hello Frank

On Wed, 23 Jul 2003, Frank Warmerdam wrote:

> Paul Kelly wrote:
> > Still, IMHO I don't think it should have been trying to load the datum
> > shift file when the output co-ordinate system did not have a datum
> > specified, i.e. whether or not the input point is inside the region
> > covered by the shift file is relevant only if we are doing a datum shift,
> > not if just doing the inverse projection from projected co-ordinates to
> > lat/long. In this case I feel the PROJ pj_transform() function should just
> > do the inverse projection (although the fact that it refused was a useful
> > warning to you that something was wrong with your input values).
> >
> > Maybe Frank can help (copied to PROJ list). I feel the call to
> > pj_datum_transform() in pj_transform() should be conditional on
> > ( srcdefn->datum_type != PJD_UNKNOWN
> >   && dstdefn->datum_type != PJD_UNKNOWN )
> > however it is far from clear what is the obvious solution and I don't
> > pretend to fully understand the code there.
>
> Paul,
>
> I am boggled.  My belief had always been that pj_transform() would fail
> if the two coordinate systems had different ellpsoids, but no explicit
> datum information on one of them (or both).  However, trying this out,
> and examining the code I can find no support for why I thought this was
> the case.

Actually no that is a different problem (I'll come to what I was talking
about below later). In this case my understanding was wrong as well. I
had always thought it was not meaningful to convert between two ellipsoids
in this way (treating the point as having the same geocentric co-ordinates on
both ellipsoids, which is how I understand what you are saying) when datum
shift parameters were not available.

So I assumed that the second ellipsoid was ignored and the projection was
done assuming the first ellipsoid for both source and destination. Come to
think of it I can see no obvious reason why I made that assumption either.

> That said, I don't want to avoid the pj_datum_transform() call if on of the
> datums is PJD_UNKNOWN, because pj_datum_transform() also takes care of the
> ellipse change.

OK, makes sense when taken with my new understanding as described in the
previous paragraph.

> Hmmm, on further review, this isn't always done right either.  The ellpsoid
> conversion is only applied in the situation where I do a conversion through
> geocentric coordinates ... which is only done if one of the two coordinate
> systems is using a 3 or 7 parameter grid shift.   So currently we have the
> odd situation that if we do a conversion between ellpsoids with both
> coordinate system marked as PJD_UNKNOWN we get no change at all:
>
> eg.
>
> vc6 at gdal2200% cs2cs +proj=latlong +ellps=clrk66 +towgs84=0,0,0 +to +proj=latlong
> +ellps=WGS84 +towgs84=0,0,0
> 10 10
> 10dE    9d59'57.36"N 62.250
>
> But if we do the very same thing, but provide a zero shift to WGS84 like
> this, we at least get the ellipsoid conversion:
>
> vc6 at gdal2200% cs2cs  +proj=latlong +ellps=clrk66 +to +proj=latlong +ellps=WGS84
> 10 10
> 10dE    10dN 0.000
>
> Clearly I need to review this area of code.

I take it the order of these two examples should be switched round but I
get the point.

>
> However, Paul, your point I think was that you don't feel either "side" of
> a datum shift should be done unless the datum information for both the
> source and destionation are available, is that right?  I am not certain why
> you feel this is necessarily the case.

This is just how I always assumed it worked. I suppose it made sense to me
for the purposes of backward compatibility between two co-ordinate systems
that were the same, except one had datum transformation parameters
included in the PROJ definition, e.g.

$ echo '-103 44' | cs2cs +proj=longlat +ellps=clrk66 +nadgrids=./conus +to
+proj=longlat +ellps=clrk66 +nadgrids=./conus
103dW   44dN 0.000

produces no shift which is correct, but if the second co-ordinate system
didn't explicitly specify datum parameters

$ echo '-103 44' | cs2cs +proj=longlat +ellps=clrk66 +nadgrids=./conus +to
+proj=longlat +ellps=clrk66
103d0'1.629"W   43d59'59.959"N 0.000

then this is the wrong result in my opinion/understanding as no shift
should have taken place. Or even now that I understand about ellipsoid
shifts, both of these systems are on the same ellipsoid so again there
should be no shift.

Back to the original problem now. Very closely related to the second
example above... in fact that is it in a nutshell. I just felt that
it shouldn't even try to load the gridshift file if the two ellipsoids
were the same, as it would not need it. But in fact it is not just being
loaded but being used erroneously. This is a real problem at the minute in
GRASS modules that convert projected co-ordinates to latitude/longitude.

One more example:
it appears what you are saying is that if datum shift parameters are
provided for one side, but not the other then the shift should go ahead
on the source side, followed by projection onto the destination ellipsoid
(assuming a zero shift on this side)? But this is not what happens. E.g.
if this were the case then the following example from the website
cs2cs +proj=latlong +ellps=GRS80 +towgs84=-199.87,74.79,246.62 \
    +to +proj=latlong +datum=WGS84
could be equivalently written as
cs2cs +proj=latlong +ellps=GRS80 +towgs84=-199.87,74.79,246.62 \
    +to +proj=latlong +ellps=WGS84
But if that is done then no transformation happens because a datum shift
has only been specified on one side.

I think that examples like these were maybe the root cause of my
understanding that datum shift parameters needed to be supplied on both
sides for anything to happen at all.

Thanks for your attention to this

Paul

(Sorry, not sure how to copy this to the bug tracker)
----------------------------------------
PROJ.4 Discussion List
See http://www.remotesensing.org/proj for subscription, unsubscription
and other information.