[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()
  • Ad

    Er zit een foutje in de tweede versie in Kpq: -0.0039 moet zijn -0.00039 (-3.9E-4 in eerste versie). Verder prima code en netjes geschreven.

    • thomasvnl

      Bedankt, goed gezien. Zal het aanpassen!

  • rob

    Bedankt, ik had de code al een keer geschreven in Golang maar sinds mijn afscheid van Github kon ik er niet meer bij 🙂 Ik meen me wel te herinneren dat ik code had geschreven die ook buiten NL goed werkte (deze code is een benadering die op (volgens mij) locaties steeds meer richting de polen te grote fouten oplevert). Ik zal eens kijken of ik die code weer kan vinden dan deel ik ‘m op mijn bitbucket account. In ieder geval bedankt voor de code (ik heb er voor mijn doeleinden classes van gemaakt).

    Groet,
    Rob ( https://bitbucket.org/breinbaas/ (deze code staat onder bbtools/gis )

    ps. leuke site, met interessante topics!

    • thomasvnl

      Hej Rob,

      Graag gedaan en bedankt voor je comment ;)! Goed om te zien dat je er wat aan hebt. Heb je ook het commentaar van Ad gezien? Hij geeft aan dat er een fout in Kpq zit. Misschien dat dit die afwijking kan verklaren (ben er zelf al een tijd niet mee bezig geweest).

    • Paul Sohier

      RD is alleen te gebruiken op het Nederlands continentaal plat, buiten dit gebied zal er bij het omrekenen altijd verschillen zijn, van een aantal centimeters tot honderden meters.

      Dit document [1] bevat uitleg wat de oorzaak hiervan is, en wat de geldigheid is van RD (Pagina 33).
      Om correcte transformaties te doen zal je altijd een correctiegrid moeten gebruiken, dit correctiegrid is echter niet beschikbaar buiten het geldigheids gebied.

      [1] http://www.kadaster.nl/web/artikel/download/De-geodetische-referentiestelsels-van-Nederland-1.htm

  • P. Altenburg

    Hey Thomas, ziet er goed uit deze code, lekker duidelijk 🙂
    Alleen die link naar de referentie is dood, gezien stichting De Koepel blijkbaar in liquidatie is…
    Dus ik vroeg me af:
    Die officiele methode van RWS: RDNaptrans08 (oid) is helemaal op de cm nauwkeurig vanwege dat correctiegrid natuurlijk, maar hoe nauwkeurig is bovenstaande methode? Iets van 1 meter ofzo? (Wat ook prima is voor de meeste toepassingen)

    • P. Altenburg

      K kan hem zelf beantwoorden 🙂
      Ik vond die transformatie.pdf nog, en heb die twee voorbeelden van de Martinitoren en de Westertoren daaruit even ingevoerd in de RDNaptrans08 functie van PCTrans 4.2.10 (van RWS), en beide punten zijn binnen de 50cm nauwkeurig.

      • thomasvnl

        Bedankt voor het beantwoorden. Jammer dat stichting De Koepel overleden is. Ik hou me de laatste tijd niet echt meer bezig met deze formules dus ik had je vraag al niet kunnen beantwoorden op dit moment.

        Ik moet echt eens zorgen dat ik updates krijg wanneer iemand hier wat plaatst, mijn excuses daarom ook dat mijn antwoord op zich liet wachten. Vinkje stond uit, verklaard een hoop.

  • Vincent

    Hi Thomas,

    Ik weet niet of je deze berichten nog leest, maar ik wilde even bedanken. Ik heb me rot gezocht naar de formule om van GPS Rijksdriehoekscoördinaten te maken en het officiële bestandje op de site van het kadaster vond ik redelijk onoverzichtelijk. Op basis van jouw Python code heb ik in Stata een programma kunnen maken om al mijn data te kunnen omzetten (160.000+ observaties, dus niet echt met de hand te doen).

    Mocht iemand in te toekomst nog geïnteresseerd zijn naar de Stata code, deze staat op http://pastebin.com/WtfuETqy. Het zal niet de beste of meest efficiënte code zijn, maar het werkt.

    Groet,
    Vincent

    • thomasvnl

      Graag gedaan. Altijd mooi om te horen dat iemand er wat aan heeft en bedankt voor het delen van je code, hopelijk hebben anderen daar ook nog wat aan!

    • Kees

      Ha Vincent, ik was op zoek naar Stata code die dit kon doen. Kan ik jouw bestandje als ado-file opslaan en gebruiken? Zou je me kunnen zeggen hoe ik dat moet doen?
      Bedankt, Kees