src/common/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 <QDir>
  16 #include <QFile>
  17 #include <QRegularExpression>
  18 #include <QLocale>
  19 #include <private/qzipreader_p.h>
  20 #include "rectc.h"
  21 #include "dem.h"
  22
  23
  24 static unsigned int isqrt(unsigned int x)
  25 {
  26         unsigned int r = 0;
  27
  28         while ((r + 1) * (r + 1) <= x)
  29                 r++;
  30
  31         return r;
  32 }
  33
  34 static double interpolate(double dx, double dy, double p0, double p1, double p2,
  35   double p3)
  36 {
  37         return p0 * (1.0 - dx) * (1.0 - dy) + p1 * dx * (1.0 - dy)
  38           + p2 * dy * (1.0 - dx) + p3 * dx * dy;
  39 }
  40
  41 static double value(int col, int row, int samples, const QByteArray &data)
  42 {
  43         int pos = ((samples - 1 - row) * samples + col) * 2;
  44         qint16 val = qFromBigEndian(*((const qint16*)(data.constData() + pos)));
  45
  46         return (val == -32768) ? NAN : val;
  47 }
  48
  49 QMutex DEM::_lock;
  50
  51 DEM::Entry::Entry(const QByteArray &data) : _data(data)
  52 {
  53         _samples = isqrt(_data.size() / 2);
  54 }
  55
  56 QString DEM::Tile::latStr() const
  57 {
  58         const char ns = (_lat >= 0) ? 'N' : 'S';
  59         return QString("%1%2").arg(ns).arg(qAbs(_lat), 2, 10, QChar('0'));
  60 }
  61
  62 QString DEM::Tile::lonStr() const
  63 {
  64         const char ew = (_lon >= 0) ? 'E' : 'W';
  65         return QString("%1%2").arg(ew).arg(qAbs(_lon), 3, 10, QChar('0'));
  66 }
  67
  68 QString DEM::Tile::baseName() const
  69 {
  70         return QString("%1%2.hgt").arg(latStr(), lonStr());
  71 }
  72
  73 QString DEM::_dir;
  74 DEM::TileCache DEM::_data;
  75
  76 void DEM::setCacheSize(int size)
  77 {
  78         _data.setMaxCost(size * 1024);
  79 }
  80
  81 void DEM::setDir(const QString &path)
  82 {
  83         _dir = path;
  84 }
  85
  86 void DEM::clearCache()
  87 {
  88         _data.clear();
  89 }
  90
  91 double DEM::height(const Coordinates &c, const Entry *e)
  92 {
  93         if (!e->samples())
  94                 return NAN;
  95
  96         double lat = (c.lat() - floor(c.lat())) * (e->samples() - 1);
  97         double lon = (c.lon() - floor(c.lon())) * (e->samples() - 1);
  98         int row = (int)lat;
  99         int col = (int)lon;
 100
 101         double p0 = value(col, row, e->samples(), e->data());
 102         double p1 = value(col + 1, row, e->samples(), e->data());
 103         double p2 = value(col, row + 1, e->samples(), e->data());
 104         double p3 = value(col + 1, row + 1, e->samples(), e->data());
 105
 106         return interpolate(lon - col, lat - row, p0, p1, p2, p3);
 107 }
 108
 109 DEM::Entry *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 Entry(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 Entry();
 124                 } else
 125                         return new Entry(file.readAll());
 126         }
 127 }
 128
 129 double DEM::elevation(const Coordinates &c)
 130 {
 131         if (_dir.isEmpty())
 132                 return NAN;
 133
 134         Tile tile(floor(c.lon()), floor(c.lat()));
 135         Entry *e = _data.object(tile);
 136         double ele;
 137
 138         if (!e) {
 139                 e = loadTile(tile);
 140                 ele = height(c, e);
 141                 _data.insert(tile, e, e->data().size());
 142         } else
 143                 ele = height(c, e);
 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                 area.setStyle(PolygonStyle(QColor(0xFF, 0, 0, 0x40),
 174                   QColor(0xFF, 0, 0, 0x80), 2));
 175
 176                 list.append(area);
 177         }
 178
 179         return list;
 180 }
 181
 182 #ifndef QT_NO_DEBUG
 183 QDebug operator<<(QDebug dbg, const DEM::Tile &tile)
 184 {
 185         dbg.nospace() << "Tile(" << tile.baseName() << ")";
 186         return dbg.space();
 187 }
 188 #endif // QT_NO_DEBUG