math - sferiche - sistemi di coordinate




Come convertire le coordinate di velocità sferiche in cartesiano (2)

Ho un vettore di velocità in altitudine, longitudine, altitudine, vorrei convertirlo in coordinate cartesiane, vx, vy, vz. Il formato è dallo standard WGS84.

ecco la formula

  //------------------------------------------------------------------------------
    template <class T> 
    TVectorXYZ<T> WGS84::ToCartesian(T latitude, T longitude, T elevation)
    //------------------------------------------------------------------------------
    {
        double sinlat, coslat;
        double sinlon, coslon;
        sincos_degree(latitude,  sinlat, coslat);
        sincos_degree(longitude, sinlon, coslon);  

        const double v = a / sqrt(1 - WGS84::ee * sinlat*sinlat);

        TVectorXYZ<T> coord
        (
            static_cast<T>((v + elevation) * coslat * sinlon),
            static_cast<T>(((1 - WGS84::ee) * v + elevation) * sinlat),
            static_cast<T>((v + elevation) * coslat * coslon)                                    
        );

        return coord;
    }

Prendi la formula che usi per convertire le posizioni da coordinate geografiche a coordinate cartesiane. Questo è un vettore p (λ, φ, h) ∈ ℝ³, ovvero si trasformano la latitudine, la longitudine e l'altitudine in un vettore a tre elementi delle coordinate x, y, z. Calcola ora le derivate parziali di questa formula rispetto ai tre parametri. Otterrete tre vettori, che dovrebbero essere ortogonali tra loro. La derivata rispetto alla longitudine λ dovrebbe puntare localmente a est, quella rispetto alla latitudine φ che punta a nord, quella rispetto all'altitudine h rivolta verso l'alto. Moltiplicare questi vettori con le velocità per ottenere un vettore di velocità cartesiano.

Osserva come le unità corrispondono: la posizione è in metri, le prime due derivate sono metri per grado e la velocità sarebbe gradi al secondo. O qualcos'altro, forse miglia e radianti.

Tutto questo è abbastanza facile per la sfera. Per l'ellissoide WGS84 la formula di posizione è un po 'più coinvolta e tale complessità porterà il calcolo.


OK in base alla domanda precedente e lungo il flusso di commenti, si assume che il proprio input sia:

lon [rad], lat [rad], alt [m] // WGS84 position
vlon [m/s], vlat [m/s], alt [m/s] // speed in WGS84  lon,lat,alt directions but in [m/s]

E vuoi uscita:

x,y,z // Cartesian position [m/s]
vx,vy,vz // Cartesian velocity [m/s]

E avere una valida trasformazione in coordinate cartesiane per le posizioni a tua disposizione questo è mio:

void WGS84toXYZ(double &x,double &y,double &z,double lon,double lat,double alt) // [rad,rad,m] -> [m,m,m]
    {
    const double  _earth_a=6378137.00000;   // [m] WGS84 equator radius
    const double  _earth_b=6356752.31414;   // [m] WGS84 epolar radius
    const double  _earth_e=8.1819190842622e-2; //  WGS84 eccentricity
    const double _aa=_earth_a*_earth_a;
    const double _ee=_earth_e*_earth_e;
    double  a,b,x,y,z,h,l,c,s;
    a=lon;
    b=lat;
    h=alt;
    c=cos(b);
    s=sin(b);
    // WGS84 from eccentricity
    l=_earth_a/sqrt(1.0-(_ee*s*s));
    x=(l+h)*c*cos(a);
    y=(l+h)*c*sin(a);
    z=(((1.0-_ee)*l)+h)*s;
    }

E routine per normalizzare il vettore in base alle dimensioni dell'unità:

void normalize(double &x,double &y,double &z)
    {
    double l=sqrt(x*x+y*y+z*z);
    if (l>1e-6) l=1.0/l;
    x*=l; y*=l; z*=l;
    }

