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 double interpolate(double dx, double dy, double p0, double p1, double p2,
  32   double p3)
  33 {
  34         return p0 * (1.0 - dx) * (1.0 - dy) + p1 * dx * (1.0 - dy)
  35           + p2 * dy * (1.0 - dx) + p3 * dx * dy;
  36 }
  37
  38 static double 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 double 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         double lat = (c.lat() - floor(c.lat())) * (samples - 1);
  60         double lon = (c.lon() - floor(c.lon())) * (samples - 1);
  61         int row = (int)lat;
  62         int col = (int)lon;
  63
  64         double p0 = value(col, row, samples, data);
  65         double p1 = value(col + 1, row, samples, data);
  66         double p2 = value(col, row + 1, samples, data);
  67         double p3 = value(col + 1, row + 1, samples, data);
  68
  69         return interpolate(lon - col, lat - row, p0, p1, p2, p3);
  70 }
  71
  72 QMutex DEM::_lock;
  73
  74 QString DEM::Tile::latStr() const
  75 {
  76         const char ns = (_lat >= 0) ? 'N' : 'S';
  77         return QString("%1%2").arg(ns).arg(qAbs(_lat), 2, 10, QChar('0'));
  78 }
  79
  80 QString DEM::Tile::lonStr() const
  81 {
  82         const char ew = (_lon >= 0) ? 'E' : 'W';
  83         return QString("%1%2").arg(ew).arg(qAbs(_lon), 3, 10, QChar('0'));
  84 }
  85
  86 QString DEM::Tile::baseName() const
  87 {
  88         return QString("%1%2.hgt").arg(latStr(), lonStr());
  89 }
  90
  91 QString DEM::_dir;
  92 DEM::TileCache DEM::_data;
  93
  94 void DEM::setCacheSize(int size)
  95 {
  96         _data.setMaxCost(size * 1024);
  97 }
  98
  99 void DEM::setDir(const QString &path)
 100 {
 101         _dir = path;
 102 }
 103
 104 void DEM::clearCache()
 105 {
 106         _data.clear();
 107 }
 108
 109 QByteArray *DEM::loadTile(const Tile &tile)
 110 {
 111         QString bn(tile.baseName());
 112         QString fn(QDir(_dir).absoluteFilePath(bn));
 113         QString zn(fn + ".zip");
 114
 115         if (QFileInfo::exists(zn)) {
 116                 QZipReader zip(zn, QIODevice::ReadOnly);
 117                 return new QByteArray(zip.fileData(bn));
 118         } else {
 119                 QFile file(fn);
 120                 if (!file.open(QIODevice::ReadOnly)) {
 121                         qWarning("%s: %s", qPrintable(file.fileName()),
 122                           qPrintable(file.errorString()));
 123                         return new QByteArray();
 124                 } else
 125                         return new QByteArray(file.readAll());
 126         }
 127 }
 128
 129 double DEM::elevation(const Coordinates &c)
 130 {
 131         if (_dir.isEmpty())
 132                 return NAN;
 133
 134         Tile tile(qFloor(c.lon()), qFloor(c.lat()));
 135         QByteArray *ba = _data.object(tile);
 136         double ele;
 137
 138         if (!ba) {
 139                 ba = loadTile(tile);
 140                 ele = height(c, ba);
 141                 _data.insert(tile, ba, ba->size());
 142         } else
 143                 ele = height(c, ba);
 144
 145         return ele;
 146 }
 147
 148 QList<Area> DEM::tiles()
 149 {
 150         static const QRegularExpression re("([NS])([0-9]{2})([EW])([0-9]{3})");
 151         QDir dir(_dir);
 152         QFileInfoList files(dir.entryInfoList(QDir::Files | QDir::Readable));
 153         QLocale l(QLocale::system());
 154         QList<Area> list;
 155
 156         for (int i = 0; i < files.size(); i++) {
 157                 QString basename(files.at(i).baseName());
 158                 QRegularExpressionMatch match(re.match(basename));
 159                 if (!match.hasMatch())
 160                         continue;
 161
 162                 int lat = match.captured(2).toInt();
 163                 int lon = match.captured(4).toInt();
 164                 if (match.captured(1) == "S")
 165                         lat = -lat;
 166                 if (match.captured(3) == "W")
 167                         lon = -lon;
 168
 169                 Area area(RectC(Coordinates(lon, lat + 1), Coordinates(lon + 1, lat)));
 170                 area.setName(basename);
 171                 area.setDescription(files.at(i).suffix().toUpper() + ", "
 172                   + l.formattedDataSize(files.at(i).size()));
 173
 174                 list.append(area);
 175         }
 176
 177         return list;
 178 }
 179
 180 #ifndef QT_NO_DEBUG
 181 QDebug operator<<(QDebug dbg, const DEM::Tile &tile)
 182 {
 183         dbg.nospace() << "Tile(" << tile.baseName() << ")";
 184         return dbg.space();
 185 }
 186 #endif // QT_NO_DEBUG