1
0
mirror of https://github.com/tumic0/GPXSee.git synced 2025-06-27 11:39:16 +02:00

Switched from hardcoded values to CSV files for datum definitions.

Non-WGS84 maps now fit correctly when defined using map grid.
This commit is contained in:
2017-04-17 18:03:04 +02:00
parent e26fa92ce6
commit 4de3ddf71b
11 changed files with 459 additions and 80 deletions

View File

@ -13,6 +13,8 @@
#define FONT_SIZE 12
#define SCREEN_DPI 96.0
#define ELLIPSOID_FILE QString("ellipsoids.csv")
#define DATUM_FILE QString("datums.csv")
#define MAP_FILE QString("maps.txt")
#define MAP_DIR QString("maps")
#define POI_DIR QString("POI")
@ -29,12 +31,16 @@
#define GLOBAL_DIR QString("/usr/share/gpxsee")
#endif
#define USER_MAP_DIR USER_DIR + QString("/") + MAP_DIR
#define USER_MAP_FILE USER_DIR + QString("/") + MAP_FILE
#define USER_POI_DIR USER_DIR + QString("/") + POI_DIR
#define GLOBAL_MAP_DIR GLOBAL_DIR + QString("/") + MAP_DIR
#define GLOBAL_MAP_FILE GLOBAL_DIR + QString("/") + MAP_FILE
#define GLOBAL_POI_DIR GLOBAL_DIR + QString("/") + POI_DIR
#define TILES_DIR USER_DIR + QString("/tiles")
#define USER_ELLIPSOID_FILE USER_DIR + QString("/") + ELLIPSOID_FILE
#define USER_DATUM_FILE USER_DIR + QString("/") + DATUM_FILE
#define USER_MAP_DIR USER_DIR + QString("/") + MAP_DIR
#define USER_MAP_FILE USER_DIR + QString("/") + MAP_FILE
#define USER_POI_DIR USER_DIR + QString("/") + POI_DIR
#define GLOBAL_ELLIPSOID_FILE GLOBAL_DIR + QString("/") + ELLIPSOID_FILE
#define GLOBAL_DATUM_FILE GLOBAL_DIR + QString("/") + DATUM_FILE
#define GLOBAL_MAP_DIR GLOBAL_DIR + QString("/") + MAP_DIR
#define GLOBAL_MAP_FILE GLOBAL_DIR + QString("/") + MAP_FILE
#define GLOBAL_POI_DIR GLOBAL_DIR + QString("/") + POI_DIR
#define TILES_DIR USER_DIR + QString("/tiles")
#endif /* CONFIG_H */

83
src/datum.cpp Normal file
View File

@ -0,0 +1,83 @@
#include <QFile>
#include "wgs84.h"
#include "datum.h"
static QMap<QString, Datum> WGS84()
{
QMap<QString, Datum> map;
map.insert("WGS 84", Datum(Ellipsoid(WGS84_RADIUS, WGS84_FLATTENING),
0, 0, 0));
return map;
}
QMap<QString, Datum> Datum::_datums = WGS84();
QString Datum::_errorString;
int Datum::_errorLine = 0;
Datum Datum::datum(const QString &name)
{
QMap<QString, Datum>::const_iterator it = _datums.find(name);
if (it == _datums.end())
return Datum();
return it.value();
}
bool Datum::loadList(const QString &path)
{
QFile file(path);
bool res;
if (!file.open(QFile::ReadOnly)) {
_errorString = qPrintable(file.errorString());
return false;
}
_errorLine = 1;
_errorString.clear();
while (!file.atEnd()) {
QByteArray line = file.readLine();
QList<QByteArray> list = line.split(',');
if (list.size() != 6) {
_errorString = "Format error";
return false;
}
int eid = list[2].trimmed().toInt(&res);
if (!res) {
_errorString = "Invalid ellipsoid id";
return false;
}
double dx = list[3].trimmed().toDouble(&res);
if (!res) {
_errorString = "Invalid dx";
return false;
}
double dy = list[4].trimmed().toDouble(&res);
if (!res) {
_errorString = "Invalid dy";
return false;
}
double dz = list[5].trimmed().toDouble(&res);
if (!res) {
_errorString = "Invalid dz";
return false;
}
Ellipsoid e = Ellipsoid::ellipsoid(eid);
if (e.isNull()) {
_errorString = "Unknown ellipsoid ID";
return false;
}
Datum d(e, dx, dy, dz);
_datums.insert(list[0].trimmed(), d);
_errorLine++;
}
return true;
}

