How to convert from one ellipsoid to another?

Question

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?

Answer 1

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.

Answer 2

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:

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.