[Proj] Lat-Lon values under different ellipsoides

[Jan Hartmann]

Thanks Noel. It seems to be a bit more complex (they use both 3D and 2D 
transformations from appendixes B and C), but essentially you are right. 

I'm sorry to not have given the URL of this Ordnance Survey publication:

http://www.ordnancesurvey.co.uk/oswebsite/gps/docs/A_Guide_to_Coordinate_Systems_in_Great_Britain.pdf

This booklet ( I would advise everyone to print it out: it's very good) 
also answered my question about how different lat-lon values are 
computed for different Datums (p. 16). Essentially, one initial point is 
chosen with reference to the the center, orientation, size and shape of 
the chosen ellipsoid,  and with that a coordinate system is defined  in 
which all other points on the earth have unique coordinates. Different 
datums give different parameters for the initial control point, and so 
for the whole coordinate system.  I still don't get the relationship 
between these values and astronomical lat-lon values though.

This also makes clear to me a discussion on this list about the 
difference between the DATUM and ELLPS parameters in PROJ. As explained 
on p. 16 of the OS book, an ellipse has three sets of parameters: origin 
(x,y,z), rotation (around the x,y,z axes), and size (minor and major 
axes). The original ELLPS parameter only takes the third one into 
account: size, aka minor and major axis, while assuming that origin and 
rotation are the same. To transform both origin and rotation of the 
ellipsoid (but not its size), a so called Helmert transformation is 
used, needing seven parameters: three for the origin, three for the 
rotation, and one additional scale parameter (p. 32). Essentially, this 
is what the seven-number +TOWGS84 parameter in PROJ does. Also, there is 
a formula to transform both the origin and the size of the ellipsoid 
(but not the rotation): the Molodensky transformation (p. 34). I guess 
that this is the three-number +TOWGS84 parameter in PROJ, used for 
origin-translation only. There seems to be no one-for-all formula for a 
direct transformation of all three features of an ellipsoid: origin, 
rotation and size.

This may be the source of confusion in PROJ: the +ELLPS parameter 
transforms only the size of the ellipsoid, the +TOWGS84 parameter 
transforms the origin and rotation of the ellipsoid, and the +DATUM 
parameter transform all three features: origin, rotation and size. The 
+DATUM parameter is not a single transformation however: it is just 
shorthand for +ELLPS=... +TOWGS84=... Although from a mathematical and 
programming point of view there is nothing wrong with this, PROJ's 
terminology is confusing: +TOWGS84 does not transform to the WGS84 
datum: it just transforms origin and rotation between two arbitrary 
ellipsoids. To transform to the WGS84 Datum you still need the 
+ELLPS=WGS84 parameter. From a generalized point of view, the +DATUM 
parameter is unnecessary. It would be much more precise to always 
specify both the +ELLPS and the +TOWGS84 ("Helmert") parameters. I can 
understand why the DATUM solution was chosen, as almost all datum 
transforms are from and to the WGS84 datum, but it obscures the more 
general issue. Same goes for the decision to disregard the +ELLPS 
parameter in PROJ when no further DATUM information is given 
(http://lists.maptools.org/pipermail/proj/2009-January/004303.html), 
although there too I understand the immediate cause for this. And 
finally, the "nadgrids=@null"  would become superfluous with the more 
general approach: always specify two transformations for a datum shift, 
and one for an ellipsoid transformation.

I agree with Gerald Evenden on most of his remarks about the separation 
of ellipsoid and datum computations: there seems to be no complete 
analytical transformation formula for ellipsoids as there is in planar 
geometry, it is always a two-step procedure. Having said that, I 
nevertheless think that the complete procedure can be easily 
implementend within PROJ. It's more a matter of terminology than of 
mathematics.

So at last I have seen the light in this matter, although after 
considerable head-scratching. If I am wrong though, please let the list 
know. There is an enormous amount of knowledge about this matter on this 
list, and I am just a simple historian (medievalist) who happens to have 
worked a lot with old maps and had to learn something about programming 
in the meantime. I still don't know much about geodesy, but I have been 
surprised how much you need to know about the exact form of the earth 
and its geoid to get an exact georeference: that was definitely not part 
of our curriculum.  Yet, and this is a fundamental point for me, this 
precision should be aimed at. Theoretically, precisions of one to five 
meters should be possible for large scale maps, and I believe that I can 
get very far in that direction with PROJ, GDAL and the answers from this 
list. A lot of Internet cartography lacks this final precision, not only 
with historical maps, but also e.g. in the case of Google maps (there 
have been extensive discussions about Google's accuracy on this list). 
It's relatively easy to put lots and lots of maps on the Internet, but 
it is enormously difficult to get them as precise as they could (and 
should) be. Fifty meters just isn't good enough. Personally, I find it 
worthwile (and terrifically interesting) to aim at this maximum 
precision, and it will certainly pay back, but only in the long run. In 
that same long run however, much of the rubbish that is now on Internet 
can be happily thrown away.

And if someone on this list needs help with his or her medieval 
charters, I would be happy to assist :-)

Jan


Noel Zinn wrote:
>
> Jan,
>
>  
>
> Not derived from *projected* coordinates (i.e. 3D => 2D) at all.  The 
> Ordnance Survey have taken you on a tour from one 3D space 
> (lat/lon/height) to another 3D space (what we're calling geocentric 
> Cartesian coordinates X/Y/Z, also known as ECEF) and back to the 
> original 3D space (lat/lon/height, but referenced to a different 
> ellipsoid).  So you've never left three dimensions even though the 
> Ordnance Survey suggested that you toss datum B ellipsoid height.
>
>  
>
> You've previously asked about an ellipsoid switch.  If nothing happens 
> in the geocentric Cartesian domain (ECEF is so much easier to write!), 
> that is, no translations at the geocenter in X, Y or Z and no 
> rotations about any of these axes and no scale change, then you've 
> **just** switched ellipsoids.  Nevertheless, the resulting change in 
> ellipsoid height is an important clue that you're not **on** the new 
> ellipsoid.  Once the ellipsoid height is tossed, the point at that 
> lat/lon on the new ellipsoid with zero height is a different point 
> than that you started with.
>
>  
>
> Noel Zinn
>
>  
>
> ------------------------------------------------------------------------
>
> *From:* proj-bounces at lists.maptools.org 
> [mailto:proj-bounces at lists.maptools.org] *On Behalf Of *Jan Hartmann
> *Sent:* Friday, February 06, 2009 8:19 AM
> *To:* PROJ.4 and general Projections Discussions
> *Subject:* Re: [Proj] Lat-Lon values under different ellipsoides
>
>  
>
> As an addition, look how the British Ordnance survey solves the 
> problem of converting Lat-Lon values from one ellipse to another:
>
> " To summarise: for a simple datum change of latitude and longitude 
> coordinates from datum A to datum B, first convert to Cartesian 
> coordinates (formulae in annexe B) taking all ellipsoid heights as 
> zero and using the ellipsoid parameters of datum A; then apply a 
> Helmert transformation from datum A to datum B using equation (3); 
> finally convert back to latitude and longitude using the ellipsoid 
> parameters of datum B (formulae in annexe C), discarding the datum B 
> ellipsoid height. "
>
> This would mean that different lat-lon values for different ellipsoids 
> are *derived* from projected coordinates, not the other way round, as 
> I thought. I still don't get the relationship between those computed 
> lat-lon values and the astronomical ones though.
>
> Jan
>
> Frank Warmerdam wrote:
>
> Jan Hartmann wrote:
>   
>> Hi,
>>  
>> This is something I have long been banging my head against. I think it 
>> is a bug, but I am not sure. If I take a lat-lon value, computed on a 
>> particular ellipsoid, and convert it to the lat-lon value on another 
>> ellipsoid, I should get a different value, right? (e.g. cs2cs 
>> +proj=longlat +ELLPS=bessel +to +proj=longlat +ELLPS=WGS84). PROJ4 
>> always gives the same value, but I have an extensive list of coordinates 
>> of church towers in the Netherlands with their lat-lon values in 1850, 
>> based an a slightly smaller ellipsoid than we use nowadays, and the 
>> lat-lon of the same towers derived from our modern RD-system, based on 
>> the Bessel-ellipsoid, and without the WGS84 correction. There is a 
>> systematic difference of about 50 meters. If I do the same computation 
>> with the projected coordinates, I get the correct answer. Moreover, in 
>> that case the transformation changes when I change the 
>> ellipse-parameter, something that does not happen with lat-lon coordinates.
>>  
>> So, is this a bug in PROJ? If so, can someone with geodetic experience 
>> here explain to me how people can get different lat-lon values for the 
>> same point, based on astronomical measurements?
>>     
>  
> Jan,
>  
> As of PROJ 4.6.0, the policy is to not attempt any conversion of lat/long
> values between coordinate systems where only an ellipsoid is given.  So, to
> get a datum shift it is now necessary to provide some sort of datum
> shift information for both the source and destination coordinate systems.
>  
> This is a deliberate change of policy to avoid lots of other complaints in
> the past.
>  
> I would suggest using something appropriate like +datum=potsdam for your
> Bessel data, and +datum=WGS84 instead of +ELLPS=WGS84.
>  
> Best regards,
>   
> ------------------------------------------------------------------------
>
> _______________________________________________
> Proj mailing list
> Proj at lists.maptools.org
> http://lists.maptools.org/mailman/listinfo/proj
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://lists.maptools.org/pipermail/proj/attachments/20090207/47d61e56/attachment-0001.htm