38
src/datum.h Normal file
View File

@ -0,0 +1,38 @@
#ifndef DATUM_H
#define DATUM_H
#include <QMap>
#include "ellipsoid.h"
class Datum
{
public:
Datum() : _ellipsoid(Ellipsoid()), _dx(0), _dy(0), _dz(0) {}
Datum(const Ellipsoid &ellipsoid, double dx, double dy, double dz)
: _ellipsoid(ellipsoid), _dx(dx), _dy(dy), _dz(dz) {}
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();}
bool isWGS84() const
{return _ellipsoid.isWGS84() && _dx == 0 && _dy == 0 && _dz == 0;}
static bool loadList(const QString &path);
static const QString &errorString() {return _errorString;}
static int errorLine() {return _errorLine;}
static Datum datum(const QString &name);
private:
Ellipsoid _ellipsoid;
double _dx, _dy, _dz;
static QMap<QString, Datum> _datums;
static QString _errorString;
static int _errorLine;
};
#endif // DATUM_H

View File

@ -1,42 +1,62 @@
#include <cmath>
#include <QString>
#include "wgs84.h"
#include <QFile>
#include "ellipsoid.h"
QMap<int, Ellipsoid> Ellipsoid::_ellipsoids;
QString Ellipsoid::_errorString;
int Ellipsoid::_errorLine = 0;
#define INTERNATIONAL_RADIUS 6378388.0
#define INTERNATIONAL_FLATTENING (1.0/297.0)
#define KRASSOVSKY_RADIUS 6378245.0
#define KRASSOVSKY_FLATTENING (1.0/298.3)
#define CLARKE1866_RADIUS 6378206.4
#define CLARKE1866_FLATTENING (1.0/294.9786982)
#define GRS80_RADIUS 6378137.0
#define GRS80_FLATTENING (1.0/298.257222101)
#define WGS72_RADIUS 6378135.0
#define WGS72_FLATTENING (1.0/298.26)
#define GRS67_RADIUS 6378160.0
#define GRS67_FLATTENING (1.0/298.25)
// Must be in Ellipsoid::Name order!
struct {double radius; double flattening;} static ellipsoids[] = {
{CLARKE1866_RADIUS, CLARKE1866_FLATTENING},
{GRS80_RADIUS, GRS80_FLATTENING},
{INTERNATIONAL_RADIUS, INTERNATIONAL_FLATTENING},
{KRASSOVSKY_RADIUS, KRASSOVSKY_FLATTENING},
{WGS84_RADIUS, WGS84_FLATTENING},
{WGS72_RADIUS, WGS72_FLATTENING},
{GRS67_RADIUS, GRS67_FLATTENING}
};
Ellipsoid::Ellipsoid()
Ellipsoid Ellipsoid::ellipsoid(int id)
{
_radius = WGS84_RADIUS;
_flattening = WGS84_FLATTENING;
QMap<int, Ellipsoid>::const_iterator it = _ellipsoids.find(id);
if (it == _ellipsoids.end())
return Ellipsoid();
return it.value();
}
Ellipsoid::Ellipsoid(Name name)
bool Ellipsoid::loadList(const QString &path)
{
_radius = ellipsoids[name].radius;
_flattening = ellipsoids[name].flattening;
QFile file(path);
bool res;
if (!file.open(QFile::ReadOnly)) {
_errorString = qPrintable(file.errorString());
return false;
}
_errorLine = 1;
_errorString.clear();
while (!file.atEnd()) {
QByteArray line = file.readLine();
QList<QByteArray> list = line.split(',');
if (list.size() != 4) {
_errorString = "Format error";
return false;
}
int id = list[0].trimmed().toInt(&res);
if (!res) {
_errorString = "Invalid ellipsoid id";
return false;
}
double radius = list[2].trimmed().toDouble(&res);
if (!res) {
_errorString = "Invalid ellipsoid radius";
return false;
}
double flattening = list[3].trimmed().toDouble(&res);
if (!res) {
_errorString = "Invalid ellipsoid flattening";
return false;
}
Ellipsoid e(radius, 1.0/flattening);
_ellipsoids.insert(id, e);
_errorLine++;
}
return true;
}

