EOV to WGS84

Algoritmusok

Sziasztok,

 

Azzal a kihívással találkoztam, hogy EOVX és EOVY koordinátákat kellene átszámolnom GPS lon és lat formátumokra.

Találtam online átszámítókat, de maga a képlet nincs meg. Szeretnék írni egy kis programot, ami egy csv fájlt alakít át lokálisan.

Nézegettem pár szakdolgozatban, ami elérhető online, de sajnos nem találtam meg a képletet.

Köszi a segítséget előre is.

 

Itt tartok:

 

using System;


namespace GPS_transformer
{
    class Position
    {
        public DateTime timeStamp { get; private set; }
        public string address { get; private set; }
        public int eovx { get; private set; }
        public int eovy { get; private set; }

        public double longitude { get; private set; }

        public double latitude { get; private set; }

        public Position(DateTime timeStamp, string address, int eovx, int eovy)
        {
            this.timeStamp = timeStamp;
            this.address = address;
            this.eovx = eovx;
            this.eovy = eovy;
            this.EOV_to_GPS84();
        }

        private void EOV_to_GPS84()
        {
            //longitude = 47;
            //latitude = 19;
        }
    }
}

Üdv,

 

Zoli

( SDN | 2021. 06. 24., cs – 07:53 )

Szerintem egyszerűbb valami proj4 wrappert, vagy DotSpatial lib-et használni.

( uid_6201 v | 2021. 06. 24., cs – 08:17 )

Nem lesz rá egzakt képlet. A két vetület közötti átjárás csak fokozatos közelítő mágia segítségével lehetséges.

( sebist v | 2021. 06. 24., cs – 08:32 )

Most nem tudtam megtalálni azt a leírást(valószínűleg a változónevek is azt követik), ami alapján a lenti jó régi C kód született, de talán használható belőle valami.

 


long double rad ( long double x)
{
    return (x/180.0*M_PI);
}    

long double deg ( long double x)
{
    return (x/M_PI*180.0);
}    

struct WGS84 {
    long double fi,la;
};
    
    
struct WGS84 eov2wgs (long double EOVY, long double EOVX)
{

struct WGS84 ret;
long double a1=6378160,b1=6356774.52,a2=6378137,b2=6356752.31;
long double dX=52.684,dY=-71.194,dZ=-13.975,eX=0.312,eY=0.1063,eZ=0.3729,k=1.0191e-6;
long double FI2,LA2,st,ct,teta,P,ev2,e22,f2,Zv,Yv,Xv;
long double h=0; // magassag
long double Z,Y,X,N,e2,f,fi,la;

long double C1=1.0007197049;
long double E1=rad(19.048571778);

long double G1=rad(47.1);

long double H1=6379296.419;

long double I1=47.0+7.0/60.0+20.0578/3600.0;
long double J1=rad(I1);

long double K1=EOVX-200000;
long double L1=EOVY-650000;

long double M1=2.0*(atan(exp(K1/H1))-M_PI/4.0);
long double N1=L1/H1;

long double O1=47+1.0/6.0;
long double P1=asin(cos(G1)*sin(M1)+sin(G1)*cos(M1)*cos(N1));

long double Q1=asin(sin(N1)*cos(M1)/cos(P1));

//R1
long double T1=rad(O1);

long double U1=6378160;
long double V1=6356774.516;

long double W1=(U1*U1-V1*V1)*cos(T1)*cos(T1)/V1/V1;    
long double X1=180*3600/M_PI;

long double Y1=sqrt(1.0+W1);
long double Z1=1.5*W1*tan(T1)/X1;

long double AA1=0.5*W1*(-1.0+tan(T1)*tan(T1)-W1+5*W1*tan(T1))/Y1/X1/X1;
long double S1=(P1-J1)*X1;

fi=(T1+S1*Y1/X1-S1*S1*Z1/X1+S1*S1*S1*AA1/X1);
la=(E1+Q1/C1);
    
f=(a1-b1)/a1;
e2=2*f-f*f;

N=a1/(sqrt(1-e2*sinl(fi)*sinl(fi)));
X=(N+h)*cos(fi)*cosl(la);

Y=(N+h)*cosl(fi)*sinl(la);
Z=(N*(1-e2)+h)*sinl(fi);
    
Xv=dX+(1+k)*( X+rad(eZ/3600)*Y-rad(eY/3600)*Z);
Yv=dY+(1+k)*(-X*rad(eZ/3600)+Y+Z*rad(eX/3600));
Zv=dZ+(1+k)*( X*rad(eY/3600)-Y*rad(eX/3600)+Z);

f2=(a2-b2)/a2;
e22=2*f2-f2*f2;
ev2=(a2*a2-b2*b2)/b2/b2;    

P=sqrt(Xv*Xv+Yv*Yv);
teta=atan2l(Zv*a2,P*b2);

ct=cosl(teta);
st=sinl(teta);
FI2=atan2l(Zv+ev2*b2*st*st*st,P-e22*a2*ct*ct*ct);
LA2=atan2l(Yv,Xv);

ret.fi=deg(FI2);
ret.la=deg(LA2);
        
return ret;

}

