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

Add support for GeoTIFF images

Browse Source
This commit is contained in:
Martin Tůma
2018-01-08 23:47:45 +01:00
parent 84f41b5aa9
commit 5fa52b0166
44 changed files with 5295 additions and 762 deletions
Download Patch File Download Diff File Expand all files Collapse all files

159
src/map/datum.cpp
View File

@ -4,26 +4,103 @@
#include "datum.h"
static QMap<QString, Datum> WGS84()
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;
};
// Abridged Molodensky transformation
static Coordinates molodensky(const Coordinates &c, const Datum &from,
const Datum &to)
{
QMap<QString, Datum> map;
map.insert("WGS 84", Datum(Ellipsoid(WGS84_RADIUS, WGS84_FLATTENING),
0.0, 0.0, 0.0));
return map;
double rlat = deg2rad(c.lat());
double rlon = deg2rad(c.lon());
double slat = sin(rlat);
double clat = cos(rlat);
double slon = sin(rlon);
double clon = cos(rlon);
double ssqlat = slat * slat;
double dx = from.dx() - to.dx();
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 df = to_f - from_f;
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);
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;
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));
}
QMap<QString, Datum> Datum::_datums = WGS84();
QList<Datum::Entry> Datum::_datums = WGS84();
QString Datum::_errorString;
int Datum::_errorLine = 0;
Datum Datum::datum(const QString &name)
QList<Datum::Entry> Datum::WGS84()
{
QMap<QString, Datum>::const_iterator it = _datums.find(name);
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;
}
if (it == _datums.end())
return Datum();
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;
}
return it.value();
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();
}
bool Datum::loadList(const QString &path)
@ -47,6 +124,15 @@ bool Datum::loadList(const QString &path)
return false;
}
QString f1 = list[1].trimmed();
int id = 0;
if (!f1.isEmpty()) {
id = f1.toInt(&res);
if (!res) {
_errorString = "Invalid datum id";
return false;
}
}
int eid = list[2].trimmed().toInt(&res);
if (!res) {
_errorString = "Invalid ellipsoid id";
@ -68,14 +154,14 @@ bool Datum::loadList(const QString &path)
return false;
}
Ellipsoid e = Ellipsoid::ellipsoid(eid);
Ellipsoid e(eid);
if (e.isNull()) {
_errorString = "Unknown ellipsoid ID";
return false;
}
Datum d(e, dx, dy, dz);
_datums.insert(list[0].trimmed(), d);
_datums.append(Entry(list[0].trimmed(), id, d));
_errorLine++;
}
@ -83,38 +169,21 @@ bool Datum::loadList(const QString &path)
return true;
}
// Abridged Molodensky transformation
Coordinates Datum::toWGS84(const Coordinates &c) const
{
if (_ellipsoid.radius() == WGS84_RADIUS
&& _ellipsoid.flattening() == WGS84_FLATTENING
&& _dx == 0.0 && _dy == 0.0 && _dz == 0.0)
return c;
double rlat = deg2rad(c.lat());
double rlon = deg2rad(c.lon());
double slat = sin(rlat);
double clat = cos(rlat);
double slon = sin(rlon);
double clon = cos(rlon);
double ssqlat = slat * slat;
double from_f = ellipsoid().flattening();
double df = WGS84_FLATTENING - from_f;
double from_a = ellipsoid().radius();
double da = WGS84_RADIUS - from_a;
double from_esq = ellipsoid().flattening()
* (2.0 - 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;
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));
return _WGS84 ? c : molodensky(c, *this, Datum(Ellipsoid(WGS84_RADIUS,
WGS84_FLATTENING), 0.0, 0.0, 0.0));
}
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);
}
QDebug operator<<(QDebug dbg, const Datum &datum)
{
dbg.nospace() << "Datum(" << datum.ellipsoid() << ", " << datum.dx() << ", "
<< datum.dy() << ", " << datum.dz() << ")";
return dbg.space();
}

27
src/map/datum.h
View File

@ -1,16 +1,20 @@
#ifndef DATUM_H
#define DATUM_H
#include <QMap>
#include <cmath>
#include <QList>
#include <QDebug>
#include "ellipsoid.h"
#include "common/coordinates.h"
class Datum
{
public:
Datum() : _ellipsoid(Ellipsoid()), _dx(0.0), _dy(0.0), _dz(0.0) {}
Datum(const Ellipsoid &ellipsoid, double dx, double dy, double dz)
: _ellipsoid(ellipsoid), _dx(dx), _dy(dy), _dz(dz) {}
Datum() : _ellipsoid(Ellipsoid()), _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);
const Ellipsoid &ellipsoid() const {return _ellipsoid;}
double dx() const {return _dx;}
@ -18,23 +22,30 @@ public:
double dz() const {return _dz;}
bool isNull() const
{return (_ellipsoid.isNull() && _dx == 0.0 && _dy == 0.0 && _dz == 0.0);}
{return (_ellipsoid.isNull() && 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;}
static Datum datum(const QString &name);
private:
class Entry;
static QList<Entry> WGS84();
Ellipsoid _ellipsoid;
double _dx, _dy, _dz;
bool _WGS84;
static QMap<QString, Datum> _datums;
static QList<Entry> _datums;
static QString _errorString;
static int _errorLine;
};
QDebug operator<<(QDebug dbg, const Datum &datum);
#endif // DATUM_H

19
src/map/ellipsoid.cpp
View File

@ -1,18 +1,19 @@
#include <QFile>
#include <QDebug>
#include "ellipsoid.h"
QMap<int, Ellipsoid> Ellipsoid::_ellipsoids;
QString Ellipsoid::_errorString;
int Ellipsoid::_errorLine = 0;
Ellipsoid Ellipsoid::ellipsoid(int id)
Ellipsoid::Ellipsoid(int id)
{
QMap<int, Ellipsoid>::const_iterator it = _ellipsoids.find(id);
if (it == _ellipsoids.end())
return Ellipsoid();
return it.value();
*this = Ellipsoid();
else
*this = it.value();
}
bool Ellipsoid::loadList(const QString &path)
@ -20,6 +21,7 @@ bool Ellipsoid::loadList(const QString &path)
QFile file(path);
bool res;
if (!file.open(QFile::ReadOnly)) {
_errorString = qPrintable(file.errorString());
return false;
@ -36,7 +38,7 @@ bool Ellipsoid::loadList(const QString &path)
return false;
}
int id = list[0].trimmed().toInt(&res);
int id = list[1].trimmed().toInt(&res);
if (!res) {
_errorString = "Invalid ellipsoid id";
return false;
@ -60,3 +62,10 @@ bool Ellipsoid::loadList(const QString &path)
return true;
}
QDebug operator<<(QDebug dbg, const Ellipsoid &ellipsoid)
{
dbg.nospace() << "Ellipsoid(" << ellipsoid.radius() << ", "
<< 1.0 / ellipsoid.flattening() << ")";
return dbg.space();
}

12
src/map/ellipsoid.h
View File

@ -1,27 +1,29 @@
#ifndef ELLIPSOID_H
#define ELLIPSOID_H
#include <cmath>
#include <QString>
#include <QMap>
#include <QDebug>
class Ellipsoid
{
public:
Ellipsoid() : _radius(0.0), _flattening(0.0) {}
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 (_radius == 0.0 && _flattening == 0.0);}
bool isNull() 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 Ellipsoid ellipsoid(int id);
private:
double _radius;
double _flattening;
@ -31,4 +33,6 @@ private:
static int _errorLine;
};
QDebug operator<<(QDebug dbg, const Ellipsoid &ellipsoid);
#endif // ELLIPSOID_H

457
src/map/geotiff.cpp Normal file
View File

