Spring naar bijdragen

Coördinaat transformatie's; RD <-> WGS

gpsgek

Door gpsgek
in Algemene vragen

 Share

Aanbevolen berichten

gpsgek

gpsgek

Geplaatst
Geplaatst (bewerkt)

Niet zozeer een newbie-vraag, geloof ik, maar toch. Kruimel, als je een betere plaats weet mag je het topic verplaatsen...

 

 

 

 

Wie kan mij helpen aan een - relatief eenvoudige - formule waarmee RD wordt omgerekend naar WGS84 of andersom? Geen tooltjes dus om dit te doen, die heb ik zat.

 

Ik wil graag de formule implementeren in een rekenprogramma (nee, geen excel...).

bewerkt door gpsgek
ReizenReizen

Biertjuh

Biertjuh

Geplaatst
Geplaatst (bewerkt)
Niet zozeer een newbie-vraag, geloof ik, maar toch. Kruimel, als je een betere plaats weet mag je het topic verplaatsen...

Wie kan mij helpen aan een - relatief eenvoudige - formule waarmee RD wordt omgerekend naar WGS84 of andersom? Geen tooltjes dus om dit te doen, die heb ik zat.

 

Ik wil graag de formule implementeren in een rekenprogramma (nee, geen excel...).

Kruimel zal zich zeer vereerd voelen dat hij het mag verplaatsen :xmas:

Eenvoudiger kan ik het niet zeggen. Kijk hier even :lol:

Ruud :rolleyes:

bewerkt door Biertjuh
deBruineBeren

deBruineBeren

Geplaatst
Geplaatst

Microsoft Visual Basic .Net :

 

Public Class Compute

Public Shared Function RD2WGS84(ByVal X As Integer, ByVal Y As Integer)

	Dim dX As Double
	Dim dY As Double
	Dim SomN As Double
	Dim SomE As Double
	Dim Latitude As String
	Dim LatitudeGraden As Integer
	Dim LatitudeMinuten As String
	Dim Longitude As Double
	Dim LongitudeGraden As Integer
	Dim LongitudeMinuten As String

	dX = (X - 155000) * 10 ^ -5
	dY = (Y - 463000) * 10 ^ -5

	SomN = (3235.65389 * dY) + (-32.58297 * dX ^ 2) + (-0.2475 * dY ^ 2) + (-0.84978 * dX ^ 2 * dY) + (-0.0655 * dY ^ 3) + (-0.01709 * dX ^ 2 * dY ^ 2) + (-0.00738 * dX) + (0.0053 * dX ^ 4) + (-0.00039 * dX ^ 2 * dY ^ 3) + (0.00033 * dX ^ 4 * dY) + (-0.00012 * dX * dY)
	SomE = (5260.52916 * dX) + (105.94684 * dX * dY) + (2.45656 * dX * dY ^ 2) + (-0.81885 * dX ^ 3) + (0.05594 * dX * dY ^ 3) + (-0.05607 * dX ^ 3 * dY) + (0.01199 * dY) + (-0.00256 * dX ^ 3 * dY ^ 2) + (0.00128 * dX * dY ^ 4) + (0.00022 * dY ^ 2) + (-0.00022 * dX ^ 2) + (0.00026 * dX ^ 5)

	Latitude = 52.15517 + (SomN / 3600)
	Longitude = 5.387206 + (SomE / 3600)

	LatitudeGraden = Int(Latitude)
	LongitudeGraden = Int(Longitude)

	LatitudeMinuten = (Latitude - LatitudeGraden) * 60.0
	LongitudeMinuten = (Longitude - LongitudeGraden) * 60.0

	Dim Coordinate As String() = New String() {LatitudeGraden, LatitudeMinuten, LongitudeGraden, LongitudeMinuten}
	Return Coordinate

End Function
End Class

Tornado Bram

Tornado Bram

Geplaatst
Geplaatst
Microsoft Visual Basic .Net :

 

Public Class Compute