Köszi! Átültetem c sharp osztályba, igazabol majdnem keszen vagyok, csak nem talalom a Csharp statikus Math osztályában a sinl meg cosl függvényeket. Ezek kihelyettesíthetőek?

 

using System;

namespace GPS_transformer
{
    struct _WGS84
    {
        double fi, la;
    };
    class WGS84
    {
        static double rad(double x)
        {
            return (x / 180.0 * Math.PI);
        }

        static double deg(double x)
        {
            return (x / Math.PI * 180.0);
        }



        public _WGS84 eov2wgs(double EOVY, double EOVX)
        {
            _WGS84 ret;
            double a1 = 6378160, b1 = 6356774.52, a2 = 6378137, b2 = 6356752.31;
            double dX = 52.684, dY = -71.194, dZ = -13.975, eX = 0.312, eY = 0.1063, eZ = 0.3729, k = 1.0191e-6;
            double FI2, LA2, st, ct, teta, P, ev2, e22, f2, Zv, Yv, Xv;
            double h = 0; // magassag
            double Z, Y, X, N, e2, f, fi, la;

            double C1 = 1.0007197049;
            double E1 = rad(19.048571778);

            double G1 = rad(47.1);

            double H1 = 6379296.419;

            double I1 = 47.0 + 7.0 / 60.0 + 20.0578 / 3600.0;
            double J1 = rad(I1);

            double K1 = EOVX - 200000;
            double L1 = EOVY - 650000;

            double M1 = 2.0 * (Math.Atan(Math.Exp(K1 / H1)) - Math.PI / 4.0);
            double N1 = L1 / H1;

            double O1 = 47 + 1.0 / 6.0;
            double P1 = Math.Asin(Math.Cos(G1) * Math.Sin(M1) + Math.Sin(G1) * Math.Cos(M1) * Math.Cos(N1));

            double Q1 = Math.Asin(Math.Sin(N1) * Math.Cos(M1) / Math.Cos(P1));

            //R1
            double T1 = rad(O1);

            double U1 = 6378160;
            double V1 = 6356774.516;

            double W1 = (U1 * U1 - V1 * V1) * Math.Cos(T1) * Math.Cos(T1) / V1 / V1;
            double X1 = 180 * 3600 / Math.PI;

            double Y1 = Math.Sqrt(1.0 + W1);
            double Z1 = 1.5 * W1 * Math.Tan(T1) / X1;

            double AA1 = 0.5 * W1 * (-1.0 + Math.Tan(T1) * Math.Tan(T1) - W1 + 5 * W1 * Math.Tan(T1)) / Y1 / X1 / X1;
            double S1 = (P1 - J1) * X1;

            fi = (T1 + S1 * Y1 / X1 - S1 * S1 * Z1 / X1 + S1 * S1 * S1 * AA1 / X1);
            la = (E1 + Q1 / C1);

            f = (a1 - b1) / a1;
            e2 = 2 * f - f * f;

            N = a1 / (Math.Sqrt(1 - e2 * sinl(fi) * sinl(fi)));
            X = (N + h) * Math.Cos(fi) * cosl(la);

            Y = (N + h) * Math.Cos(fi) * sinl(la);
            Z = (N * (1 - e2) + h) * sinl(fi);

            Xv = dX + (1 + k) * (X + rad(eZ / 3600) * Y - rad(eY / 3600) * Z);
            Yv = dY + (1 + k) * (-X * rad(eZ / 3600) + Y + Z * rad(eX / 3600));
            Zv = dZ + (1 + k) * (X * rad(eY / 3600) - Y * rad(eX / 3600) + Z);

            f2 = (a2 - b2) / a2;
            e22 = 2 * f2 - f2 * f2;
            ev2 = (a2 * a2 - b2 * b2) / b2 / b2;

            P = sqrt(Xv * Xv + Yv * Yv);
            teta = atan2l(Zv * a2, P * b2);

            ct = cosl(teta);
            st = sinl(teta);
            FI2 = atan2l(Zv + ev2 * b2 * st * st * st, P - e22 * a2 * ct * ct * ct);
            LA2 = atan2l(Yv, Xv);

            ret.fi = deg(FI2);
            ret.la = deg(LA2);

            return ret;
        }
    }
}