Sì, puoi provare a ricavare la formula che suggerisci @MvG ma dai tuoi errori da principiante dubito fortemente che possa portare a risultati positivi. Invece puoi farlo:

  1. ottieni i vettori di lon,lat,alt per la tua posizione (x,y,z)

    è facile usare solo qualche piccolo incremento di gradino nella posizione WGS84 convertire in substract cartesiano e normalizzare in vettori unitari. Chiamiamo questi vettori di base della direzione U,V,W

    double Ux,Uy,Uz;    // [m]
    double Vx,Vy,Vz;    // [m]
    double Wx,Wy,Wz;    // [m]
    double da=1.567e-7; // [rad] angular step ~ 1.0 m in lon direction
    double dl=1.0;      // [m] altitide step 1.0 m
    WGS84toXYZ( x, y, z,lon   ,lat,alt   ); // actual position
    WGS84toXYZ(Ux,Uy,Uz,lon+da,lat,alt   ); // lon direction Nort
    WGS84toXYZ(Vx,Vy,Vz,lon,lat+da,alt   ); // lat direction East
    WGS84toXYZ(Wx,Wy,Wz,lon,lat   ,alt+dl); // alt direction High/Up
    Ux-=x; Uy-=y; Uz-=z;
    Vx-=x; Vy-=y; Vz-=z;
    Wx-=x; Wy-=y; Wz-=z;
    normalize(Ux,Uy,Uz);
    normalize(Vx,Vy,Vz);
    normalize(Wx,Wy,Wz);
  2. convertire la velocità da lon,lat,alt a vx,vy,vz

    vx = vlon*Ux + vlat*Vx + valt*Wx;
    vy = vlon*Uy + vlat*Vy + valt*Wy;
    vz = vlon*Uz + vlat*Vz + valt*Wz;

Spero sia abbastanza chiaro. Come al solito fai attenzione alle unità deg/rad e m/ft/km perché le unità contano molto.

Btw I vettori di base U,V,W formano il frame di riferimento NEH e nello stesso tempo sono le derivate di direzione che MvG sta citando.

[Modifica1] conversioni più precise

//---------------------------------------------------------------------------
//--- WGS84 transformations ver: 1.00 ---------------------------------------
//---------------------------------------------------------------------------
#ifndef _WGS84_h
#define _WGS84_h
//---------------------------------------------------------------------------
// http://www.navipedia.net/index.php/Ellipsoidal_and_Cartesian_Coordinates_Conversion
//---------------------------------------------------------------------------
// WGS84(a,b,h) = (long,lat,alt) [rad,rad,m]
// XYZ(x,y,z) [m]
//---------------------------------------------------------------------------
const double  _earth_a=6378137.00000;   // [m] WGS84 equator radius
const double  _earth_b=6356752.31414;   // [m] WGS84 epolar radius
const double  _earth_e=8.1819190842622e-2; //  WGS84 eccentricity
//const double  _earth_e=sqrt(1.0-((_earth_b/_earth_a)*(_earth_b/_earth_a)));
const double  _earth_ee=_earth_e*_earth_e;
//---------------------------------------------------------------------------
const double kmh=1.0/3.6;               // [km/h] -> [m/s]
//---------------------------------------------------------------------------
void XYZtoWGS84       (double *abh                  ,double *xyz                  ); // [m,m,m] -> [rad,rad,m]
void WGS84toXYZ       (double *xyz                  ,double *abh                  ); // [rad,rad,m] -> [m,m,m]
void WGS84toXYZ_posvel(double *xyzpos,double *xyzvel,double *abhpos,double *abhvel); // [rad,rad,m],[m/s,m/s,m/s] -> [m,m,m],[m/s,m/s,m/s]
void WGS84toNEH       (reper &neh                   ,double *abh                  ); // [rad,rad,m] -> NEH [m]
void WGS84_m2rad      (double &da,double &db,double *abh);                           // [rad,rad,m] -> [rad],[rad] representing 1m angle step
void XYZ_interpolate  (double *pt,double *p0,double *p1,double t);                   // [m,m,m] pt = p0 + (p1-p0)*t in ellipsoid space t = <0,1>
//---------------------------------------------------------------------------
void XYZtoWGS84(double *abh,double *xyz)
    {
    int i;
    double  a,b,h,l,n,db,s;
    a=atanxy(xyz[0],xyz[1]);
    l=sqrt((xyz[0]*xyz[0])+(xyz[1]*xyz[1]));
    // estimate lat
    b=atanxy((1.0-_earth_ee)*l,xyz[2]);
    // iterate to improve accuracy
    for (i=0;i<100;i++)
        {
        s=sin(b); db=b;
        n=divide(_earth_a,sqrt(1.0-(_earth_ee*s*s)));
        h=divide(l,cos(b))-n;
        b=atanxy((1.0-divide(_earth_ee*n,n+h))*l,xyz[2]);
        db=fabs(db-b);
        if (db<1e-12) break;
        }
    if (b>0.5*pi) b-=pi2;
    abh[0]=a;
    abh[1]=b;
    abh[2]=h;
    }
