[Update] RD coördinaten omzetten naar WGS84 coördinaten

Het volgende stuk code hieronder leek me wel handig te delen, deels omdat ik het vast nog wel vaker nodig ga hebben maar ook omdat anderen er misschien ook wel baat bij hebben uiteraard. De code is geschreven in python maar elke gangbare programmeertaal zal zich er wel voor lenen om iets vergelijkbaars ervan te maken.

import math
class Converter:
    RDOriginX = 155E3
    RDOriginY = 463E3
    GpsOriginLat = 52.1551744
    GpsOriginLon = 5.38720621
    def __init__(self):
        self.lat = []
        self.lon = []
        for i in range(0,11):
            self.lat.insert(i, [])
            self.lon.insert(i, [])
            self.lon.insert(i+1,[])

        ## Latitude calculation
        self.lat[0] = [0,1,3235.65389]
        self.lat[1] = [2,0,-32.58297]
        self.lat[2] = [0,2,-0.2475]
        self.lat[3] = [2,1,-0.84978]
        self.lat[4] = [0,3,-0.0665]
        self.lat[5] = [2,2,-0.01709]
        self.lat[6] = [1,0,-0.00738]
        self.lat[7] = [4,0,0.0053]
        self.lat[8] = [2,3,-3.9E-4]
        self.lat[9] = [4,1,3.3E-4]
        self.lat[10] = [1,1,-1.2E-4]
        
        ## Longitude calculation
        self.lon[0] = [1,0,5260.52916]
        self.lon[1] = [1,1,105.94684]
        self.lon[2] = [1,2,2.45656]
        self.lon[3] = [3,0,-0.81885]
        self.lon[4] = [1,3,0.05594]
        self.lon[5] = [3,1,-0.05607]
        self.lon[6] = [0,1,0.01199]
        self.lon[7] = [3,2,-0.00256]
        self.lon[8] = [1,4,0.00128]
        self.lon[9] = [0,0,2.2E-4]
        self.lon[10] = [2,0,-2.2E-4]
        self.lon[11] = [5,0,2.6E-4]

    def toLat(self,rdX,rdY):
        a = 0
        dX = 1E-5 * (rdX - self.RDOriginX)
        dY = 1E-5 * (rdY - self.RDOriginY)

        for i in range(0,11):
            a = a + ( self.lat[i][2] * math.pow(dX, self.lat[i][0]) * math.pow(dY, self.lat[i][1]) )

        return round(self.GpsOriginLat + ( a / 3600 ), 9)

    def toLon(self,rdX,rdY):
        a = 0
        dX = 1E-5 * (rdX - self.RDOriginX)
        dY = 1E-5 * (rdY - self.RDOriginY)
        
        for i in range(0,12):
            a = a + ( self.lon[i][2] * math.pow(dX,self.lon[i][0]) * math.pow(dY, self.lon[i][1]) )

        return round(self.GpsOriginLon + ( a / 3600 ), 9)

    def main():
        converter = Converter()

        ## Test (Rotterdam Centraal)
        print(converter.toLat(91819, 437802))
        print(converter.toLon(91819, 437802))

if __name__ == "__main__":
    main()

Update

Iets nettere en uitgebreidere versie met twee weg conversie. Inclusief bronvermelding.

# Formules voor benadering zijn gebaseerd op http://www.dekoepel.nl/pdf/Transformatieformules.pdf
# Bovenstaande link werkt helaas niet meer, daar Stiching de Koepel opgeheven is. Backup link: http://media.thomasv.nl/2015/07/Transformatieformules.pdf

class RDWGSConverter:

    X0      = 155000
    Y0      = 463000
    phi0    = 52.15517440
    lam0    = 5.38720621

    def fromRdToWgs( self, coords ):

        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 * ( coords[0] - self.X0 )
        dY = 1E-5 * ( coords[1] - self.Y0 )
        
        phi = 0
        lam = 0

        for k in range(len(Kpq)):
            phi = phi + ( Kpq[k] * dX**Kp[k] * dY**Kq[k] )
        phi = self.phi0 + phi / 3600

        for l in range(len(Lpq)):
            lam = lam + ( Lpq[l] * dX**Lp[l] * dY**Lq[l] )
        lam = self.lam0 + lam / 3600

        return [phi,lam]

    def fromWgsToRd( self, coords ):
        
        Rp = [0,1,2,0,1,3,1,0,2]
        Rq = [1,1,1,3,0,1,3,2,3]
        Rpq = [190094.945,-11832.228,-114.221,-32.391,-0.705,-2.340,-0.608,-0.008,0.148]

        Sp = [1,0,2,1,3,0,2,1,0,1]
        Sq = [0,2,0,2,0,1,2,1,4,4]
        Spq = [309056.544,3638.893,73.077,-157.984,59.788,0.433,-6.439,-0.032,0.092,-0.054]

        dPhi = 0.36 * ( coords[0] - self.phi0 )
        dLam = 0.36 * ( coords[1] - self.lam0 )

        X = 0
        Y = 0

        for r in range( len( Rpq ) ):
            X = X + ( Rpq[r] * dPhi**Rp[r] * dLam**Rq[r] ) 
        X = self.X0 + X

        for s in range( len( Spq ) ):
            Y = Y + ( Spq[s] * dPhi**Sp[s] * dLam**Sq[s] )
        Y = self.Y0 + Y

        return [X,Y]

def main():
    # Rotterdam Centraal Station
    coords = [91819, 437802]
    conv = RDWGSConverter()
    wgsCoords = conv.fromRdToWgs( coords )
    newCoords = conv.fromWgsToRd( wgsCoords )
    # Laat de coordinaten voor conversie en de coordinaten na dubbele conversie op het scherm zien
    # om te bewijzen dat het een benadering betreft. Tot slot worden ook nog de WGS coordinaten weergegeven.
    print( coords, newCoords, wgsCoords )

if __name__ == "__main__":
    main()