Public Shared Function RD2WGS84(ByVal X As Integer, ByVal Y As Integer)

	Dim dX As Double
	Dim dY As Double
	Dim SomN As Double
	Dim SomE As Double
	Dim Latitude As String
	Dim LatitudeGraden As Integer
	Dim LatitudeMinuten As String
	Dim Longitude As Double
	Dim LongitudeGraden As Integer
	Dim LongitudeMinuten As String

	dX = (X - 155000) * 10 ^ -5
	dY = (Y - 463000) * 10 ^ -5

	SomN = (3235.65389 * dY) + (-32.58297 * dX ^ 2) + (-0.2475 * dY ^ 2) + (-0.84978 * dX ^ 2 * dY) + (-0.0655 * dY ^ 3) + (-0.01709 * dX ^ 2 * dY ^ 2) + (-0.00738 * dX) + (0.0053 * dX ^ 4) + (-0.00039 * dX ^ 2 * dY ^ 3) + (0.00033 * dX ^ 4 * dY) + (-0.00012 * dX * dY)
	SomE = (5260.52916 * dX) + (105.94684 * dX * dY) + (2.45656 * dX * dY ^ 2) + (-0.81885 * dX ^ 3) + (0.05594 * dX * dY ^ 3) + (-0.05607 * dX ^ 3 * dY) + (0.01199 * dY) + (-0.00256 * dX ^ 3 * dY ^ 2) + (0.00128 * dX * dY ^ 4) + (0.00022 * dY ^ 2) + (-0.00022 * dX ^ 2) + (0.00026 * dX ^ 5)

	Latitude = 52.15517 + (SomN / 3600)
	Longitude = 5.387206 + (SomE / 3600)

	LatitudeGraden = Int(Latitude)
	LongitudeGraden = Int(Longitude)

	LatitudeMinuten = (Latitude - LatitudeGraden) * 60.0
	LongitudeMinuten = (Longitude - LongitudeGraden) * 60.0

	Dim Coordinate As String() = New String() {LatitudeGraden, LatitudeMinuten, LongitudeGraden, LongitudeMinuten}
	Return Coordinate

End Function
End Class

Ik weet het niet hoor of je dit op je rekenmachine kunt intoetsen? Volgens mij lukt dat niet?

Biertjuh

Biertjuh

Geplaatst
Geplaatst
Ik weet het niet hoor of je dit op je rekenmachine kunt intoetsen? Volgens mij lukt dat niet?

Werk jij nog met een telraam dan :xmas:

Ruud :rolleyes:

gpsgek

gpsgek

Geplaatst
  • Auteur
Geplaatst
Microsoft Visual Basic .Net :

Bedankt Bruine Beren, was precies wat ik zocht! :rolleyes:

 

Ik weet het niet hoor of je dit op je rekenmachine kunt intoetsen? Volgens mij lukt dat niet?

Ja hoor, dat gaat prima...

Biertjuh

Biertjuh

Geplaatst
Geplaatst

Microsoft Visual Basic .Net :

Bedankt Bruine Beren, was precies wat ik zocht! :xmas:

Ik weet het niet hoor of je dit op je rekenmachine kunt intoetsen? Volgens mij lukt dat niet?

Ja hoor, dat gaat prima...

Antwoord gegeven: Draadje kan gesloten worden.

Ruud :rolleyes:

gpsgek

gpsgek

Geplaatst
  • Auteur
Geplaatst

Microsoft Visual Basic .Net :

Bedankt Bruine Beren, was precies wat ik zocht! :xmas:

Ik weet het niet hoor of je dit op je rekenmachine kunt intoetsen? Volgens mij lukt dat niet?

Ja hoor, dat gaat prima...

Antwoord gegeven: Draadje kan gesloten worden.

Ruud :rolleyes:

Genoeg te doen voor Kruimel... :lol:

Desondanks was het draadje zeker nuttig :yes:

gpsgek

gpsgek

Geplaatst
  • Auteur
Geplaatst
Wil je ook nog van WGS84 naar RD ?

 

Want dan moet ik ff zoeken.

Ik weet dat ik het heb, in eMbedded Visual Basic.