//---------------------------------------------------------------------------
void WGS84toXYZ(double *xyz,double *abh)
    {
    double  a,b,h,l,c,s;
    a=abh[0];
    b=abh[1];
    h=abh[2];
    c=cos(b);
    s=sin(b);
    // WGS84 from eccentricity
    l=_earth_a/sqrt(1.0-(_earth_ee*s*s));
    xyz[0]=(l+h)*c*cos(a);
    xyz[1]=(l+h)*c*sin(a);
    xyz[2]=(((1.0-_earth_ee)*l)+h)*s;
    }
//---------------------------------------------------------------------------
void WGS84toNEH(reper &neh,double *abh)
    {
    double N[3],E[3],H[3];                  // [m]
    double p[3],xyzpos[3];
    const double da=1.567e-7;               // [rad] angular step ~ 1.0 m in lon direction
    const double dl=1.0;                    // [m] altitide step 1.0 m
    vector_copy(p,abh);
    // actual position
    WGS84toXYZ(xyzpos,abh);
    // NEH
    p[0]+=da; WGS84toXYZ(N,p); p[0]-=da;
    p[1]+=da; WGS84toXYZ(E,p); p[1]-=da;
    p[2]+=dl; WGS84toXYZ(H,p); p[2]-=dl;
    vector_sub(N,N,xyzpos);
    vector_sub(E,E,xyzpos);
    vector_sub(H,H,xyzpos);
    vector_one(N,N);
    vector_one(E,E);
    vector_one(H,H);
    neh._rep=1;
    neh._inv=0;
    // axis X
    neh.rep[ 0]=N[0];
    neh.rep[ 1]=N[1];
    neh.rep[ 2]=N[2];
    // axis Y
    neh.rep[ 4]=E[0];
    neh.rep[ 5]=E[1];
    neh.rep[ 6]=E[2];
    // axis Z
    neh.rep[ 8]=H[0];
    neh.rep[ 9]=H[1];
    neh.rep[10]=H[2];
    // gpos
    neh.rep[12]=xyzpos[0];
    neh.rep[13]=xyzpos[1];
    neh.rep[14]=xyzpos[2];
    neh.orto(1);
    }
//---------------------------------------------------------------------------
void WGS84toXYZ_posvel(double *xyzpos,double *xyzvel,double *abhpos,double *abhvel)
    {
    reper neh;
    WGS84toNEH(neh,abhpos);
    neh.gpos_get(xyzpos);
    neh.l2g_dir(xyzvel,abhvel);
    }
//---------------------------------------------------------------------------
void WGS84_m2rad(double &da,double &db,double *abh)
    {
    // WGS84 from eccentricity
    double p[3],rr;
    WGS84toXYZ(p,abh);
    rr=(p[0]*p[0])+(p[1]*p[1]);
    da=divide(1.0,sqrt(rr));
    rr+=p[2]*p[2];
    db=divide(1.0,sqrt(rr));
    }
//---------------------------------------------------------------------------
void XYZ_interpolate(double *pt,double *p0,double *p1,double t)
    {
    const double  mz=_earth_a/_earth_b;
    const double _mz=_earth_b/_earth_a;
    double p[3],r,r0,r1;
    // compute spherical radiuses of input points
    r0=sqrt((p0[0]*p0[0])+(p0[1]*p0[1])+(p0[2]*p0[2]*mz*mz));
    r1=sqrt((p1[0]*p1[0])+(p1[1]*p1[1])+(p1[2]*p1[2]*mz*mz));
    // linear interpolation
    r   = r0   +(r1   -r0   )*t;
    p[0]= p0[0]+(p1[0]-p0[0])*t;
    p[1]= p0[1]+(p1[1]-p0[1])*t;
    p[2]=(p0[2]+(p1[2]-p0[2])*t)*mz;
    // correct radius and rescale back
    r/=sqrt((p[0]*p[0])+(p[1]*p[1])+(p[2]*p[2]));
    pt[0]=p[0]*r;
    pt[1]=p[1]*r;
    pt[2]=p[2]*r*_mz;
    }
//---------------------------------------------------------------------------
#endif
//---------------------------------------------------------------------------

Tuttavia essi richiedono la matematica di base del vettore 3D vedere qui per le equazioni:

  • Comprendere le matrici di trasformazione omogenea 4x4