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; /