Zelf heb ik het niet nodig, maar wellicht is het handig het voor de volledigheid aan het draadje toe te voegen. Voor mensen die later iets soortgelijks zoeken!

Prof. Y. Lupardi

Prof. Y. Lupardi

Geplaatst
Geplaatst

Hou wel effe in de gaten dat wat hier gebruikt wordt een polynome omzetting is.

Dit zal voor gewoon gebruik wel voldoende nauwkeurig zijn.

Als er behoefte is dan heb ik wel een echte omzetting via een Molodensky transformatie maar die staat in M-code voor Matlab. (M-code lijkt wat op een blabla-Basic)

gpsgek

gpsgek

Geplaatst
  • Auteur
Geplaatst (bewerkt)
Als er behoefte is dan heb ik wel een echte omzetting via een Molodensky transformatie maar die staat in M-code voor Matlab. (M-code lijkt wat op een blabla-Basic)

Das helemaal mooi :D Gpsgek is enthousiast Matlab gebruiker en heeft bovenstaande code ook daarin toegepast :P

bewerkt door gpsgek
Prof. Y. Lupardi

Prof. Y. Lupardi

Geplaatst
Geplaatst

Eerst van RD meters naar graden op de Bessel ellipsoïde:

 

 

% rd2bessel

%

% Omrekenen geo coordinaten (tov Bessel) naar RD-grid coordinaten

% Formules van Geodetic Datum of Australia website

% 1 december 2001

 

function geo = rd2bessel(rd);

 

RDx = rd(1);

RDy = rd(2);

 

deg = pi/180; % omrekenfactor naar radialen

 

% Bessel 1841 datum:

a = 6377397.155;

f = 1/299.1528128;

 

% RD grid parameters:

lab0 = 5.3876333*deg; % origin

Efalse = 155000; % false easting

K0 = 0.9999079; % scale factor

Nfalse = -5316592; % false northing

 

% omrekenen:

 

e2 = f*(2-f);

e = sqrt(e2);

 

a0 = 1 - 1/4*e^2 - 3/64*e^4 - 5/256*e^6;

a2 = 3/8*(e^2 + 1/4*e^4 + 15/128*e^6);

a4 = 15/256*(e^4 + 3/4*e^6);

a6 = 35/3072*e^6;

 

m = (RDy-Nfalse)/K0;

 

n = f/(2-f);

G = a*(1-n)*(1-n^2)*(1+9/4*n^2 + 225/64*n^4)*pi/180;

sig = m/G*pi/180;

 

% phic = phi-accent.

phic = sig +((3*n/2)-(27*n^3/32))*sin(2*sig) + ((21*n^2/16)-(55*n^4/32))*sin(4*sig)...

+(151*n^3/96)*sin(6*sig) + (1097*n^4/512)*sin(8*sig);

 

rho = a*(1-e^2)/(1-e^2*(sin(phic))^2)^(3/2);

nu = a/sqrt(1-e^2*(sin(phic))^2);

psi = nu/rho;

 

% Grid to Geographical:

 

t = tan(phic);

Ec = RDx - Efalse;

x = Ec/(K0*nu);

 

term1 = (t/(K0*rho))*(x*Ec/2);

term2 = (t/(K0*rho))*(Ec*x^3/24)*(-4*psi^2+9*psi*(1-t^2)+12*t^2);

term3 = (t/(K0*rho))*((Ec*x^5)/720)*(8*psi^4*(11-24*t^2)-12*psi^3*(21-71*t^2)+15*(psi^2*...

15-98*t^2+15*t^4)+180*psi*(5*t^2-3*t^4)+360*t^4);

term4 = (t/(K0*rho))*(Ec*x^7/40320)*(1385+3633*t^2+4095*t^4+1575*t^6);

 

phi = phic - term1 + term2 - term3 + term4;

 

term1 = x*sec(phi);

term2 = (x^3/6)*sec(phi)*(psi+2*t^2);

term3 = (x^5/120)*sec(phi)*(-4*psi^3*(1-6*t^2)+psi^2*(9-68*t^2)+72*psi*t^2+24*t^4);

