src/data/dem.cpp

   1 /*
   2         WARNING: This code uses internal Qt API - the QZipReader class for reading
   3         ZIP files - and things may break if Qt changes the API. For Qt5 this is not
   4         a problem as we can "see the future" now and there are no changes in all
   5         the supported Qt5 versions up to the last one (5.15). In Qt6 the class
   6         might change or even disappear in the future, but this is very unlikely
   7         as there were no changes for several years and The Qt Company's policy
   8         is: "do not invest any resources into any desktop related stuff unless
   9         absolutely necessary". There is an issue (QTBUG-3897) since the year 2009 to
  10         include the ZIP reader into the public API, which aptly illustrates the
  11         effort The Qt Company is willing to make about anything desktop related...
  12 */
  13
  14 #include <QtEndian>
  15 #include <QtMath>
  16 #include <QDir>
  17 #include <QFile>
  18 #include <QRegularExpression>
  19 #include <QLocale>
  20 #include <private/qzipreader_p.h>
  21 #include "common/rectc.h"
  22 #include "dem.h"
  23
  24 #define SRTM3_SAMPLES  1201
  25 #define SRTM1_SAMPLES  3601
  26 #define SRTM05_SAMPLES 7201
  27
  28 #define SRTM_SIZE(samples) \
  29         ((samples) * (samples) * 2)
  30
  31 static qreal interpolate(qreal dx, qreal dy, qreal p0, qreal p1, qreal p2,
  32   qreal p3)
  33 {
  34         return p0 * (1 - dx) * (1 - dy) + p1 * dx * (1 - dy) + p2 * dy * (1 - dx)
  35           + p3 * dx * dy;
  36 }
  37
  38 static qreal value(int col, int row, int samples, const QByteArray *data)
  39 {
  40         int pos = ((samples - 1 - row) * samples + col) * 2;
  41         qint16 val = qFromBigEndian(*((const qint16*)(data->constData() + pos)));
  42
  43         return (val == -32768) ? NAN : val;
  44 }
  45
  46 static qreal height(const Coordinates &c, const QByteArray *data)
  47 {
  48         int samples;
  49
  50         if (data->size() == SRTM_SIZE(SRTM3_SAMPLES))
  51                 samples = SRTM3_SAMPLES;
  52         else if (data->size() == SRTM_SIZE(SRTM1_SAMPLES))
  53                 samples = SRTM1_SAMPLES;
  54         else if (data->size() == SRTM_SIZE(SRTM05_SAMPLES))
  55                 samples = SRTM05_SAMPLES;
  56         else
  57                 return NAN;
  58
  59         int row = (int)((c.lat() - qFloor(c.lat())) * (samples - 1));
  60         int col = (int)((c.lon() - qFloor(c.lon())) * (samples - 1));
  61         qreal dx = ((c.lon() - qFloor(c.lon())) * (samples - 1)) - col;
  62         qreal dy = ((c.lat() - qFloor(c.lat())) * (samples - 1)) - row;
  63
  64         qreal p0 = value(col, row, samples, data);
  65         qreal p1 = value(col + 1, row, samples, data);
  66         qreal p2 = value(col, row + 1, samples, data);
  67         qreal p3 = value(col + 1, row + 1, samples, data);
  68
  69         return interpolate(dx, dy, p0, p1, p2, p3);
  70 }
  71
  72 QString DEM::Tile::latStr() const
  73 {
  74         const char ns = (_lat >= 0) ? 'N' : 'S';
  75         return QString("%1%2").arg(ns).arg(qAbs(_lat), 2, 10, QChar('0'));
  76 }
  77
  78 QString DEM::Tile::lonStr() const
  79 {
  80         const char ew = (_lon >= 0) ? 'E' : 'W';
  81         return QString("%1%2").arg(ew).arg(qAbs(_lon), 3, 10, QChar('0'));
  82 }
  83
  84 QString DEM::Tile::baseName() const
  85 {
  86         return QString("%1%2.hgt").arg(latStr(), lonStr());
  87 }
  88
  89 QString DEM::_dir;
  90 DEM::TileCache DEM::_data;
  91
  92 QString DEM::fileName(const QString &baseName)
  93 {
  94         return QDir(_dir).absoluteFilePath(baseName);
  95 }
  96
  97 void DEM::setCacheSize(int size)
  98 {
  99         _data.setMaxCost(size);
 100 }
 101
 102 void DEM::setDir(const QString &path)
 103 {
 104         _dir = path;
 105 }
 106
 107 void DEM::clearCache()
 108 {
 109         _data.clear();
 110 }
 111
 112 qreal DEM::elevation(const Coordinates &c)
 113 {
 114         if (_dir.isEmpty())
 115                 return NAN;
 116
 117         Tile tile(qFloor(c.lon()), qFloor(c.lat()));
 118
 119         QByteArray *ba = _data[tile];
 120         if (!ba) {
 121                 QString bn(tile.baseName());
 122                 QString fn(fileName(bn));
 123                 QString zn(fn + ".zip");
 124
 125                 if (QFileInfo::exists(zn)) {
 126                         QZipReader zip(zn, QIODevice::ReadOnly);
 127                         ba = new QByteArray(zip.fileData(bn));
 128                         qreal ele = height(c, ba);
 129                         _data.insert(tile, ba, ba->size());
 130                         return ele;
 131                 } else {
 132                         QFile file(fn);
 133                         if (!file.open(QIODevice::ReadOnly)) {
 134                                 qWarning("%s: %s", qPrintable(file.fileName()),
 135                                   qPrintable(file.errorString()));
 136                                 _data.insert(tile, new QByteArray());
 137                                 return NAN;
 138                         } else {
 139                                 ba = new QByteArray(file.readAll());
 140                                 qreal ele = height(c, ba);
 141                                 _data.insert(tile, ba, ba->size());
 142                                 return ele;
 143                         }
 144                 }
 145         } else
 146                 return height(c, ba);
 147 }
 148
 149 QList<Area> DEM::tiles()
 150 {
 151         QDir dir(_dir);
 152         QFileInfoList files(dir.entryInfoList(QDir::Files | QDir::Readable));
 153         QRegularExpression re("([NS])([0-9]{2})([EW])([0-9]{3})");
 154         QLocale l(QLocale::system());
 155         QList<Area> list;
 156
 157         for (int i = 0; i < files.size(); i++) {
 158                 QString basename(files.at(i).baseName());
 159                 QRegularExpressionMatch match(re.match(basename));
 160                 if (!match.hasMatch())
 161                         continue;
 162
 163                 int lat = match.captured(2).toInt();
 164                 int lon = match.captured(4).toInt();
 165                 if (match.captured(1) == "S")
 166                         lat = -lat;
 167                 if (match.captured(3) == "W")
 168                         lon = -lon;
 169
 170                 Area area(RectC(Coordinates(lon, lat + 1), Coordinates(lon + 1, lat)));
 171                 area.setName(basename);
 172                 area.setDescription(files.at(i).suffix().toUpper() + ", "
 173                   + l.formattedDataSize(files.at(i).size()));
 174
 175                 list.append(area);
 176         }
 177
 178         return list;
 179 }
 180
 181 #ifndef QT_NO_DEBUG
 182 QDebug operator<<(QDebug dbg, const DEM::Tile &tile)
 183 {
 184         dbg.nospace() << "Tile(" << tile.baseName() << ")";
 185         return dbg.space();
 186 }
 187 #endif // QT_NO_DEBUG