Accuracy of EPSG_AMERSFOORT_RD_NEW ⇒ EPSG_WGS_84 conversion

[ameenross]

I must say first of all I'm quite impressed by this library. Great job on it.

Now, I've noticed that I'm not the first one who has a question about the accuracy of the conversion, and I appreciate how hard of a problem this is. So with that out of the way.

I've got the following code:

$point = ProjectedPoint::createFromEastingNorthing(
    new Metre(86058.2205),
    new Metre(445832.8025),
    Projected::fromSRID(Projected::EPSG_AMERSFOORT_RD_NEW)
);

$gpsPoint = $point->convert(Geographic2D::fromSRID(Geographic2D::EPSG_WGS_84));

// Result:
(51.997030373929, 4.383346111941)

// Expected (taken from Google maps):
(51.996591, 4.3831851)

This turns out to be almost 100m off, about 60/70m to the east and 10/20m to the north. The location comes from official government sources, but they only provide the coordinates in EPSG:28992. Though I do expect the data to be quite accurate.

Edit: made a mistake with google map coords.

[ameenross]

Made a couple of quick&dirty calculations and it's about 50m off (11m east and 49m north).

[ameenross]

epsg.io gives these coordinates: 51.99659194444445 4.383328611111111

If I check that on Google Maps it's still not perfect, but it is more accurate. It's about 10m east.

[ameenross]

For reference, I looked at another location (Borkumweg 25, 9979 XJ Eemshaven). It's at the other end of the country.

X: 250074.42 Y: 607478.61

EPSG_WGS_84: 53.445467237338 6.8181373251081

Found via MAPS: 53.4449774 6.8177501

According to epsg.io: 53.44503777777778 6.818100555555556

In both cases discrepancy is about 25m east and epsg.io is about 11m north, this library is 55m north.

[ameenross]

That seems to suggest there's a slight "drift" for lack of a better word. Another location, now we have the south-west corner and the north-east corner of the country.

Tol 1, 4529 PT Eede

X: 17612.8175 Y: 365058.794

EPSG_WGS_84: 51.258629611511 3.4189012309712

Found via MAPS: 51.2582437 3.418674

According to epsg.io: 51.25818694444445 3.418891111111111

Discrepancy: in both cases 15m east, epsg is about 6m south, this library is 43m north.

[dvdoug]

Hi @ameenross

Thanks for this bug report, I can confirm this isn't working quite right. According to PCTrans 5.2 version 20210630 which is the tool recommended by NSGI for desktop use, the correct result for your original point is Latitude 51.996591992 Longitude 4.383328547

(edit - technically, the PCTrans numbers are to ETRS89)

[ameenross]

Those are almost the exact same coordinates as given by epsg.io. Only a 7mm difference. Those coordinates are still not perfect, but much closer to reality (that is, assuming Google Maps is as close to reality as practically possible — which I assume to be the case since they're used for navigation).

Is there anything I can do to get a higher accuracy? Or do we need to fix the library? I'm even considering adding a transform to the result, with a simple linear function on both axes. Is that generally an acceptable idea?

[ameenross]

Only noticed the update after commenting. It's pretty much identical to epsg.io and PCTrans now:

51.996591926354 4.3833285740489

Thanks a lot for that!

Still would like to know your opinion on this (in this case just to get within a couple meters accuracy):

I'm even considering adding a transform to the result, with a simple linear function on both axes. Is that generally an acceptable idea?

[dvdoug]

On master which I don't think has any relevant changes from the most recent release, I get very slightly different numbers: 51.996591989427, 4.3833285482488 which is a 0.0003m difference from PCTrans.

Are you using the Europe datapack (php-coord/datapack-europe), or just the base library? The Amersfoort to ETRS89 transformation is one that does benefit from greater accuracy using a gridfile.

[ameenross]

I have tried the datapack last friday and saw no difference. Trying again now, I do see a difference

no datapack	with datapack	difference
51.996591926354	51.996591989427	0.000000063073
4.3833285740489	4.3833285482488	0.000000258001

This is only a couple mm difference. I'm seeing about a 10m difference (to the east) with Google Maps though. Any thoughts on that? Is my assumption that Google Maps is accurate wrong? Or do they not use WGS_84?

[ameenross]

Never mind. I made a mistake in that the center of the map in Google Maps is not actually the center of the map. When I hide the left sidebar it is spot on. Also OSM with the same coordinates is spot on.

[dvdoug]

Awesome 👍

[ameenross]

The accuracy of that location is great, but further north-east the accuracy drops somewhat:

Borkumweg 25, 9979 XJ Eemshaven (NE)

	Easting/northing	PHPCoord w/ datapack	Google Maps	difference	in meter
Y	607478.61	53.445037557349	53.44504906505906	0.00001150771006	~1
X	250074.42	6.8181002319938	6.818043901005739	0.000056330988061	~4

Tol 1, 4529 PT Eede (SW)

	Easting/northing	PHPCoord w/ datapack	Google Maps	difference	in meter
Y	365058.794	51.258187380977	51.25818486313902	0.0000251783798	~0.2
X	17612.8175	3.4188942815227	3.4188915993136852	0.0000026822090148	~0.3

Still rather accurate though, I'm impressed! Good enough for my use.

[dvdoug]

Happy to check into those, do you have the projection eastings and northings for those addresses?

[ameenross]

I'll edit the comment and add them. I took the centerpoint of the buildings using a bounding box (based on easting northing). I took the Google maps coordinates by hand... should be accurate to about a meter (but definitely not sub-meter precision). I suspect that the SW point is perfectly accurate because of that. Only NE is slightly off.

[ameenross]

Looking again, I actually suspect that in this case the Google Maps data is inaccurate (since the building polygon looks completely different from the satellite image)... let's look at another location/building close by.

Some building near Borkumkade 1, 9979 XX Eemshaven

	Easting/northing	PHPCoord w/ datapack	Google Maps	difference	in meter
Y	608823.625	53.456987023216	53.457000404196215	0.000013380980215	~1.5
X	250810.449	6.8295801118956	6.829541596650915	0.000038515244685	~2.5

Still a bit off. However looking at the maps widget from ESRG with the coordinates from PHPCoord is absolutely dead on. So Google maps probably a bit off after all (even though the building polygon at least matches satellite image and everything else — as opposed to Borkumweg 25)?

[ameenross]

Or could it be that Google Maps is more accurate because it's constantly calibrated by GPS data of its users, and OSM is based on converted easting/northing government data? This is where I'm really losing it.

[dvdoug]

I'll check the others later, but PHPCoord and PCTrans are matching to within a mm for Borkumweg 25 on the conversion from RD New -> Amersfoort -> ETRS89 (it wouldn't have before your bug report, so thanks for that)