-- Zoli ---

Lenovo T400 @ Crunchbang "Waldorf @ SSD

( KAMI911 v | 2021. 06. 24., cs – 08:57 )

Szerkesztve: 2021. 06. 24., cs – 08:58

Proj4?

Esetleg ez hátha segítség:

https://github.com/KAMI911/lactransformer

Ez csináltam vele Pythonban. Nem nagyon néztem rá mostanában, de ha hibát találsz akkor szólj és próbálom javítani. LAS, CSV és még egy pár formátumon működik. VITEL2014-es korrekciót is, ha deciméter alatti pontosság kell.

http://eht.gnssnet.hu/images/hdif.jpg

Ha pedig lakossági RTK-ra van szükséged, akkor erről itt egy cikkem:

https://linuxmint.hu/blog/2020/10/gps-utan-pontosabb-meres-nyomvonal-ke…

KAMI | 神
Linux Mint: https://linuxmint.hu LibreOffice: http://libreoffice.hu SeaMonkey : http://is.gd/smhun
Legyél Te is Mozilla és Linux Mint önkéntes

( knzoltan | 2021. 06. 24., cs – 09:43 )

elvileg a sinl csak annyiban különbözik hogy long double ben adja vissza...

-- Zoli ---

Lenovo T400 @ Crunchbang "Waldorf @ SSD

( harlequin | 2021. 06. 24., cs – 12:07 )

Szerkesztve: 2021. 06. 24., cs – 12:07

pattints föl egy PostGIS-t, és queryvel szépen és gyorsan átszámolja, rákapcsolódni meg bármilyen nyelvből rátudsz.

Valami ilyesmi: SELECT ST_Transform(ST_SetSRID( ST_Point(x, y), 23700), 4326) (de lehet hogy nem pont így, reg hasznaltam :) )

“Any book worth banning is a book worth reading.”

( falu v | 2021. 06. 24., cs – 17:37 )

Nos az a baj, hogy a két rendszer között (mivel nagyon nem azonos az alapfelület) nem létezik egzakt átszámítási módszer.

A lényeg az, hogy először az EOV koordinátákból kell GRS67 geocentrikus (X, Y, Z) koordinátákat számolni, majd valamilyen trükkel átszámítani a WGS84 ellipszoidra. Az átszámítás általában egy térbeli hasonlósági transzformáció, melynek paramétereit olyan pontokból számítjuk, amelyeknek mind a két rendszerben ismert a koordinátája.

A szakdolgozatomban egy olyan megoldást mutattam be, amely először közelítő (jellemzően 1 m alatti) pontossággal számítja át a pontokat, majd egy javítórács segítségével pontosítja a számítási eredményeket. A javítórács úgy készül, hogy a pl 1 km * 1 km-es rácsháló pontjait kiszámítjuk az közelítő és az előző  pontos módszerrel, és az eltéréseket tároljuk egy adatbázisban. Innentől kezdve nincs más dolgunk, mint lekérdezni az adatbázisból az átszámítandó pont közelébe eső pl 4 rácspontot és interpolációval megbecsülni, hogy mekkora lehet az eltérés a kérdéses ponton, majd a becsült javítási étékekkel módosítani az előzetes koordinátát.

[Falu.me](Tárhely|Domain)