term4 = (x^7/5040)*sec(phi)*(61+662*t^2+1320*t^4+720*t^6);

 

omega = term1-term2+term3-term4;

 

lab = lab0 + omega;

 

geo(1) = phi/deg; % lengte [graden]

geo(2) = lab/deg; % breedte [graden]

%

%-------------------------------------------------------------------

 

En dan van Bessel naar de WGS84 ellipsoïde:

 

 

% molodensky datum omzetting van Bessel1841 naar WGS84

% Formules volgens website Peter Dana.

% Gecontroleerd met rekenvoorbeeld en GPS.

% In/uitvoer in graden en meters

 

% verbeterd nav check met EPS: fout in berekening Rm gehaald.

% 8 feb 2001

% 1 dec 2001: hoogte niet *deg!

 

function [to_geo]=bessel2wgs(from_geo);

 

deg = pi/180; % omrekenfactor naar radialen

 

from_phi = from_geo(1)*deg; % from_hoogte boven evenaar

from_lab = from_geo(2)*deg; % from_oosterlengte

from_h = from_geo(3); % from_hoogte boven zeespiegel

 

% parameters voor WGS 84:

to_a = 6378137.000; % semi-major axis of ellipsoid

to_f = 1/298.257223563; % 1/flattening

 

% parameters voor Bessel 1841:

from_a = 6377397.155;

from_f = 1/299.1528128;

 

% correctiefactoren uitrekenen:

da = to_a - from_a;

df = to_f - from_f;

bda = 1-from_f;

from_es = 2*from_f - from_f^2;

 

dx = 593;

dy = 26;

dz = 478;

 

% Nu coordinaten omrekenen van Bessel 1841 naar WGS84

 

sphi = sin(from_phi);

cphi = cos(from_phi);

slab = sin(from_lab);

clab = cos(from_lab);

 

Rn = from_a/sqrt(1.0 - from_es*(sin(from_phi))^2);

 

% oud/fout: Rm = from_a*(1-from_es)/(1.0-from_es*(sin(from_phi)^2)^(3/2));

Rm = from_a*(1-from_es)/((1.0-from_es*(sin(from_phi)^2))^(3/2));

 

term1 = (((-dx*sphi*clab-dy*sphi*slab) + dz*cphi) + da*Rn*from_es*sphi*cphi/from_a);

term2 = df*(Rm/bda + Rn*bda)*sphi*cphi;

 

dphi = (term1 + term2)/(Rm+from_h);

dlab = (-dx*slab + dy*clab)/((Rn+from_h)*cphi);

dh = dx*cphi*clab + dy*cphi*slab + dz*sphi - da*from_a/Rn + df*bda*Rn*sphi*sphi;

 

to_phi = from_phi+dphi;

to_lab = from_lab + dlab;

to_h = from_h + dh;

 

to_phi = to_phi/deg;

to_lab = to_lab/deg;

to_h = from_h + dh;

 

to_geo = [to_phi,to_lab,to_h];

%

%-------------------------------------------

 

 

Zo doe je dat dus :D

LTeunis

LTeunis

Geplaatst
Geplaatst
Eerst van RD meters naar graden op de Bessel ellipsoïde:

 

 

% rd2bessel

%

% Omrekenen geo coordinaten (tov Bessel) naar RD-grid coordinaten

% Formules van Geodetic Datum of Australia website

% 1 december 2001

 

function geo = rd2bessel(rd);

 

RDx = rd(1);

RDy = rd(2);

 

deg = pi/180; % omrekenfactor naar radialen

 

% Bessel 1841 datum:

a = 6377397.155;

f = 1/299.1528128;

 

% RD grid parameters:

lab0 = 5.3876333*deg; % origin

Efalse = 155000; % false easting

K0 = 0.9999079; % scale factor

Nfalse = -5316592; % false northing

 

% omrekenen:

 

e2 = f*(2-f);

e = sqrt(e2);

 

a0 = 1 - 1/4*e^2 - 3/64*e^4 - 5/256*e^6;