The bit that gets complicated after that is the transformation from ETRS89 to WGS84. Although they're often treated interchangeably, that's only true with within ~1m of accuracy. They were the same in 1989, but the effects of continental drift mean they get a few cm further apart every year. To do a fully accurate ETRS89->WGS84 conversion involves additional time-based calculations to compensate for this (I'll add an example later).

So it wouldn't surprise me if Google Maps and OSM are both a combination (in different ways) of recent WGS84 (accurate to ~1cm), slightly older WGS84 (accurate to ~1cm at the time, ~10-20cm now) and government data that's actually ETRS89 but is treated as WGS84 (accurate ~1m). This blog post from OSM is informative: https://www.openstreetmap.org/user/StephaneP/diary/390290

[dvdoug]

Just to check that you see what I see, for "Some building near Borkumkade 1", using the above calculated PHPCoord w/datapack coordinates (53.456987023216, 6.8295801118956). That's the centre-ish of the building, but clearly not the actual centre. Good enough for most purposes though.

Using a more sophisticated conversion from ETRS89 to WGS84, we can get the centre of the roofline, but still slightly off-centre in the other direction. I suspect that's partially perspective (this is not a direct overhead) , partially the inherent uncertainness when obtaining coordinates to sub-metre accuracy in a moving object (a plane). To check, lets look at Bing Maps.

Using the exact same coordinates, it's still the centre of the roofline but slightly offset from centre but on the other side 🤣

[dvdoug]

Code for a more sophisticated conversion. If ~1m accuracy is good enough, don't bother implementing this!

// what you already have
$from = ProjectedPoint::createFromEastingNorthing(new Metre(250810.449), new Metre(608823.625), Projected::fromSRID(Projected::EPSG_AMERSFOORT_RD_NEW));
$etrs89 = Geographic2D::fromSRID(Geographic2D::EPSG_ETRS89);
$asETRS89 = $from->convert($etrs89);

//ETRS89 is a broad "ensemble" (see https://phpcoord.net/en/stable/builtin_crs_geographic2d.html#etrs89), from here it's necessary to pick one. Note that PHPCoord does not do any direct conversions from the generic ensemble to a more specific realisation as which one is applicable requires user-knowledge

//here I am constructing a 3D point. A 2D point should work but requires one-too-many conversions to be "chained" for the current limits (needs to be converted to a 3D point internally) . I'll look at tweaking the limit
$asETRF2014 = GeographicPoint::create(
    Geographic3D::fromSRID(Geographic3D::EPSG_ETRF2014),
    $asETRS89->getLatitude(),
    $asETRS89->getLongitude(),
    new Metre(0),
    new DateTime() // convert coordinates to "today" WGS84, different dates will give different results
);
$toCRS = Geographic2D::fromSRID(Geographic2D::EPSG_WGS_84_G1762); // note the specific WGS84 realisation
$asWGS84 = $asETRF2014->convert($toCRS);

Internally that would makes the following happen: ~ETRF2014 (2D)~ -> ETRF2014 (3D) -> ETRF2014 (geocentric) -> ITRF2014 (geocentric) -> WGS84 G1762 (geocentric) -> WGS84 G1762 (3D) -> WGS84 G1762 (2D)

[ameenross]

So it wouldn't surprise me if Google Maps and OSM are both a combination (in different ways) of recent WGS84 (accurate to ~1cm), slightly older WGS84 (accurate to ~1cm at the time, ~10-20cm now) and government data that's actually ETRS89 but is treated as WGS84 (accurate ~1m).

Yeah I suspected something like that would be the case. Interesting.

Just to check that you see what I see

That's the building all right, but I didn't wanna look at the satellite image (well I have, but not for pinpointing the center) exactly for the reasons you mentioned. I assumed/hoped that the polygon would not be based on a single satellite image.

Google Maps is a bit off. I added crosshairs to indicate the coords of PHPCoord (red) and the actual center of the object (green).

OSM is dead on the money.

[ameenross]

Good enough for most purposes though.

Yes it's good enough for my application as well. However still interested in knowing what the cause of that 3m discrepancy is. I guess I'll assume it's Google Maps.

[dvdoug]

I've done a cross-check of the date-based transformation from ETRS89 to WGS84 for that building against the EUREF Permanent GNSS Network http://www.epncb.oma.be/_productsservices/coord_trans/index.php assuming that WGS84 and ITRF are the same (they basically are). That website and PHPCoord are matching to ~0.0001m.

The "dumb" calculation where ETRS89 coordinates are simply treated as interchangeable with WGS84 give a point ~0.7663m away from that.

So 3m is ¯_(ツ)_/¯

[ameenross]

This has been quite an interesting learning experience. Thanks for your fix and hand-holding, very much appreciated.