6

Introduction

I have the following problem: I want to convert from WGS84 to "Anguilla 1957 / British West Indies Grid" (EPSG 2000). This seems to be straight-forward, as the transformation is a simple Transverse Mercator and the projection uses the Clarke 1880 Ellipsoid without any datum shift. My result differs significantly (some 400m) from the one i.e. epsg.io is returning.

Problem

The reason seems to be the datum transformation chain. Since the code is implemented such that it can convert from any cs to any cs, it also in this case does

  1. WGS84 geographic -> WGS84 Geocentric (using the WGS ellipsoid)
  2. WGS84 Geocentric -> Clarke 1880 Geocentric (nothing to do here)
  3. Clarke 1880 Geocentric -> Clarke 1880 Geographic (now using the Clarke Ellipsoid)

The last step is giving me a headache now, because the output of step 3 is not equal to the input at step 1, while the same conversion returns the exact input on reference pages I've checked.

Question

Should and can a "Geographic Ellipsoid 1->Geocentric->Geographic Ellipsoid 2" Conversion change latitude and/or longitude?

Update

Here's the steps I take: I basically want to convert from WGS84 to Anguilla, British West Indies Grid (EPSG 2000). According to http://epsg.io/2000 this uses Datum 4600. This datum consists of the Clarke 1880 RGS ellipsoid 7012 and no datum transformation.

I'm only considering the conversion from WGS84 to Datum 4600 here, because the rest is not relevant here.

I'm getting -63.000000000000014, 18.001851635270203, -78.860902368091047, which is the exact same result as mkennedy below using a 0 0 0 transformation. However, If I input the same coordinates at http://epsg.io/4600/map, The output is precisely -63, 18. That's where the confusion comes from. May I assume that the web site is wrong and whatever coordinate system library they're using is incorrectly assuming that a conversion whose datum parameters are 0 can be skipped altogether?

  • So you've implemented the transformation yourself? Not using some extant tool (like cs2cs) or library? – BradHards Mar 4 '14 at 10:07
  • Yes, I'm basically writing a transformation library. So far I've got a series of test cases which return the correct result, and also this returns (almost) the same values as the reference website when I'm directly feeding WGS lat/lon to the projection. – PMF Mar 4 '14 at 10:19
  • You know that there are well tested ones already (see metacrs)? – BradHards Mar 4 '14 at 10:56
  • Yes, I do, and I'm partially using existing code, but there are reasons I need to do some things myself. I'm extending a library that has been used for quite some time now and stumbled over this problem while performing some tests. – PMF Mar 4 '14 at 11:40
  • 1
    I think epsg.io is not performing the ellipsoid change if no transformation is chosen. It's an equally valid choice--both results are WRONG. In some ways the epsg.io result is better. It's more clear that the result isn't correct. – mkennedy Mar 6 '14 at 14:46
3

Note: simplified explanation

The relationship between two datums (geographic coordinate reference systems) is not just that they may have different ellipsoids. Differing ellipsoids mean different sizes and shapes. If you treat both datums as being earth-centered (center at 0,0,0 in 3D Cartesian coordinates), and convert between them, you may see latitude and height differences.

However, particularly older GeoCRS aren't earth-centered. It's easiest for me to think of them as attached at the surface and their centers as offset from the earth's center. See this diagram as an example.

When converting between two GeoCRS, you also have to take into account the offsets/rotations/scaling. The EPSG online registry lists one transformation between Anguilla 1957 and WGS 1984. It's applied as arc-second shifts:

latitude: -18"  
longitude: +4.4"

This is from Anguilla 1957 to WGS 1984.

One place to find more information on transformation including the equations and worked examples is EPSG's Guidance Note 7-2.

Edit: Including Esri projection engine test values

From WGS84 to Anguilla 1957 using the EPSG:1447 transformation
Input: 18.0N 63.0W
Output: 18.005N 63.0012222222222W

From Anguilla 1957 to the ProjCRS EPSG:2000:
Output: 294014.7725 1990652.6476

If I use the original coordinates instead (no transformation):
Output: 294141.1785 1990098.7974

If I set up a 0,0,0 (same as +wgs84=0,0,0,0,0,0,0) transformation between Anguilla 1957 and WGS84:
Output: 18.00185163527N 63.0W -78.86 m
and this value to EPSG:2000:
Output: 294142.2835 1990303.6444

The first and last test cases differ because the first just changes the CRS info on the data--literally nothing happens. On the last, because we're converting to/from XYZ space, the differences in flattening and size between the two ellipsoids are handled. If you can try to picture it, one ellipsoid is within the other. Even though their semimajor radii are different, there's no way to model it, so the longitude is unchanged. The semiminor axes (or flattening) are also different but the differences it reflected in a different latitude and new ellipsoid height value.

The Y/Northing difference on the no-transformation example may be to different Clarke 1880 definitions. EPSG has Clarke 1880 (RGS). I don't remember which version PROJ4 uses for clrk80.

