[Proj] Transformations in small areas - was Stereo 1970 (EPSG 31700)

[Noel Zinn (cc)]

Thanks, Mikael.  I appreciate the exchange.

It may be useful to identify our different perspectives on this topic and to 
put some error bars around the term "accuracy".  Perspective and accuracy 
are "apples and oranges" issues, but they inextricably impact the 
discussion.

Although you are a software developer, your perspective seems to be that of 
the consumer seeking a "black-box mechanism that can transform long/lat from 
one datum to another, as accurately as possible".  You provide that to your 
customer (the consumer).  My perspective has been that of a deriver of datum 
transformations in (small) oil exploration frontiers throughout the world, 
e.g. the offshore lease blocks of some West African country.  Having the 
means to derive 3-, 7- or 10-parameter transformations, my choice was always 
a 3-parameter translation.  Why?  Because software support for M-B is 
regrettably scarce (a situation that you and Proj.4 can improve upon), and 
because 7-parameter derivations in small frontiers are ill conditioned (my 
thesis), and also because users all too frequently get the signs of the 
rotations in a 7-parameter transformation wrong, and, finally, because 
7-parameter transformations are no more "accurate" in a small area than a 
3-parameter translation derived from the same data set.

I gracefully accept that, from your consumer perspective, my advice to use 
3-parameter translations in a small area may seem "bad".  After all, why not 
use a 7-parameter transformation in a small area if you've got one, 
especially if it's from an authoritative source?  But from my deriver's 
perspective - and as an advocate of sound geodetic practice - it's the only 
responsible position to take and I don't recant.

Regarding accuracy options, you cite a couple of Great Britain OSGB 1936 to 
WGS84 transformations, namely tfm code 1195 (3 parameters) by the DMA good 
to 21 meters over 38 points, a remarkably large number for the DMA, but 
dating to 1987 which is effectively pre-GPS (lots of Doppler error there), 
and 1314 (7 parameters) of uncertain derivation for "oil exploration" but 
cited by UKOOA and Eurographics, and reportedly good to 2 or 4 or 5 meters. 
You've not cited tfm code 1039, derived by the Ordnance Survey of Great 
Britain itself, providing 0.1 meters between OSGB 1936 and ETRS89, which is 
WGS84 to less than a meter and the transformation one should be using if one 
is interested in accuracy, but it's a grid interpolation.

I assert that a 7-parameter transformation derived in a small area is no 
more accurate than a 3-parameter translation derived from the same data set. 
Because the conventional wisdom is just the opposite we have too many bad 
7-parameter transformations and too few good 3-parameter translations.  Not 
good geodetic practice in my opinion.  Tfm Code 1039 provides us the 
opportunity to test my assertion if you're interested.  "Small" can be 
defined quantitatively based on the dilution of precision of the 7-parameter 
adjustment, but let's just accept that the United Kingdom is small for our 
purposes.  Using 1039 we can populate some number of randomly distributed 
points throughout the UK in both OSGB 1936 and a WGS84 surrogate (ETRS89) to 
within 0.1 meters per point.  Then we need a bit of random error for both 
OSGB 1936 and WGS84 in amounts to be determined, 1039 having removed the 
systematic errors (datum distortions).  This "noise" is important for the 
propagation of random error to the derived parameters.  These points become 
our truth data.  Given this design it's then just mechanics to derive 3- and 
7-parameter transformations, back compute one datum from the other, and test 
my assertion.  We can also compare our derived transformations with Tfms 
Codes 1195 and 1314.  The exercise can be repeated with different random 
noise (in the same quantities) to test the stability of the derived 
parameters.

More than an evening's work, but an interesting thought experiment and a 
quantitative way to test whether small-area 7-parameter transformations are 
more accurate than 3-parameter transformations.

Regards,
Noel

Noel Zinn, Principal, Hydrometronics LLC
+1-832-539-1472 (office), +1-281-221-0051 (cell)
noel.zinn at hydrometronics.com (email)
http://www.hydrometronics.com (website)

--------------------------------------------------
From: "Mikael Rittri" <Mikael.Rittri at carmenta.com>
Sent: Thursday, October 07, 2010 9:46 AM
To: "PROJ.4 and general Projections Discussions" <proj at lists.maptools.org>
Subject: Re: [Proj] Stereo 1970 (EPSG 31700)