@ -0,0 +1,457 @@
#include <QFileInfo>
#include <QtEndian>
#include "pcs.h"
#include "latlon.h"
#include "geotiff.h"
#define TIFF_DOUBLE 12
#define TIFF_SHORT 3
#define TIFF_LONG 4
#define ModelPixelScaleTag 33550
#define ModelTiepointTag 33922
#define ModelTransformationTag 34264
#define GeoKeyDirectoryTag 34735
#define GeoDoubleParamsTag 34736
#define ImageWidth 256
#define ImageHeight 257
#define GTModelTypeGeoKey 1024
#define GTRasterTypeGeoKey 1025
#define GeographicTypeGeoKey 2048
#define GeogGeodeticDatumGeoKey 2050
#define GeogEllipsoidGeoKey 2056
#define ProjectedCSTypeGeoKey 3072
#define ProjectionGeoKey 3074
#define ProjCoordTransGeoKey 3075
#define ProjStdParallel1GeoKey 3078
#define ProjStdParallel2GeoKey 3079
#define ProjNatOriginLongGeoKey 3080
#define ProjNatOriginLatGeoKey 3081
#define ProjFalseEastingGeoKey 3082
#define ProjFalseNorthingGeoKey 3083
#define ProjScaleAtNatOriginGeoKey 3092
#define ModelTypeProjected 1
#define ModelTypeGeographic 2
#define ModelTypeGeocentric 3
#define IS_SET(map, key) \
((map).contains(key) && (map).value(key).SHORT != 32767)
#define ARRAY_SIZE(a) \
(sizeof(a) / sizeof(*a))
typedef struct {
quint16 KeyDirectoryVersion;
quint16 KeyRevision;
quint16 MinorRevision;
quint16 NumberOfKeys;
} Header;
typedef struct {
quint16 KeyID;
quint16 TIFFTagLocation;
quint16 Count;
quint16 ValueOffset;
} KeyEntry;
bool GeoTIFF::readEntry(TIFFFile &file, Ctx &ctx)
{
quint16 tag;
quint16 type;
quint32 count, offset;
if (!file.readValue(tag))
return false;
if (!file.readValue(type))
return false;
if (!file.readValue(count))
return false;
if (!file.readValue(offset))
return false;
switch (tag) {
case ImageWidth:
if (!(type == TIFF_SHORT || type == TIFF_LONG))
return false;
_size.setWidth(offset);
break;
case ImageHeight:
if (!(type == TIFF_SHORT || type == TIFF_LONG))
return false;
_size.setHeight(offset);
break;
case ModelPixelScaleTag:
if (type != TIFF_DOUBLE || count != 3)
return false;
ctx.scale = offset;
break;
case ModelTiepointTag:
if (type != TIFF_DOUBLE || count < 6 || count % 6)
return false;
ctx.tiepoints = offset;
ctx.tiepointCount = count / 6;
break;
case GeoKeyDirectoryTag:
if (type != TIFF_SHORT)
return false;
ctx.keys = offset;
break;
case GeoDoubleParamsTag:
if (type != TIFF_DOUBLE)
return false;
ctx.values = offset;
break;
case ModelTransformationTag:
if (type != TIFF_DOUBLE || count != 16)
return false;
ctx.matrix = offset;
break;
}
return true;
}
bool GeoTIFF::readIFD(TIFFFile &file, quint32 &offset, Ctx &ctx)
{
quint16 count;
if (!file.seek(offset))
return false;
if (!file.readValue(count))
return false;
for (quint16 i = 0; i < count; i++)
if (!readEntry(file, ctx))
return false;
if (!file.readValue(offset))
return false;
return true;
}
bool GeoTIFF::readScale(TIFFFile &file, quint32 offset, QPointF &scale)
{
if (!file.seek(offset))
return false;
if (!file.readValue(scale.rx()))
return false;
if (!file.readValue(scale.ry()))
return false;
return true;
}
bool GeoTIFF::readTiepoints(TIFFFile &file, quint32 offset, quint32 count,
QList<ReferencePoint> &points)
{
double z, pz;
QPointF xy, pp;
if (!file.seek(offset))
return false;
for (quint32 i = 0; i < count; i++) {
if (!file.readValue(xy.rx()))
return false;
if (!file.readValue(xy.ry()))
return false;
if (!file.readValue(z))
return false;
if (!file.readValue(pp.rx()))
return false;
if (!file.readValue(pp.ry()))
return false;
if (!file.readValue(pz))
return false;
ReferencePoint p;
p.xy = xy.toPoint();
p.pp = pp;
points.append(p);
}
return true;
}
bool GeoTIFF::readKeys(TIFFFile &file, Ctx &ctx, QMap<quint16, Value> &kv)
{
Header header;
KeyEntry entry;
Value value;
if (!file.seek(ctx.keys))
return false;
if (!file.readValue(header.KeyDirectoryVersion))
return false;
if (!file.readValue(header.KeyRevision))
return false;
if (!file.readValue(header.MinorRevision))
return false;
if (!file.readValue(header.NumberOfKeys))
return false;
for (int i = 0; i < header.NumberOfKeys; i++) {
if (!file.readValue(entry.KeyID))
return false;
if (!file.readValue(entry.TIFFTagLocation))
return false;
if (!file.readValue(entry.Count))
return false;
if (!file.readValue(entry.ValueOffset))
return false;
switch (entry.KeyID) {
case GeographicTypeGeoKey:
case ProjectedCSTypeGeoKey:
case ProjCoordTransGeoKey:
case GTModelTypeGeoKey:
case GTRasterTypeGeoKey:
case GeogGeodeticDatumGeoKey:
case ProjectionGeoKey:
case GeogEllipsoidGeoKey:
if (entry.TIFFTagLocation != 0 || entry.Count != 1)
return false;
value.SHORT = entry.ValueOffset;
kv.insert(entry.KeyID, value);
break;
case ProjScaleAtNatOriginGeoKey:
case ProjNatOriginLongGeoKey:
case ProjNatOriginLatGeoKey:
case ProjFalseEastingGeoKey:
case ProjFalseNorthingGeoKey:
case ProjStdParallel1GeoKey:
case ProjStdParallel2GeoKey:
if (!readGeoValue(file, ctx.values, entry.ValueOffset,
value.DOUBLE))
return false;
kv.insert(entry.KeyID, value);
break;
default:
break;
}
}
return true;
}
bool GeoTIFF::readGeoValue(TIFFFile &file, quint32 offset, quint16 index,
double &val)
{
qint64 pos = file.pos();
if (!file.seek(offset + index * sizeof(double)))
return false;
if (!file.readValue(val))
return false;
if (!file.seek(pos))
return false;
return true;
}
Datum GeoTIFF::datum(QMap<quint16, Value> &kv)
{
Datum datum;
if (IS_SET(kv, GeographicTypeGeoKey)) {
datum = Datum(kv.value(GeographicTypeGeoKey).SHORT);
if (datum.isNull())
_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";
return datum;
}
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 false;
}
quint16 index = kv.value(ProjCoordTransGeoKey).SHORT;
epsg = (index >= ARRAY_SIZE(table)) ? 0 : table[index];
if (!epsg) {
_errorString = QString("Unknown coordinate transformation method");
return false;
}
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())
_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")
.arg(kv.value(GeographicTypeGeoKey).SHORT)
.arg(kv.value(ProjectionGeoKey).SHORT);
} else {
Datum d = datum(kv);
if (d.isNull())
return false;
Projection::Method m = method(kv);
if (m.isNull())
return false;
Projection::Setup setup(
kv.value(ProjNatOriginLatGeoKey).DOUBLE,
kv.value(ProjNatOriginLongGeoKey).DOUBLE,
kv.value(ProjScaleAtNatOriginGeoKey).DOUBLE,
kv.value(ProjFalseEastingGeoKey).DOUBLE,
kv.value(ProjFalseNorthingGeoKey).DOUBLE,
kv.value(ProjStdParallel1GeoKey).DOUBLE,
kv.value(ProjStdParallel2GeoKey).DOUBLE
);
pcs = PCS(d, 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())
return false;
_projection = new LatLon();
return true;
}
bool GeoTIFF::load(const QString &path)
{
quint32 ifd;
QList<ReferencePoint> points;
QPointF scale;
QMap<quint16, Value> kv;
Ctx ctx;
TIFFFile file;
file.setFileName(path);
if (!file.open(QIODevice::ReadOnly)) {
_errorString = QString("Error opening TIFF file: %1")
.arg(file.errorString());
return false;
}
if (!file.readHeader(ifd)) {
_errorString = "Invalid TIFF header";
return false;
}
while (ifd) {
if (!readIFD(file, ifd, ctx)) {
_errorString = "Invalid IFD";
return false;
}
}
if (!_size.isValid()) {
_errorString = "Invalid/missing image size";
return false;
}
if (!ctx.keys) {
_errorString = "Not a GeoTIFF file";
return false;
}
if (ctx.scale) {
if (!readScale(file, ctx.scale, scale)) {
_errorString = "Error reading model pixel scale";
return false;
}
}
if (ctx.tiepoints) {
if (!readTiepoints(file, ctx.tiepoints, ctx.tiepointCount, points)) {
_errorString = "Error reading raster->model tiepoint pairs";
return false;
}
}
if (!readKeys(file, ctx, kv)) {
_errorString = "Error reading Geo key/value";
return false;
}
switch (kv.value(GTModelTypeGeoKey).SHORT) {
case ModelTypeProjected:
if (!projectedModel(kv))
return false;
break;
case ModelTypeGeographic:
if (!geographicModel(kv))
return false;
break;
case ModelTypeGeocentric:
_errorString = "TODO ModelTypeGeocentric";
return false;
default:
_errorString = "Unknown/missing model type";
return false;
}
if (ctx.scale && ctx.tiepoints)
_transform = QTransform(scale.x(), 0, 0, -scale.y(),
points.first().pp.x(), points.first().pp.y()).inverted();
else if (ctx.tiepointCount > 1) {
Transform t(points);
if (t.isNull()) {
_errorString = t.errorString();
return false;
}
_transform = t.transform();
} else if (ctx.matrix) {
_errorString = "TODO ModelTransformationTag";
return false;
} else {
_errorString = "Incomplete/missing model transformation info";
return false;
}
QFileInfo fi(path);
_name = fi.fileName();
return true;
}

69
src/map/geotiff.h Normal file
View File

