root/GPXSee
1
0
Fork 0
mirror of https://github.com/tumic0/GPXSee.git synced 2024-11-28 13:41:16 +01:00
Code Issues Projects Releases Activity

Optimized projection computation

Browse Source
This commit is contained in:
Martin Tůma 2017-04-02 22:18:03 +02:00
parent a56aa4e706
commit e46bba18f2
3 changed files with 51 additions and 61 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

103
src/transversemercator.cpp
View File

@ -4,14 +4,6 @@
#include "transversemercator.h" #include "transversemercator.h"
TransverseMercator::TransverseMercator()
{
_centralMeridian = 0;
_scale = 1.0;
_falseEasting = 0;
_falseNorthing = 0;
}
TransverseMercator::TransverseMercator(double centralMeridian, double scale, TransverseMercator::TransverseMercator(double centralMeridian, double scale,
double falseEasting, double falseNorthing) double falseEasting, double falseNorthing)
{ {
@ -19,48 +11,58 @@ TransverseMercator::TransverseMercator(double centralMeridian, double scale,
_scale = scale; _scale = scale;
_falseEasting = falseEasting; _falseEasting = falseEasting;
_falseNorthing = falseNorthing; _falseNorthing = falseNorthing;
const double e2 = WGS84_FLATTENING * (2 - WGS84_FLATTENING);
const double n = WGS84_FLATTENING / (2 - WGS84_FLATTENING);
_rectifyingRadius = WGS84_RADIUS / (1 + n)
* (1 + 0.25*pow(n, 2) + 0.015625*pow(n, 4));
_A = e2;
_B = (5 * pow(e2, 2) - pow(e2, 3)) / 6.0;
_C = (104 * pow(e2, 3) - 45 * pow(e2, 4)) / 120.0;
_D = (1237 * pow(e2, 4)) / 1260.0;
_beta1 = 1/2.0 * n - 2/3.0 * pow(n, 2) + 5/16.0 * pow(n, 3) + 41/180.0
* pow(n, 4);
_beta2 = 13/48.0 * pow(n, 2) - 3/5.0 * pow(n, 3) + 557/1440.0 * pow(n, 4);
_beta3 = 61/240.0 * pow(n, 3) - 103/140.0 * pow(n, 4);
_beta4 = 49561/161280.0 * pow(n, 4);
_delta1 = 1/2.0 * n - 2/3.0 * pow(n, 2) + 37/96.0 * pow(n, 3) - 1/360.0
* pow(n, 4);
_delta2 = 1/48.0 * pow(n, 2) + 1/15.0 * pow(n, 3) - 437/1440.0 * pow(n, 4);
_delta3 = 17/480.0 * pow(n, 3) - 37/840.0 * pow(n, 4);
_delta4 = 4397/161280.0 * pow(n, 4);
_AStar = e2 + pow(e2, 2) + pow(e2, 3) + pow(e2, 4);
_BStar = (7 * pow(e2, 2) + 17 * pow(e2, 3) + 30 * pow(e2, 4)) / -6;
_CStar = (224 * pow(e2, 3) + 889 * pow(e2, 4)) / 120;
_DStar = (4279 * pow(e2, 4)) / -1260;
} }
QPointF TransverseMercator::ll2xy(const Coordinates &c) const QPointF TransverseMercator::ll2xy(const Coordinates &c) const
{ {
QPointF p; QPointF p;
const double e2 = WGS84_FLATTENING * (2 - WGS84_FLATTENING);
const double n = WGS84_FLATTENING / (2 - WGS84_FLATTENING);
const double rectifyingRadius = WGS84_RADIUS / (1 + n)
* (1 + 0.25*pow(n, 2) + 0.015625*pow(n, 4));
double A = e2;
double B = (5 * pow(e2, 2) - pow(e2, 3)) / 6.0;
double C = (104 * pow(e2, 3) - 45 * pow(e2, 4)) / 120.0;
double D = (1237 * pow(e2, 4)) / 1260.0;
double phi = deg2rad(c.lat()); double phi = deg2rad(c.lat());
double lambda = deg2rad(c.lon()); double lambda = deg2rad(c.lon());
double lambda0 = deg2rad(_centralMeridian); double lambda0 = deg2rad(_centralMeridian);
double deltaLambda = lambda - lambda0; double deltaLambda = lambda - lambda0;
double phiStar = phi - sin(phi) * cos(phi) * (A + B*pow(sin(phi), 2) double phiStar = phi - sin(phi) * cos(phi) * (_A + _B*pow(sin(phi), 2)
+ C*pow(sin(phi), 4) + D*pow(sin(phi), 6)); + _C*pow(sin(phi), 4) + _D*pow(sin(phi), 6));
double xiPrim = atan(tan(phiStar) / cos(deltaLambda)); double xiPrim = atan(tan(phiStar) / cos(deltaLambda));
double etaPrim = atanh(cos(phiStar) * sin(deltaLambda)); double etaPrim = atanh(cos(phiStar) * sin(deltaLambda));
double beta1 = 1/2.0 * n - 2/3.0 * pow(n, 2) + 5/16.0 * pow(n, 3) p.ry() = _falseNorthing + _scale * _rectifyingRadius * (xiPrim + _beta1
+ 41/180.0 * pow(n, 4); * sin(2*xiPrim) * cosh(2*etaPrim) + _beta2 * sin(4*xiPrim)
double beta2 = 13/48.0 * pow(n, 2) - 3/5.0 * pow(n, 3) + 557/1440.0 * cosh(4*etaPrim) + _beta3 * sin(6*xiPrim) * cosh(6*etaPrim) + _beta4
* pow(n, 4);
double beta3 = 61/240.0 * pow(n, 3) - 103/140.0 * pow(n, 4);
double beta4 = 49561/161280.0 * pow(n, 4);
p.ry() = _falseNorthing + _scale * rectifyingRadius * (xiPrim + beta1
* sin(2*xiPrim) * cosh(2*etaPrim) + beta2 * sin(4*xiPrim)
* cosh(4*etaPrim) + beta3 * sin(6*xiPrim) * cosh(6*etaPrim) + beta4
* sin(8*xiPrim) * cosh(8*etaPrim)); * sin(8*xiPrim) * cosh(8*etaPrim));
p.rx() = _falseEasting + _scale * rectifyingRadius * (etaPrim + beta1 p.rx() = _falseEasting + _scale * _rectifyingRadius * (etaPrim + _beta1
* cos(2*xiPrim) * sinh(2*etaPrim) + beta2 * cos(4*xiPrim) * cos(2*xiPrim) * sinh(2*etaPrim) + _beta2 * cos(4*xiPrim)
* sinh(4*etaPrim) + beta3 * cos(6*xiPrim) * sinh(6*etaPrim) + beta4 * sinh(4*etaPrim) + _beta3 * cos(6*xiPrim) * sinh(6*etaPrim) + _beta4
* cos(8*xiPrim) * sinh(8*etaPrim)); * cos(8*xiPrim) * sinh(8*etaPrim));
return p; return p;
@ -68,38 +70,21 @@ QPointF TransverseMercator::ll2xy(const Coordinates &c) const
Coordinates TransverseMercator::xy2ll(const QPointF &p) const Coordinates TransverseMercator::xy2ll(const QPointF &p) const
{ {
const double e2 = WGS84_FLATTENING * (2 - WGS84_FLATTENING); double xi = (p.y() - _falseNorthing) / (_scale * _rectifyingRadius);
const double n = WGS84_FLATTENING / (2 - WGS84_FLATTENING); double eta = (p.x() - _falseEasting) / (_scale * _rectifyingRadius);
const double rectifyingRadius = WGS84_RADIUS / (1 + n)
* (1 + 0.25*pow(n, 2) + 0.015625*pow(n, 4));
double xi = (p.y() - _falseNorthing) / (_scale * rectifyingRadius); double xiPrim = xi - _delta1 * sin(2*xi) * cosh(2*eta) - _delta2 * sin(4*xi)
double eta = (p.x() - _falseEasting) / (_scale * rectifyingRadius); * cosh(4*eta) - _delta3 * sin(6*xi) * cosh(6*eta) - _delta4 * sin(8*xi)
double delta1 = 1/2.0 * n - 2/3.0 * pow(n, 2) + 37/96.0 * pow(n, 3)
- 1/360.0 * pow(n, 4);
double delta2 = 1/48.0 * pow(n, 2) + 1/15.0 * pow(n, 3) - 437/1440.0
* pow(n, 4);
double delta3 = 17/480.0 * pow(n, 3) - 37/840.0 * pow(n, 4);
double delta4 = 4397/161280.0 * pow(n, 4);
double xiPrim = xi - delta1 * sin(2*xi) * cosh(2*eta) - delta2 * sin(4*xi)
* cosh(4*eta) - delta3 * sin(6*xi) * cosh(6*eta) - delta4 * sin(8*xi)
* cosh(8*eta); * cosh(8*eta);
double etaPrim = eta - delta1 * cos(2*xi) * sinh(2*eta) - delta2 * cos(4*xi) double etaPrim = eta - _delta1 * cos(2*xi) * sinh(2*eta) - _delta2
* sinh(4*eta) - delta3 * cos(6*xi) * sinh(6*eta) - delta4 * cos(8*xi) * cos(4*xi) * sinh(4*eta) - _delta3 * cos(6*xi) * sinh(6*eta) - _delta4
* sinh(8*eta); * cos(8*xi) * sinh(8*eta);
double phiStar = asin(sin(xiPrim) / cosh(etaPrim)); double phiStar = asin(sin(xiPrim) / cosh(etaPrim));
double deltaLambda = atan(sinh(etaPrim) / cos(xiPrim)); double deltaLambda = atan(sinh(etaPrim) / cos(xiPrim));
double AStar = e2 + pow(e2, 2) + pow(e2, 3) + pow(e2, 4); double phi = phiStar + sin(phiStar) * cos(phiStar) * (_AStar + _BStar
double BStar = (7 * pow(e2, 2) + 17 * pow(e2, 3) + 30 * pow(e2, 4)) / -6; * pow(sin(phiStar), 2) + _CStar * pow(sin(phiStar), 4) + _DStar
double CStar = (224 * pow(e2, 3) + 889 * pow(e2, 4)) / 120;
double DStar = (4279 * pow(e2, 4)) / -1260;
double phi = phiStar + sin(phiStar) * cos(phiStar) * (AStar + BStar
* pow(sin(phiStar), 2) + CStar * pow(sin(phiStar), 4) + DStar
* pow(sin(phiStar), 6)); * pow(sin(phiStar), 6));
return Coordinates(_centralMeridian + rad2deg(deltaLambda), rad2deg(phi)); return Coordinates(_centralMeridian + rad2deg(deltaLambda), rad2deg(phi));

7
src/transversemercator.h
View File

@ -6,7 +6,6 @@
class TransverseMercator : public Projection class TransverseMercator : public Projection
{ {
public: public:
TransverseMercator();
TransverseMercator(double centralMeridian, double scale, TransverseMercator(double centralMeridian, double scale,
double falseEasting, double falseNorthing); double falseEasting, double falseNorthing);
@ -18,6 +17,12 @@ private:
double _scale; double _scale;
double _falseEasting; double _falseEasting;
double _falseNorthing; double _falseNorthing;
double _rectifyingRadius;
double _A, _B, _C, _D;
double _beta1, _beta2, _beta3, _beta4;
double _delta1, _delta2, _delta3, _delta4;
double _AStar, _BStar, _CStar, _DStar;
}; };
#endif // TRANSVERSEMERCATOR_H #endif // TRANSVERSEMERCATOR_H

2
src/utm.cpp
View File

@ -1,6 +1,6 @@
#include "utm.h" #include "utm.h"
UTM::UTM(const Coordinates &c) UTM::UTM(const Coordinates &c) : _tm(0, 1.0, 0, 0)
{ {
int zone = int((c.lon() + 180)/6) + 1; int zone = int((c.lon() + 180)/6) + 1;