> Hello Noel,
>
> I think you are right in everything you say - except the final conclusion
> that 3-parameter transformations are to be preferred if the software
> does not support Molodensky-Badekas.
>
> Some specific comments:
>
>> (Perhaps you can advise me if this indictment of geodetic software
>> applies to Proj4.)
>
> Yes, as far as I know, Proj4 (or cs2cs) does not support 10 parameters
> (M-B parameters)
> in the +towgs84 parameter.
>
>> It may be true (not positive as I haven't tested this) that the
>> coordinate
>> results of a given M-B transformation can be reproduced by some
>> 7-parameter
>> transformation as you suggest.
>
> It was Melita Kennedy who taught me how to translate an M-B transformation
> into an equivalent 7-parameter transformation. (Strictly speaking, I
> haven't
> been able to test if I understood everything, since I don't have access to
> independent software that supports M-B. But I am pretty confident.)
>
>> ... observe the first 3 parameters of the 7-parameter transformations ...
>> Shocking, isn't it?  They're scattered meaninglessly all over the place,
>> ...
>
> Well, let's say it is unfortunate and disappointing, if it is
> important to be able to interpret and understand the parameters.
>
> On the other hand, most of my customers just need some kind of black-box
> mechanism that can transform long/lat from one datum to another, as
> accurately as possible.  (Or as accurately as necessary, perhaps.)
> So they don't care if the first 3 parameters are widely different
> between 7-parameter transformations for neighbouring regions.
>
>> Until all of (y)our software supports M-B, use a 3-parameter translation
>> instead!
>
> No, that's bad advice. There are plenty of examples for regions that
> you say are too small (like medium-sized European countries), where
> the best available 3-parameter transformation can give errors greater
> than 15 meters, while there are published 7-parameter transformations
> that give errors of just a few meters.
>
> Example: transforming the British OSGB 1936 to WGS84:
>
> Best published 3-parameter transformation: EPSG:1195, "OSGB 1936 to WGS 84
> (1)".
> EPSG Accuracy entry: 21 metre.
> EPSG Scope says: "Accuracy 10m, 10m and 15m in X, Y and Z axes."
>
> Best published 7-parameter transformation: EPSG:1314, "OSGB 1936 to WGS 84
> (6)".
> EPSG Accuracy entry: 2 metre.
> EPSG Scope says: "Accuracy better than 4m and generally better than 2m."
> (British sources say the worst error is 5 meters.)
>
> There are many applications where a 5-meter error can be acceptable,
> while 21 meters or more wouldn't be. (Fleet management, ambulance
> dispatch...)
>
> So, if my software supports 3-parameter and 7-parameter datum shifts,
> but not Molodensky-Badekas (10-parameter) or grid shift files (or if
> you don't need the accuracy of a grid shift file), then what?
> Surely it is better to get 2 - 4 meter accuracy instead of 21 meters,
> even if I cannot interpret or understand the seven parameters at all
> (just as long as they give good results).
>
> Best regards,
>
> Mikael Rittri
> Carmenta AB
> Sweden
> www.carmenta.com
>
> -----Original Message-----
> From: proj-bounces at lists.maptools.org
> [mailto:proj-bounces at lists.maptools.org] On Behalf Of Noel Zinn (cc)
> Sent: den 7 oktober 2010 12:17
> To: PROJ.4 and general Projections Discussions
> Subject: Re: [Proj] Stereo 1970 (EPSG 31700)
>
> Hi Mikael,
>
> It is true that there is no numerical problem using a 7-parameter
> transformation derived from observations in a small area, which means, I
> will add, derived poorly.
>
> It is an unfortunate reality that more geodetic software supports
> 7-parameter transformations than supports Molodensky-Badekas (M-B)
> transformations, sometimes called 10-parameter transformations due to the
> additional 3 coordinates of the topographic rotation center, but these are
> givens of the adjustment and not additional parameters to be solved for.
> (Perhaps you can advise me if this indictment of geodetic software applies
> to Proj4.)  Anyway, in the example you provide (EPSG Tfm Codes 4827 and
> 4829) the Slovakian authorities had this deficiency in mind in supplying a
> 7-parameter transformation to replace a M-B transformation.
>
> It may be true (not positive as I haven't tested this) that the coordinate
> results of a given M-B transformation can be reproduced by some
> 7-parameter transformation as you suggest.
>
> The problem with 7-parameter transformations derived in small area is more
> fundamental, however.  The 7 parameters of a 7-parameter (Helmert,
> Bursa-Wolfe, or similarity) transformation have well-known and
> well-documented physical interpretations.  The 3 translations represent
> the relative displacement of the geocenters of the two ellipsoids involved
> (in the ECEF X, Y and Z directions); the 3 rotations represent rotations
> about the X, Y and Z axes; and the scale represents an inflation (or
> deflation) of an ellipsoid (i.e. the linear unit).  These 7 parameters are
> not like the parameters of a multiple regression equation (MRE) of the
> kind that used to be published by the Defense Mapping Agency (DMA,
> predecessor to NIMA and
> NGA) for transforming datums with more "accuracy" than a 3-parameter
> translation.  Those parameters can be pencil whipped into submission
> because they "mean" nothing.  The 7 parameters of a Helmert (et al)
> transformation mean something and that meaning is badly corrupted by the
> ill conditioning of an adjustment over a small area.
>
> Let me illustrate this point with my own EPSG tfm code examples.  Codes
> 15714, 15716, 15718, 15720, 15722, 15724, 15716 and 15728 offer
> 7-parameter transformations from Bogota 1975 to Magna-Sirgas for regions I
> through VIII in Colombia and codes 15730 through 15737 offer M-B
> parameters for the same regions respectively.  Code 1125 offers a
> 3-parameter translation between Bogota 1975 and WGS84, which EPSG regard
> to be the same as Magna-Sirgas at the 1-meter level.  So, the 3 parameters
> of tfm 1125 provide a 1-meter assessment of the geocentric displacement
> between the International 1924 ellipsoid of Bogota 1975 and GRS80 of
> Magna-Sirgas.  Now, observe that these same 3 translations (307,
> 304, -318) are pretty much the same as the first 3 parameters of all the
> M-B transformations for the eight regions (tfm codes 15730 through 15737).
> Finally, observe the first 3 parameters of the 7-parameter transformations
> for the eight regions (tfm codes 15714, 15716, 15718, 15720, 15722, 15
> 724, 15716 and 15728).  Shocking, isn't it?  They're scattered
> meaninglessly all over the place, a consequence of the ill conditioning of
> deriving a 7-parameter transformation in an area an eighth the size of
> Colombia.
>
> The derivation of a 7-parameter Helmert (et al) transformation over a
> small area is pencil whipping, not geodesy.  Until all of (y)our software
> supports M-B, use a 3-parameter translation instead!
>
> Regards,
> Noel
>
> Noel Zinn, Principal, Hydrometronics LLC
> +1-832-539-1472 (office), +1-281-221-0051 (cell)
> noel.zinn at hydrometronics.com (email)
> http://www.hydrometronics.com (website)
>
> --------------------------------------------------
> From: "Mikael Rittri" <Mikael.Rittri at carmenta.com>
> Sent: Thursday, October 07, 2010 2:13 AM
> To: "PROJ.4 and general Projections Discussions" <proj at lists.maptools.org>
> Subject: Re: [Proj] Stereo 1970 (EPSG 31700)
>
>> Hello Noel,
>>
>>> The least-squares adjustment of 7 transformation parameters
>>> (3 translations, 3 rotations and a scale) for an area as small as the
>>> country of Romania is badly ill-conditioned.
>>
>> Now I am confused.
>>
>> What you say is that it is a difficult statistical problem to derive a
>> 7-parameter transformation from observations in a small area. But if
>> someone I can trust has already derived a 7-parameter transformation,
>> then there is no numerical problems when using it, is there?
>>
>>> In these cases use a 3-parameter transformation instead.
>>> Or a Molodensky-Badekas transformation ... as has become common in
>>> some smaller countries (e.g. Luxembourg, Venezuela) ... if you can
>>> find (or derive) one.
>>
>> Well, but if a Molodensky-Badekas transformation is okay, then there
>> is always a 7-parameter transformation that gives the same result.
>> For all I know, the Romanian authorities may have started out by
>> deriving a Molodensky-Badekas transformation, and then decided to
>> publish it as the equivalent 7-parameter transformation.  If that is
>> what they did, there would be no reason to suspect that their
>> transformation is less accurate than the 3-parameter transformations,
>> is there?
>>
>> (An example: the datum shift EPSG:4829, "S-JTSK to ETRS89 (3)", is a
>> Molodensky-Badekas transformation, and the datum shift EPSG:4827,
>> "S-JTSK to ETRS89 (4)", is the same transformation expressed as a
>> 7-parameter transformation.)
>>
>> Best regards,
>>
>> Mikael Rittri
>> Carmenta AB
>> Sweden
>> www.carmenta.com
>>
>> -----Original Message-----
>> From: proj-bounces at lists.maptools.org
>> [mailto:proj-bounces at lists.maptools.org] On Behalf Of Noel Zinn (cc)
>> Sent: den 6 oktober 2010 16:47
>> To: PROJ.4 and general Projections Discussions
>> Subject: Re: [Proj] Stereo 1970 (EPSG 31700)
>>
>> Although the third of the three EPSG transformations cited herein has
>> the
>> (Romanian) National Agency for Cadastre and Land Registration as its
>> source (and is, therefore, authoritative) and it has more parameters
>> (7 versus 3 for the Shell and NIMA transformations), it should not be
>> regarded as more accurate.  The least-squares adjustment of 7
>> transformation parameters (3 translations, 3 rotations and a scale)
>> for an area as small as the country of Romania is badly
>> ill-conditioned.  The 7 parameters are highly correlated and,
>> therefore, are driven by small survey errors, not physical reality.
>> Simply stated, the translations pull in one direction and the
>> rotations in another.  Australia is big enough, but not Romania.  Not
>> even Germany.  In these cases use a 3-parameter transformation
>> instead.  Or a Molodensky-Badekas transformation ... as has become
>> common in some smaller countries (e.g. Luxembourg, Venezuela) ... if
>> you can find (or derive) one.
>>
>> Noel Zinn, Principal, Hydrometronics LLC
>> +1-832-539-1472 (office), +1-281-221-0051 (cell)
>> noel.zinn at hydrometronics.com (email)
>> http://www.hydrometronics.com (website)
>>
>> --------------------------------------------------
>> From: "Mikael Rittri" <Mikael.Rittri at carmenta.com>
>> Sent: Wednesday, October 06, 2010 2:38 AM
>> To: "PROJ.4 and general Projections Discussions"
>> <proj at lists.maptools.org>
>> Subject: Re: [Proj] Stereo 1970 (EPSG 31700)
>>
>>> Hello, Thibaut.
>>>
>>> For the record, www.epsg-registry.org gives three datum shifts from
>>> Pulkovo 1942(58) to WGS84 that are valid in Romania.
>>>
>>> Code Accuracy Source Proj.4 syntax
>>> 15496 10 meters Shell SIEP +towgs84=44.107,-116.147,-54.648
>>> 15497 7 meters NIMA TR8350.2 +towgs84=28,-121,-77
>>> 15995 3 meters www.ancpi.ro (*)
>>> +towgs84=2.329,-147.042,-92.08,-0.309,0.325,0.497,5.69
>>> (*) National Agency for Cadastre and Land Registration
>>>
>>> The third one seems to be more accurate than the two you mentioned.
>>> Of course, this may not matter if your main concern is to get ArcMap
>>> and Qgis to agree.
>>>
>>> Actually, it is possible that the Shell datum shift is more accurate
>>> than the NIMA one, since accuracies that EPSG quotes from different
>>> sources are not comparable.  I think accuracies from NIMA are usually
>>> one-sigma (about 67 percent confidence), but I don't know about Shell.
>>> But my impression is that the ANCPI datum shift is the most accurate.
>>>
>>> By the way, I noticed that the code EPSG:31700 is deprecated, because
>>> it used the the geodetic datum Dealul Piscului 1970, EPSG:6317.  This
>>> was deprecated in September 2008 because
>>>
>>>   "Datum does not exist but is an alias for S-42 in Romania" (Change
>>> ID 2008.011).
>>>
>>> And S-42 is (in this context) the same as Pulkovo 1942(58).
>>> So, EPSG:31700 has been replaced by EPSG:3844, "Pulkovo 1942(58) /
>>> Stereo70".
>>>
>>> Best regards,
>>>
>>> Mikael Rittri
>>> Carmenta AB
>>> Sweden
>>> www.carmenta.com
>>>
>>> ________________________________
>>>
>>> From: proj-bounces at lists.maptools.org
>>> [mailto:proj-bounces at lists.maptools.org] On Behalf Of Thibaut Gheysen
>>> Sent: den 4 oktober 2010 17:24
>>> To: PROJ.4 and general Projections Discussions
>>> Subject: Re: [Proj] Stereo 1970 (EPSG 31700)
>>>
>>>
>>> I checked in Arcmap, there are two "geographic coordinate system
>>> transformations" for stereo 1970 to wgs84. I used the first with
>>> those parameters : +towgs84=44.107,-116.147,-54.648,0,0,0,0. The
>>> suggested
>>> proj.4 definition correspond to the second "geographic coordinate
>>> system transformations" in ArcMap.
>>>
>>> All is fine now.
>>>
>>> Thanks for your help,
>>> Thibaut.
>>>
>>>
>>> 2010/10/4 Thibaut Gheysen <gheysen.t at gmail.com>
>>>
>>>
>>> Thanks Frank.
>>> I use QGis 1.5 from osgeo4w (I don't know exactly which version of
>>> Proj.4 is used).
>>> I just try with your suggested definition. It's better but I always
>>> have a deviation of about 5 meters.
>>>
>>> Thibaut.
>>>
>>>
>>> 2010/10/4 Frank Warmerdam <warmerdam at pobox.com>
>>>
>>>
>>> Thibaut Gheysen wrote:
>>> > Hi,
>>> >
>>> > I have a problem in QGIS with the Stereo 1970 coordinate system
>>> > (EPSG 31700).
>>> > I have a shapefile in the WGS 84 coordinate system (from a GPS
>>> > device).
>>> > I projected this shapefile in Stereo 1970 using ArcMap.
>>> > In ArcMap, the source and projected points are in the same place.
>>> > In Qgis, I have a deviation of about 100 meters when displaying
>>> > this 2 shapefiles and activating the "on the fly projection"
>>> > (stereo.jpg).
>>> > This deviation is the same that the one I obtain in ArcMap without
>>> > "geographic coordinate system transformations".
>>>
>>>
>>> Thibault,
>>>
>>> I'm not sure exactly what version of PROJ.4 you are using, but after
>>> my "datum shift identification upgrade" the EPSG dictionary for proj
>>> now has this as the suggested definition of 31700:
>>>
>>>
>>> +proj=sterea +lat_0=46 +lon_0=25 +k=0.99975 +x_0=500000 +y_0=500000
>>> +ellps=krass +towgs84=28,-121,-77,0,0,0,0 +units=m +no_defs
>>>
>>> Perhaps you could try with this?
>>>
>>> Best regards,
>>> --
>>> ---------------------------------------+-----------------------------
>>> ---------------------------------------+-
>>> ---------------------------------------+--------
>>> I set the clouds in motion - turn up   | Frank Warmerdam,
>>> warmerdam at pobox.com
>>> light and sound - activate the windows | http://pobox.com/~warmerdam
>>> <http://pobox.com/%7Ewarmerdam>
>>> and watch the world go round - Rush    | Geospatial Programmer for Rent
>>>
>>> _______________________________________________
>>> Proj mailing list
>>> Proj at lists.maptools.org
>>> http://lists.maptools.org/mailman/listinfo/proj
>>>
>>>
>>>
>>>
>>> _______________________________________________
>>> Proj mailing list
>>> Proj at lists.maptools.org
>>> http://lists.maptools.org/mailman/listinfo/proj
>>>
>>
>> _______________________________________________
>> Proj mailing list
>> Proj at lists.maptools.org
>> http://lists.maptools.org/mailman/listinfo/proj
>> _______________________________________________
>> Proj mailing list
>> Proj at lists.maptools.org
>> http://lists.maptools.org/mailman/listinfo/proj
>>
>
> _______________________________________________
> Proj mailing list
> Proj at lists.maptools.org
> http://lists.maptools.org/mailman/listinfo/proj
> _______________________________________________
> Proj mailing list
> Proj at lists.maptools.org
> http://lists.maptools.org/mailman/listinfo/proj
>