@ -0,0 +1,69 @@
#ifndef GEOTIFF_H
#define GEOTIFF_H
#include <QTransform>
#include <QFile>
#include <QMap>
#include "datum.h"
#include "projection.h"
#include "tifffile.h"
#include "transform.h"
class GeoTIFF
{
public:
GeoTIFF() : _projection(0) {}
bool load(const QString &path);
const Datum &datum() const {return _datum;}
Projection *projection() const {return _projection;}
const QTransform &transform() const {return _transform;}
const QString name() const {return _name;}
const QSize &size() const {return _size;}
const QString &errorString() const {return _errorString;}
private:
union Value {
quint16 SHORT;
double DOUBLE;
};
struct Ctx {
quint32 scale;
quint32 tiepoints;
quint32 tiepointCount;
quint32 matrix;
quint32 keys;
quint32 values;
Ctx() : scale(0), tiepoints(0), tiepointCount(0), matrix(0), keys(0),
values(0) {}
};
bool readEntry(TIFFFile &file, Ctx &ctx);
bool readIFD(TIFFFile &file, quint32 &offset, Ctx &ctx);
bool readScale(TIFFFile &file, quint32 offset, QPointF &scale);
bool readTiepoints(TIFFFile &file, quint32 offset, quint32 count,
QList<ReferencePoint> &points);
bool readKeys(TIFFFile &file, Ctx &ctx, QMap<quint16, Value> &kv);
bool readGeoValue(TIFFFile &file, quint32 offset, quint16 index,
double &val);
Datum datum(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;
Projection *_projection;
QSize _size;
QString _name;
QString _errorString;
};
#endif // GEOTIFF_H

8
src/map/lambertazimuthal.cpp
View File

@ -54,13 +54,13 @@ LambertAzimuthal::LambertAzimuthal(const Ellipsoid &ellipsoid,
{
double es2, es4, es6;
_e = ellipsoid;
es2 = 2 * _e.flattening() - _e.flattening() * _e.flattening();
es2 = 2 * ellipsoid.flattening() - ellipsoid.flattening()
* ellipsoid.flattening();
es4 = es2 * es2;
es6 = es4 * es2;
_ra = _e.radius() * (1.0 - es2 / 6.0 - 17.0 * es4 / 360.0 - 67.0 * es6
/ 3024.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);
_cosLatOrigin = cos(_latOrigin);

2
src/map/lambertazimuthal.h
View File

@ -14,8 +14,6 @@ public:
virtual Coordinates xy2ll(const QPointF &p) const;
private:
Ellipsoid _e;
double _ra;
double _sinLatOrigin;
double _cosLatOrigin;

294
src/map/lambertconic.cpp
View File

@ -1,3 +1,46 @@
/*
* Based on libgeotrans with the following Source Code Disclaimer:
1. The GEOTRANS source code ("the software") is provided free of charge by
the National Imagery and Mapping Agency (NIMA) of the United States
Department of Defense. Although NIMA makes no copyright claim under Title 17
U.S.C., NIMA claims copyrights in the source code under other legal regimes.
NIMA hereby grants to each user of the software a license to use and
distribute the software, and develop derivative works.
2. Warranty Disclaimer: The software was developed to meet only the internal
requirements of the U.S. National Imagery and Mapping Agency. The software
is provided "as is," and no warranty, express or implied, including but not
limited to the implied warranties of merchantability and fitness for
particular purpose or arising by statute or otherwise in law or from a
course of dealing or usage in trade, is made by NIMA as to the accuracy and
functioning of the software.
3. NIMA and its personnel are not required to provide technical support or
general assistance with respect to the software.
4. Neither NIMA nor its personnel will be liable for any claims, losses, or
damages arising from or connected with the use of the software. The user
agrees to hold harmless the United States National Imagery and Mapping
Agency. The user's sole and exclusive remedy is to stop using the software.
5. NIMA requests that products developed using the software credit the
source of the software with the following statement, "The product was
developed using GEOTRANS, a product of the National Imagery and Mapping
Agency and U.S. Army Engineering Research and Development Center."
6. For any products developed using the software, NIMA requires a disclaimer
that use of the software does not indicate endorsement or approval of the
product by the Secretary of Defense or the National Imagery and Mapping
Agency. Pursuant to the United States Code, 10 U.S.C. Sec. 2797, the name of
the National Imagery and Mapping Agency, the initials "NIMA", the seal of
the National Imagery and Mapping Agency, or any colorable imitation thereof
shall not be used to imply approval, endorsement, or authorization of a
product without prior written permission from United States Secretary of
Defense.
*/
#include <cmath>
#include "lambertconic.h"
@ -5,76 +48,229 @@
#define M_PI_2 1.57079632679489661923
#endif // M_PI_2
#ifndef M_PI_4
#define M_PI_4 0.78539816339744830962
#endif // M_PI_4
#define M_PI_4 0.785398163397448309616
#endif /* M_PI_4 */
static double q(const Ellipsoid &el, double b)
#define LAMBERT_m(clat, essin) (clat / sqrt(1.0 - essin * essin))
#define LAMBERT2_t(lat, essin, es_over_2) \
(tan(M_PI_4 - lat / 2) * pow((1.0 + essin) / (1.0 - essin), es_over_2))
#define LAMBERT1_t(lat, essin, es_over_2) \
(tan(M_PI_4 - lat / 2) / pow((1.0 - essin) / (1.0 + essin), es_over_2))
LambertConic1::LambertConic1(const Ellipsoid &ellipsoid, double latitudeOrigin,
double longitudeOrigin, double scale, double falseEasting,
double falseNorthing)
{
double e = sqrt(el.flattening() * (2. - el.flattening()));
double esb = e * sin(b);
return log(tan(M_PI_4 + b / 2.) * pow((1. - esb) / (1. + esb), e / 2.));
double es2;
double es_sin;
double m0;
double lat_orig;
lat_orig = deg2rad(latitudeOrigin);
_longitudeOrigin = deg2rad(longitudeOrigin);
if (_longitudeOrigin > M_PI)
_longitudeOrigin -= 2 * M_PI;
_falseEasting = falseEasting;
_falseNorthing = falseNorthing;
es2 = 2.0 * ellipsoid.flattening() - ellipsoid.flattening()
* ellipsoid.flattening();
_es = sqrt(es2);
_es_over_2 = _es / 2.0;
_n = sin(lat_orig);
es_sin = _es * sin(lat_orig);
m0 = LAMBERT_m(cos(lat_orig), es_sin);
_t0 = LAMBERT1_t(lat_orig, es_sin, _es_over_2);
_rho0 = ellipsoid.radius() * scale * m0 / _n;
_rho_olat = _rho0;
}
static double iq(const Ellipsoid &el, double q)
QPointF LambertConic1::ll2xy(const Coordinates &c) const
{
double e = sqrt(el.flattening() * (2. - el.flattening()));
double b0 = 0.;
double b = 2. * atan(exp(q)) - M_PI_2;
double t;
double rho;
double dlam;
double theta;
double lat = deg2rad(c.lat());
do {
b0 = b;
double esb = e * sin(b);
b = 2. * atan(exp(q) * pow((1. - esb) / (1. + esb), -e / 2.)) - M_PI_2;
} while (fabs(b - b0) > 1e-10);
return b;
if (fabs(fabs(lat) - M_PI_2) > 1.0e-10) {
t = LAMBERT1_t(lat, _es * sin(lat), _es_over_2);
rho = _rho0 * pow(t / _t0, _n);
} else
rho = 0.0;
dlam = deg2rad(c.lon()) - _longitudeOrigin;
if (dlam > M_PI)
dlam -= 2 * M_PI;
if (dlam < -M_PI)
dlam += 2 * M_PI;
theta = _n * dlam;
return QPointF(rho * sin(theta) + _falseEasting, _rho_olat - rho
* cos(theta) + _falseNorthing);
}
static double nradius(const Ellipsoid &el, double phi)
Coordinates LambertConic1::xy2ll(const QPointF &p) const
{
double sin_phi = sin(phi);
return (el.radius() / sqrt(1. - (el.flattening() * (2. - el.flattening()))
* sin_phi * sin_phi));
double dx;
double dy;
double rho;
double rho_olat_minus_dy;
double t;
double PHI;
double es_sin;
double tempPHI = 0.0;
double theta = 0.0;
double tolerance = 4.85e-10;
int count = 30;
double lat, lon;
dy = p.y() - _falseNorthing;
dx = p.x() - _falseEasting;
rho_olat_minus_dy = _rho_olat - dy;
rho = sqrt(dx * dx + (rho_olat_minus_dy) * (rho_olat_minus_dy));
if (_n < 0.0) {
rho *= -1.0;
dx *= -1.0;
rho_olat_minus_dy *= -1.0;
}
if (rho != 0.0) {
theta = atan2(dx, rho_olat_minus_dy) / _n;
t = _t0 * pow(rho / _rho0, 1 / _n);
PHI = M_PI_2 - 2.0 * atan(t);
while (fabs(PHI - tempPHI) > tolerance && count) {
tempPHI = PHI;
es_sin = _es * sin(PHI);
PHI = M_PI_2 - 2.0 * atan(t * pow((1.0 - es_sin) / (1.0 + es_sin),
_es_over_2));
count--;
}
if (!count)
return Coordinates();
lat = PHI;
lon = theta + _longitudeOrigin;
if (fabs(lat) < 2.0e-7)
lat = 0.0;
if (lat > M_PI_2)
lat = M_PI_2;
else if (lat < -M_PI_2)
lat = -M_PI_2;
if (lon > M_PI) {
if (lon - M_PI < 3.5e-6)
lon = M_PI;
else
lon -= 2 * M_PI;
}
if (lon < -M_PI) {
if (fabs(lon + M_PI) < 3.5e-6)
lon = -M_PI;
else
lon += 2 * M_PI;
}
if (fabs(lon) < 2.0e-7)
lon = 0.0;
if (lon > M_PI)
lon = M_PI;
else if (lon < -M_PI)
lon = -M_PI;
} else {
if (_n > 0.0)
lat = M_PI_2;
else
lat = -M_PI_2;
lon = _longitudeOrigin;
}
return Coordinates(rad2deg(lon), rad2deg(lat));
}
LambertConic::LambertConic(const Ellipsoid &ellipsoid, double standardParallel1,
double standardParallel2, double latitudeOrigin, double longitudeOrigin,
double scale, double falseEasting, double falseNorthing) : _e(ellipsoid)
LambertConic2::LambertConic2(const Ellipsoid &ellipsoid,
double standardParallel1, double standardParallel2, double latitudeOrigin,
double longitudeOrigin, double falseEasting, double falseNorthing)
{
_cm = deg2rad(longitudeOrigin);
_fe = falseEasting;
_fn = falseNorthing;
double es, es_over_2, es2, es_sin;
double lat0;
double k0;
double t0;
double t1, t2;
double t_olat;
double m0;
double m1;
double m2;
double n;
double const_value;
double sp1, sp2;
double lat_orig;
double sp1 = deg2rad(standardParallel1);
double sp2 = deg2rad(standardParallel2);
double N1 = nradius(_e, sp1);
double N2 = nradius(_e, sp2);
lat_orig = deg2rad(latitudeOrigin);
sp1 = deg2rad(standardParallel1);
sp2 = deg2rad(standardParallel2);
_q0 = q(_e, deg2rad(latitudeOrigin));
double q1 = q(_e, sp1);
double q2 = q(_e, sp2);
if (fabs(sp1 - sp2) > 1.0e-10) {
es2 = 2 * ellipsoid.flattening() - ellipsoid.flattening()
* ellipsoid.flattening();
es = sqrt(es2);
es_over_2 = es / 2.0;
_n = log((N1 * cos(sp1)) / (N2 * cos(sp2))) / (q2 - q1);
double R1 = N1 * cos(sp1) / _n;
_R0 = scale * R1 * exp(_n * (q1 - _q0));
es_sin = es * sin(lat_orig);
t_olat = LAMBERT2_t(lat_orig, es_sin, es_over_2);
es_sin = es * sin(sp1);
m1 = LAMBERT_m(cos(sp1), es_sin);
t1 = LAMBERT2_t(sp1, es_sin, es_over_2);
es_sin = es * sin(sp2);
m2 = LAMBERT_m(cos(sp2), es_sin);
t2 = LAMBERT2_t(sp2, es_sin, es_over_2);
n = log(m1 / m2) / log(t1 / t2);
lat0 = asin(n);
es_sin = es * sin(lat0);
m0 = LAMBERT_m(cos(lat0), es_sin);
t0 = LAMBERT2_t(lat0, es_sin, es_over_2);
k0 = (m1 / m0) * (pow(t0 / t1, n));
const_value = ((ellipsoid.radius() * m2) / (n * pow(t2, n)));
falseNorthing += (const_value * pow(t_olat, n)) - (const_value
* pow(t0, n));
} else {
lat0 = sp1;
k0 = 1.0;
}
_lc1 = LambertConic1(ellipsoid, rad2deg(lat0), longitudeOrigin, k0,
falseEasting, falseNorthing);
}
QPointF LambertConic::ll2xy(const Coordinates &c) const
QPointF LambertConic2::ll2xy(const Coordinates &c) const
{
double dl = _n * (deg2rad(c.lon()) - _cm);
double R = _R0 * exp(_n * (_q0 - q(_e, deg2rad(c.lat()))));
return QPointF(_fe + R * sin(dl), _fn + _R0 - R * cos(dl));
return _lc1.ll2xy(c);
}
Coordinates LambertConic::xy2ll(const QPointF &p) const
Coordinates LambertConic2::xy2ll(const QPointF &p) const
{
double dl = atan((p.x() - _fe) / (_fn + _R0 - p.y()));
double dx = p.x() - _fe;
double dy = p.y() - _fn - _R0;
double R = sqrt(dx * dx + dy * dy);
double q = _q0 - log(fabs(R / _R0)) / _n;
return Coordinates(rad2deg(_cm + dl / _n), rad2deg(iq(_e, q)));
return _lc1.xy2ll(p);
}

38
src/map/lambertconic.h
View File

@ -4,26 +4,42 @@
#include "ellipsoid.h"
#include "projection.h"
class LambertConic : public Projection
class LambertConic1 : public Projection
{
public:
LambertConic(const Ellipsoid &ellipsoid, double standardParallel1,
double standardParallel2, double latitudeOrigin, double longitudeOrigin,
double scale, double falseEasting, double falseNorthing);
LambertConic1() {}
LambertConic1(const Ellipsoid &ellipsoid, double latitudeOrigin,
double longitudeOrigin, double scale, double falseEasting,
double falseNorthing);
virtual QPointF ll2xy(const Coordinates &c) const;
virtual Coordinates xy2ll(const QPointF &p) const;
private:
Ellipsoid _e;
double _longitudeOrigin;
double _falseEasting;
double _falseNorthing;
double _cm;
double _fe;
double _fn;
double _q0;
double _R0;
double _es;
double _es_over_2;
double _n;
double _t0;
double _rho0;
double _rho_olat;
};
class LambertConic2 : public Projection
{
public:
LambertConic2(const Ellipsoid &ellipsoid, double standardParallel1,
double standardParallel2, double latitudeOrigin, double longitudeOrigin,
double falseEasting, double falseNorthing);
virtual QPointF ll2xy(const Coordinates &c) const;
virtual Coordinates xy2ll(const QPointF &p) const;
private:
LambertConic1 _lc1;
};
#endif // LAMBERTCONIC_H

265
src/map/mapfile.cpp Normal file
View File

@ -0,0 +1,265 @@
#include "latlon.h"
#include "utm.h"
#include "mapfile.h"
static double parameter(const QString &str, bool *res)
{
QString field = str.trimmed();
if (field.isEmpty()) {
*res = true;
return NAN;
}
return field.toDouble(res);
}
int MapFile::parse(QIODevice &device, QList<CalibrationPoint> &points,
QString &projection, Projection::Setup &setup, QString &datum)
{
bool res, r[8];
int ln = 1;
while (!device.atEnd()) {
QByteArray line = device.readLine();
if (ln == 1) {
if (!line.trimmed().startsWith("OziExplorer Map Data File"))
return ln;
} else if (ln == 2)
_name = line.trimmed();
else if (ln == 3)
_image = line.trimmed();
else if (ln == 5)
datum = line.split(',').at(0).trimmed();
else {
QList<QByteArray> list = line.split(',');
QString key(list.at(0).trimmed());
bool ll = true; bool pp = true;
if (key.startsWith("Point") && list.count() == 17
&& !list.at(2).trimmed().isEmpty()) {
CalibrationPoint p;
int x = list.at(2).trimmed().toInt(&res);
if (!res)
return ln;
int y = list.at(3).trimmed().toInt(&res);
if (!res)
return ln;
int latd = list.at(6).trimmed().toInt(&res);
if (!res)
ll = false;
qreal latm = list.at(7).trimmed().toFloat(&res);
if (!res)
ll = false;
int lond = list.at(9).trimmed().toInt(&res);
if (!res)
ll = false;
qreal lonm = list.at(10).trimmed().toFloat(&res);
if (!res)
ll = false;
if (ll && list.at(8).trimmed() == "S") {
latd = -latd;
latm = -latm;
}
if (ll && list.at(11).trimmed() == "W") {
lond = -lond;
lonm = -lonm;
}
p.zone = list.at(13).trimmed().toInt(&res);
if (!res)
p.zone = 0;
qreal ppx = list.at(14).trimmed().toFloat(&res);
if (!res)
pp = false;
qreal ppy = list.at(15).trimmed().toFloat(&res);
if (!res)
pp = false;
if (list.at(16).trimmed() == "S")
p.zone = -p.zone;
p.rp.xy = QPoint(x, y);
if (ll) {
p.ll = Coordinates(lond + lonm/60.0, latd + latm/60.0);
if (p.ll.isValid())
points.append(p);
else
return ln;
} else if (pp) {
p.rp.pp = QPointF(ppx, ppy);
points.append(p);
} else
return ln;
} else if (key == "IWH") {
if (list.count() < 4)
return ln;
int w = list.at(2).trimmed().toInt(&res);
if (!res)
return ln;
int h = list.at(3).trimmed().toInt(&res);
if (!res)
return ln;
_size = QSize(w, h);
} else if (key == "Map Projection") {
if (list.count() < 2)
return ln;
projection = list.at(1);
} else if (key == "Projection Setup") {
if (list.count() < 8)
return ln;
setup = Projection::Setup(
parameter(list[1], &r[1]), parameter(list[2], &r[2]),
parameter(list[3], &r[3]), parameter(list[4], &r[4]),
parameter(list[5], &r[5]), parameter(list[6], &r[6]),
parameter(list[7], &r[7]));
for (int i = 1; i < 8; i++)
if (!r[i])
return ln;
}
}
ln++;
}
return 0;
}
bool MapFile::parseMapFile(QIODevice &device, QList<CalibrationPoint> &points,
QString &projection, Projection::Setup &setup, QString &datum)
{
int el;
if (!device.open(QIODevice::ReadOnly)) {
_errorString = QString("Error opening map file: %1")
.arg(device.errorString());
return false;
}
if ((el = parse(device, points, projection, setup, datum))) {
_errorString = QString("Map file parse error on line %1").arg(el);
return false;
}
return true;
}
bool MapFile::createDatum(const QString &datum)
{
_datum = Datum(datum);
if (_datum.isNull()) {
_errorString = QString("%1: Unknown datum").arg(datum);
return false;
}
return true;
}
bool MapFile::createProjection(const QString &name,
const Projection::Setup &setup, QList<CalibrationPoint> &points)
{
if (name == "Mercator")
_projection = Projection::projection(_datum, 1024, setup);
else if (name == "Transverse Mercator")
_projection = Projection::projection(_datum, 9807, setup);
else if (name == "Latitude/Longitude")
_projection = new LatLon();
else if (name == "Lambert Conformal Conic")
_projection = Projection::projection(_datum, 9802, setup);
else if (name == "Albers Equal Area")
_projection = Projection::projection(_datum, 9822, setup);
else if (name == "(A)Lambert Azimuthual Equal Area")
_projection = Projection::projection(_datum, 9820, setup);
else if (name == "(UTM) Universal Transverse Mercator") {
Projection::Setup s;
if (points.first().zone)
s = UTM::setup(points.first().zone);
else if (!points.first().ll.isNull())
s = UTM::setup(UTM::zone(points.first().ll));
else {
_errorString = "Can not determine UTM zone";
return 0;
}
_projection = Projection::projection(_datum, 9807, s);
} else if (name == "(NZTM2) New Zealand TM 2000")
_projection = Projection::projection(_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));
else if (name == "(IG) Irish Grid")
_projection = Projection::projection(_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));
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));
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));
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));
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));
else if (name == "(VICGRID) Victoria Australia")
_projection = Projection::projection(_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));
else {
_errorString = QString("%1: Unknown map projection").arg(name);
return false;
}
return true;
}
bool MapFile::computeTransformation(QList<CalibrationPoint> &points)
{
QList<ReferencePoint> rp;
for (int i = 0; i < points.size(); i++) {
if (points.at(i).rp.pp.isNull())
points[i].rp.pp = _projection->ll2xy(points.at(i).ll);
rp.append(points.at(i).rp);
}
Transform t(rp);
if (t.isNull()) {
_errorString = t.errorString();
return false;
}
_transform = t.transform();
return true;
}
bool MapFile::load(QIODevice &file)
{
QList<CalibrationPoint> points;
QString projection, datum;
Projection::Setup setup;
if (!parseMapFile(file, points, projection, setup, datum))
return false;
if (!createDatum(datum))
return false;
if (!createProjection(projection, setup, points))
return false;
if (!computeTransformation(points))
return false;
return true;
}

