#include "generator/srtm_parser.hpp" #include "platform/platform.hpp" #include "coding/endianness.hpp" #include "coding/zip_reader.hpp" #include "base/file_name_utils.hpp" #include "base/logging.hpp" #include #include #include #include namespace generator { namespace { int constexpr kArcSecondsInDegree = 60 * 60; size_t constexpr kSrtmTileSize = (kArcSecondsInDegree + 1) * (kArcSecondsInDegree + 1) * 2; struct UnzipMemDelegate : public ZipFileReader::Delegate { explicit UnzipMemDelegate(std::vector & buffer) : m_buffer(buffer), m_completed(false) { m_buffer.reserve(kSrtmTileSize); } // ZipFileReader::Delegate overrides: void OnBlockUnzipped(size_t size, char const * data) override { m_buffer.insert(m_buffer.end(), data, data + size); } void OnStarted() override { m_buffer.clear(); m_completed = false; } void OnCompleted() override { m_completed = true; } std::vector & m_buffer; bool m_completed; }; std::string GetSrtmContFileName(std::string const & dir, std::string const & base) { return base::JoinPath(dir, base + ".SRTMGL1.hgt.zip"); } } // namespace // SrtmTile ---------------------------------------------------------------------------------------- SrtmTile::SrtmTile(SrtmTile && rhs) : m_data(std::move(rhs.m_data)), m_valid(rhs.m_valid) { rhs.Invalidate(); } void SrtmTile::Init(std::string const & dir, ms::LatLon const & coord) { Invalidate(); std::string const base = GetBase(coord); std::string const cont = GetSrtmContFileName(dir, base); std::string file = base + ".hgt"; // Original files are stored in zip archives. Alternatively, they can be loaded // from uncompressed files like "N34E012.hgt". static bool const loadFromZip = static_cast(std::ifstream(cont)); if (loadFromZip) { UnzipMemDelegate delegate(m_data); try { ZipFileReader::UnzipFile(cont, file, delegate); } catch (ZipFileReader::LocateZipException const &) { // Sometimes packed file has different name. See N39E051 measure. file = base + ".SRTMGL1.hgt"; ZipFileReader::UnzipFile(cont, file, delegate); } if (!delegate.m_completed) { LOG(LWARNING, ("Can't decompress SRTM file:", cont)); Invalidate(); return; } } else { m_data = base::ReadFile(base::JoinPath(dir, file)); } if (m_data.size() != kSrtmTileSize) { LOG(LWARNING, ("Bad decompressed SRTM file size:", cont, m_data.size())); Invalidate(); return; } m_valid = true; } // static ms::LatLon SrtmTile::GetCoordInSeconds(ms::LatLon const & coord) { double ln = coord.m_lon - static_cast(coord.m_lon); if (ln < 0) ln += 1; double lt = coord.m_lat - static_cast(coord.m_lat); if (lt < 0) lt += 1; lt = 1 - lt; // from North to South return {kArcSecondsInDegree * lt, kArcSecondsInDegree * ln}; } geometry::Altitude SrtmTile::GetHeight(ms::LatLon const & coord) const { if (!IsValid()) return geometry::kInvalidAltitude; auto const ll = GetCoordInSeconds(coord); return GetHeightRC(std::lround(ll.m_lat), std::lround(ll.m_lon)); } geometry::Altitude SrtmTile::GetHeightRC(size_t row, size_t col) const { size_t const ix = row * (kArcSecondsInDegree + 1) + col; CHECK_LESS(ix, Size(), (row, col)); return ReverseByteOrder(Data()[ix]); } double SrtmTile::GetTriangleHeight(ms::LatLon const & coord) const { if (!IsValid()) return geometry::kInvalidAltitude; auto const ll = GetCoordInSeconds(coord); m2::Point const p1(std::lround(ll.m_lon), std::lround(ll.m_lat)); auto p2 = p1; if (p2.x > ll.m_lon) { if (p2.x > 0) --p2.x; } else if (p2.x < ll.m_lon) { if (p2.x < kArcSecondsInDegree) ++p2.x; } auto p3 = p1; if (p3.y > ll.m_lat) { if (p3.y > 0) --p3.y; } else if (p3.y < ll.m_lat) { if (p3.y < kArcSecondsInDegree) ++p3.y; } // Approximate height from triangle p1, p2, p3. // p1.y == p2.y; p1.x == p3.x // https://stackoverflow.com/questions/36090269/finding-height-of-point-on-height-map-triangles int const det = (p2.y - p3.y) * (p1.x - p3.x) + (p3.x - p2.x) * (p1.y - p3.y); if (det == 0) return GetHeightRC(p1.y, p1.x); double const a1 = ((p2.y - p3.y) * (ll.m_lon - p3.x) + (p3.x - p2.x) * (ll.m_lat - p3.y)) / det; double const a2 = ((p3.y - p1.y) * (ll.m_lon - p3.x) + (p1.x - p3.x) * (ll.m_lat - p3.y)) / det; double const a3 = 1 - a1 - a2; return a1 * GetHeightRC(p1.y, p1.x) + a2 * GetHeightRC(p2.y, p2.x) + a3 * GetHeightRC(p3.y, p3.x); } double SrtmTile::GetBilinearHeight(ms::LatLon const & coord) const { if (!IsValid()) return geometry::kInvalidAltitude; auto const ll = GetCoordInSeconds(coord); m2::Point const p1(static_cast(ll.m_lon), static_cast(ll.m_lat)); auto p2 = p1; if (p2.x < kArcSecondsInDegree) ++p2.x; if (p2.y < kArcSecondsInDegree) ++p2.y; // https://en.wikipedia.org/wiki/Bilinear_interpolation double const denom = (p2.x - p1.x) * (p2.y - p1.y); if (denom == 0) return GetHeightRC(p1.y, p1.x); return (GetHeightRC(p1.y, p1.x) * (p2.x - ll.m_lon) * (p2.y - ll.m_lat) + GetHeightRC(p1.y, p2.x) * (ll.m_lon - p1.x) * (p2.y - ll.m_lat) + GetHeightRC(p2.y, p1.x) * (p2.x - ll.m_lon) * (ll.m_lat - p1.y) + GetHeightRC(p2.y, p2.x) * (ll.m_lon - p1.x) * (ll.m_lat - p1.y)) / denom; } // static std::string SrtmTile::GetPath(std::string const & dir, std::string const & base) { return GetSrtmContFileName(dir, base); } // static SrtmTile::LatLonKey SrtmTile::GetKey(ms::LatLon const & coord) { ms::LatLon center{floor(coord.m_lat) + 0.5, floor(coord.m_lon) + 0.5}; if (coord.m_lat < 0) center.m_lat -= 1.0; if (coord.m_lon < 0) center.m_lon -= 1.0; return {static_cast(center.m_lat), static_cast(center.m_lon)}; } // static std::string SrtmTile::GetBase(ms::LatLon const & coord) { auto key = GetKey(coord); std::ostringstream ss; if (coord.m_lat < 0) { ss << "S"; key.first = -key.first; } else ss << "N"; ss << std::setw(2) << std::setfill('0') << key.first; if (coord.m_lon < 0) { ss << "W"; key.second = -key.second; } else ss << "E"; ss << std::setw(3) << key.second; return ss.str(); } geometry::Altitude * SrtmTile::DataForTests(size_t & sz) { m_valid = true; sz = kArcSecondsInDegree + 1; m_data.resize(kSrtmTileSize, 0); return reinterpret_cast(m_data.data()); } void SrtmTile::Invalidate() { m_data.clear(); m_data.shrink_to_fit(); m_valid = false; } // SrtmTileManager --------------------------------------------------------------------------------- SrtmTile const & SrtmTileManager::GetTile(ms::LatLon const & coord) { auto res = m_tiles.emplace(SrtmTile::GetKey(coord), SrtmTile()); if (res.second) { try { res.first->second.Init(m_dir, coord); } catch (RootException const & e) { std::string const base = SrtmTile::GetBase(coord); LOG(LINFO, ("Can't init SRTM tile:", base, "reason:", e.Msg())); } } return res.first->second; } void SrtmTileManager::Purge() { MapT().swap(m_tiles); } } // namespace generator