View File

@ -1,28 +1,37 @@
#ifndef ELLIPSOID_H
#define ELLIPSOID_H
#include <QString>
#include <QMap>
#include "wgs84.h"
class Ellipsoid
{
public:
enum Name {
Clarke1866,
GRS80,
International1924,
Krassowsky1940,
WGS84,
WGS72,
GRS67
};
Ellipsoid();
Ellipsoid(Name name);
Ellipsoid() : _radius(-1.0), _flattening(-1.0) {}
Ellipsoid(double radius, double flattening)
: _radius(radius), _flattening(flattening) {}
double radius() const {return _radius;}
double flattening() const {return _flattening;}
bool isNull() const {return _radius < 0 || _flattening < 0;}
bool isWGS84() const
{return _radius == WGS84_RADIUS && _flattening == WGS84_FLATTENING;}
static bool loadList(const QString &path);
static const QString &errorString() {return _errorString;}
static int errorLine() {return _errorLine;}
static Ellipsoid ellipsoid(int id);
private:
double _radius;
double _flattening;
static QMap<int, Ellipsoid> _ellipsoids;
static QString _errorString;
static int _errorLine;
};
#endif // ELLIPSOID_H

View File

@ -26,6 +26,8 @@
#include "keys.h"
#include "settings.h"
#include "data.h"
#include "ellipsoid.h"
#include "datum.h"
#include "map.h"
#include "maplist.h"
#include "mapdir.h"
@ -47,6 +49,7 @@
GUI::GUI()
{
loadDatums();
loadMaps();
loadPOIs();
@ -117,6 +120,45 @@ void GUI::createBrowser()
_browser->setFilter(filter);
}
void GUI::loadDatums()
{
QString ef, df;
if (QFile::exists(USER_ELLIPSOID_FILE))
ef = USER_ELLIPSOID_FILE;
else if (QFile::exists(GLOBAL_ELLIPSOID_FILE))
ef = GLOBAL_ELLIPSOID_FILE;
else
qWarning("No ellipsoids file found.");
if (QFile::exists(USER_DATUM_FILE))
df = USER_DATUM_FILE;
else if (QFile::exists(GLOBAL_DATUM_FILE))
df = GLOBAL_DATUM_FILE;
else
qWarning("No datums file found.");
if (!ef.isNull() && !df.isNull()) {
if (!Ellipsoid::loadList(ef)) {
if (Ellipsoid::errorLine())
qWarning("%s: parse error on line %d: %s", qPrintable(ef),
Ellipsoid::errorLine(), qPrintable(Ellipsoid::errorString()));
else
qWarning("%s: %s", qPrintable(ef), qPrintable(
Ellipsoid::errorString()));
} else {
if (!Datum::loadList(df)) {
if (Datum::errorLine())
qWarning("%s: parse error on line %d: %s", qPrintable(ef),
Datum::errorLine(), qPrintable(Datum::errorString()));
else
qWarning("%s: %s", qPrintable(ef), qPrintable(
Datum::errorString()));
}
}
}
}
void GUI::loadMaps()
{
QList<Map*> online, offline;

View File

@ -79,6 +79,7 @@ private slots:
private:
typedef QPair<QDate, QDate> DateRange;
void loadDatums();
void loadMaps();
void loadPOIs();
void closeFiles();

View File

@ -12,7 +12,7 @@
#include "wgs84.h"
#include "coordinates.h"
#include "matrix.h"
#include "ellipsoid.h"
#include "datum.h"
#include "latlon.h"
#include "mercator.h"
#include "transversemercator.h"
@ -22,17 +22,39 @@
#include "offlinemap.h"
struct {const char *name; Ellipsoid::Name ellipsoid;} static datums[] = {
{"S42", Ellipsoid::Krassowsky1940},
{"Pulkovo 1942", Ellipsoid::Krassowsky1940},
{"European 1950", Ellipsoid::International1924},
{"European 1979", Ellipsoid::International1924},
{"NZGD1949", Ellipsoid::International1924},
{"NAD27", Ellipsoid::Clarke1866},
{"NAD83", Ellipsoid::GRS80},
{"WGS 72", Ellipsoid::WGS72},
{"South American 1969", Ellipsoid::GRS67}
};
// Abridged Molodensky transformation
static Coordinates toWGS84(Coordinates c, const Datum &datum)
{
Ellipsoid WGS84(WGS84_RADIUS, WGS84_FLATTENING);
double dX = datum.dx();
double dY = datum.dy();
double dZ = datum.dz();
double slat = sin(deg2rad(c.lat()));
double clat = cos(deg2rad(c.lat()));
double slon = sin(deg2rad(c.lon()));
double clon = cos(deg2rad(c.lon()));
double ssqlat = slat * slat;
double from_f = datum.ellipsoid().flattening();
double df = WGS84.flattening() - from_f;
double from_a = datum.ellipsoid().radius();
double da = WGS84.radius() - from_a;
double from_esq = datum.ellipsoid().flattening()
* (2.0 - datum.ellipsoid().flattening());
double adb = 1.0 / (1.0 - from_f);
double rn = from_a / sqrt(1 - from_esq * ssqlat);
double rm = from_a * (1 - from_esq) / pow((1 - from_esq * ssqlat), 1.5);
double from_h = 0.0; // we're flat!
double dlat = (-dX * slat * clon - dY * slat * slon + dZ * clat + da * rn
* from_esq * slat * clat / from_a + +df * (rm * adb + rn / adb) * slat
* clat) / (rm + from_h);
double dlon = (-dX * slon + dY * clon) / ((rn + from_h) * clat);
return Coordinates(c.lon() + rad2deg(dlon), c.lat() + rad2deg(dlat));
}
int OfflineMap::parseMapFile(QIODevice &device, QList<ReferencePoint> &points,
QString &projection, ProjectionSetup &setup, QString &datum)
@ -161,32 +183,30 @@ bool OfflineMap::createProjection(const QString &datum,
return false;
}
Ellipsoid::Name ellipsoid = Ellipsoid::WGS84;
for (size_t i = 0; i < ARRAY_SIZE(datums); i++) {
if (datum.startsWith(datums[i].name)) {
ellipsoid = datums[i].ellipsoid;
break;
}
Datum d = Datum::datum(datum);
if (d.isNull()) {
qWarning("%s: %s: unknown datum", qPrintable(_name), qPrintable(datum));
return false;
}
if (projection == "Mercator")
_projection = new Mercator();
else if (projection == "Transverse Mercator")
_projection = new TransverseMercator(Ellipsoid(ellipsoid),
_projection = new TransverseMercator(d.ellipsoid(),
setup.longitudeOrigin, setup.scale, setup.falseEasting,
setup.falseNorthing);
else if (projection == "Latitude/Longitude")
_projection = new LatLon();
else if (projection == "Lambert Conformal Conic")
_projection = new LambertConic(Ellipsoid(ellipsoid),
_projection = new LambertConic(d.ellipsoid(),
setup.standardParallel1, setup.standardParallel2,
setup.latitudeOrigin, setup.longitudeOrigin, setup.scale,
setup.falseEasting, setup.falseNorthing);
else if (projection == "(UTM) Universal Transverse Mercator") {
if (setup.zone)
_projection = new UTM(Ellipsoid(ellipsoid), setup.zone);
_projection = new UTM(d.ellipsoid(), setup.zone);
else if (!points.first().ll.isNull())
_projection = new UTM(Ellipsoid(ellipsoid), points.first().ll);
_projection = new UTM(d.ellipsoid(), points.first().ll);
else {
qWarning("%s: Can not determine UTM zone", qPrintable(_name));
return false;
@ -197,9 +217,14 @@ bool OfflineMap::createProjection(const QString &datum,
return false;
}
for (int i = 0; i < points.size(); i++)
if (points.at(i).ll.isNull())
points[i].ll = _projection->xy2ll(points.at(i).pp);
for (int i = 0; i < points.size(); i++) {
if (points.at(i).ll.isNull()) {
if (d.isWGS84())
points[i].ll = _projection->xy2ll(points.at(i).pp);
else
points[i].ll = toWGS84(_projection->xy2ll(points.at(i).pp), d);
}
}
return true;
}