52
src/map/mapfile.h Normal file
View File

@ -0,0 +1,52 @@
#ifndef MAPFILE_H
#define MAPFILE_H
#include <QIODevice>
#include <QTransform>
#include "datum.h"
#include "transform.h"
#include "projection.h"
class MapFile
{
public:
bool load(QIODevice &file);
const Datum &datum() const {return _datum;}
Projection *projection() const {return _projection;}
const QTransform &transform() const {return _transform;}
const QString &name() const {return _name;}
const QString &image() const {return _image;}
const QSize &size() const {return _size;}
const QString &errorString() const {return _errorString;}
private:
struct CalibrationPoint {
ReferencePoint rp;
Coordinates ll;
int zone;
};
int parse(QIODevice &device, QList<CalibrationPoint> &points,
QString &projection, Projection::Setup &setup, QString &datum);
bool parseMapFile(QIODevice &device, QList<CalibrationPoint> &points,
QString &projection, Projection::Setup &setup, QString &datum);
bool createDatum(const QString &datum);
bool createProjection(const QString &projection,
const Projection::Setup &setup, QList<CalibrationPoint> &points);
bool computeTransformation(QList<CalibrationPoint> &points);
QString _name;
QString _image;
QSize _size;
Datum _datum;
QTransform _transform;
Projection *_projection;
QString _errorString;
};
#endif // MAPFILE_H