a2 = 3/8*(e^2 + 1/4*e^4 + 15/128*e^6);

a4 = 15/256*(e^4 + 3/4*e^6);

a6 = 35/3072*e^6;

 

m = (RDy-Nfalse)/K0;

 

n = f/(2-f);

G = a*(1-n)*(1-n^2)*(1+9/4*n^2 + 225/64*n^4)*pi/180;

sig = m/G*pi/180;

 

% phic = phi-accent.

phic = sig +((3*n/2)-(27*n^3/32))*sin(2*sig) + ((21*n^2/16)-(55*n^4/32))*sin(4*sig)...

+(151*n^3/96)*sin(6*sig) + (1097*n^4/512)*sin(8*sig);

 

rho = a*(1-e^2)/(1-e^2*(sin(phic))^2)^(3/2);

nu = a/sqrt(1-e^2*(sin(phic))^2);

psi = nu/rho;

 

% Grid to Geographical:

 

t = tan(phic);

Ec = RDx - Efalse;

x = Ec/(K0*nu);

 

term1 = (t/(K0*rho))*(x*Ec/2);

term2 = (t/(K0*rho))*(Ec*x^3/24)*(-4*psi^2+9*psi*(1-t^2)+12*t^2);

term3 = (t/(K0*rho))*((Ec*x^5)/720)*(8*psi^4*(11-24*t^2)-12*psi^3*(21-71*t^2)+15*(psi^2*...

15-98*t^2+15*t^4)+180*psi*(5*t^2-3*t^4)+360*t^4);

term4 = (t/(K0*rho))*(Ec*x^7/40320)*(1385+3633*t^2+4095*t^4+1575*t^6);

 

phi = phic - term1 + term2 - term3 + term4;

 

term1 = x*sec(phi);

term2 = (x^3/6)*sec(phi)*(psi+2*t^2);

term3 = (x^5/120)*sec(phi)*(-4*psi^3*(1-6*t^2)+psi^2*(9-68*t^2)+72*psi*t^2+24*t^4);

term4 = (x^7/5040)*sec(phi)*(61+662*t^2+1320*t^4+720*t^6);

 

omega = term1-term2+term3-term4;

 

lab = lab0 + omega;

 

geo(1) = phi/deg; % lengte [graden]

geo(2) = lab/deg; % breedte [graden]

%

%-------------------------------------------------------------------

 

En dan van Bessel naar de WGS84 ellipsoïde:

 

 

% molodensky datum omzetting van Bessel1841 naar WGS84

% Formules volgens website Peter Dana.

% Gecontroleerd met rekenvoorbeeld en GPS.

% In/uitvoer in graden en meters

 

% verbeterd nav check met EPS: fout in berekening Rm gehaald.

% 8 feb 2001

% 1 dec 2001: hoogte niet *deg!

 

function [to_geo]=bessel2wgs(from_geo);

 

deg = pi/180; % omrekenfactor naar radialen

 

from_phi = from_geo(1)*deg; % from_hoogte boven evenaar

from_lab = from_geo(2)*deg; % from_oosterlengte

from_h = from_geo(3); % from_hoogte boven zeespiegel

 

% parameters voor WGS 84:

to_a = 6378137.000; % semi-major axis of ellipsoid

to_f = 1/298.257223563; % 1/flattening

 

% parameters voor Bessel 1841:

from_a = 6377397.155;

from_f = 1/299.1528128;

 

% correctiefactoren uitrekenen:

da = to_a - from_a;

df = to_f - from_f;

bda = 1-from_f;

from_es = 2*from_f - from_f^2;

 

dx = 593;

dy = 26;

dz = 478;

 

% Nu coordinaten omrekenen van Bessel 1841 naar WGS84

 

sphi = sin(from_phi);

cphi = cos(from_phi);

slab = sin(from_lab);

clab = cos(from_lab);

 

Rn = from_a/sqrt(1.0 - from_es*(sin(from_phi))^2);

 

% oud/fout: Rm = from_a*(1-from_es)/(1.0-from_es*(sin(from_phi)^2)^(3/2));

