Converting RD coordinates to WGS84 in PL/SQL

I am currently working on an application which uses the BAG (Basis Administratie Gebouwen) of the Dutch Kadaster (Dutch Cadastre, Land Registry and national mapping agency). When I was trying to integrate Google Maps in my OBIEE report I noticed the the chart would not work. Looking into the problem I noticed the coordinates looked a bit different than normal. It turns out that the Kadaster uses a coordinate system specific to the Netherlands. Searching on the internet I found several websites capable of converting the coordinates, but I wanted to add them to my table as a normal longitude-latitude pair. I eventually found an algorithm in Python (Dutch). I converted the code to PL/SQL and put it in a package. The functions take 1 parameter which contains the coordinates from the BAG, for example: ‘163186.737 425229.796 0.0’.

Many thanks to Thomas V. for his original post.

CREATE OR REPLACE PACKAGE bag_rd_conv
AS
   FUNCTION conv_long(p_pos VARCHAR2)
      RETURN NUMBER
      DETERMINISTIC;

   FUNCTION conv_lat(p_pos VARCHAR2)
      RETURN NUMBER
      DETERMINISTIC;
END;
/

CREATE OR REPLACE PACKAGE BODY bag_rd_conv
AS
   -- Convert RD (Rijksdriehoeksysteem) coordinate to WGS84
   -- Latitude
   FUNCTION conv_lat(p_pos VARCHAR2)
      RETURN NUMBER
      DETERMINISTIC
   AS
      rdx   NUMBER;
      rdy   NUMBER;

      x0    NUMBER := 155000;
      y0    NUMBER := 463000;
      phi0  NUMBER := 52.15517440;
      lam0  NUMBER := 5.38720621;

      phi   NUMBER;
      lam   NUMBER;

      dx    NUMBER;
      dy    NUMBER;

      TYPE kt IS TABLE OF NUMBER;

      kp    kt := kt(0, 2, 0, 2, 0, 2, 1, 4, 2, 4, 10);
      kq    kt := kt(1, 0, 2, 1, 3, 2, 0, 0, 3, 1, 1);
      kpq   kt := kt(3235.65389, -32.58297, -0.24750, -0.84978, -0.06550, -0.01709, -0.00738, 0.00530, -0.00039, 0.00033, -0.00012);
   BEGIN
      rdx   := TO_NUMBER(SUBSTR(p_pos, 1, INSTR(p_pos, ' ') - 1));
      rdy   := TO_NUMBER(SUBSTR(p_pos, INSTR(p_pos, ' ') + 1, INSTR(p_pos, ' ', 1, 2) - INSTR(p_pos, ' ') - 1));

      dx    := 1E-5 * (rdx - x0);
      dy    := 1E-5 * (rdy - y0);

      phi   := 0;
      lam   := 0;

      FOR k IN 1 .. kpq.COUNT LOOP
         phi   := phi + (kpq(k) * dx ** kp(k) * dy ** kq(k));

      END LOOP;

      phi   := phi0 + phi / 3600;
      RETURN (ROUND(phi, 5));
   EXCEPTION
      WHEN OTHERS THEN
         RETURN (NULL);
   END;

   -- Convert RD (Rijksdriehoeksysteem) coordinate to WGS84
   -- Longitude
   FUNCTION conv_long(p_pos VARCHAR2)
      RETURN NUMBER
      DETERMINISTIC
   AS
      rdx   NUMBER;
      rdy   NUMBER;

      x0    NUMBER := 155000;
      y0    NUMBER := 463000;
      phi0  NUMBER := 52.15517440;
      lam0  NUMBER := 5.38720621;

      phi   NUMBER;
      lam   NUMBER;

      dx    NUMBER;
      dy    NUMBER;

      TYPE kt IS TABLE OF NUMBER;

      kp    kt := kt(0, 2, 0, 2, 0, 2, 1, 4, 2, 4, 10);
      kq    kt := kt(1, 0, 2, 1, 3, 2, 0, 0, 3, 1, 1);

      kpq   kt := kt(3235.65389, -32.58297, -0.24750, -0.84978, -0.06550, -0.01709, -0.00738, 0.00530, -0.00039, 0.00033, -0.00012);

      lp    kt := kt(1, 1, 1, 3, 1, 3, 0, 3, 1, 0, 2, 5);
      lq    kt := kt(0, 1, 2, 0, 3, 1, 1, 2, 4, 2, 0, 0);
      lpq   kt := kt(5260.52916, 105.94684, 2.45656, -0.81885, 0.05594, -0.05607, 0.01199, -0.00256, 0.00128, 0.00022, -0.00022, 0.00026);
   BEGIN
      rdx   := TO_NUMBER(SUBSTR(p_pos, 1, INSTR(p_pos, ' ') - 1));
      rdy   := TO_NUMBER(SUBSTR(p_pos, INSTR(p_pos, ' ') + 1, INSTR(p_pos, ' ', 1, 2) - INSTR(p_pos, ' ') - 1));

      dx    := 1E-5 * (rdx - x0);
      dy    := 1E-5 * (rdy - y0);

      lam   := 0;

      FOR l IN 1 .. lpq.COUNT LOOP
         lam   := lam + (lpq(l) * dx ** lp(l) * dy ** lq(l));

      END LOOP;

      lam   := lam0 + lam / 3600;

      RETURN (ROUND(lam, 5));
   EXCEPTION
      WHEN OTHERS THEN
         RETURN (NULL);
   END;

END;
/