10
src/map/maplist.cpp
View File

@ -80,7 +80,7 @@ bool MapList::loadMap(const QString &path)
}
}
bool MapList::loadTba(const QString &path)
bool MapList::loadAtlas(const QString &path)
{
Atlas *atlas = new Atlas(path, this);
if (atlas->isValid()) {
@ -124,10 +124,10 @@ bool MapList::loadFile(const QString &path)
if (suffix == "txt")
return loadList(path);
else if (suffix == "map")
else if (suffix == "map" || suffix == "tif" || suffix == "tiff")
return loadMap(path);
else if (suffix == "tba")
return loadTba(path);
return loadAtlas(path);
else if (suffix == "tar")
return loadTar(path);
else {
@ -138,13 +138,13 @@ bool MapList::loadFile(const QString &path)
QString MapList::formats()
{
return tr("Map files (*.map *.tba *.tar)") + ";;"
return tr("Map files (*.map *.tba *.tar *.tif *.tiff)") + ";;"
+ tr("URL list files (*.txt)");
}
QStringList MapList::filter()
{
QStringList filter;
filter << "*.map" << "*.tba" << "*.tar" << "*.txt";
filter << "*.map" << "*.tba" << "*.tar" << "*.txt" << "*.tif" << "*.tiff";
return filter;
}

2
src/map/maplist.h
View File

@ -25,7 +25,7 @@ private:
bool loadList(const QString &path);
bool loadMap(const QString &path);
bool loadTba(const QString &path);
bool loadAtlas(const QString &path);
bool loadTar(const QString &path);
QList<Map*> _maps;

403
src/map/offlinemap.cpp
View File

@ -9,297 +9,20 @@
#include <QPixmapCache>
#include "common/coordinates.h"
#include "common/rectc.h"
#include "datum.h"
#include "utm.h"
#include "matrix.h"
#include "ozf.h"
#include "mapfile.h"
#include "geotiff.h"
#include "offlinemap.h"
int OfflineMap::parse(QIODevice &device, QList<ReferencePoint> &points,
QString &projection, Projection::Setup &setup, QString &datum)
void OfflineMap::computeResolution()
{
bool res;
int ln = 1;
Coordinates tl = xy2ll((bounds().topLeft()));
Coordinates br = xy2ll(bounds().bottomRight());
while (!device.atEnd()) {
QByteArray line = device.readLine();
qreal ds = tl.distanceTo(br);
qreal ps = QLineF(bounds().topLeft(), bounds().bottomRight()).length();
if (ln == 1) {
if (!line.trimmed().startsWith("OziExplorer Map Data File"))
return ln;
} else if (ln == 2)
_name = line.trimmed();
else if (ln == 3)
_imgPath = line.trimmed();
else if (ln == 5)
datum = line.split(',').at(0).trimmed();
else {
QList<QByteArray> list = line.split(',');
QString key(list.at(0).trimmed());
bool ll = true; bool pp = true;
if (key.startsWith("Point") && list.count() == 17
&& !list.at(2).trimmed().isEmpty()) {
int x = list.at(2).trimmed().toInt(&res);
if (!res)
return ln;
int y = list.at(3).trimmed().toInt(&res);
if (!res)
return ln;
int latd = list.at(6).trimmed().toInt(&res);
if (!res)
ll = false;
qreal latm = list.at(7).trimmed().toFloat(&res);
if (!res)
ll = false;
int lond = list.at(9).trimmed().toInt(&res);
if (!res)
ll = false;
qreal lonm = list.at(10).trimmed().toFloat(&res);
if (!res)
ll = false;
if (ll && list.at(8).trimmed() == "S") {
latd = -latd;
latm = -latm;
}
if (ll && list.at(11).trimmed() == "W") {
lond = -lond;
lonm = -lonm;
}
setup.zone = list.at(13).trimmed().toInt(&res);
if (!res)
setup.zone = 0;
qreal ppx = list.at(14).trimmed().toFloat(&res);
if (!res)
pp = false;
qreal ppy = list.at(15).trimmed().toFloat(&res);
if (!res)
pp = false;
if (list.at(16).trimmed() == "S")
setup.zone = -setup.zone;
ReferencePoint p;
p.xy = QPoint(x, y);
if (ll) {
p.ll = Coordinates(lond + lonm/60.0, latd + latm/60.0);
if (p.ll.isValid())
points.append(p);
else
return ln;
} else if (pp) {
p.pp = QPointF(ppx, ppy);
points.append(p);
} else
return ln;
} else if (key == "IWH") {
if (list.count() < 4)
return ln;
int w = list.at(2).trimmed().toInt(&res);
if (!res)
return ln;
int h = list.at(3).trimmed().toInt(&res);
if (!res)
return ln;
_size = QSize(w, h);
} else if (key == "Map Projection") {
if (list.count() < 2)
return ln;
projection = list.at(1);
} else if (key == "Projection Setup") {
if (list.count() < 8)
return ln;
setup.latitudeOrigin = list.at(1).trimmed().toFloat(&res);
if (!res)
setup.latitudeOrigin = 0;
setup.longitudeOrigin = list.at(2).trimmed().toFloat(&res);
if (!res)
setup.longitudeOrigin = 0;
setup.scale = list.at(3).trimmed().toFloat(&res);
if (!res)
setup.scale = 1.0;
setup.falseEasting = list.at(4).trimmed().toFloat(&res);
if (!res)
setup.falseEasting = 0;
setup.falseNorthing = list.at(5).trimmed().toFloat(&res);
if (!res)
setup.falseNorthing = 0;
setup.standardParallel1 = list.at(6).trimmed().toFloat(&res);
if (!res)
setup.standardParallel1 = 0;
setup.standardParallel2 = list.at(7).trimmed().toFloat(&res);
if (!res)
setup.standardParallel2 = 0;
}
}
ln++;
}
if (!setup.zone && !points.first().ll.isNull())
setup.zone = UTM::zone(points.first().ll);
return 0;
}
bool OfflineMap::parseMapFile(QIODevice &device, QList<ReferencePoint> &points,
QString &projection, Projection::Setup &setup, QString &datum)
{
int el;
if (!device.open(QIODevice::ReadOnly)) {
_errorString = QString("Error opening map file: %1")
.arg(device.errorString());
return false;
}
if ((el = parse(device, points, projection, setup, datum))) {
_errorString = QString("Map file parse error on line %1").arg(el);
return false;
}
return true;
}
bool OfflineMap::createProjection(const QString &datum,
const QString &projection, const Projection::Setup &setup,
QList<ReferencePoint> &points)
{
if (points.count() < 2) {
_errorString = "Insufficient number of reference points";
return false;
}
Datum d = Datum::datum(datum);
if (d.isNull()) {
_errorString = QString("%1: Unknown datum").arg(datum);
return false;
}
_projection = Projection::projection(projection, d.ellipsoid(), setup);
if (!_projection) {
_errorString = Projection::errorString();
return false;
}
for (int i = 0; i < points.size(); i++) {
if (points.at(i).ll.isNull())
points[i].ll = d.toWGS84(_projection->xy2ll(points.at(i).pp));
else
points[i].ll = d.toWGS84(points[i].ll);
}
return true;
}
bool OfflineMap::simpleTransformation(const QList<ReferencePoint> &points)
{
if (points.at(0).xy.x() == points.at(1).xy.x()
|| points.at(0).xy.y() == points.at(1).xy.y()) {
_errorString = "Invalid reference points tuple";
return false;
}
QPointF p0(_projection->ll2xy(points.at(0).ll));
QPointF p1(_projection->ll2xy(points.at(1).ll));
qreal dX, dY, lat0, lon0;
dX = (p0.x() - p1.x()) / (points.at(0).xy.x() - points.at(1).xy.x());
dY = (p1.y() - p0.y()) / (points.at(1).xy.y() - points.at(0).xy.y());
lat0 = p0.y() - points.at(0).xy.y() * dY;
lon0 = p1.x() - points.at(1).xy.x() * dX;
_transform = QTransform(1.0/dX, 0, 0, 1.0/dY, -lon0/dX, -lat0/dY);
_inverted = _transform.inverted();
return true;
}
bool OfflineMap::affineTransformation(const QList<ReferencePoint> &points)
{
Matrix c(3, 2);
c.zeroize();
for (size_t i = 0; i < c.h(); i++) {
for (size_t j = 0; j < c.w(); j++) {
for (int k = 0; k < points.size(); k++) {
double f[3], t[2];
QPointF p = _projection->ll2xy(points.at(k).ll);
f[0] = p.x();
f[1] = p.y();
f[2] = 1.0;
t[0] = points.at(k).xy.x();
t[1] = points.at(k).xy.y();
c.m(i,j) += f[i] * t[j];
}
}
}
Matrix Q(3, 3);
Q.zeroize();
for (int qi = 0; qi < points.size(); qi++) {
double v[3];
QPointF p = _projection->ll2xy(points.at(qi).ll);
v[0] = p.x();
v[1] = p.y();
v[2] = 1.0;
for (size_t i = 0; i < Q.h(); i++)
for (size_t j = 0; j < Q.w(); j++)
Q.m(i,j) += v[i] * v[j];
}
Matrix M = Q.augemented(c);
if (!M.eliminate()) {
_errorString = "Singular transformation matrix";
return false;
}
_transform = QTransform(M.m(0,3), M.m(0,4), M.m(1,3), M.m(1,4), M.m(2,3),
M.m(2,4));
_inverted = _transform.inverted();
return true;
}
bool OfflineMap::computeTransformation(const QList<ReferencePoint> &points)
{
Q_ASSERT(points.size() >= 2);
if (points.size() == 2)
return simpleTransformation(points);
else
return affineTransformation(points);
}
bool OfflineMap::computeResolution(QList<ReferencePoint> &points)
{
Q_ASSERT(points.count() >= 2);
int maxLon = 0, minLon = 0, maxLat = 0, minLat = 0;
qreal dLon, pLon, dLat, pLat;
for (int i = 1; i < points.size(); i++) {
if (points.at(i).ll.lon() < points.at(minLon).ll.lon())
minLon = i;
if (points.at(i).ll.lon() > points.at(maxLon).ll.lon())
maxLon = i;
if (points.at(i).ll.lat() < points.at(minLat).ll.lat())
minLat = i;
if (points.at(i).ll.lat() > points.at(maxLat).ll.lat())
maxLat = i;
}
dLon = points.at(minLon).ll.distanceTo(points.at(maxLon).ll);
pLon = QLineF(points.at(minLon).xy, points.at(maxLon).xy).length();
dLat = points.at(minLat).ll.distanceTo(points.at(maxLat).ll);
pLat = QLineF(points.at(minLat).xy, points.at(maxLat).xy).length();
_resolution = (dLon/pLon + dLat/pLat) / 2.0;
return true;
_resolution = ds/ps;
}
bool OfflineMap::getImageInfo(const QString &path)
@ -384,13 +107,9 @@ bool OfflineMap::totalSizeSet()
OfflineMap::OfflineMap(const QString &fileName, QObject *parent)
: Map(parent)
{
QList<ReferencePoint> points;
QString proj, datum;
Projection::Setup setup;
QFileInfo fi(fileName);
QString suffix = fi.suffix().toLower();
_valid = false;
_img = 0;
_projection = 0;
@ -411,22 +130,52 @@ OfflineMap::OfflineMap(const QString &fileName, QObject *parent)
return;
}
QBuffer mapFile(&ba);
if (!parseMapFile(mapFile, points, proj, setup, datum))
MapFile mf;
if (!mf.load(mapFile)) {
_errorString = mf.errorString();
return;
} else if (suffix =="map") {
} else {
_name = mf.name();
_size = mf.size();
_imgPath = mf.image();
_datum = mf.datum();
_projection = mf.projection();
_transform = mf.transform();
}
} else if (suffix == "map") {
MapFile mf;
QFile mapFile(fileName);
if (!parseMapFile(mapFile, points, proj, setup, datum))
if (!mf.load(mapFile)) {
_errorString = mf.errorString();
return;
} else {
_name = mf.name();
_size = mf.size();
_imgPath = mf.image();
_datum = mf.datum();
_projection = mf.projection();
_transform = mf.transform();
}
} else if (suffix == "tif" || suffix == "tiff") {
GeoTIFF gt;
if (!gt.load(fileName)) {
_errorString = gt.errorString();
return;
} else {
_name = gt.name();
_size = gt.size();
_imgPath = fileName;
_datum = gt.datum();
_projection = gt.projection();
_transform = gt.transform();
}
} else {
_errorString = "Not a map file";
return;
}
if (!createProjection(datum, proj, setup, points))
return;
if (!computeTransformation(points))
return;
computeResolution(points);
_inverted = _transform.inverted();
computeResolution();
if (_tar.isOpen()) {
if (!totalSizeSet())
@ -454,10 +203,8 @@ OfflineMap::OfflineMap(const QString &fileName, QObject *parent)
OfflineMap::OfflineMap(const QString &fileName, Tar &tar, QObject *parent)
: Map(parent)
{
QList<ReferencePoint> points;
QString proj, datum;
Projection::Setup setup;
QFileInfo fi(fileName);
MapFile mf;
_valid = false;
_img = 0;
@ -477,28 +224,28 @@ OfflineMap::OfflineMap(const QString &fileName, Tar &tar, QObject *parent)
}
QBuffer buffer(&ba);
if (!parseMapFile(buffer, points, proj, setup, datum))
if (!mf.load(buffer)) {
_errorString = mf.errorString();
return;
if (!createProjection(datum, proj, setup, points))
return;
if (!totalSizeSet())
return;
if (!computeTransformation(points))
return;
computeResolution(points);
}
_name = mf.name();
_size = mf.size();
_datum = mf.datum();
_projection = mf.projection();
_transform = mf.transform();
_inverted = _transform.inverted();
computeResolution();
_tarPath = fi.absolutePath() + "/" + fi.completeBaseName() + ".tar";
_imgPath = QString();
_valid = true;
}
OfflineMap::~OfflineMap()
{
if (_img)
delete _img;
if (_projection)
delete _projection;
delete _img;
delete _projection;
}
void OfflineMap::load()
@ -522,10 +269,8 @@ void OfflineMap::load()
void OfflineMap::unload()
{
if (_img) {
delete _img;
_img = 0;
}
delete _img;
_img = 0;
}
void OfflineMap::drawTiled(QPainter *painter, const QRectF &rect)
@ -599,13 +344,9 @@ void OfflineMap::drawOZF(QPainter *painter, const QRectF &rect)
void OfflineMap::drawImage(QPainter *painter, const QRectF &rect)
{
if (!_img || _img->isNull())
painter->fillRect(rect, _backgroundColor);
else {
QRect r(rect.toRect());
painter->drawImage(r.left(), r.top(), *_img, r.left(), r.top(),
r.width(), r.height());
}
QRect r(rect.toRect());
painter->drawImage(r.left(), r.top(), *_img, r.left(), r.top(),
r.width(), r.height());
}
void OfflineMap::draw(QPainter *painter, const QRectF &rect)
@ -614,26 +355,28 @@ void OfflineMap::draw(QPainter *painter, const QRectF &rect)
drawOZF(painter, rect);
else if (_tileSize.isValid())
drawTiled(painter, rect);
else
else if (_img && !_img->isNull())
drawImage(painter, rect);
else
painter->fillRect(rect, _backgroundColor);
}
QPointF OfflineMap::ll2xy(const Coordinates &c)
{
if (_ozf.isOpen()) {
QPointF p(_transform.map(_projection->ll2xy(c)));
QPointF p(_transform.map(_projection->ll2xy(_datum.fromWGS84(c))));
return QPointF(p.x() * _scale.x(), p.y() * _scale.y());
} else
return _transform.map(_projection->ll2xy(c));
return _transform.map(_projection->ll2xy(_datum.fromWGS84(c)));
}
Coordinates OfflineMap::xy2ll(const QPointF &p)
{
if (_ozf.isOpen()) {
return _projection->xy2ll(_inverted.map(QPointF(p.x() / _scale.x(),
p.y() / _scale.y())));
return _datum.toWGS84(_projection->xy2ll(_inverted.map(QPointF(p.x()
/ _scale.x(), p.y() / _scale.y()))));
} else
return _projection->xy2ll(_inverted.map(p));
return _datum.toWGS84(_projection->xy2ll(_inverted.map(p)));
}
QRectF OfflineMap::bounds() const

29
src/map/offlinemap.h
View File

@ -3,7 +3,9 @@
#include <QTransform>
#include "common/coordinates.h"
#include "datum.h"
#include "projection.h"
#include "transform.h"
#include "map.h"
#include "tar.h"
#include "ozf.h"
@ -50,40 +52,23 @@ public:
{return _transform.map(p);}
private:
struct ReferencePoint {
QPoint xy;
Coordinates ll;
QPointF pp;
};
int parse(QIODevice &device, QList<ReferencePoint> &points,
QString &projection, Projection::Setup &setup, QString &datum);
bool parseMapFile(QIODevice &device, QList<ReferencePoint> &points,
QString &projection, Projection::Setup &setup, QString &datum);
bool totalSizeSet();
bool createProjection(const QString &datum, const QString &projection,
const Projection::Setup &setup, QList<ReferencePoint> &points);
bool simpleTransformation(const QList<ReferencePoint> &points);
bool affineTransformation(const QList<ReferencePoint> &points);
bool computeTransformation(const QList<ReferencePoint> &points);
bool computeResolution(QList<ReferencePoint> &points);
bool getTileInfo(const QStringList &tiles, const QString &path = QString());
bool getImageInfo(const QString &path);
bool totalSizeSet();
void drawTiled(QPainter *painter, const QRectF &rect);
void drawOZF(QPainter *painter, const QRectF &rect);
void drawImage(QPainter *painter, const QRectF &rect);
void computeResolution();
void rescale(int zoom);
QString _name;
bool _valid;
QString _errorString;
QSize _size;
Datum _datum;
Projection *_projection;
QTransform _transform, _inverted;
qreal _resolution;
OZF _ozf;
Tar _tar;
@ -94,7 +79,11 @@ private:
QString _tileName;
int _zoom;
qreal _resolution;
QPointF _scale;
bool _valid;
QString _errorString;
};
#endif // OFFLINEMAP_H

