root/GPXSee
1
0
Fork 0
mirror of https://github.com/tumic0/GPXSee.git synced 2025-06-25 18:49:16 +02:00
Code Issues Projects Releases Activity

Added support for prime meridians and angular units to GCS

Browse Source
This commit is contained in:
Martin Tůma
2018-01-20 20:13:56 +01:00
parent e9ef68a81c
commit 1739625896
38 changed files with 886 additions and 773 deletions
Download Patch File Download Diff File Expand all files Collapse all files

8
src/map/albersequal.cpp
View File

@ -56,7 +56,7 @@ Defense.
((clat) / sqrt(one_minus_sqr_es_sin))
AlbersEqual::AlbersEqual(const Ellipsoid &ellipsoid, double standardParallel1,
AlbersEqual::AlbersEqual(const Ellipsoid *ellipsoid, double standardParallel1,
double standardParallel2, double latitudeOrigin, double longitudeOrigin,
double falseEasting, double falseNorthing)
{
@ -78,7 +78,7 @@ AlbersEqual::AlbersEqual(const Ellipsoid &ellipsoid, double standardParallel1,
sp1 = deg2rad(standardParallel1);
sp2 = deg2rad(standardParallel2);
_es2 = 2 * _e.flattening() - _e.flattening() * _e.flattening();
_es2 = 2 * _e->flattening() - _e->flattening() * _e->flattening();
_es = sqrt(_es2);
_one_minus_es2 = 1 - _es2;
_two_es = 2 * _es;
@ -108,7 +108,7 @@ AlbersEqual::AlbersEqual(const Ellipsoid &ellipsoid, double standardParallel1,
_n = sin_lat1;
_C = sqr_m1 + _n * q1;
_a_over_n = _e.radius() / _n;
_a_over_n = _e->radius() / _n;
nq0 = _n * q0;
_rho0 = (_C < nq0) ? 0 : _a_over_n * sqrt(_C - nq0);
}
@ -171,7 +171,7 @@ Coordinates AlbersEqual::xy2ll(const QPointF &p) const
if (rho != 0.0)
theta = atan2(dx, rho0_minus_dy);
rho_n = rho * _n;
q = (_C - (rho_n * rho_n) / (_e.radius() * _e.radius())) / _n;
q = (_C - (rho_n * rho_n) / (_e->radius() * _e->radius())) / _n;
qc = 1 - ((_one_minus_es2) / (_two_es)) * log((1.0 - _es) / (1.0 + _es));
if (fabs(fabs(qc) - fabs(q)) > 1.0e-6) {
q_over_2 = q / 2.0;

4
src/map/albersequal.h
View File

@ -7,7 +7,7 @@
class AlbersEqual : public Projection
{
public:
AlbersEqual(const Ellipsoid &ellipsoid, double standardParallel1,
AlbersEqual(const Ellipsoid *ellipsoid, double standardParallel1,
double standardParallel2, double latitudeOrigin, double longitudeOrigin,
double falseEasting, double falseNorthing);
@ -15,7 +15,7 @@ public:
virtual Coordinates xy2ll(const QPointF &p) const;
private:
Ellipsoid _e;
const Ellipsoid *_e;
double _latitudeOrigin;
double _longitudeOrigin;

41
src/map/angularunits.cpp Normal file
View File

@ -0,0 +1,41 @@
#include "common/coordinates.h"
#include "angularunits.h"
AngularUnits::AngularUnits(int code)
{
switch (code) {
case 9101:
_f = 180.0 / M_PI;
break;
case 9102:
case 9107:
case 9108:
case 9110:
case 9122:
_f = 1.0;
break;
case 9103:
_f = 1 / 60.0;
break;
case 9104:
_f = 1 / 3600.0;
break;
case 9105:
_f = 180.0 / 200.0;
break;
case 9106:
_f = 180.0 / 200.0;
break;
case 9109:
_f = 180.0 / (M_PI * 1000000.0);
break;
default:
_f = NAN;
}
}
QDebug operator<<(QDebug dbg, const AngularUnits &au)
{
dbg.nospace() << "AngularUnits(" << deg2rad(au._f) << ")";
return dbg.maybeSpace();
}

31
src/map/angularunits.h Normal file
View File

@ -0,0 +1,31 @@
#ifndef ANGULARUNITS_H
#define ANGULARUNITS_H
#include <cmath>
#include <QDebug>
class AngularUnits
{
public:
AngularUnits() : _f(NAN) {}
AngularUnits(int code);
bool isNull() const {return std::isnan(_f);}
bool isValid() const {return !std::isnan(_f);}
double toDegrees(double val) const {return val * _f;}
double fromDegrees(double val) const {return val / _f;}
friend bool operator==(const AngularUnits &au1, const AngularUnits &au2);
friend QDebug operator<<(QDebug dbg, const AngularUnits &au);
private:
double _f;
};
inline bool operator==(const AngularUnits &au1, const AngularUnits &au2)
{return (au1._f == au2._f);}
QDebug operator<<(QDebug dbg, const AngularUnits &au);
#endif // ANGULARUNITS_H

152
src/map/datum.cpp
View File

@ -1,24 +1,9 @@
#include <cmath>
#include <QFile>
#include "common/wgs84.h"
#include "datum.h"
class Datum::Entry {
public:
Entry(const QString &name, int epsg, const Datum &datum)
: _name(name), _epsg(epsg), _datum(datum) {}
const QString &name() const {return _name;}
int epsg() const {return _epsg;}
const Datum &datum() const {return _datum;}
private:
QString _name;
int _epsg;
Datum _datum;
};
static Ellipsoid WGS84e = Ellipsoid(WGS84_RADIUS, WGS84_FLATTENING);
static Datum WGS84 = Datum(&WGS84e, 0.0, 0.0, 0.0);
// Abridged Molodensky transformation
static Coordinates molodensky(const Coordinates &c, const Datum &from,
@ -37,11 +22,11 @@ static Coordinates molodensky(const Coordinates &c, const Datum &from,
double dy = from.dy() - to.dy();
double dz = from.dz() - to.dz();
double from_f = from.ellipsoid().flattening();
double to_f = to.ellipsoid().flattening();
double from_f = from.ellipsoid()->flattening();
double to_f = to.ellipsoid()->flattening();
double df = to_f - from_f;
double from_a = from.ellipsoid().radius();
double to_a = to.ellipsoid().radius();
double from_a = from.ellipsoid()->radius();
double to_a = to.ellipsoid()->radius();
double da = to_a - from_a;
double from_esq = from_f * (2.0 - from_f);
double adb = 1.0 / (1.0 - from_f);
@ -58,138 +43,25 @@ static Coordinates molodensky(const Coordinates &c, const Datum &from,
return Coordinates(c.lon() + rad2deg(dlon), c.lat() + rad2deg(dlat));
}
QList<Datum::Entry> Datum::_datums = WGS84();
QString Datum::_errorString;
int Datum::_errorLine = 0;
QList<Datum::Entry> Datum::WGS84()
Datum::Datum(const Ellipsoid *ellipsoid, double dx, double dy,
double dz) : _ellipsoid(ellipsoid), _dx(dx), _dy(dy), _dz(dz)
{
QList<Datum::Entry> list;
list.append(Datum::Entry("WGS 84", 4326,
Datum(Ellipsoid(WGS84_RADIUS, WGS84_FLATTENING), 0.0, 0.0, 0.0)));
return list;
}
Datum::Datum(const Ellipsoid &ellipsoid, double dx, double dy, double dz)
: _ellipsoid(ellipsoid), _dx(dx), _dy(dy), _dz(dz)
{
_WGS84 = (_ellipsoid.radius() == WGS84_RADIUS
&& _ellipsoid.flattening() == WGS84_FLATTENING
&& _dx == 0.0 && _dy == 0.0 && _dz == 0.0)
? true : false;
}
Datum::Datum(int id)
{
for (int i = 0; i < _datums.size(); i++) {
if (_datums.at(i).epsg() == id) {
*this = _datums.at(i).datum();
return;
}
}
*this = Datum();
}
Datum::Datum(const QString &name)
{
for (int i = 0; i < _datums.size(); i++) {
if (_datums.at(i).name() == name) {
*this = _datums.at(i).datum();
return;
}
}
*this = Datum();
}
void Datum::error(const QString &str)
{
_errorString = str;
_datums = WGS84();
}
bool Datum::loadList(const QString &path)
{
QFile file(path);
bool res;
if (!file.open(QFile::ReadOnly)) {
error(file.errorString());
return false;
}
_errorLine = 1;
_errorString.clear();
while (!file.atEnd()) {
QByteArray line = file.readLine();
QList<QByteArray> list = line.split(',');
if (list.size() != 6) {
error("Format error");
return false;
}
QString f1 = list[1].trimmed();
int id = 0;
if (!f1.isEmpty()) {
id = f1.toInt(&res);
if (!res) {
error("Invalid datum id");
return false;
}
}
int eid = list[2].trimmed().toInt(&res);
if (!res) {
error("Invalid ellipsoid id");
return false;
}
double dx = list[3].trimmed().toDouble(&res);
if (!res) {
error("Invalid dx");
return false;
}
double dy = list[4].trimmed().toDouble(&res);
if (!res) {
error("Invalid dy");
return false;
}
double dz = list[5].trimmed().toDouble(&res);
if (!res) {
error("Invalid dz");
return false;
}
Ellipsoid e(eid);
if (e.isNull()) {
error("Unknown ellipsoid ID");
return false;
}
Datum d(e, dx, dy, dz);
_datums.append(Entry(list[0].trimmed(), id, d));
_errorLine++;
}
return true;
_WGS84 = (*this == WGS84) ? true : false;
}
Coordinates Datum::toWGS84(const Coordinates &c) const
{
return _WGS84 ? c : molodensky(c, *this, Datum(Ellipsoid(WGS84_RADIUS,
WGS84_FLATTENING), 0.0, 0.0, 0.0));
return _WGS84 ? c : molodensky(c, *this, WGS84);
}
Coordinates Datum::fromWGS84(const Coordinates &c) const
{
return _WGS84 ? c : molodensky(c, Datum(Ellipsoid(WGS84_RADIUS,
WGS84_FLATTENING), 0.0, 0.0, 0.0), *this);
return _WGS84 ? c : molodensky(c, WGS84, *this);
}
QDebug operator<<(QDebug dbg, const Datum &datum)
{
dbg.nospace() << "Datum(" << datum.ellipsoid() << ", " << datum.dx() << ", "
<< datum.dy() << ", " << datum.dz() << ")";
return dbg.space();
return dbg.maybeSpace();
}

32
src/map/datum.h
View File

@ -10,43 +10,35 @@
class Datum
{
public:
Datum() : _ellipsoid(Ellipsoid()), _dx(NAN), _dy(NAN), _dz(NAN),
Datum() : _ellipsoid(0), _dx(NAN), _dy(NAN), _dz(NAN),
_WGS84(false) {}
Datum(const Ellipsoid &ellipsoid, double dx, double dy, double dz);
Datum(const QString &name);
Datum(int id);
Datum(const Ellipsoid *ellipsoid, double dx, double dy, double dz);
const Ellipsoid &ellipsoid() const {return _ellipsoid;}
const Ellipsoid *ellipsoid() const {return _ellipsoid;}
double dx() const {return _dx;}
double dy() const {return _dy;}
double dz() const {return _dz;}
bool isNull() const
{return (_ellipsoid.isNull() && std::isnan(_dx) && std::isnan(_dy)
{return (!_ellipsoid && std::isnan(_dx) && std::isnan(_dy)
&& std::isnan(_dz));}
bool isValid() const
{return (_ellipsoid && !std::isnan(_dx) && !std::isnan(_dy)
&& !std::isnan(_dz));}
Coordinates toWGS84(const Coordinates &c) const;
Coordinates fromWGS84(const Coordinates &c) const;
static bool loadList(const QString &path);
static const QString &errorString() {return _errorString;}
static int errorLine() {return _errorLine;}
private:
class Entry;
static QList<Entry> WGS84();
static void error(const QString &str);
Ellipsoid _ellipsoid;
const Ellipsoid *_ellipsoid;
double _dx, _dy, _dz;
bool _WGS84;
static QList<Entry> _datums;
static QString _errorString;
static int _errorLine;
};
inline bool operator==(const Datum &d1, const Datum &d2)
{return (d1.ellipsoid() == d2.ellipsoid() && d1.dx() == d2.dx()
&& d1.dy() == d2.dy() && d1.dz() == d2.dz());}
QDebug operator<<(QDebug dbg, const Datum &datum);
#endif // DATUM_H

59
src/map/ellipsoid.cpp
View File

@ -1,77 +1,74 @@
#include <QFile>
#include <QDebug>
#include "common/wgs84.h"
#include "ellipsoid.h"
QMap<int, Ellipsoid> Ellipsoid::_ellipsoids;
QString Ellipsoid::_errorString;
int Ellipsoid::_errorLine = 0;
QMap<int, Ellipsoid> Ellipsoid::_ellipsoids = WGS84();
Ellipsoid::Ellipsoid(int id)
QMap<int, Ellipsoid> Ellipsoid::WGS84()
{
QMap<int, Ellipsoid> map;
map.insert(7030, Ellipsoid(WGS84_RADIUS, WGS84_FLATTENING));
return map;
}
const Ellipsoid *Ellipsoid::ellipsoid(int id)
{
QMap<int, Ellipsoid>::const_iterator it = _ellipsoids.find(id);
if (it == _ellipsoids.end())
*this = Ellipsoid();
return 0;
else
*this = it.value();
return &(it.value());
}
void Ellipsoid::error(const QString &str)
{
_errorString = str;
_ellipsoids.clear();
}
bool Ellipsoid::loadList(const QString &path)
void Ellipsoid::loadList(const QString &path)
{
QFile file(path);
bool res;
int ln = 0;
if (!file.open(QFile::ReadOnly)) {
error(file.errorString());
return false;
qWarning("Error opening ellipsoids file: %s: %s", qPrintable(path),
qPrintable(file.errorString()));
return;
}
_errorLine = 1;
_errorString.clear();
while (!file.atEnd()) {
ln++;
QByteArray line = file.readLine();
QList<QByteArray> list = line.split(',');
if (list.size() != 4) {
error("Format error");
return false;
qWarning("%s: %d: Format error", qPrintable(path), ln);
continue;
}
int id = list[1].trimmed().toInt(&res);
if (!res) {
error("Invalid ellipsoid id");
return false;
qWarning("%s: %d: Invalid ellipsoid code", qPrintable(path), ln);
continue;
}
double radius = list[2].trimmed().toDouble(&res);
if (!res) {
error("Invalid ellipsoid radius");
return false;
qWarning("%s: %d: Invalid radius", qPrintable(path), ln);
continue;
}
double flattening = list[3].trimmed().toDouble(&res);
if (!res) {
error("Invalid ellipsoid flattening");
return false;
qWarning("%s: %d: Invalid flattening", qPrintable(path), ln);
continue;
}
Ellipsoid e(radius, 1.0/flattening);
_ellipsoids.insert(id, e);
_errorLine++;
}
return true;
}
QDebug operator<<(QDebug dbg, const Ellipsoid &ellipsoid)
{
dbg.nospace() << "Ellipsoid(" << ellipsoid.radius() << ", "
<< 1.0 / ellipsoid.flattening() << ")";
return dbg.space();
return dbg.maybeSpace();
}

16
src/map/ellipsoid.h
View File

@ -12,29 +12,29 @@ public:
Ellipsoid() : _radius(NAN), _flattening(NAN) {}
Ellipsoid(double radius, double flattening)
: _radius(radius), _flattening(flattening) {}
Ellipsoid(int id);
double radius() const {return _radius;}
double flattening() const {return _flattening;}
bool isNull() const
{return (std::isnan(_radius) && std::isnan(_flattening));}
bool isValid() const
{return !(std::isnan(_radius) || std::isnan(_flattening));}
static bool loadList(const QString &path);
static const QString &errorString() {return _errorString;}
static int errorLine() {return _errorLine;}
static const Ellipsoid *ellipsoid(int id);
static void loadList(const QString &path);
private:
static void error(const QString &str);
double _radius;
double _flattening;
static QMap<int, Ellipsoid> WGS84();
static QMap<int, Ellipsoid> _ellipsoids;
static QString _errorString;
static int _errorLine;
};
inline bool operator==(const Ellipsoid &e1, const Ellipsoid &e2)
{return (e1.radius() == e2.radius() && e1.flattening() == e2.flattening());}
QDebug operator<<(QDebug dbg, const Ellipsoid &ellipsoid);
#endif // ELLIPSOID_H

191
src/map/gcs.cpp Normal file
View File

@ -0,0 +1,191 @@
#include <QFile>
#include "common/wgs84.h"
#include "gcs.h"
class GCS::Entry {
public:
Entry(int id, int gd, const QString &name, const GCS &gcs)
: _id(id), _gd(gd), _name(name), _gcs(gcs) {}
int id() const {return _id;}
int gd() const {return _gd;}
const QString &name() const {return _name;}
const GCS &gcs() const {return _gcs;}
private:
int _id, _gd;
QString _name;
GCS _gcs;
};
static Ellipsoid WGS84e = Ellipsoid(WGS84_RADIUS, WGS84_FLATTENING);
static Datum WGS84d = Datum(&WGS84e, 0.0, 0.0, 0.0);
static int parameter(const QString &str, bool *res)
{
QString field = str.trimmed();
if (field.isEmpty()) {
*res = true;
return 0;
}
return field.toInt(res);
}
QList<GCS::Entry> GCS::_gcss = WGS84();
QList<GCS::Entry> GCS::WGS84()
{
QList<GCS::Entry> list;
list.append(GCS::Entry(4326, 6326, "WGS 84", GCS(WGS84d, 8901, 9122)));
return list;
}
const GCS *GCS::gcs(int id)
{
for (int i = 0; i < _gcss.size(); i++)
if (_gcss.at(i).id() == id)
return &(_gcss.at(i).gcs());
return 0;
}
const GCS *GCS::gcs(int geodeticDatum, int primeMeridian, int angularUnits)
{
for (int i = 0; i < _gcss.size(); i++) {
const Entry &e = _gcss.at(i);
if (e.gd() == geodeticDatum && e.gcs().primeMeridian() == primeMeridian
&& e.gcs().angularUnits() == angularUnits)
return &(e.gcs());
}
return 0;
}
const GCS *GCS::gcs(const QString &name)
{
for (int i = 0; i < _gcss.size(); i++)
if (_gcss.at(i).name() == name)
return &(_gcss.at(i).gcs());
return 0;
}
void GCS::loadList(const QString &path)
{
QFile file(path);
bool res;
int ln = 0;
const Ellipsoid *e;
if (!file.open(QFile::ReadOnly)) {
qWarning("Error opening PCS file: %s: %s", qPrintable(path),
qPrintable(file.errorString()));
return;
}
while (!file.atEnd()) {
ln++;
QByteArray line = file.readLine();
QList<QByteArray> list = line.split(',');
if (list.size() < 10) {
qWarning("%s: %d: Format error", qPrintable(path), ln);
continue;
}
int id = parameter(list[1], &res);
if (!res) {
qWarning("%s: %d: Invalid GCS code", qPrintable(path), ln);
continue;
}
int gd = parameter(list[2], &res);
if (!res) {
qWarning("%s: %d: Invalid geodetic datum code", qPrintable(path),
ln);
continue;
}
int au = list[3].trimmed().toInt(&res);
if (!res) {
qWarning("%s: %d: Invalid angular units code", qPrintable(path),
ln);
continue;
}
int el = list[4].trimmed().toInt(&res);
if (!res) {
qWarning("%s: %d: Invalid ellipsoid code", qPrintable(path), ln);
continue;
}
int pm = list[5].trimmed().toInt(&res);
if (!res) {
qWarning("%s: %d: Invalid prime meridian code", qPrintable(path),
ln);
continue;
}
int ct = list[6].trimmed().toInt(&res);
if (!res) {
qWarning("%s: %d: Invalid coordinates transformation code",
qPrintable(path), ln);
continue;
}
double dx = list[7].trimmed().toDouble(&res);
if (!res) {
qWarning("%s: %d: Invalid dx", qPrintable(path), ln);
continue;
}
double dy = list[8].trimmed().toDouble(&res);
if (!res) {
qWarning("%s: %d: Invalid dy", qPrintable(path), ln);
continue;
}
double dz = list[9].trimmed().toDouble(&res);
if (!res) {
qWarning("%s: %d: Invalid dz", qPrintable(path), ln);
continue;
}
if (!(e = Ellipsoid::ellipsoid(el))) {
qWarning("%s: %d: Unknown ellipsoid code", qPrintable(path), ln);
continue;
}
Datum datum;
switch (ct) {
case 9603:
datum = Datum(e, dx, dy, dz);
break;
default:
qWarning("%s: %d: Unknown coordinates transformation method",
qPrintable(path), ln);
continue;
}
GCS gcs(datum, pm, au);
if (gcs.isValid())
_gcss.append(Entry(id, gd, list[0].trimmed(), gcs));
else
qWarning("%s: %d: Unknown prime meridian/angular units code",
qPrintable(path), ln);
}
}
Coordinates GCS::toWGS84(const Coordinates &c) const
{
return datum().toWGS84(Coordinates(_primeMeridian.toGreenwich(c.lon()),
c.lat()));
}
Coordinates GCS::fromWGS84(const Coordinates &c) const
{
Coordinates ds(datum().fromWGS84(c));
return Coordinates(_primeMeridian.fromGreenwich(ds.lon()), ds.lat());
}
QDebug operator<<(QDebug dbg, const GCS &gcs)
{
dbg.nospace() << "GCS(" << gcs.datum() << ", " << gcs.primeMeridian()
<< ", " << gcs.angularUnits() << ")";
return dbg.maybeSpace();
}

55
src/map/gcs.h Normal file
View File

@ -0,0 +1,55 @@
#ifndef GCS_H
#define GCS_H
#include "datum.h"
#include "angularunits.h"
#include "primemeridian.h"
class GCS
{
public:
GCS() {}
GCS(const Datum &datum, const PrimeMeridian &primeMeridian,
const AngularUnits &angularUnits) : _datum(datum),
_primeMeridian(primeMeridian), _angularUnits(angularUnits) {}
const PrimeMeridian &primeMeridian() const {return _primeMeridian;}
const AngularUnits &angularUnits() const {return _angularUnits;}
const Datum &datum() const {return _datum;}
bool isNull() const {return _datum.isNull() && _angularUnits.isNull()
&& _primeMeridian.isNull();}
bool isValid() const {return _datum.isValid() && _angularUnits.isValid()
&& _primeMeridian.isValid();}
Coordinates toWGS84(const Coordinates &c) const;
Coordinates fromWGS84(const Coordinates &c) const;
static const GCS *gcs(int id);
static const GCS *gcs(int geodeticDatum, int primeMeridian,
int angularUnits);
static const GCS *gcs(const QString &name);
static void loadList(const QString &path);
private:
class Entry;
static QList<Entry> WGS84();
Datum _datum;
PrimeMeridian _primeMeridian;
AngularUnits _angularUnits;
static QList<Entry> _gcss;
};
inline bool operator==(const GCS &gcs1, const GCS &gcs2)
{return (gcs1.datum() == gcs2.datum()
&& gcs1.primeMeridian() == gcs2.primeMeridian()
&& gcs1.angularUnits() == gcs2.angularUnits());}
QDebug operator<<(QDebug dbg, const GCS &gcs);
#endif // GCS_H

116
src/map/geotiff.cpp
View File

@ -21,7 +21,8 @@
#define GTRasterTypeGeoKey 1025
#define GeographicTypeGeoKey 2048
#define GeogGeodeticDatumGeoKey 2050
#define GeogEllipsoidGeoKey 2056
#define GeogPrimeMeridianGeoKey 2051
#define GeogAngularUnitsGeoKey 2054
#define ProjectedCSTypeGeoKey 3072
#define ProjectionGeoKey 3074
#define ProjCoordTransGeoKey 3075
@ -218,7 +219,8 @@ bool GeoTIFF::readKeys(TIFFFile &file, Ctx &ctx, QMap<quint16, Value> &kv) const
case GTRasterTypeGeoKey:
case GeogGeodeticDatumGeoKey:
case ProjectionGeoKey:
case GeogEllipsoidGeoKey:
case GeogPrimeMeridianGeoKey:
case GeogAngularUnitsGeoKey:
if (entry.TIFFTagLocation != 0 || entry.Count != 1)
return false;
value.SHORT = entry.ValueOffset;
@ -260,80 +262,84 @@ bool GeoTIFF::readGeoValue(TIFFFile &file, quint32 offset, quint16 index,
return true;
}
Datum GeoTIFF::datum(QMap<quint16, Value> &kv)
const GCS *GeoTIFF::gcs(QMap<quint16, Value> &kv)
{
Datum datum;
const GCS *gcs;
if (IS_SET(kv, GeographicTypeGeoKey)) {
datum = Datum(kv.value(GeographicTypeGeoKey).SHORT);
if (datum.isNull())
if (!(gcs = GCS::gcs(kv.value(GeographicTypeGeoKey).SHORT)))
_errorString = QString("%1: unknown GCS")
.arg(kv.value(GeographicTypeGeoKey).SHORT);
} else if (IS_SET(kv, GeogGeodeticDatumGeoKey)) {
datum = Datum(kv.value(GeogGeodeticDatumGeoKey).SHORT - 2000);
if (datum.isNull())
_errorString = QString("%1: unknown geodetic datum")
.arg(kv.value(GeogGeodeticDatumGeoKey).SHORT);
} else if (IS_SET(kv, GeogEllipsoidGeoKey)
&& kv.value(GeogEllipsoidGeoKey).SHORT == 7019) {
datum = Datum(4326);
} else
_errorString = "Missing datum";
int pm = IS_SET(kv, GeogPrimeMeridianGeoKey)
? kv.value(GeogPrimeMeridianGeoKey).SHORT : 8901;
int au = IS_SET(kv, GeogAngularUnitsGeoKey)
? kv.value(GeogAngularUnitsGeoKey).SHORT : 9122;
return datum;
if (!(gcs = GCS::gcs(kv.value(GeogGeodeticDatumGeoKey).SHORT, pm, au)))
_errorString = QString("%1+%2: unknown geodetic datum + prime"
" meridian combination")
.arg(kv.value(GeogGeodeticDatumGeoKey).SHORT)
.arg(kv.value(GeogPrimeMeridianGeoKey).SHORT);
} else
_errorString = "Can not determine GCS";
return gcs;
}
Projection::Method GeoTIFF::method(QMap<quint16, Value> &kv)
{
int epsg;
int table[] = {0, 9807, 0, 0, 0, 0, 0, 1024, 9801, 9802, 9820,
9822};
if (!IS_SET(kv, ProjCoordTransGeoKey)) {
_errorString = "Missing coordinate transformation method";
return Projection::Method();
}
quint16 index = kv.value(ProjCoordTransGeoKey).SHORT;
epsg = (index >= ARRAY_SIZE(table)) ? 0 : table[index];
if (!epsg) {
_errorString = QString("Unknown coordinate transformation method");
return Projection::Method();
switch (kv.value(ProjCoordTransGeoKey).SHORT) {
case 1:
return Projection::Method(9807);
case 7:
return Projection::Method(1024);
case 8:
return Projection::Method(9801);
case 9:
return Projection::Method(9802);
case 10:
return Projection::Method(9820);
case 11:
return Projection::Method(9822);
default:
_errorString = QString("%1: unknown coordinate transformation method")
.arg(kv.value(ProjCoordTransGeoKey).SHORT);
return Projection::Method();
}
return Projection::Method(epsg);
}
bool GeoTIFF::projectedModel(QMap<quint16, Value> &kv)
{
PCS pcs;
if (IS_SET(kv, ProjectedCSTypeGeoKey)) {
pcs = PCS(kv.value(ProjectedCSTypeGeoKey).SHORT);
if (pcs.isNull())
const PCS *pcs;
if (!(pcs = PCS::pcs(kv.value(ProjectedCSTypeGeoKey).SHORT))) {
_errorString = QString("%1: unknown PCS")
.arg(kv.value(ProjectedCSTypeGeoKey).SHORT);
} else if (IS_SET(kv, GeographicTypeGeoKey)
&& IS_SET(kv, ProjectionGeoKey)) {
pcs = PCS(kv.value(GeographicTypeGeoKey).SHORT,
kv.value(ProjectionGeoKey).SHORT);
if (pcs.isNull())
_errorString = QString("%1+%2: unknown GCS+projection combination")
return false;
}
_gcs = pcs->gcs();
_projection = Projection::projection(_gcs->datum(), pcs->method(),
pcs->setup());
} else if (IS_SET(kv, ProjectionGeoKey)) {
if (!(_gcs = gcs(kv)))
return false;
const PCS *pcs;
if (!(pcs = PCS::pcs(_gcs, kv.value(ProjectionGeoKey).SHORT))) {
_errorString = QString("%1: unknown projection code")
.arg(kv.value(GeographicTypeGeoKey).SHORT)
.arg(kv.value(ProjectionGeoKey).SHORT);
} else if (IS_SET(kv, GeogGeodeticDatumGeoKey)
&& IS_SET(kv, ProjectionGeoKey)) {
pcs = PCS(kv.value(GeogGeodeticDatumGeoKey).SHORT - 2000,
kv.value(ProjectionGeoKey).SHORT);
if (pcs.isNull())
_errorString =
QString("%1+%2: unknown geodetic datum+projection combination")
.arg(kv.value(GeogGeodeticDatumGeoKey).SHORT)
.arg(kv.value(ProjectionGeoKey).SHORT);
return false;
}
_projection = Projection::projection(_gcs->datum(), pcs->method(),
pcs->setup());
} else {
Datum d = datum(kv);
if (d.isNull())
if (!(_gcs = gcs(kv)))
return false;
Projection::Method m = method(kv);
if (m.isNull())
@ -348,24 +354,18 @@ bool GeoTIFF::projectedModel(QMap<quint16, Value> &kv)
kv.value(ProjStdParallel1GeoKey).DOUBLE,
kv.value(ProjStdParallel2GeoKey).DOUBLE
);
pcs = PCS(d, m, setup);
_projection = Projection::projection(_gcs->datum(), m, setup);
}
_datum = pcs.datum();
_projection = Projection::projection(pcs.datum(), pcs.method(),
pcs.setup());
return true;
}
bool GeoTIFF::geographicModel(QMap<quint16, Value> &kv)
{
_datum = datum(kv);
if (_datum.isNull())
if (!(_gcs = gcs(kv)))
return false;
_projection = new LatLon();
_projection = new LatLon(_gcs->angularUnits());
return true;
}

8
src/map/geotiff.h
View File

@ -4,7 +4,7 @@
#include <QTransform>
#include <QFile>
#include <QMap>
#include "datum.h"
#include "gcs.h"
#include "projection.h"
#include "tifffile.h"
#include "transform.h"
@ -16,7 +16,7 @@ public:
bool load(const QString &path);
const Datum &datum() const {return _datum;}
const GCS *gcs() const {return _gcs;}
Projection *projection() const {return _projection;}
const QTransform &transform() const {return _transform;}
@ -50,13 +50,13 @@ private:
bool readGeoValue(TIFFFile &file, quint32 offset, quint16 index,
double &val) const;
Datum datum(QMap<quint16, Value> &kv);
const GCS *gcs(QMap<quint16, Value> &kv);
Projection::Method method(QMap<quint16, Value> &kv);
bool geographicModel(QMap<quint16, Value> &kv);
bool projectedModel(QMap<quint16, Value> &kv);
QTransform _transform;
Datum _datum;
const GCS *_gcs;
Projection *_projection;
QString _errorString;

8
src/map/lambertazimuthal.cpp
View File

@ -48,18 +48,18 @@ Defense.
#define M_PI_2 1.57079632679489661923
#endif // M_PI_2
LambertAzimuthal::LambertAzimuthal(const Ellipsoid &ellipsoid,
LambertAzimuthal::LambertAzimuthal(const Ellipsoid *ellipsoid,
double latitudeOrigin, double longitudeOrigin, double falseEasting,
double falseNorthing)
{
double es2, es4, es6;
es2 = 2 * ellipsoid.flattening() - ellipsoid.flattening()
* ellipsoid.flattening();
es2 = 2 * ellipsoid->flattening() - ellipsoid->flattening()
* ellipsoid->flattening();
es4 = es2 * es2;
es6 = es4 * es2;
_ra = ellipsoid.radius() * (1.0 - es2 / 6.0 - 17.0 * es4 / 360.0 - 67.0
_ra = ellipsoid->radius() * (1.0 - es2 / 6.0 - 17.0 * es4 / 360.0 - 67.0
* es6 / 3024.0);
_latOrigin = deg2rad(latitudeOrigin);
_sinLatOrigin = sin(_latOrigin);

2
src/map/lambertazimuthal.h
View File

@ -7,7 +7,7 @@
class LambertAzimuthal : public Projection
{
public:
LambertAzimuthal(const Ellipsoid &ellipsoid, double latitudeOrigin,
LambertAzimuthal(const Ellipsoid *ellipsoid, double latitudeOrigin,
double longitudeOrigin, double falseEasting, double falseNorthing);
virtual QPointF ll2xy(const Coordinates &c) const;

16
src/map/lambertconic.cpp
View File

@ -58,7 +58,7 @@ Defense.
(tan(M_PI_4 - lat / 2) / pow((1.0 - essin) / (1.0 + essin), es_over_2))
LambertConic1::LambertConic1(const Ellipsoid &ellipsoid, double latitudeOrigin,
LambertConic1::LambertConic1(const Ellipsoid *ellipsoid, double latitudeOrigin,
double longitudeOrigin, double scale, double falseEasting,
double falseNorthing)
{
@ -76,8 +76,8 @@ LambertConic1::LambertConic1(const Ellipsoid &ellipsoid, double latitudeOrigin,
_falseEasting = falseEasting;
_falseNorthing = falseNorthing;
es2 = 2.0 * ellipsoid.flattening() - ellipsoid.flattening()
* ellipsoid.flattening();
es2 = 2.0 * ellipsoid->flattening() - ellipsoid->flattening()
* ellipsoid->flattening();
_es = sqrt(es2);
_es_over_2 = _es / 2.0;
@ -87,7 +87,7 @@ LambertConic1::LambertConic1(const Ellipsoid &ellipsoid, double latitudeOrigin,
m0 = LAMBERT_m(cos(lat_orig), es_sin);
_t0 = LAMBERT1_t(lat_orig, es_sin, _es_over_2);
_rho0 = ellipsoid.radius() * scale * m0 / _n;
_rho0 = ellipsoid->radius() * scale * m0 / _n;
_rho_olat = _rho0;
}
@ -202,7 +202,7 @@ Coordinates LambertConic1::xy2ll(const QPointF &p) const
return Coordinates(rad2deg(lon), rad2deg(lat));
}
LambertConic2::LambertConic2(const Ellipsoid &ellipsoid,
LambertConic2::LambertConic2(const Ellipsoid *ellipsoid,
double standardParallel1, double standardParallel2, double latitudeOrigin,
double longitudeOrigin, double falseEasting, double falseNorthing)
{
@ -226,8 +226,8 @@ LambertConic2::LambertConic2(const Ellipsoid &ellipsoid,
sp2 = deg2rad(standardParallel2);
if (fabs(sp1 - sp2) > 1.0e-10) {
es2 = 2 * ellipsoid.flattening() - ellipsoid.flattening()
* ellipsoid.flattening();
es2 = 2 * ellipsoid->flattening() - ellipsoid->flattening()
* ellipsoid->flattening();
es = sqrt(es2);
es_over_2 = es / 2.0;
@ -252,7 +252,7 @@ LambertConic2::LambertConic2(const Ellipsoid &ellipsoid,
k0 = (m1 / m0) * (pow(t0 / t1, n));
const_value = ((ellipsoid.radius() * m2) / (n * pow(t2, n)));
const_value = ((ellipsoid->radius() * m2) / (n * pow(t2, n)));
falseNorthing += (const_value * pow(t_olat, n)) - (const_value
* pow(t0, n));

4
src/map/lambertconic.h
View File

@ -8,7 +8,7 @@ class LambertConic1 : public Projection
{
public:
LambertConic1() {}
LambertConic1(const Ellipsoid &ellipsoid, double latitudeOrigin,
LambertConic1(const Ellipsoid *ellipsoid, double latitudeOrigin,
double longitudeOrigin, double scale, double falseEasting,
double falseNorthing);
@ -31,7 +31,7 @@ private:
class LambertConic2 : public Projection
{
public:
LambertConic2(const Ellipsoid &ellipsoid, double standardParallel1,
LambertConic2(const Ellipsoid *ellipsoid, double standardParallel1,
double standardParallel2, double latitudeOrigin, double longitudeOrigin,
double falseEasting, double falseNorthing);

10
src/map/latlon.h
View File

@ -2,14 +2,20 @@
#define LATLON_H
#include "projection.h"
#include "angularunits.h"
class LatLon : public Projection
{
public:
LatLon(const AngularUnits &au) : _au(au) {}
virtual QPointF ll2xy(const Coordinates &c) const
{return QPointF(c.lon(), c.lat());}
{return QPointF(_au.fromDegrees(c.lon()), _au.fromDegrees(c.lat()));}
virtual Coordinates xy2ll(const QPointF &p) const
{return Coordinates(p.x(), p.y());}
{return Coordinates(_au.toDegrees(p.x()), _au.toDegrees(p.y()));}
private:
AngularUnits _au;
};
#endif // LATLON_H

63
src/map/mapfile.cpp
View File

@ -1,5 +1,6 @@
#include "latlon.h"
#include "utm.h"
#include "gcs.h"
#include "mapfile.h"
@ -151,8 +152,7 @@ bool MapFile::parseMapFile(QIODevice &device, QList<CalibrationPoint> &points,
bool MapFile::createDatum(const QString &datum)
{
_datum = Datum(datum);
if (_datum.isNull()) {
if (!(_gcs = GCS::gcs(datum))) {
_errorString = QString("%1: Unknown datum").arg(datum);
return false;
}
@ -164,17 +164,17 @@ bool MapFile::createProjection(const QString &name,
const Projection::Setup &setup, QList<CalibrationPoint> &points)
{
if (name == "Mercator")
_projection = Projection::projection(_datum, 1024, setup);
_projection = Projection::projection(_gcs->datum(), 1024, setup);
else if (name == "Transverse Mercator")
_projection = Projection::projection(_datum, 9807, setup);
_projection = Projection::projection(_gcs->datum(), 9807, setup);
else if (name == "Latitude/Longitude")
_projection = new LatLon();
_projection = new LatLon(9102);
else if (name == "Lambert Conformal Conic")
_projection = Projection::projection(_datum, 9802, setup);
_projection = Projection::projection(_gcs->datum(), 9802, setup);
else if (name == "Albers Equal Area")
_projection = Projection::projection(_datum, 9822, setup);
_projection = Projection::projection(_gcs->datum(), 9822, setup);
else if (name == "(A)Lambert Azimuthual Equal Area")
_projection = Projection::projection(_datum, 9820, setup);
_projection = Projection::projection(_gcs->datum(), 9820, setup);
else if (name == "(UTM) Universal Transverse Mercator") {
int zone;
if (points.first().zone)
@ -185,37 +185,42 @@ bool MapFile::createProjection(const QString &name,
_errorString = "Can not determine UTM zone";
return 0;
}
_projection = Projection::projection(_datum, 9807, UTM::setup(zone));
_projection = Projection::projection(_gcs->datum(), 9807,
UTM::setup(zone));
} else if (name == "(NZTM2) New Zealand TM 2000")
_projection = Projection::projection(_datum, 9807, Projection::Setup(
0, 173.0, 0.9996, 1600000, 10000000, NAN, NAN));
_projection = Projection::projection(_gcs->datum(), 9807,
Projection::Setup(0, 173.0, 0.9996, 1600000, 10000000, NAN, NAN));
else if (name == "(BNG) British National Grid")
_projection = Projection::projection(_datum, 9807, Projection::Setup(
49, -2, 0.999601, 400000, -100000, NAN, NAN));
_projection = Projection::projection(_gcs->datum(), 9807,
Projection::Setup(49, -2, 0.999601, 400000, -100000, NAN, NAN));
else if (name == "(IG) Irish Grid")
_projection = Projection::projection(_datum, 9807, Projection::Setup(
53.5, -8, 1.000035, 200000, 250000, NAN, NAN));
_projection = Projection::projection(_gcs->datum(), 9807,
Projection::Setup(53.5, -8, 1.000035, 200000, 250000, NAN, NAN));
else if (name == "(SG) Swedish Grid")
_projection = Projection::projection(_datum, 9807, Projection::Setup(
0, 15.808278, 1, 1500000, 0, NAN, NAN));
_projection = Projection::projection(_gcs->datum(), 9807,
Projection::Setup(0, 15.808278, 1, 1500000, 0, NAN, NAN));
else if (name == "(I) France Zone I")
_projection = Projection::projection(_datum, 9802, Projection::Setup(
49.5, 2.337229, NAN, 600000, 1200000, 48.598523, 50.395912));
_projection = Projection::projection(_gcs->datum(), 9802,
Projection::Setup(49.5, 2.337229, NAN, 600000, 1200000, 48.598523,
50.395912));
else if (name == "(II) France Zone II")
_projection = Projection::projection(_datum, 9802, Projection::Setup(
46.8, 2.337229, NAN, 600000, 2200000, 45.898919, 47.696014));
_projection = Projection::projection(_gcs->datum(), 9802,
Projection::Setup(46.8, 2.337229, NAN, 600000, 2200000, 45.898919,
47.696014));
else if (name == "(III) France Zone III")
_projection = Projection::projection(_datum, 9802, Projection::Setup(
44.1, 2.337229, NAN, 600000, 3200000, 43.199291, 44.996094));
_projection = Projection::projection(_gcs->datum(), 9802,
Projection::Setup(44.1, 2.337229, NAN, 600000, 3200000, 43.199291,
44.996094));
else if (name == "(IV) France Zone IV")
_projection = Projection::projection(_datum, 9802, Projection::Setup(
42.165, 2.337229, NAN, 234.358, 4185861.369, 41.560388, 42.767663));
_projection = Projection::projection(_gcs->datum(), 9802,
Projection::Setup(42.165, 2.337229, NAN, 234.358, 4185861.369,
41.560388, 42.767663));
else if (name == "(VICGRID) Victoria Australia")
_projection = Projection::projection(_datum, 9802, Projection::Setup(
-37, 145, NAN, 2500000, 4500000, -36, -38));
_projection = Projection::projection(_gcs->datum(), 9802,
Projection::Setup(-37, 145, NAN, 2500000, 4500000, -36, -38));
else if (name == "(VG94) VICGRID94 Victoria Australia")
_projection = Projection::projection(_datum, 9802, Projection::Setup(
-37, 145, NAN, 2500000, 2500000, -36, -38));
_projection = Projection::projection(_gcs->datum(), 9802,
Projection::Setup(-37, 145, NAN, 2500000, 2500000, -36, -38));
else {
_errorString = QString("%1: Unknown map projection").arg(name);
return false;

8
src/map/mapfile.h
View File

@ -3,7 +3,7 @@
#include <QIODevice>
#include <QTransform>
#include "datum.h"
#include "gcs.h"
#include "transform.h"
#include "projection.h"
@ -12,7 +12,7 @@ class MapFile
public:
bool load(QIODevice &file);
const Datum &datum() const {return _datum;}
const GCS *gcs() const {return _gcs;}
Projection *projection() const {return _projection;}
const QTransform &transform() const {return _transform;}
@ -42,9 +42,9 @@ private:
QString _image;
QSize _size;
Datum _datum;
QTransform _transform;
const GCS *_gcs;
Projection *_projection;
QTransform _transform;
QString _errorString;
};

16
src/map/offlinemap.cpp
View File

@ -138,7 +138,7 @@ OfflineMap::OfflineMap(const QString &fileName, QObject *parent)
_name = mf.name();
_size = mf.size();
_imgPath = mf.image();
_datum = mf.datum();
_gcs = mf.gcs();
_projection = mf.projection();
_transform = mf.transform();
}
@ -152,7 +152,7 @@ OfflineMap::OfflineMap(const QString &fileName, QObject *parent)
_name = mf.name();
_size = mf.size();
_imgPath = mf.image();
_datum = mf.datum();
_gcs = mf.gcs();
_projection = mf.projection();
_transform = mf.transform();
}
@ -164,7 +164,7 @@ OfflineMap::OfflineMap(const QString &fileName, QObject *parent)
} else {
_name = fi.fileName();
_imgPath = fileName;
_datum = gt.datum();
_gcs = gt.gcs();
_projection = gt.projection();
_transform = gt.transform();
}
@ -230,7 +230,7 @@ OfflineMap::OfflineMap(const QString &fileName, Tar &tar, QObject *parent)
_name = mf.name();
_size = mf.size();
_datum = mf.datum();
_gcs = mf.gcs();
_projection = mf.projection();
_transform = mf.transform();
@ -363,19 +363,19 @@ void OfflineMap::draw(QPainter *painter, const QRectF &rect)
QPointF OfflineMap::ll2xy(const Coordinates &c)
{
if (_ozf.isOpen()) {
QPointF p(_transform.map(_projection->ll2xy(_datum.fromWGS84(c))));
QPointF p(_transform.map(_projection->ll2xy(_gcs->fromWGS84(c))));
return QPointF(p.x() * _scale.x(), p.y() * _scale.y());
} else
return _transform.map(_projection->ll2xy(_datum.fromWGS84(c)));
return _transform.map(_projection->ll2xy(_gcs->fromWGS84(c)));
}
Coordinates OfflineMap::xy2ll(const QPointF &p)
{
if (_ozf.isOpen()) {
return _datum.toWGS84(_projection->xy2ll(_inverted.map(QPointF(p.x()
return _gcs->toWGS84(_projection->xy2ll(_inverted.map(QPointF(p.x()
/ _scale.x(), p.y() / _scale.y()))));
} else
return _datum.toWGS84(_projection->xy2ll(_inverted.map(p)));
return _gcs->toWGS84(_projection->xy2ll(_inverted.map(p)));
}
QRectF OfflineMap::bounds() const

4
src/map/offlinemap.h
View File

@ -3,7 +3,7 @@
#include <QTransform>
#include "common/coordinates.h"
#include "datum.h"
#include "gcs.h"
#include "projection.h"
#include "transform.h"
#include "map.h"
@ -66,7 +66,7 @@ private:
QString _name;
QSize _size;
Datum _datum;
const GCS *_gcs;
Projection *_projection;
QTransform _transform, _inverted;

161
src/map/pcs.cpp
View File

@ -4,24 +4,18 @@
class PCS::Entry {
public:
Entry(int id, int gcs, int proj, const PCS &pcs)
: _id(id), _gcs(gcs), _proj(proj), _pcs(pcs) {}
Entry(int id, int proj, const PCS &pcs) : _id(id), _proj(proj), _pcs(pcs) {}
int id() const {return _id;}
int gcs() const {return _gcs;}
int proj() const {return _proj;}
const PCS &pcs() const {return _pcs;}
private:
int _id;
int _gcs;
int _proj;
int _id, _proj;
PCS _pcs;
};
QList<PCS::Entry> PCS::_pcss;
QString PCS::_errorString;
int PCS::_errorLine;
static double parameter(const QString &str, bool *res)
{
@ -34,99 +28,114 @@ static double parameter(const QString &str, bool *res)
return field.toDouble(res);
}
PCS::PCS(int id)
{
for (int i = 0; i < _pcss.size(); i++) {
if (_pcss.at(i).id() == id) {
*this = _pcss.at(i).pcs();
return;
}
}
*this = PCS();
const PCS *PCS::pcs(int id)
{
for (int i = 0; i < _pcss.size(); i++)
if (_pcss.at(i).id() == id)
return &(_pcss.at(i).pcs());
return 0;
}
PCS::PCS(int gcs, int proj)
const PCS *PCS::pcs(const GCS *gcs, int proj)
{
for (int i = 0; i < _pcss.size(); i++) {
if (_pcss.at(i).gcs() == gcs && _pcss.at(i).proj() == proj) {
*this = _pcss.at(i).pcs();
return;
}
}
for (int i = 0; i < _pcss.size(); i++)
if (_pcss.at(i).proj() == proj && *(_pcss.at(i).pcs().gcs()) == *gcs)
return &(_pcss.at(i).pcs());
*this = PCS();
return 0;
}
void PCS::error(const QString &str)
{
_errorString = str;
_pcss.clear();
}
bool PCS::loadList(const QString &path)
void PCS::loadList(const QString &path)
{
QFile file(path);
bool res[12];
int id, gcs, proj, transform;
bool res;
int ln = 0;
if (!file.open(QFile::ReadOnly)) {
error(file.errorString());
return false;
qWarning("Error opening PCS file: %s: %s", qPrintable(path),
qPrintable(file.errorString()));
return;
}
_errorLine = 1;
_errorString.clear();
while (!file.atEnd()) {
ln++;
QByteArray line = file.readLine();
QList<QByteArray> list = line.split(',');
if (list.size() != 12) {
error("Format error");
return false;
qWarning("%s: %d: Format error", qPrintable(path), ln);
continue;
}
id = list[1].trimmed().toInt(&res[1]);;
gcs = list[2].trimmed().toInt(&res[2]);
proj = list[3].trimmed().toInt(&res[3]);
transform = list[4].trimmed().toInt(&res[4]);
Projection::Setup setup(
parameter(list[5], &res[5]), parameter(list[6], &res[6]),
parameter(list[7], &res[7]), parameter(list[8], &res[8]),
parameter(list[9], &res[9]), parameter(list[10], &res[10]),
parameter(list[11], &res[11]));
for (int i = 1; i < 12; i++) {
if (!res[i]) {
error("Parse error");
return false;
}
int id = list[1].trimmed().toInt(&res);
if (!res) {
qWarning("%s: %d: Invalid PCS code", qPrintable(path), ln);
continue;
}
int gcs = list[2].trimmed().toInt(&res);
if (!res) {
qWarning("%s: %d: Invalid GCS code", qPrintable(path), ln);
continue;
}
int proj = list[3].trimmed().toInt(&res);
if (!res) {
qWarning("%s: %d: Invalid projection code", qPrintable(path), ln);
continue;
}
int transform = list[4].trimmed().toInt(&res);
if (!res) {
qWarning("%s: %d: Invalid coordinate transformation code",
qPrintable(path), ln);
continue;
}
double lat0 = parameter(list[5], &res);
if (!res) {
qWarning("%s: %d: Invalid latitude origin", qPrintable(path), ln);
continue;
}
double lon0 = parameter(list[6], &res);
if (!res) {
qWarning("%s: %d: Invalid longitude origin", qPrintable(path), ln);
continue;
}
double scale = parameter(list[7], &res);
if (!res) {
qWarning("%s: %d: Invalid scale", qPrintable(path), ln);
continue;
}
double fe = parameter(list[8], &res);
if (!res) {
qWarning("%s: %d: Invalid false easting", qPrintable(path), ln);
continue;
}
double fn = parameter(list[9], &res);
if (!res) {
qWarning("%s: %d: Invalid false northing", qPrintable(path), ln);
continue;
}
double sp1 = parameter(list[10], &res);
if (!res) {
qWarning("%s: %d: Invalid standard parallel #1", qPrintable(path),
ln);
continue;
}
double sp2 = parameter(list[11], &res);
if (!res) {
qWarning("%s: %d: Invalid standard parallel #2", qPrintable(path),
ln);
continue;
}
Datum datum(gcs);
if (datum.isNull()) {
error("Unknown datum");
return false;
}
Projection::Method method(transform);
if (method.isNull()) {
error("Unknown coordinates transformation method");
return false;
}
_pcss.append(Entry(id, gcs, proj, PCS(datum, method, setup)));
_errorLine++;
_pcss.append(Entry(id, proj, PCS(GCS::gcs(gcs), transform,
Projection::Setup(lat0, lon0, scale, fe, fn, sp1, sp2))));
}
return true;
}
QDebug operator<<(QDebug dbg, const PCS &pcs)
{
dbg.nospace() << "PCS(" << pcs.datum() << ", " << pcs.method()
dbg.nospace() << "PCS(" << *pcs.gcs() << ", " << pcs.method()
<< ", " << pcs.setup() << ")";
return dbg.space();
return dbg.maybeSpace();
}

31
src/map/pcs.h
View File

@ -3,44 +3,37 @@
#include <QDebug>
#include <QList>
#include "datum.h"
#include "gcs.h"
#include "projection.h"
class PCS
{
public:
PCS() {}
PCS(const Datum &datum, const Projection::Method &method,
const Projection::Setup &setup)
: _datum(datum), _method(method), _setup(setup) {}
PCS(int id);
PCS(int gcs, int proj);
PCS() : _gcs(0) {}
PCS(const GCS *gcs, const Projection::Method &m, const Projection::Setup &s)
: _gcs(gcs), _method(m), _setup(s) {}
PCS(const GCS *gcs, int proj);
const Datum &datum() const {return _datum;}
const GCS *gcs() const {return _gcs;}
const Projection::Method &method() const {return _method;}
const Projection::Setup &setup() const {return _setup;}
bool isNull() const
{return (_datum.isNull() && _method.isNull());}
bool isValid() const
{return !(_datum.isNull() || _method.isNull());}
{return (_gcs->isNull() && _method.isNull() && _setup.isNull());}
static bool loadList(const QString &path);
static const QString &errorString() {return _errorString;}
static int errorLine() {return _errorLine;}
static void loadList(const QString &path);
static const PCS *pcs(int id);
static const PCS *pcs(const GCS *gcs, int proj);
private:
class Entry;
static void error(const QString &str);
Datum _datum;
const GCS *_gcs;
Projection::Method _method;
Projection::Setup _setup;
static QList<Entry> _pcss;
static QString _errorString;
static int _errorLine;
static GCS _nullGCS;
};
QDebug operator<<(QDebug dbg, const PCS &pcs);

73
src/map/primemeridian.cpp Normal file
View File

@ -0,0 +1,73 @@
#include "primemeridian.h"
static double shift(double lon, double offset)
{
double ret = lon + offset;
if (ret > 180.0)
ret -= 360.0;
if (ret < -180)
ret += 360.0;
return ret;
}
PrimeMeridian::PrimeMeridian(int code)
{
switch (code) {
case 8901:
_pm = 0.0;
break;
case 8902:
_pm = -9.1319061111;
break;
case 8903:
_pm = 2.3372291666;
break;
case 8904:
_pm = -74.0809166666;
break;
case 8905:
_pm = -3.6879388888;
break;
case 8906:
_pm = 12.4523333333;
break;
case 8907:
_pm = 7.4395833333;
break;
case 8908:
_pm = 106.8077194444;
break;
case 8909:
_pm = -17.6666666666;
break;
case 8910:
_pm = 4.3679750000;
break;
case 8911:
_pm = 18.0582777777;
break;
case 8913:
_pm = 10.7229166666;
break;
default:
_pm = NAN;
}
}
double PrimeMeridian::toGreenwich(double val) const
{
return shift(val, _pm);
}
double PrimeMeridian::fromGreenwich(double val) const
{
return shift(val, -_pm);
}
QDebug operator<<(QDebug dbg, const PrimeMeridian &pm)
{
dbg.nospace() << "PrimeMeridian(" << pm._pm << ")";
return dbg.maybeSpace();
}

31
src/map/primemeridian.h Normal file
View File

@ -0,0 +1,31 @@
#ifndef PRIMEMERIDIAN_H
#define PRIMEMERIDIAN_H
#include <cmath>
#include <QDebug>
class PrimeMeridian
{
public:
PrimeMeridian() : _pm(NAN) {}
PrimeMeridian(int code);
bool isNull() const {return std::isnan(_pm);}
bool isValid() const {return !std::isnan(_pm);}
double toGreenwich(double val) const;
double fromGreenwich(double val) const;
friend bool operator==(const PrimeMeridian &pm1, const PrimeMeridian &pm2);
friend QDebug operator<<(QDebug dbg, const PrimeMeridian &pm);
private:
double _pm;
};
inline bool operator==(const PrimeMeridian &pm1, const PrimeMeridian &pm2)
{return (pm1._pm == pm2._pm);}
QDebug operator<<(QDebug dbg, const PrimeMeridian &pm);
#endif // PRIMEMERIDIAN_H

6
src/map/projection.cpp
View File

@ -27,7 +27,7 @@ Projection::Method::Method(int id)
Projection *Projection::projection(const Datum &datum, const Method &method,
const Setup &setup)
{
const Ellipsoid &ellipsoid = datum.ellipsoid();
const Ellipsoid *ellipsoid = datum.ellipsoid();
switch (method.id()) {
case 9807:
@ -66,11 +66,11 @@ QDebug operator<<(QDebug dbg, const Projection::Setup &setup)
<< setup.longitudeOrigin() << ", " << setup.scale() << ", "
<< setup.falseEasting() << ", " << setup.falseNorthing() << ", "
<< setup.standardParallel1() << ", " << setup.standardParallel2() << ")";
return dbg.space();
return dbg.maybeSpace();
}
QDebug operator<<(QDebug dbg, const Projection::Method &method)
{
dbg.nospace() << "Method(" << method.id() << ")";
return dbg.space();
return dbg.maybeSpace();
}

5
src/map/projection.h
View File

@ -30,6 +30,11 @@ public:
double standardParallel1() const {return _standardParallel1;}
double standardParallel2() const {return _standardParallel2;}
bool isNull() const {return std::isnan(_latitudeOrigin)
&& std::isnan(_longitudeOrigin) && std::isnan(_scale)
&& std::isnan(_falseEasting) && std::isnan(_falseNorthing)
&& std::isnan(_standardParallel1) && std::isnan(_standardParallel2);}
private:
double _latitudeOrigin;
double _longitudeOrigin;

22
src/map/transversemercator.cpp
View File

@ -47,17 +47,17 @@ Defense.
#define SPHSN(lat) \
((double)(_e.radius() / sqrt(1.e0 - _es * pow(sin(lat), 2))))
((double)(_e->radius() / sqrt(1.e0 - _es * pow(sin(lat), 2))))
#define SPHTMD(lat) \
((double)(_ap * lat - _bp * sin(2.e0 * lat) + _cp * sin(4.e0 * lat) \
- _dp * sin(6.e0 * lat) + _ep * sin(8.e0 * lat)))
#define DENOM(lat) \
((double)(sqrt(1.e0 - _es * pow(sin(lat),2))))
#define SPHSR(lat) \
((double)(_e.radius() * (1.e0 - _es) / pow(DENOM(lat), 3)))
((double)(_e->radius() * (1.e0 - _es) / pow(DENOM(lat), 3)))
TransverseMercator::TransverseMercator(const Ellipsoid &ellipsoid,
TransverseMercator::TransverseMercator(const Ellipsoid *ellipsoid,
double latitudeOrigin, double longitudeOrigin, double scale,
double falseEasting, double falseNorthing)
{
@ -72,24 +72,24 @@ TransverseMercator::TransverseMercator(const Ellipsoid &ellipsoid,
_falseEasting = falseEasting;
_falseNorthing = falseNorthing;
_es = 2 * _e.flattening() - _e.flattening() * _e.flattening();
_es = 2 * _e->flattening() - _e->flattening() * _e->flattening();
_ebs = (1 / (1 - _es)) - 1;
b = _e.radius() * (1 - _e.flattening());
b = _e->radius() * (1 - _e->flattening());
tn = (_e.radius() - b) / (_e.radius() + b);
tn = (_e->radius() - b) / (_e->radius() + b);
tn2 = tn * tn;
tn3 = tn2 * tn;
tn4 = tn3 * tn;
tn5 = tn4 * tn;
_ap = _e.radius() * (1.e0 - tn + 5.e0 * (tn2 - tn3) / 4.e0 + 81.e0
_ap = _e->radius() * (1.e0 - tn + 5.e0 * (tn2 - tn3) / 4.e0 + 81.e0
* (tn4 - tn5) / 64.e0);
_bp = 3.e0 * _e.radius() * (tn - tn2 + 7.e0 * (tn3 - tn4) / 8.e0 + 55.e0
_bp = 3.e0 * _e->radius() * (tn - tn2 + 7.e0 * (tn3 - tn4) / 8.e0 + 55.e0
* tn5 / 64.e0 ) / 2.e0;
_cp = 15.e0 * _e.radius() * (tn2 - tn3 + 3.e0 * (tn4 - tn5 ) / 4.e0) / 16.0;
_dp = 35.e0 * _e.radius() * (tn3 - tn4 + 11.e0 * tn5 / 16.e0) / 48.e0;
_ep = 315.e0 * _e.radius() * (tn4 - tn5) / 512.e0;
_cp = 15.e0 * _e->radius() * (tn2 - tn3 + 3.e0 * (tn4 - tn5 ) / 4.e0) / 16.0;
_dp = 35.e0 * _e->radius() * (tn3 - tn4 + 11.e0 * tn5 / 16.e0) / 48.e0;
_ep = 315.e0 * _e->radius() * (tn4 - tn5) / 512.e0;
}
QPointF TransverseMercator::ll2xy(const Coordinates &c) const

4
src/map/transversemercator.h
View File

@ -7,7 +7,7 @@
class TransverseMercator : public Projection
{
public:
TransverseMercator(const Ellipsoid &ellipsoid, double latitudeOrigin,
TransverseMercator(const Ellipsoid *ellipsoid, double latitudeOrigin,
double longitudeOrigin, double scale, double falseEasting,
double falseNorthing);
@ -15,7 +15,7 @@ public:
virtual Coordinates xy2ll(const QPointF &p) const;
private:
Ellipsoid _e;
const Ellipsoid *_e;
double _longitudeOrigin;
double _latitudeOrigin;
double _scale;