Rm = from_a*(1-from_es)/((1.0-from_es*(sin(from_phi)^2))^(3/2));

 

term1 = (((-dx*sphi*clab-dy*sphi*slab) + dz*cphi) + da*Rn*from_es*sphi*cphi/from_a);

term2 = df*(Rm/bda + Rn*bda)*sphi*cphi;

 

dphi = (term1 + term2)/(Rm+from_h);

dlab = (-dx*slab + dy*clab)/((Rn+from_h)*cphi);

dh = dx*cphi*clab + dy*cphi*slab + dz*sphi - da*from_a/Rn + df*bda*Rn*sphi*sphi;

 

to_phi = from_phi+dphi;

to_lab = from_lab + dlab;

to_h = from_h + dh;

 

to_phi = to_phi/deg;

to_lab = to_lab/deg;

to_h = from_h + dh;

 

to_geo = [to_phi,to_lab,to_h];

%

%-------------------------------------------

 

 

Zo doe je dat dus :P

 

makkie :D

gpsgek

gpsgek

Geplaatst
  • Auteur
Geplaatst (bewerkt)

Prof Lupardi, ik ben ermee aan het stoeien gegaan maar loop een klein beetje vast. Voor de laatste transformatie is een hoogte nodig. Ik heb daarvoor 0 ingevoerd, maar krijg aan het einde van de berekening wèl een hoogte in het WGS formaat :P Bovendien klopt de positie in WGS ook niet helemaal. :D

 

Ik heb de code zo aangepast:

 

% rd2bessel
%
% Omrekenen geo coordinaten (tov Bessel) naar RD-grid coordinaten
% Formules van Geodetic Datum of Australia website
% 1 december 2001

%function geo = rd2bessel(rd);

%INVOER
rd(1)=84000;
rd(2)=448000;

RDx = rd(1);
RDy = rd(2);

deg = pi/180; % omrekenfactor naar radialen

% Bessel 1841 datum:
a = 6377397.155;
f = 1/299.1528128;

% RD grid parameters:
lab0 = 5.3876333*deg; % origin
Efalse = 155000; % false easting
K0 = 0.9999079; % scale factor
Nfalse = -5316592; % false northing

% omrekenen:

e2 = f*(2-f);
e = sqrt(e2);

a0 = 1 - 1/4*e^2 - 3/64*e^4 - 5/256*e^6;
a2 = 3/8*(e^2 + 1/4*e^4 + 15/128*e^6);
a4 = 15/256*(e^4 + 3/4*e^6);
a6 = 35/3072*e^6;

m = (RDy-Nfalse)/K0;

n = f/(2-f);
G = a*(1-n)*(1-n^2)*(1+9/4*n^2 + 225/64*n^4)*pi/180;
sig = m/G*pi/180;

% phic = phi-accent.
phic = sig +((3*n/2)-(27*n^3/32))*sin(2*sig) + ((21*n^2/16)-(55*n^4/32))*sin(4*sig)...
+(151*n^3/96)*sin(6*sig) + (1097*n^4/512)*sin(8*sig);

rho = a*(1-e^2)/(1-e^2*(sin(phic))^2)^(3/2);
nu = a/sqrt(1-e^2*(sin(phic))^2);
psi = nu/rho;

% Grid to Geographical:

t = tan(phic);
Ec = RDx - Efalse;
x = Ec/(K0*nu);

term1 = (t/(K0*rho))*(x*Ec/2);
term2 = (t/(K0*rho))*(Ec*x^3/24)*(-4*psi^2+9*psi*(1-t^2)+12*t^2);
term3 = (t/(K0*rho))*((Ec*x^5)/720)*(8*psi^4*(11-24*t^2)-12*psi^3*(21-71*t^2)+15*(psi^2*15-98*t^2+15*t^4)+180*psi*(5*t^2-3*t^4)+360*t^4);
term4 = (t/(K0*rho))*(Ec*x^7/40320)*(1385+3633*t^2+4095*t^4+1575*t^6);

phi = phic - term1 + term2 - term3 + term4;