125
src/map/pcs.cpp Normal file
View File

@ -0,0 +1,125 @@
#include "pcs.h"
class PCS::Entry {
public:
Entry(int id, int gcs, int proj, const PCS &pcs)
: _id(id), _gcs(gcs), _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;
PCS _pcs;
};
QList<PCS::Entry> PCS::_pcss;
QString PCS::_errorString;
int PCS::_errorLine;
static double parameter(const QString &str, bool *res)
{
QString field = str.trimmed();
if (field.isEmpty()) {
*res = true;
return NAN;
}
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();
}
PCS::PCS(int 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;
}
}
*this = PCS();
}
bool PCS::loadList(const QString &path)
{
QFile file(path);
bool res[12];
int id, gcs, proj, transform;
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() != 12) {
_errorString = "Format error";
return false;
}
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]) {
_errorString = "Parse error";
return false;
}
}
Datum datum(gcs);
if (datum.isNull()) {
_errorString = "Unknown datum";
return false;
}
Projection::Method method(transform);
if (method.isNull()) {
_errorString = "Unknown coordinates transformation method";
return false;
}
_pcss.append(Entry(id, gcs, proj, PCS(datum, method, setup)));
_errorLine++;
}
return true;
}
QDebug operator<<(QDebug dbg, const PCS &pcs)
{
dbg.nospace() << "PCS(" << pcs.datum() << ", " << pcs.method()
<< ", " << pcs.setup() << ")";
return dbg.space();
}

