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 QString DEM::fileName(const QString &baseName)
  95 {
  96         return QDir(_dir).absoluteFilePath(baseName);
  97 }
  98
  99 void DEM::setCacheSize(int size)
 100 {
 101         _data.setMaxCost(size * 1024);
 102 }
 103
 104 void DEM::setDir(const QString &path)
 105 {
 106         _dir = path;
 107 }
 108
 109 void DEM::clearCache()
 110 {
 111         _data.clear();
 112 }
 113
 114 double DEM::elevation(const Coordinates &c)
 115 {
 116         if (_dir.isEmpty())
 117                 return NAN;
 118
 119         Tile tile(qFloor(c.lon()), qFloor(c.lat()));
 120
 121         QByteArray *ba = _data[tile];
 122         if (!ba) {
 123                 QString bn(tile.baseName());
 124                 QString fn(fileName(bn));
 125                 QString zn(fn + ".zip");
 126
 127                 if (QFileInfo::exists(zn)) {
 128                         QZipReader zip(zn, QIODevice::ReadOnly);
 129                         ba = new QByteArray(zip.fileData(bn));
 130                         double ele = height(c, ba);
 131                         _data.insert(tile, ba, ba->size());
 132                         return ele;
 133                 } else {
 134                         QFile file(fn);
 135                         if (!file.open(QIODevice::ReadOnly)) {
 136                                 qWarning("%s: %s", qPrintable(file.fileName()),
 137                                   qPrintable(file.errorString()));
 138                                 _data.insert(tile, new QByteArray());
 139                                 return NAN;
 140                         } else {
 141                                 ba = new QByteArray(file.readAll());
 142                                 double ele = height(c, ba);
 143                                 _data.insert(tile, ba, ba->size());
 144                                 return ele;
 145                         }
 146                 }
 147         } else
 148                 return height(c, ba);
 149 }
 150
 151 QList<Area> DEM::tiles()
 152 {
 153         static const QRegularExpression re("([NS])([0-9]{2})([EW])([0-9]{3})");
 154         QDir dir(_dir);
 155         QFileInfoList files(dir.entryInfoList(QDir::Files | QDir::Readable));
 156         QLocale l(QLocale::system());
 157         QList<Area> list;
 158
 159         for (int i = 0; i < files.size(); i++) {
 160                 QString basename(files.at(i).baseName());
 161                 QRegularExpressionMatch match(re.match(basename));
 162                 if (!match.hasMatch())
 163                         continue;
 164
 165                 int lat = match.captured(2).toInt();
 166                 int lon = match.captured(4).toInt();
 167                 if (match.captured(1) == "S")
 168                         lat = -lat;
 169                 if (match.captured(3) == "W")
 170                         lon = -lon;
 171
 172                 Area area(RectC(Coordinates(lon, lat + 1), Coordinates(lon + 1, lat)));
 173                 area.setName(basename);
 174                 area.setDescription(files.at(i).suffix().toUpper() + ", "
 175                   + l.formattedDataSize(files.at(i).size()));
 176
 177                 list.append(area);
 178         }
 179
 180         return list;
 181 }
 182
 183 #ifndef QT_NO_DEBUG
 184 QDebug operator<<(QDebug dbg, const DEM::Tile &tile)
 185 {
 186         dbg.nospace() << "Tile(" << tile.baseName() << ")";
 187         return dbg.space();
 188 }
 189 #endif // QT_NO_DEBUG