term1 = x*sec(phi);
term2 = (x^3/6)*sec(phi)*(psi+2*t^2);
term3 = (x^5/120)*sec(phi)*(-4*psi^3*(1-6*t^2)+psi^2*(9-68*t^2)+72*psi*t^2+24*t^4);
term4 = (x^7/5040)*sec(phi)*(61+662*t^2+1320*t^4+720*t^6);

omega = term1-term2+term3-term4;

lab = lab0 + omega;

geo(1) = phi/deg; % lengte [graden]
geo(2) = lab/deg; % breedte [graden]
%
%-------------------------------------------------------------------



% molodensky datum omzetting van Bessel1841 naar WGS84
% Formules volgens website Peter Dana.
% Gecontroleerd met rekenvoorbeeld en GPS.
% In/uitvoer in graden en meters

% verbeterd nav check met EPS: fout in berekening Rm gehaald.
% 8 feb 2001
% 1 dec 2001: hoogte niet *deg!

%function [to_geo]=bessel2wgs(from_geo);


%INVOER
from_geo(1)=geo(1);
from_geo(2)=geo(2);
from_geo(3)=0;

deg = pi/180; % omrekenfactor naar radialen

from_phi = from_geo(1)*deg; % from_hoogte boven evenaar
from_lab = from_geo(2)*deg; % from_oosterlengte
from_h = from_geo(3); % from_hoogte boven zeespiegel

% parameters voor WGS 84:
to_a = 6378137.000; % semi-major axis of ellipsoid
to_f = 1/298.257223563; % 1/flattening

% parameters voor Bessel 1841:
from_a = 6377397.155;
from_f = 1/299.1528128;

% correctiefactoren uitrekenen:
da = to_a - from_a;
df = to_f - from_f;
bda = 1-from_f;
from_es = 2*from_f - from_f^2;

dx = 593;
dy = 26;
dz = 478;

% Nu coordinaten omrekenen van Bessel 1841 naar WGS84

sphi = sin(from_phi);
cphi = cos(from_phi);
slab = sin(from_lab);
clab = cos(from_lab);

Rn = from_a/sqrt(1.0 - from_es*(sin(from_phi))^2);

% oud/fout: Rm = from_a*(1-from_es)/(1.0-from_es*(sin(from_phi)^2)^(3/2));
Rm = from_a*(1-from_es)/((1.0-from_es*(sin(from_phi)^2))^(3/2));

term1 = (((-dx*sphi*clab-dy*sphi*slab) + dz*cphi) + da*Rn*from_es*sphi*cphi/from_a);
term2 = df*(Rm/bda + Rn*bda)*sphi*cphi;

dphi = (term1 + term2)/(Rm+from_h);
dlab = (-dx*slab + dy*clab)/((Rn+from_h)*cphi);
dh = dx*cphi*clab + dy*cphi*slab + dz*sphi - da*from_a/Rn + df*bda*Rn*sphi*sphi;

to_phi = from_phi+dphi;
to_lab = from_lab + dlab;
to_h = from_h + dh;

to_phi = to_phi/deg;
to_lab = to_lab/deg;
to_h = from_h + dh;


LatitudeGraden = floor(to_phi)
LongitudeGraden = floor(to_lab)

LatitudeMinuten = (to_phi - LatitudeGraden) * 60.0
LongitudeMinuten = (to_lab - LongitudeGraden) * 60.0

%to_geo = [to_phi,to_lab,to_h]
%
%-------------------------------------------

bewerkt door gpsgek

Maak een account aan of meld je aan om een opmerking te plaatsen

Je moet lid zijn om een opmerking achter te kunnen laten

Account aanmaken

Maak een account aan in onze gemeenschap. Het is makkelijk!

Registreer een nieuw account

Aanmelden

Ben je al lid? Meld je hier aan.

Nu aanmelden
 Share
Ga naar discussieoverzicht

  • Onlangs hier   0 leden

    • Er kijken geen geregistreerde gebruikers naar deze pagina.
×
×
  • Nieuwe aanmaken...