46
src/map/pcs.h Normal file
View File

@ -0,0 +1,46 @@
#ifndef PCS_H
#define PCS_H
#include <QDebug>
#include <QList>
#include "datum.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);
const Datum &datum() const {return _datum;}
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());}
static bool loadList(const QString &path);
static const QString &errorString() {return _errorString;}
static int errorLine() {return _errorLine;}
private:
class Entry;
Datum _datum;
Projection::Method _method;
Projection::Setup _setup;
static QList<Entry> _pcss;
static QString _errorString;
static int _errorLine;
};
QDebug operator<<(QDebug dbg, const PCS &pcs);
#endif // PCS_H

133
src/map/projection.cpp
View File

@ -1,89 +1,76 @@
#include "latlon.h"
#include "datum.h"
#include "mercator.h"
#include "transversemercator.h"
#include "utm.h"
#include "lambertconic.h"
#include "albersequal.h"
#include "lambertazimuthal.h"
#include "projection.h"
QString Projection::_errorString;
Projection *Projection::projection(const QString &name,
const Ellipsoid &ellipsoid, const Setup &setup)
Projection::Method::Method(int id)
{
if (setup.latitudeOrigin < -90.0 || setup.latitudeOrigin > 90.0
|| setup.longitudeOrigin < -180.0 || setup.longitudeOrigin > 180.0
|| setup.standardParallel1 < -90.0 || setup.standardParallel1 > 90.0
|| setup.standardParallel2 < -90.0 || setup.standardParallel2 > 90.0) {
_errorString = "Invalid projection setup";
return 0;
switch (id) {
case 1024:
case 9801:
case 9802:
case 9807:
case 9820:
case 9822:
case 9841:
_id = id;
break;
default:
_id = 0;
}
}
if (name == "Mercator")
return new Mercator();
else if (name == "Transverse Mercator")
return new TransverseMercator(ellipsoid,
setup.latitudeOrigin, setup.longitudeOrigin, setup.scale,
setup.falseEasting, setup.falseNorthing);
else if (name == "Latitude/Longitude")
return new LatLon();
else if (name == "Lambert Conformal Conic")
return new LambertConic(ellipsoid,
setup.standardParallel1, setup.standardParallel2,
setup.latitudeOrigin, setup.longitudeOrigin, setup.scale,
setup.falseEasting, setup.falseNorthing);
else if (name == "Albers Equal Area")
return new AlbersEqual(ellipsoid, setup.standardParallel1,
setup.standardParallel2, setup.latitudeOrigin, setup.longitudeOrigin,
setup.falseEasting, setup.falseNorthing);
else if (name == "(A)Lambert Azimuthual Equal Area")
return new LambertAzimuthal(ellipsoid, setup.latitudeOrigin,
setup.longitudeOrigin, setup.falseEasting, setup.falseNorthing);
else if (name == "(UTM) Universal Transverse Mercator") {
if (setup.zone)
return new UTM(ellipsoid, setup.zone);
else {
_errorString = "Can not determine UTM zone";
Projection *Projection::projection(const Datum &datum, const Method &method,
const Setup &setup)
{
const Ellipsoid &ellipsoid = datum.ellipsoid();
switch (method.id()) {
case 9807:
return new TransverseMercator(ellipsoid, setup.latitudeOrigin(),
setup.longitudeOrigin(), setup.scale(), setup.falseEasting(),
setup.falseNorthing());
case 1024:
case 9841:
return new Mercator();
case 9802:
return new LambertConic2(ellipsoid, setup.standardParallel1(),
setup.standardParallel2(), setup.latitudeOrigin(),
setup.longitudeOrigin(), setup.falseEasting(),
setup.falseNorthing());
case 9801:
return new LambertConic1(ellipsoid, setup.latitudeOrigin(),
setup.longitudeOrigin(), setup.scale(), setup.falseEasting(),
setup.falseNorthing());
case 9820:
return new LambertAzimuthal(ellipsoid, setup.latitudeOrigin(),
setup.longitudeOrigin(), setup.falseEasting(),
setup.falseNorthing());
case 9822:
return new AlbersEqual(ellipsoid, setup.standardParallel1(),
setup.standardParallel2(), setup.latitudeOrigin(),
setup.longitudeOrigin(), setup.falseEasting(),
setup.falseNorthing());
default:
return 0;
}
} else if (name == "(NZTM2) New Zealand TM 2000")
return new TransverseMercator(ellipsoid, 0, 173.0, 0.9996,
1600000, 10000000);
else if (name == "(BNG) British National Grid")
return new TransverseMercator(ellipsoid, 49, -2, 0.999601,
400000, -100000);
else if (name == "(IG) Irish Grid")
return new TransverseMercator(ellipsoid, 53.5, -8, 1.000035,
200000, 250000);
else if (name == "(SG) Swedish Grid")
return new TransverseMercator(ellipsoid, 0, 15.808278, 1,
1500000, 0);
else if (name == "(I) France Zone I")
return new LambertConic(ellipsoid, 48.598523, 50.395912,
49.5, 2.337229, 1 /*0.99987734*/, 600000, 1200000);
else if (name == "(II) France Zone II")
return new LambertConic(ellipsoid, 45.898919, 47.696014,
46.8, 2.337229, 1 /*0.99987742*/, 600000, 2200000);
else if (name == "(III) France Zone III")
return new LambertConic(ellipsoid, 43.199291, 44.996094,
44.1, 2.337229, 1 /*0.99987750*/, 600000, 3200000);
else if (name == "(IV) France Zone IV")
return new LambertConic(ellipsoid, 41.560388, 42.767663,
42.165, 2.337229, 1 /*0.99994471*/, 234.358, 4185861.369);
else if (name == "(VICGRID) Victoria Australia")
return new LambertConic(ellipsoid, -36, -38, -37, 145, 1,
2500000, 4500000);
else if (name == "(VG94) VICGRID94 Victoria Australia")
return new LambertConic(ellipsoid, -36, -38, -37, 145, 1,
2500000, 2500000);
else {
_errorString = QString("%1: Unknown map projection").arg(name);
return 0;
}
}
const QString &Projection::errorString()
QDebug operator<<(QDebug dbg, const Projection::Setup &setup)
{
return _errorString;
dbg.nospace() << "Setup(" << setup.latitudeOrigin() << ", "
<< setup.longitudeOrigin() << ", " << setup.scale() << ", "
<< setup.falseEasting() << ", " << setup.falseNorthing() << ", "
<< setup.standardParallel1() << ", " << setup.standardParallel2() << ")";
return dbg.space();
}
QDebug operator<<(QDebug dbg, const Projection::Method &method)
{
dbg.nospace() << "Method(" << method.id() << ")";
return dbg.space();
}

64
src/map/projection.h
View File

@ -2,22 +2,53 @@
#define PROJECTION_H
#include <QPointF>
#include <QString>
#include <QDebug>
#include "common/coordinates.h"
class Ellipsoid;
class Datum;
class Projection {
public:
struct Setup {
double latitudeOrigin;
double longitudeOrigin;
double scale;
double falseEasting;
double falseNorthing;
double standardParallel1;
double standardParallel2;
int zone;
class Setup {
public:
Setup() : _latitudeOrigin(NAN), _longitudeOrigin(NAN), _scale(NAN),
_falseEasting(NAN), _falseNorthing(NAN), _standardParallel1(NAN),
_standardParallel2(NAN) {}
Setup(double latitudeOrigin, double longitudeOrigin, double scale,
double falseEasting, double falseNorthing, double standardParallel1,
double standardParallel2) : _latitudeOrigin(latitudeOrigin),
_longitudeOrigin(longitudeOrigin), _scale(scale),
_falseEasting(falseEasting), _falseNorthing(falseNorthing),
_standardParallel1(standardParallel1),
_standardParallel2(standardParallel2) {}
double latitudeOrigin() const {return _latitudeOrigin;}
double longitudeOrigin() const {return _longitudeOrigin;}
double scale() const {return _scale;}
double falseEasting() const {return _falseEasting;}
double falseNorthing() const {return _falseNorthing;}
double standardParallel1() const {return _standardParallel1;}
double standardParallel2() const {return _standardParallel2;}
private:
double _latitudeOrigin;
double _longitudeOrigin;
double _scale;
double _falseEasting;
double _falseNorthing;
double _standardParallel1;
double _standardParallel2;
};
class Method {
public:
Method() : _id(0) {}
Method(int id);
int id() const {return _id;}
bool isNull() const {return (_id == 0);}
private:
int _id;
};
virtual ~Projection() {}
@ -25,12 +56,11 @@ public:
virtual QPointF ll2xy(const Coordinates &c) const = 0;
virtual Coordinates xy2ll(const QPointF &p) const = 0;
static Projection *projection(const QString &name,
const Ellipsoid &ellipsoid, const Setup &setup);
static const QString &errorString();
private:
static QString _errorString;
static Projection *projection(const Datum &datum, const Method &method,
const Setup &setup);
};
QDebug operator<<(QDebug dbg, const Projection::Setup &setup);
QDebug operator<<(QDebug dbg, const Projection::Method &method);
#endif // PROJECTION_H

28
src/map/tifffile.cpp Normal file
View File

@ -0,0 +1,28 @@
#include "tifffile.h"
#define TIFF_II 0x4949
#define TIFF_MM 0x4D4D
#define TIFF_MAGIC 42
bool TIFFFile::readHeader(quint32 &ifd)
{
quint16 endian, magic;
if (QFile::read((char*)&endian, sizeof(endian)) < (qint64)sizeof(endian))
return false;
if (endian == TIFF_II)
_be = false;
else if (endian == TIFF_MM)
_be = true;
else
return false;
if (!readValue(magic))
return false;
if (magic != TIFF_MAGIC)
return false;
if (!readValue(ifd))
return false;
return true;
}

51
src/map/tifffile.h Normal file
View File

@ -0,0 +1,51 @@
#ifndef TIFFFILE_H
#define TIFFFILE_H
#include <QFile>
#include <QtEndian>
class TIFFFile : public QFile
{
public:
TIFFFile() : QFile(), _be(false) {}
TIFFFile(const QString &path) : QFile(path), _be(false) {}
bool readHeader(quint32 &ifd);
template<class T> bool readValue(T &val)
{
T data;
if (QFile::read((char*)&data, sizeof(T)) < (qint64)sizeof(T))
return false;
if (sizeof(T) > 4) {
#if Q_BYTE_ORDER == Q_BIG_ENDIAN
if (_be)
val = data;
else {
for (size_t i = 0; i < sizeof(T); i++)
*((char *)&val + i) = *((char*)&data + sizeof(T) - 1 - i);
}
#else
if (_be) {
for (size_t i = 0; i < sizeof(T); i++)
*((char *)&val + i) = *((char*)&data + sizeof(T) - 1 - i);
} else
val = data;
#endif
} else if (sizeof(T) > 1) {
if (_be)
val = qFromBigEndian(data);
else
val = qFromLittleEndian(data);
} else
val = data;
return true;
}
private:
bool _be;
};
#endif // TIFFFILE_H

74
src/map/transform.cpp Normal file
View File

@ -0,0 +1,74 @@
#include "projection.h"
#include "matrix.h"
#include "transform.h"
void Transform::simple(const QList<ReferencePoint> &points)
{
if (points.at(0).xy.x() == points.at(1).xy.x()
|| points.at(0).xy.y() == points.at(1).xy.y()) {
_errorString = "Invalid reference points tuple";
return;
}
qreal sX, sY, dX, dY;
sX = (points.at(0).xy.x() - points.at(1).xy.x())
/ (points.at(0).pp.x() - points.at(1).pp.x());
sY = (points.at(1).xy.y() - points.at(0).xy.y())
/ (points.at(1).pp.y() - points.at(0).pp.y());
dX = points.at(1).xy.x() - points.at(1).pp.x() * sX;
dY = points.at(0).xy.y() - points.at(0).pp.y() * sY;
_transform = QTransform(sX, 0, 0, sY, dX, dY);
}
void Transform::affine(const QList<ReferencePoint> &points)
{
Matrix c(3, 2);
c.zeroize();
for (size_t i = 0; i < c.h(); i++) {
for (size_t j = 0; j < c.w(); j++) {
for (int k = 0; k < points.size(); k++) {
double f[3], t[2];
f[0] = points.at(k).pp.x();
f[1] = points.at(k).pp.y();
f[2] = 1.0;
t[0] = points.at(k).xy.x();
t[1] = points.at(k).xy.y();
c.m(i,j) += f[i] * t[j];
}
}
}
Matrix Q(3, 3);
Q.zeroize();
for (int qi = 0; qi < points.size(); qi++) {
double v[3];
v[0] = points.at(qi).pp.x();
v[1] = points.at(qi).pp.y();
v[2] = 1.0;
for (size_t i = 0; i < Q.h(); i++)
for (size_t j = 0; j < Q.w(); j++)
Q.m(i,j) += v[i] * v[j];
}
Matrix M = Q.augemented(c);
if (!M.eliminate()) {
_errorString = "Singular transformation matrix";
return;
}
_transform = QTransform(M.m(0,3), M.m(0,4), M.m(1,3), M.m(1,4), M.m(2,3),
M.m(2,4));
}
Transform::Transform(const QList<ReferencePoint> &points)
{
if (points.count() < 2)
_errorString = "Insufficient number of reference points";
else if (points.size() == 2)
simple(points);
else
affine(points);
}

30
src/map/transform.h Normal file
View File

@ -0,0 +1,30 @@
#ifndef TRANSFORM_H
#define TRANSFORM_H
#include <QTransform>
#include <QList>
#include "common/coordinates.h"
struct ReferencePoint {
QPoint xy;
QPointF pp;
};
class Transform
{
public:
Transform(const QList<ReferencePoint> &points);
bool isNull() {return _transform.type() == QTransform::TxNone;}
const QString &errorString() const {return _errorString;}
const QTransform &transform() const {return _transform;}
private:
void simple(const QList<ReferencePoint> &points);
void affine(const QList<ReferencePoint> &points);
QTransform _transform;
QString _errorString;
};
#endif // TRANSFORM_H

6
src/map/utm.cpp
View File

@ -1,10 +1,10 @@
#include "ellipsoid.h"
#include "utm.h"
UTM::UTM(const Ellipsoid &ellipsoid, int zone)
: _tm(ellipsoid, 0, (qAbs(zone) - 1)*6 - 180 + 3, 0.9996, 500000,
zone < 0 ? 10000000 : 0)
Projection::Setup UTM::setup(int zone)
{
return Projection::Setup(0, (qAbs(zone) - 1)*6 - 180 + 3, 0.9996, 500000,
zone < 0 ? 10000000 : 0, 0, 0);
}
int UTM::zone(const Coordinates &c)

14
src/map/utm.h
View File

@ -2,22 +2,12 @@
#define UTM_H
#include "projection.h"
#include "transversemercator.h"
class UTM : public Projection
class UTM
{
public:
UTM(const Ellipsoid &ellipsoid, int zone);
virtual QPointF ll2xy(const Coordinates &c) const
{return _tm.ll2xy(c);}
virtual Coordinates xy2ll(const QPointF &p) const
{return _tm.xy2ll(p);}
static int zone(const Coordinates &c);
private:
TransverseMercator _tm;
static Projection::Setup setup(int zone);
};
#endif // UTM_H