mirror of
https://github.com/tumic0/GPXSee.git
synced 2025-06-27 19:49:15 +02:00
Show the map elevation from the map DEM if available
This commit is contained in:
110
src/map/IMG/dem.cpp
Normal file
110
src/map/IMG/dem.cpp
Normal file
@ -0,0 +1,110 @@
|
||||
#include "dem.h"
|
||||
|
||||
using namespace IMG;
|
||||
|
||||
static double interpolate(double dx, double dy, double p0, double p1, double p2,
|
||||
double p3)
|
||||
{
|
||||
return p0 * (1.0 - dx) * (1.0 - dy) + p1 * dx * (1.0 - dy)
|
||||
+ p2 * dy * (1.0 - dx) + p3 * dx * dy;
|
||||
}
|
||||
|
||||
static double val(const Matrix<qint16> &m, int row, int col)
|
||||
{
|
||||
qint16 v = m.at(row, col);
|
||||
return (v == -32768) ? NAN : (double)v;
|
||||
}
|
||||
|
||||
bool DEM::edgeCb(const MapData::Elevation *e, void *context)
|
||||
{
|
||||
EdgeCTX *ctx = (EdgeCTX*)context;
|
||||
|
||||
double x = (ctx->c.lon() - e->rect.left()) / e->xr;
|
||||
double y = (e->rect.top() - ctx->c.lat()) / e->yr;
|
||||
int row = qMin((int)y, e->m.h() - 1);
|
||||
int col = qMin((int)x, e->m.w() - 1);
|
||||
|
||||
ctx->ele = val(e->m, row, col);
|
||||
|
||||
return std::isnan(ctx->ele);
|
||||
}
|
||||
|
||||
double DEM::edge(const DEMTRee &tree, const Coordinates &c)
|
||||
{
|
||||
double min[2], max[2];
|
||||
double ele = NAN;
|
||||
EdgeCTX ctx(c, ele);
|
||||
|
||||
min[0] = c.lon();
|
||||
min[1] = c.lat();
|
||||
max[0] = c.lon();
|
||||
max[1] = c.lat();
|
||||
|
||||
tree.Search(min, max, edgeCb, &ctx);
|
||||
|
||||
return ele;
|
||||
}
|
||||
|
||||
double DEM::elevation(const DEMTRee &tree, const MapData::Elevation *e,
|
||||
const Coordinates &c)
|
||||
{
|
||||
double x = (c.lon() - e->rect.left()) / e->xr;
|
||||
double y = (e->rect.top() - c.lat()) / e->yr;
|
||||
int row = qMin((int)y, e->m.h() - 1);
|
||||
int col = qMin((int)x, e->m.w() - 1);
|
||||
|
||||
double p0 = val(e->m, row, col);
|
||||
double p1 = (col == e->m.w() - 1)
|
||||
? edge(tree, Coordinates(e->rect.left() + (col + 1) * e->xr + e->xr/2,
|
||||
e->rect.top() - row * e->yr))
|
||||
: val(e->m, row, col + 1);
|
||||
double p2 = (row == e->m.h() - 1)
|
||||
? edge(tree, Coordinates(e->rect.left() + col * e->xr,
|
||||
e->rect.top() - (row + 1) * e->yr - e->yr/2))
|
||||
: val(e->m, row + 1, col);
|
||||
double p3 = ((col == e->m.w() - 1) || (row == e->m.h() - 1))
|
||||
? edge(tree, Coordinates(e->rect.left() + (col + 1) * e->xr + e->xr/2,
|
||||
e->rect.top() - (row + 1) * e->yr - e->yr/2))
|
||||
: val(e->m, row + 1, col + 1);
|
||||
|
||||
return interpolate(x - col, y - row, p0, p1, p2, p3);
|
||||
}
|
||||
|
||||
void DEM::buildTree(const QList<MapData::Elevation> &tiles, DEMTRee &tree)
|
||||
{
|
||||
double min[2], max[2];
|
||||
|
||||
for (int i = 0; i < tiles.size(); i++) {
|
||||
const MapData::Elevation &e = tiles.at(i);
|
||||
|
||||
min[0] = e.rect.left();
|
||||
min[1] = e.rect.bottom();
|
||||
max[0] = e.rect.right();
|
||||
max[1] = e.rect.top();
|
||||
|
||||
tree.Insert(min, max, &e);
|
||||
}
|
||||
}
|
||||
|
||||
bool DEM::elevationCb(const MapData::Elevation *e, void *context)
|
||||
{
|
||||
ElevationCTX *ctx = (ElevationCTX*)context;
|
||||
|
||||
ctx->ele = elevation(ctx->tree, e, ctx->c);
|
||||
return std::isnan(ctx->ele);
|
||||
}
|
||||
|
||||
void DEM::searchTree(const DEMTRee &tree, const Coordinates &c,
|
||||
double &ele)
|
||||
{
|
||||
double min[2], max[2];
|
||||
ElevationCTX ctx(tree, c, ele);
|
||||
|
||||
min[0] = c.lon();
|
||||
min[1] = c.lat();
|
||||
max[0] = c.lon();
|
||||
max[1] = c.lat();
|
||||
|
||||
ele = NAN;
|
||||
tree.Search(min, max, elevationCb, &ctx);
|
||||
}
|
||||
46
src/map/IMG/dem.h
Normal file
46
src/map/IMG/dem.h
Normal file
@ -0,0 +1,46 @@
|
||||
#ifndef IMG_DEM_H
|
||||
#define IMG_DEM_H
|
||||
|
||||
#include "common/rtree.h"
|
||||
#include "mapdata.h"
|
||||
|
||||
namespace IMG {
|
||||
|
||||
class DEM {
|
||||
public:
|
||||
|
||||
typedef RTree<const MapData::Elevation*, double, 2> DEMTRee;
|
||||
|
||||
static void buildTree(const QList<MapData::Elevation> &tiles, DEMTRee &tree);
|
||||
static void searchTree(const DEMTRee &tree, const Coordinates &c,
|
||||
double &ele);
|
||||
|
||||
private:
|
||||
struct ElevationCTX {
|
||||
ElevationCTX(const DEMTRee &tree, const Coordinates &c, double &ele)
|
||||
: tree(tree), c(c), ele(ele) {}
|
||||
|
||||
const DEMTRee &tree;
|
||||
const Coordinates &c;
|
||||
double &ele;
|
||||
};
|
||||
|
||||
struct EdgeCTX {
|
||||
EdgeCTX(const Coordinates &c, double &ele) : c(c), ele(ele) {}
|
||||
|
||||
const Coordinates &c;
|
||||
double &ele;
|
||||
};
|
||||
|
||||
MatrixD elevation() const;
|
||||
|
||||
static double edge(const DEMTRee &tree, const Coordinates &c);
|
||||
static double elevation(const DEMTRee &tree, const MapData::Elevation *e,
|
||||
const Coordinates &c);
|
||||
static bool elevationCb(const MapData::Elevation *e, void *context);
|
||||
static bool edgeCb(const MapData::Elevation *e, void *context);
|
||||
};
|
||||
|
||||
}
|
||||
|
||||
#endif // IMG_ELEVATIONTREE_H
|
||||
@ -9,6 +9,7 @@
|
||||
#include "map/hillshading.h"
|
||||
#include "style.h"
|
||||
#include "lblfile.h"
|
||||
#include "dem.h"
|
||||
#include "rastertile.h"
|
||||
|
||||
using namespace IMG;
|
||||
@ -108,20 +109,6 @@ static bool rectNearPolygon(const QPolygonF &polygon, const QRectF &rect)
|
||||
|| polygon.containsPoint(rect.bottomRight(), Qt::OddEvenFill)));
|
||||
}
|
||||
|
||||
static double interpolate(double dx, double dy, double p0, double p1, double p2,
|
||||
double p3)
|
||||
{
|
||||
return p0 * (1.0 - dx) * (1.0 - dy) + p1 * dx * (1.0 - dy)
|
||||
+ p2 * dy * (1.0 - dx) + p3 * dx * dy;
|
||||
}
|
||||
|
||||
static double val(const Matrix<qint16> &m, int row, int col)
|
||||
{
|
||||
qint16 v = m.at(row, col);
|
||||
return (v == -32768) ? NAN : (double)v;
|
||||
}
|
||||
|
||||
|
||||
const QFont *RasterTile::poiFont(Style::FontSize size, int zoom,
|
||||
bool extended) const
|
||||
{
|
||||
@ -465,100 +452,6 @@ void RasterTile::fetchData(QList<MapData::Poly> &polygons,
|
||||
_data->points(pointRectD.toRectC(_proj, 20), _zoom, &points);
|
||||
}
|
||||
|
||||
bool RasterTile::edgeCb(const MapData::Elevation *e, void *context)
|
||||
{
|
||||
EdgeCTX *ctx = (EdgeCTX*)context;
|
||||
|
||||
double x = (ctx->c.lon() - e->rect.left()) / e->xr;
|
||||
double y = (e->rect.top() - ctx->c.lat()) / e->yr;
|
||||
int row = qMin((int)y, e->m.h() - 1);
|
||||
int col = qMin((int)x, e->m.w() - 1);
|
||||
|
||||
ctx->ele = val(e->m, row, col);
|
||||
|
||||
return std::isnan(ctx->ele);
|
||||
}
|
||||
|
||||
double RasterTile::edge(const DEMTRee &tree, const Coordinates &c)
|
||||
{
|
||||
double min[2], max[2];
|
||||
double ele = NAN;
|
||||
EdgeCTX ctx(c, ele);
|
||||
|
||||
min[0] = c.lon();
|
||||
min[1] = c.lat();
|
||||
max[0] = c.lon();
|
||||
max[1] = c.lat();
|
||||
|
||||
tree.Search(min, max, edgeCb, &ctx);
|
||||
|
||||
return ele;
|
||||
}
|
||||
|
||||
double RasterTile::elevation(const DEMTRee &tree, const MapData::Elevation *e,
|
||||
const Coordinates &c)
|
||||
{
|
||||
double x = (c.lon() - e->rect.left()) / e->xr;
|
||||
double y = (e->rect.top() - c.lat()) / e->yr;
|
||||
int row = qMin((int)y, e->m.h() - 1);
|
||||
int col = qMin((int)x, e->m.w() - 1);
|
||||
|
||||
double p0 = val(e->m, row, col);
|
||||
double p1 = (col == e->m.w() - 1)
|
||||
? edge(tree, Coordinates(e->rect.left() + (col + 1) * e->xr + e->xr/2,
|
||||
e->rect.top() - row * e->yr))
|
||||
: val(e->m, row, col + 1);
|
||||
double p2 = (row == e->m.h() - 1)
|
||||
? edge(tree, Coordinates(e->rect.left() + col * e->xr,
|
||||
e->rect.top() - (row + 1) * e->yr - e->yr/2))
|
||||
: val(e->m, row + 1, col);
|
||||
double p3 = ((col == e->m.w() - 1) || (row == e->m.h() - 1))
|
||||
? edge(tree, Coordinates(e->rect.left() + (col + 1) * e->xr + e->xr/2,
|
||||
e->rect.top() - (row + 1) * e->yr - e->yr/2))
|
||||
: val(e->m, row + 1, col + 1);
|
||||
|
||||
return interpolate(x - col, y - row, p0, p1, p2, p3);
|
||||
}
|
||||
|
||||
void RasterTile::buildTree(const QList<MapData::Elevation> &tiles, DEMTRee &tree)
|
||||
{
|
||||
double min[2], max[2];
|
||||
|
||||
for (int i = 0; i < tiles.size(); i++) {
|
||||
const MapData::Elevation &e = tiles.at(i);
|
||||
|
||||
min[0] = e.rect.left();
|
||||
min[1] = e.rect.bottom();
|
||||
max[0] = e.rect.right();
|
||||
max[1] = e.rect.top();
|
||||
|
||||
tree.Insert(min, max, &e);
|
||||
}
|
||||
}
|
||||
|
||||
bool RasterTile::elevationCb(const MapData::Elevation *e, void *context)
|
||||
{
|
||||
ElevationCTX *ctx = (ElevationCTX*)context;
|
||||
|
||||
ctx->ele = elevation(ctx->tree, e, ctx->c);
|
||||
return std::isnan(ctx->ele);
|
||||
}
|
||||
|
||||
void RasterTile::searchTree(const DEMTRee &tree, const Coordinates &c,
|
||||
double &ele)
|
||||
{
|
||||
double min[2], max[2];
|
||||
ElevationCTX ctx(tree, c, ele);
|
||||
|
||||
min[0] = c.lon();
|
||||
min[1] = c.lat();
|
||||
max[0] = c.lon();
|
||||
max[1] = c.lat();
|
||||
|
||||
ele = NAN;
|
||||
tree.Search(min, max, elevationCb, &ctx);
|
||||
}
|
||||
|
||||
MatrixD RasterTile::elevation() const
|
||||
{
|
||||
MatrixD m(_rect.height() + 2, _rect.width() + 2);
|
||||
@ -586,14 +479,14 @@ MatrixD RasterTile::elevation() const
|
||||
|
||||
_data->elevations(rect, _zoom, &tiles);
|
||||
|
||||
buildTree(tiles, tree);
|
||||
DEM::buildTree(tiles, tree);
|
||||
for (int i = 0; i < ll.size(); i++)
|
||||
searchTree(tree, ll.at(i), m.at(i));
|
||||
DEM::searchTree(tree, ll.at(i), m.at(i));
|
||||
} else {
|
||||
DEM::lock();
|
||||
::DEM::lock();
|
||||
for (int i = 0; i < ll.size(); i++)
|
||||
m.at(i) = DEM::elevation(ll.at(i));
|
||||
DEM::unlock();
|
||||
m.at(i) = ::DEM::elevation(ll.at(i));
|
||||
::DEM::unlock();
|
||||
}
|
||||
|
||||
return m;
|
||||
|
||||
@ -77,15 +77,6 @@ private:
|
||||
|
||||
MatrixD elevation() const;
|
||||
|
||||
static double edge(const DEMTRee &tree, const Coordinates &c);
|
||||
static double elevation(const DEMTRee &tree, const MapData::Elevation *e,
|
||||
const Coordinates &c);
|
||||
static void buildTree(const QList<MapData::Elevation> &tiles, DEMTRee &tree);
|
||||
static void searchTree(const DEMTRee &tree, const Coordinates &c,
|
||||
double &ele);
|
||||
static bool elevationCb(const MapData::Elevation *e, void *context);
|
||||
static bool edgeCb(const MapData::Elevation *e, void *context);
|
||||
|
||||
Projection _proj;
|
||||
Transform _transform;
|
||||
MapData *_data;
|
||||
|
||||
Reference in New Issue
Block a user