Julia function to convert RD coordinates (rijksdriehoekscoordinaten) into WGS 84 coordinates.

rd_to_wss84.jl

"""
Converts RD coordinates (rijksdriehoekscoordinaten) into WGS84 coordinates

More about these systems:
- [RD](https://nl.wikipedia.org/wiki/Rijksdriehoeksco%C3%B6rdinaten)
- [WGS84](https://en.wikipedia.org/wiki/World_Geodetic_System#WGS84)

Implementation based on [this python script](https://github.com/thomasvnl/rd-to-wgs84/blob/master/conversions.py)
"""
function rd_to_wgs88(x::Number, y::Number):: Tuple{Number, Number}
	x0 = 155000
    y0 = 463000
    phi0 = 52.15517440
    lam0 = 5.38720621

	# Coefficients or the conversion from RD to WGS84
    Kp = [0, 2, 0, 2, 0, 2, 1, 4, 2, 4, 1]
    Kq = [1, 0, 2, 1, 3, 2, 0, 0, 3, 1, 1]
    Kpq = [3235.65389, -32.58297, -0.24750, -0.84978, -0.06550, -0.01709, -0.00738, 0.00530, -0.00039,
           0.00033, -0.00012]

    Lp = [1, 1, 1, 3, 1, 3, 0, 3, 1, 0, 2, 5]
    Lq = [0, 1, 2, 0, 3, 1, 1, 2, 4, 2, 0, 0]
    Lpq = [5260.52916, 105.94684, 2.45656, -0.81885, 0.05594, -0.05607, 0.01199, -0.00256, 0.00128, 0.00022,
           -0.00022, 0.00026]


	dx = 1E-5 * (x - x0)
	dy = 1E-5 * (y - y0)
	latitude = phi0 + sum(Kpq .* dx .^ Kp .* dy .^ Kq) / 3600
	longitude = lam0 + sum(Lpq .* dx .^ Lp .* dy .^ Lq) / 3600

	latitude, longitude
end