Disclosure: I'm on the subcommittee that maintains the EPSG registry.

  • Thanks for this info, I was actually trying to follow that guidance note you mention. I'll add more details about my computation parameters tomorrow when I'm back in the office. It's soon bedtime here. – PMF Mar 4 '14 at 20:13
  • I've updated my question. Seems I'm getting the same result as you using the 0,0,0 transformation. Question is why epsg.io doesn't. How do I know I have to use the epsg:1447 transformation? Neither the gml nor the wkt definitions I've seen mention it. – PMF Mar 5 '14 at 8:16
  • When I checked epsg.io, and selected 4600, there was no default transformation selected, but I could choose to use tfm:1447. – mkennedy Mar 5 '14 at 19:41
  • Do you have test points to convert from WGS84->Clarke 1880? I guess that would confirm our assumption that the epsg.io implementation is broken (our mine and esri...). – PMF Mar 6 '14 at 9:30
  • I added them to my example. – AndreJ Mar 6 '14 at 9:58
3

I ran some tests with sample coordinates -63E 18N with GDAL cs2cs:

cs2cs +init=epsg:4326 +to +init=epsg:4600 Anguilla.txt >out.txt
cs2cs +init=epsg:4326 +to +proj=latlon +ellps=clrk80 +no_defs  Anguilla.txt >>out.txt
cs2cs +init=epsg:4326 +to +proj=latlon +ellps=clrk80 +no_defs +towgs84=0,0,0,0,0,0,0  Anguilla.txt >>out.txt
cs2cs +proj=latlong +datum=WGS84 +to +proj=latlon +ellps=clrk80 +no_defs +towgs84=0,0,0,0,0,0,0  Anguilla.txt >>out.txt
cs2cs +init=epsg:4326 +to +init=epsg:2000 Anguilla.txt >>out.txt
cs2cs +init=epsg:4326 +to +proj=tmerc +lat_0=0 +lon_0=-62 +k=0.9995000000000001 +x_0=400000 +y_0=0 +ellps=clrk80 +units=m +no_defs  Anguilla.txt >>out.txt
cs2cs +init=epsg:4326 +to +proj=tmerc +lat_0=0 +lon_0=-62 +k=0.9995000000000001 +x_0=400000 +y_0=0 +ellps=WGS84 +units=m +no_defs  Anguilla.txt >>out.txt
cs2cs +init=epsg:4326 +to +proj=tmerc +lat_0=0 +lon_0=-62 +k=0.9995000000000001 +x_0=400000 +y_0=0 +ellps=clrk80 +units=m +no_defs +towgs84=0,0,0,0,0,0,0  Anguilla.txt >>out.txt

The first line converts WGS84 lat/lon to Anguilla, the second and third to Anguilla tmerc according to EPSG:2000.

Then I changed the Ellipsoid to WGS84, and in the last step I explicitely added +towgs84=0,0,0,0,0,0,0 to the clarke ellipsoid, that is no datum shift should be made.

The results are:

63dW    18dN 0.000
63dW    18dN 0.000
63dW    18d0'6.664"N -78.870
63dW    18d0'6.664"N -78.870
294141.18   1990098.85 0.00
294141.18   1990098.85 0.00
294143.59   1990271.93 0.00
294142.28   1990303.65 -78.87

strangely, number 2 and 3 and the last three do not coincide. I don't know if and how the datum shift mkennedy mentioned is applied.

@PMF: which result do you get?

@mkennedy: which result do you get with ARCGIS?


EDIT

Since PROJ4 does not apply the official datum shift, I have set up my own towgs84 parameters:

cs2cs +init=epsg:4326 +to +proj=latlong +ellps=clrk80 +no_defs +towgs84=0,0,-364.8,0,0,4.4,0  -f %%f Anguilla.txt >out.txt

which fits to longitudes (rot z is the offical arc seconds shift), and has minimal distortions on the latitudes. Here you can see the WGS84 degree grid in blue and the anguilla degree grid in green:

enter image description here

The strange thing about the degree shift is that it does not obey the different ellipsoid flattenings, which towgs84 does. Both methods would coincide only on spheres.

A ntv2 grid with constant values would resemble the degree shift in a better way than towgs84 can supply.

The basic difference between towgs84 datum shift and degree or ntv2 datum shift is that the first is done on the ellipsoid center (obeying the ellipsoids flattening) and the others on the surface (regardless of shape and position of the ellipsoids).

But if the surveying authority (British OS in this case) say a degree shift is official, it has to be taken that way.

Anyway, EPSG Guidance Notes state a low accuracy for this kind of coordinate transformation, and it might apply well for a small island in the Carribean.

  • My answer edited with test points. – mkennedy Mar 4 '14 at 21:45
  • The difference you see is exactly the problem I was observing. Actually, the output of your first line is, according to @mkennedy wrong. And it seems that the error is widely spread. – PMF Mar 5 '14 at 8:22
  • GDAL only knows the towgs84 shift, that is moving the center of the ellisoid. Shifting in degrees is used rather seldom. Maybe a constant nadgrid could solve the problem. If you have found official reference coordinates online, please add them to your question to solve the issue. – AndreJ Mar 5 '14 at 8:31
  • I have official reference coordinates for other systems, where they fit very well. I will be looking for data on this system. – PMF Mar 5 '14 at 8:33
  • Unfortunately, I couldn't find anything reliable. Maybe somebody got some reference data for another coordinate system with the same behavior (0,0,0 shift)? NAD83 could do, for instance. I may be able to find that. – PMF Mar 5 '14 at 15:50

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.