#!/usr/bin/env python3
"""
RD New / EPSG:28992 to WGS84 latitude/longitude converter.
Purpose
-------
Practical conversion of Dutch Rijksdriehoeksstelsel coordinates
(RD New / EPSG:28992) to WGS84-style geographic coordinates
(latitude / longitude, EPSG:4326-style output).
Intended use
------------
This script is suitable for:
- operational checks;
- map previews;
- DATEX II point-coordinate snippets;
- GeoJSON point creation;
- sanity checks for Dutch road, asset or traffic datasets.
Accuracy indication
-------------------
This script uses a common polynomial approximation around the Amersfoort
reference point.
In normal Dutch RD coordinate ranges, the practical accuracy is generally
around metre-level or better when the input RD coordinates are correct.
Important nuance:
If the source positions were placed manually using high-quality tooling,
satellite imagery, road-reference data and other accurate sources, the source
placement may indeed be decimetre-level. However, this script itself should
not be documented as a certified decimetre-accurate transformation until it
has been validated against the official Dutch RDNAPTRANS2018 transformation.
Disclaimer
----------
This is NOT a certified geodetic transformation.
For official Dutch geodetic, cadastral, legal, engineering, survey-grade,
or high-accuracy production use, use the official RDNAPTRANS2018 workflow
and validate the complete CRS metadata.
Coordinate conventions
----------------------
Input:
RD X, RD Y in metres.
Output:
Latitude, longitude in decimal degrees.
GeoJSON warning:
GeoJSON uses [longitude, latitude], not [latitude, longitude].
DATEX II note:
DATEX II point coordinates are typically represented as latitude and
longitude fields, not Dutch RD X/Y.
"""
from __future__ import annotations
import argparse
import json
from dataclasses import dataclass
@dataclass(frozen=True)
class GeographicPoint:
"""Geographic coordinate in decimal degrees."""
latitude: float
longitude: float
def lat_lon_text(self, decimals: int = 8) -> str:
return f"{self.latitude:.{decimals}f}, {self.longitude:.{decimals}f}"
def lon_lat_text(self, decimals: int = 8) -> str:
return f"{self.longitude:.{decimals}f}, {self.latitude:.{decimals}f}"
def to_geojson_point(self, decimals: int = 8) -> dict:
return {
"type": "Point",
"coordinates": [
round(self.longitude, decimals),
round(self.latitude, decimals),
],
}
def to_datex_ii_snippet(self, decimals: int = 8) -> str:
return (
"\n"
f" {self.latitude:.{decimals}f}\n"
f" {self.longitude:.{decimals}f}\n"
""
)
def validate_rd(x: float, y: float) -> None:
"""
Validate whether values look like plausible Dutch RD New coordinates.
This is intentionally a sanity check, not a strict geodetic boundary test.
"""
if 3 <= x <= 8 and 50 <= y <= 54:
raise ValueError(
"Input looks like longitude/latitude, not RD X/Y. "
"RD coordinates are metre values, e.g. X=155000, Y=463000."
)
if not (0 <= x <= 300000):
raise ValueError(
f"RD X={x} is outside the expected practical range 0..300000."
)
if not (300000 <= y <= 620000):
raise ValueError(
f"RD Y={y} is outside the expected practical range 300000..620000."
)
def rd_to_wgs84(x: float, y: float, *, validate: bool = True) -> GeographicPoint:
"""
Convert RD New / EPSG:28992 X/Y to WGS84-style latitude/longitude.
Parameters
----------
x:
RD X coordinate in metres.
y:
RD Y coordinate in metres.
validate:
If True, performs a practical Netherlands range check.
Returns
-------
GeographicPoint
Latitude and longitude in decimal degrees.
Notes
-----
This is a practical polynomial approximation, not RDNAPTRANS2018.
"""
if validate:
validate_rd(x, y)
# Amersfoort reference point.
x0 = 155000.0
y0 = 463000.0
lat0 = 52.15517440
lon0 = 5.38720621
dx = (x - x0) * 1e-5
dy = (y - y0) * 1e-5
# Coefficients for latitude correction in arcseconds.
lat_terms = [
(0, 1, 3235.65389),
(2, 0, -32.58297),
(0, 2, -0.24750),
(2, 1, -0.84978),
(0, 3, -0.06550),
(2, 2, -0.01709),
(1, 0, -0.00738),
(4, 0, 0.00530),
(2, 3, -0.00039),
(4, 1, 0.00033),
(1, 1, -0.00012),
]
# Coefficients for longitude correction in arcseconds.
lon_terms = [
(1, 0, 5260.52916),
(1, 1, 105.94684),
(1, 2, 2.45656),
(3, 0, -0.81885),
(1, 3, 0.05594),
(3, 1, -0.05607),
(0, 1, 0.01199),
(3, 2, -0.00256),
(1, 4, 0.00128),
(0, 2, 0.00022),
(2, 0, -0.00022),
(5, 0, 0.00026),
]
lat_seconds = sum(coef * dx**p * dy**q for p, q, coef in lat_terms)
lon_seconds = sum(coef * dx**p * dy**q for p, q, coef in lon_terms)
latitude = lat0 + lat_seconds / 3600.0
longitude = lon0 + lon_seconds / 3600.0
return GeographicPoint(latitude=latitude, longitude=longitude)
def convert_one(x: float, y: float, decimals: int = 8) -> None:
point = rd_to_wgs84(x, y)
print("Input")
print("-----")
print(f"RD X: {x:.3f}")
print(f"RD Y: {y:.3f}")
print()
print("Latitude, Longitude")
print("-------------------")
print(point.lat_lon_text(decimals))
print()
print("Longitude, Latitude")
print("-------------------")
print(point.lon_lat_text(decimals))
print()
print("GeoJSON Point")
print("-------------")
print(json.dumps(point.to_geojson_point(decimals), indent=2))
print()
print("DATEX II-style snippet")
print("----------------------")
print(point.to_datex_ii_snippet(decimals))
print()
print("CRS explanation")
print("---------------")
print("Source: RD New / EPSG:28992")
print("Output for practical GNSS/web use: WGS84 / EPSG:4326-style latitude/longitude")
print("For official Dutch high-accuracy work: use RDNAPTRANS2018")
print()
print("Accuracy note")
print("-------------")
print(
"Practical approximation, generally around metre-level or better inside "
"normal Dutch RD ranges when input coordinates are correct. Not certified. "
"Do not claim survey/legal accuracy without RDNAPTRANS2018 validation."
)
def main() -> None:
parser = argparse.ArgumentParser(
description=(
"Convert Dutch RD New / EPSG:28992 coordinates to WGS84-style "
"latitude/longitude using a practical approximation."
)
)
parser.add_argument(
"x",
type=float,
help="RD X coordinate in metres, e.g. 155000",
)
parser.add_argument(
"y",
type=float,
help="RD Y coordinate in metres, e.g. 463000",
)
parser.add_argument(
"--decimals",
type=int,
default=8,
help="Number of decimal places for latitude/longitude output. Default: 8",
)
args = parser.parse_args()
convert_one(args.x, args.y, args.decimals)
if __name__ == "__main__":
main()