raster/geotiff/geotiff.cpp
| Line | Branch | Exec | Source |
|---|---|---|---|
| 1 | #include "terraingen/raster/geotiff/geotiff.h" | ||
| 2 | |||
| 3 | #include "terraingen/raster/geotiff/tiff_format.h" | ||
| 4 | |||
| 5 | #include "terraingen/util/checks.h" | ||
| 6 | |||
| 7 | #include <cmath> | ||
| 8 | #include <cstdlib> | ||
| 9 | #include <initializer_list> | ||
| 10 | |||
| 11 | namespace terraingen { | ||
| 12 | |||
| 13 | namespace { | ||
| 14 | |||
| 15 | /// One key of a GeoKeyDirectory: its value inline (`location` 0) or at index `value` of the tag | ||
| 16 | /// `location`. | ||
| 17 | struct GeoKeyEntry { | ||
| 18 | std::uint16_t key; | ||
| 19 | std::uint16_t location; | ||
| 20 | std::uint16_t value; | ||
| 21 | }; | ||
| 22 | |||
| 23 | /// A GeoKeyDirectory of version 1.1.0 holding `entries`, each of count 1. | ||
| 24 | 20 | std::vector<std::uint16_t> GeoKeyDirectory(std::initializer_list<GeoKeyEntry> entries) { | |
| 25 | 40 | std::vector<std::uint16_t> directory = {1, 1, 0, static_cast<std::uint16_t>(entries.size())}; | |
| 26 | 180 | for (const GeoKeyEntry& entry : entries) { | |
| 27 | 160 | directory.insert(directory.end(), {entry.key, entry.location, 1, entry.value}); | |
| 28 | } | ||
| 29 | 20 | return directory; | |
| 30 | } | ||
| 31 | |||
| 32 | /// The value of attribute `name` in the XML start tag `tag`, if present. | ||
| 33 | 118 | std::optional<std::string> XmlAttribute(const std::string& tag, const std::string& name) { | |
| 34 | 118 | const std::string key = " " + name + "=\""; | |
| 35 | 118 | const std::size_t start = tag.find(key); | |
| 36 | 118 | if (start == std::string::npos) { | |
| 37 | 3 | return std::nullopt; | |
| 38 | } | ||
| 39 | 115 | const std::size_t first = start + key.size(); | |
| 40 | 115 | const std::size_t last = tag.find('"', first); | |
| 41 | 115 | if (last == std::string::npos) { | |
| 42 | ✗ | return std::nullopt; | |
| 43 | } | ||
| 44 | 115 | return tag.substr(first, last - first); | |
| 45 | 118 | } | |
| 46 | |||
| 47 | } // namespace | ||
| 48 | |||
| 49 | 628 | std::optional<std::uint16_t> GeoKeys::ShortKey(std::uint16_t key) const { | |
| 50 | 628 | if (directory.size() < 4) { | |
| 51 | 253 | return std::nullopt; | |
| 52 | } | ||
| 53 | 375 | const std::size_t keys = directory[3]; | |
| 54 | 693 | for (std::size_t k = 0; k < keys && 4 + 4 * k + 3 < directory.size(); ++k) { | |
| 55 | 657 | const std::size_t at = 4 + 4 * k; | |
| 56 | 657 | if (directory[at] == key && directory[at + 1] == 0) { | |
| 57 | 339 | return directory[at + 3]; | |
| 58 | } | ||
| 59 | } | ||
| 60 | 36 | return std::nullopt; | |
| 61 | } | ||
| 62 | |||
| 63 | 537 | bool GeoKeys::PixelIsPoint() const { | |
| 64 | 537 | return ShortKey(detail::kGtRasterTypeGeoKey) == detail::kRasterPixelIsPoint; | |
| 65 | } | ||
| 66 | |||
| 67 | 88 | bool GeoKeys::IsGeographic() const { | |
| 68 | 88 | return ShortKey(detail::kGtModelTypeGeoKey) == detail::kModelTypeGeographic; | |
| 69 | } | ||
| 70 | |||
| 71 | 233 | std::vector<BandScaling> ParseGdalBandScaling(const std::string& xml, std::size_t bands) { | |
| 72 | 233 | std::vector<BandScaling> scaling(bands); | |
| 73 | 233 | std::size_t position = 0; | |
| 74 | 292 | while ((position = xml.find("<Item", position)) != std::string::npos) { | |
| 75 | 60 | const std::size_t close = xml.find('>', position); | |
| 76 | 60 | const std::size_t end = xml.find("</Item>", close); | |
| 77 | 60 | if (close == std::string::npos || end == std::string::npos) { | |
| 78 | break; | ||
| 79 | } | ||
| 80 | 59 | const std::string tag = xml.substr(position, close - position); | |
| 81 | 59 | const std::string text = xml.substr(close + 1, end - close - 1); | |
| 82 | 59 | position = end; | |
| 83 | 118 | const auto sample = XmlAttribute(tag, "sample"); | |
| 84 | 59 | const auto role = XmlAttribute(tag, "role"); | |
| 85 | 59 | if (!sample || !role || (*role != "scale" && *role != "offset")) { | |
| 86 | 3 | continue; | |
| 87 | } | ||
| 88 | 56 | char* parsed = nullptr; | |
| 89 | 56 | const unsigned long long band = std::strtoull(sample->c_str(), &parsed, 10); | |
| 90 | 56 | if (parsed == sample->c_str() || band >= bands) { | |
| 91 | 2 | continue; | |
| 92 | } | ||
| 93 | 54 | const double value = std::strtod(text.c_str(), &parsed); | |
| 94 | 54 | if (parsed == text.c_str() || !std::isfinite(value)) { | |
| 95 | 1 | continue; | |
| 96 | } | ||
| 97 | 53 | (*role == "scale" ? scaling[band].scale : scaling[band].offset) = value; | |
| 98 | 77 | } | |
| 99 | 233 | return scaling; | |
| 100 | } | ||
| 101 | |||
| 102 | 21 | GeoKeys GeographicSphereKeys(double radius) { | |
| 103 | using namespace detail; | ||
| 104 | 21 | RequirePositive(radius, "A sphere's radius"); | |
| 105 | 20 | GeoKeys keys; | |
| 106 | 40 | keys.directory = GeoKeyDirectory({ | |
| 107 | {kGtModelTypeGeoKey, 0, kModelTypeGeographic}, | ||
| 108 | {kGtRasterTypeGeoKey, 0, kRasterPixelIsArea}, | ||
| 109 | {kGeographicTypeGeoKey, 0, kUserDefined}, | ||
| 110 | {kGeogGeodeticDatumGeoKey, 0, kUserDefined}, | ||
| 111 | {kGeogAngularUnitsGeoKey, 0, kAngularDegree}, | ||
| 112 | {kGeogEllipsoidGeoKey, 0, kUserDefined}, | ||
| 113 | {kGeogSemiMajorAxisGeoKey, kGeoDoubleParams, 0}, | ||
| 114 | {kGeogSemiMinorAxisGeoKey, kGeoDoubleParams, 1}, | ||
| 115 | 20 | }); | |
| 116 | 20 | keys.doubleParams = {radius, radius}; | |
| 117 | 20 | return keys; | |
| 118 | } | ||
| 119 | |||
| 120 | } // namespace terraingen | ||
| 121 |