RD-coordinates to WGS84

20 Ansichten (letzte 30 Tage)
Marijn
Marijn am 30 Mai 2012
Does anyone have a script to convert RD-coordinates (Rijksdriehoekstelsel, applied in The Netherlands) to WGS84? Thanks, Marijn

Akzeptierte Antwort

Marijn
Marijn am 17 Apr. 2013
Hi Roland,
My professor has given me this script, also available through modflow, so I think it should be no problem to upload it here.
Good luck. Cheers, Marijn
function [lamwgs,phiwgs,LL]=rd2wgs(x,y)
% RD2WGS: Converts Dutch rd-coordinates to GE wgs lat(Easting) long(Northing) coordinates % % USAGE: % [long,lat,LL]=rd2wgs(x,y) % % fortran 90 routine received from Peter Vermeulen % converted to Matlab by TO 090916 % % SEE ALSO: wgs2rd, kmlpath kmlpath2rd getDinoXSec % % TO 090916
if nargin<2, [lamwgs,phiwgs,LL]=selftest; return end
[phibes, lambes]=rd2bessel(x, y); [phiwgs, lamwgs]=bessel2wgs84(phibes,lambes);
N=length(x(:)); LL=cell(N,1); % allocate for i=1:length(phibes(:)) LL{i,1}=sprintf('N %.7g, E %.7g',phiwgs(i),lamwgs(i)); end
function [phi,lambda]=rd2bessel(x,y) %convert xy to Bessel
x0 = 1.55e5; y0 = 4.63e5; k=.9999079; bigr = 6382644.571; m = .003773953832; n = 1.00047585668; e = .08169683122;
lambda0 = pi * .029931327161111111; b0 = pi * .28956165138333334;
d__1 = x - x0; d__2 = y - y0;
r = sqrt(d__1 .* d__1 + d__2 .* d__2);
warning off sa = (x - x0) ./ r; ca = (y - y0) ./ r; warning on sa(r==0)=0.0; ca(r==0)=0.0;
psi = atan2(r, k * 2. * bigr) * 2.;
cpsi = cos(psi); spsi = sin(psi);
sb = ca * cos(b0) .* spsi + sin(b0) * cpsi; d__1 = sb; cb = sqrt(1. - d__1 .* d__1); b = acos(cb); sdl = sa .* spsi ./ cb; dl = asin(sdl); lambda = dl / n + lambda0; w = log(tan(b / 2.0 + pi / 4.0)); q = (w - m) / n; phiprime = atan(exp(q)) * 2. - pi / 2.;
for i = 1:4 dq = e / 2.0 * log((e * sin(phiprime) + 1.) ./ (1. - e * sin(phiprime))); phi = atan(exp(q + dq)) * 2. - pi / 2.; phiprime = phi; end
lambda = lambda / pi * 180.; phi = phi / pi * 180.;
function [phiwgs,lamwgs]=bessel2wgs84(phibes, lambes) % convert Bessel2 WGS84
a = 52.0; b = 5.0; c = -96.862; d = 11.714; e = 0.125; f = 1e-5; g = 0.329; h = 37.902; i = 14.667;
dphi = phibes - a; dlam = lambes - b; phicor = (c - dphi * d - dlam * e) * f; lamcor = (dphi * g - h - dlam * i) * f; phiwgs = phibes + phicor; lamwgs = lambes + lamcor;
function [long,lat,LL]=selftest
%Amersfoort en andere punten x = [155000; 244000; 93000; 98000; 177000]; y = [463000; 601000; 464000; 471000; 439000];
[long,lat,LL]=rd2wgs(x,y);

Weitere Antworten (0)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by