raster/geotiff/geotiff_write.cpp
| Line | Branch | Exec | Source |
|---|---|---|---|
| 1 | #include "terraingen/raster/geotiff/geotiff.h" | ||
| 2 | |||
| 3 | #include "terraingen/raster/geotiff/tiff_codecs.h" | ||
| 4 | #include "terraingen/raster/geotiff/tiff_format.h" | ||
| 5 | |||
| 6 | #include "terraingen/util/log.h" | ||
| 7 | |||
| 8 | #include <array> | ||
| 9 | #include <cstdio> | ||
| 10 | #include <cstring> | ||
| 11 | #include <fstream> | ||
| 12 | #include <map> | ||
| 13 | |||
| 14 | namespace terraingen { | ||
| 15 | |||
| 16 | namespace { | ||
| 17 | |||
| 18 | using namespace detail; | ||
| 19 | |||
| 20 | /// Builds an IFD with values stored inline or in a block after it. | ||
| 21 | class IfdWriter { | ||
| 22 | public: | ||
| 23 | 306 | explicit IfdWriter(bool bigTiff) : format_(IfdFormat::For(bigTiff)) {} | |
| 24 | |||
| 25 | template <typename T> | ||
| 26 | 4747 | void Add(std::uint16_t tag, FieldType type, const std::vector<T>& values) { | |
| 27 | 4747 | std::vector<std::uint8_t> bytes(values.size() * sizeof(T)); | |
| 28 | 4747 | std::memcpy(bytes.data(), values.data(), bytes.size()); | |
| 29 | 4747 | entries_[tag] = {type, values.size(), std::move(bytes)}; | |
| 30 | 9494 | } | |
| 31 | |||
| 32 | 36 | void AddAscii(std::uint16_t tag, const std::string& text) { | |
| 33 | 36 | std::vector<std::uint8_t> bytes(text.begin(), text.end()); | |
| 34 | 36 | bytes.push_back(0); | |
| 35 | 36 | entries_[tag] = {kAscii, bytes.size(), std::move(bytes)}; | |
| 36 | 72 | } | |
| 37 | |||
| 38 | /// The IFD and its out-of-line values, to be written at `offset`. | ||
| 39 | 306 | std::vector<std::uint8_t> Serialize(std::uint64_t offset) const { | |
| 40 | const std::size_t tableSize = | ||
| 41 | 306 | format_.countSize + entries_.size() * format_.EntrySize() + format_.wordSize; | |
| 42 | 306 | std::vector<std::uint8_t> table(tableSize, 0); | |
| 43 | 306 | std::vector<std::uint8_t> extra; | |
| 44 | 15763 | const auto put = [&table](std::size_t at, std::uint64_t value, std::size_t size) { | |
| 45 | 15763 | std::memcpy(table.data() + at, &value, size); | |
| 46 | 16069 | }; | |
| 47 | 306 | put(0, entries_.size(), format_.countSize); | |
| 48 | 306 | std::size_t at = format_.countSize; | |
| 49 | 5089 | for (const auto& [tag, entry] : entries_) { | |
| 50 | 4783 | put(at, tag, 2); | |
| 51 | 4783 | put(at + 2, entry.type, 2); | |
| 52 | 4783 | put(at + 4, entry.count, format_.wordSize); | |
| 53 | 4783 | const std::size_t valueAt = at + format_.ValueAt(); | |
| 54 | 4783 | if (entry.bytes.size() <= format_.wordSize) { | |
| 55 | 3675 | std::memcpy(table.data() + valueAt, entry.bytes.data(), entry.bytes.size()); | |
| 56 | } else { | ||
| 57 | 1108 | if (extra.size() % 2 != 0) { | |
| 58 | 3 | extra.push_back(0); | |
| 59 | } | ||
| 60 | 1108 | put(valueAt, offset + tableSize + extra.size(), format_.wordSize); | |
| 61 | 1108 | extra.insert(extra.end(), entry.bytes.begin(), entry.bytes.end()); | |
| 62 | } | ||
| 63 | 4783 | at += format_.EntrySize(); | |
| 64 | } | ||
| 65 | 306 | table.insert(table.end(), extra.begin(), extra.end()); | |
| 66 | 612 | return table; | |
| 67 | 306 | } | |
| 68 | |||
| 69 | private: | ||
| 70 | IfdFormat format_; | ||
| 71 | std::map<std::uint16_t, Entry> entries_; | ||
| 72 | }; | ||
| 73 | |||
| 74 | /// Where the tiles were written, and where the file ends. | ||
| 75 | struct TileIndex { | ||
| 76 | std::vector<std::uint64_t> offsets; | ||
| 77 | std::vector<std::uint64_t> counts; | ||
| 78 | std::uint64_t end = 0; | ||
| 79 | }; | ||
| 80 | |||
| 81 | 314 | void CheckBands(const std::filesystem::path& path, const std::vector<const Raster<float>*>& bands, | |
| 82 | const GeoTiffWriteOptions& options) { | ||
| 83 | 314 | if (bands.empty() || bands.front()->Empty()) { | |
| 84 | 2 | throw GeoTiffError("Cannot write an empty raster to " + path.string()); | |
| 85 | } | ||
| 86 | 312 | const Raster<float>& first = *bands.front(); | |
| 87 | 669 | for (const Raster<float>* band : bands) { | |
| 88 | 717 | if (band->Width() != first.Width() || band->Height() != first.Height() || | |
| 89 | 358 | !(band->Transform() == first.Transform())) { | |
| 90 | 2 | throw GeoTiffError("The bands written to " + path.string() + " must share one grid"); | |
| 91 | } | ||
| 92 | } | ||
| 93 | 310 | if (options.compression == TiffCompression::PackBits) { | |
| 94 | 1 | throw GeoTiffError("Writing PackBits is not supported; use None, Lzw or Deflate"); | |
| 95 | } | ||
| 96 | 309 | if (options.tileSize == 0 || options.tileSize % 16 != 0) { | |
| 97 | 2 | throw GeoTiffError("Tile size must be a positive multiple of 16"); | |
| 98 | } | ||
| 99 | 307 | } | |
| 100 | |||
| 101 | /// Whether `tiles` uncompressed float32 tiles of `tile` pixels square could pass 4 GiB. | ||
| 102 | 304 | bool NeedsBigTiff(std::size_t tiles, std::size_t tile) { | |
| 103 | 304 | const std::uint64_t estimate = static_cast<std::uint64_t>(tiles) * tile * tile * 4 + (1u << 20); | |
| 104 | 304 | return estimate >= (std::uint64_t{1} << 32); | |
| 105 | } | ||
| 106 | |||
| 107 | /// The bytes stored for the tile at column `tx` and row `ty` of tiles of `band`, padded with | ||
| 108 | /// `fill`. | ||
| 109 | 1684 | std::vector<std::uint8_t> EncodeTile(const Raster<float>& band, std::size_t tx, std::size_t ty, | |
| 110 | std::size_t tile, float fill, TiffCompression compression) { | ||
| 111 | 1684 | std::vector<float> samples(tile * tile); | |
| 112 | 410948 | for (std::size_t y = 0; y < tile; ++y) { | |
| 113 | 104831408 | for (std::size_t x = 0; x < tile; ++x) { | |
| 114 | 104422144 | const std::size_t col = tx * tile + x; | |
| 115 | 104422144 | const std::size_t row = ty * tile + y; | |
| 116 | 104422144 | samples[y * tile + x] = | |
| 117 | 104422144 | col < band.Width() && row < band.Height() ? band(col, row) : fill; | |
| 118 | } | ||
| 119 | } | ||
| 120 | 1684 | std::vector<std::uint8_t> bytes(samples.size() * 4); | |
| 121 | 1684 | std::memcpy(bytes.data(), samples.data(), bytes.size()); | |
| 122 | 1684 | if (compression == TiffCompression::None) { | |
| 123 | 24 | return bytes; | |
| 124 | } | ||
| 125 | 410540 | for (std::size_t row = 0; row < tile; ++row) { | |
| 126 | 408880 | EncodeFloatPredictorRow(bytes.data() + row * tile * 4, 4, tile, 1); | |
| 127 | } | ||
| 128 | 1660 | return compression == TiffCompression::Lzw ? EncodeLzw(bytes) : EncodeDeflate(bytes); | |
| 129 | 1684 | } | |
| 130 | |||
| 131 | /// Writes every tile of every band from `start`, each at an even offset. | ||
| 132 | 306 | TileIndex WriteTiles(std::ofstream& file, const std::vector<const Raster<float>*>& bands, | |
| 133 | std::size_t tile, TiffCompression compression, std::uint64_t start) { | ||
| 134 | 306 | const Raster<float>& first = *bands.front(); | |
| 135 | 306 | const std::size_t across = (first.Width() + tile - 1) / tile; | |
| 136 | 306 | const std::size_t down = (first.Height() + tile - 1) / tile; | |
| 137 | 306 | const float fill = first.NoData().value_or(0.0f); | |
| 138 | 306 | TileIndex index; | |
| 139 | 306 | index.end = start; | |
| 140 | 657 | for (const Raster<float>* band : bands) { | |
| 141 | 1009 | for (std::size_t ty = 0; ty < down; ++ty) { | |
| 142 | 2342 | for (std::size_t tx = 0; tx < across; ++tx) { | |
| 143 | const std::vector<std::uint8_t> bytes = | ||
| 144 | 1684 | EncodeTile(*band, tx, ty, tile, fill, compression); | |
| 145 | 1684 | index.offsets.push_back(index.end); | |
| 146 | 1684 | index.counts.push_back(bytes.size()); | |
| 147 | 1684 | file.write(reinterpret_cast<const char*>(bytes.data()), | |
| 148 | 1684 | static_cast<std::streamsize>(bytes.size())); | |
| 149 | 1684 | index.end += bytes.size(); | |
| 150 | 1684 | if (index.end % 2 != 0) { | |
| 151 | 702 | file.put(0); | |
| 152 | 702 | ++index.end; | |
| 153 | } | ||
| 154 | 1684 | } | |
| 155 | } | ||
| 156 | } | ||
| 157 | 306 | return index; | |
| 158 | } | ||
| 159 | |||
| 160 | /// Adds the tags describing the image, its samples and its tiles. | ||
| 161 | 306 | void AddImageTags(IfdWriter& ifd, const Raster<float>& raster, std::size_t bandCount, | |
| 162 | const GeoTiffWriteOptions& options, const TileIndex& tiles, bool bigTiff) { | ||
| 163 | 306 | const auto bands = static_cast<std::uint16_t>(bandCount); | |
| 164 | 612 | ifd.Add<std::uint32_t>(kImageWidth, kLong, {static_cast<std::uint32_t>(raster.Width())}); | |
| 165 | 918 | ifd.Add<std::uint32_t>(kImageLength, kLong, {static_cast<std::uint32_t>(raster.Height())}); | |
| 166 | 306 | ifd.Add<std::uint16_t>(kBitsPerSample, kShort, std::vector<std::uint16_t>(bands, 32)); | |
| 167 | 612 | ifd.Add<std::uint16_t>(kCompression, kShort, {CompressionCodeOf(options.compression)}); | |
| 168 | 612 | ifd.Add<std::uint16_t>(kPhotometric, kShort, {1}); | |
| 169 | 612 | ifd.Add<std::uint16_t>(kSamplesPerPixel, kShort, {bands}); | |
| 170 | 612 | ifd.Add<std::uint16_t>(kPlanarConfiguration, kShort, | |
| 171 | {static_cast<std::uint16_t>(bands == 1 ? 1 : 2)}); | ||
| 172 | 306 | if (options.compression != TiffCompression::None) { | |
| 173 | 912 | ifd.Add<std::uint16_t>(kPredictor, kShort, {3}); | |
| 174 | } | ||
| 175 | 612 | ifd.Add<std::uint32_t>(kTileWidth, kLong, {options.tileSize}); | |
| 176 | 612 | ifd.Add<std::uint32_t>(kTileLength, kLong, {options.tileSize}); | |
| 177 | 306 | if (bigTiff) { | |
| 178 | 3 | ifd.Add<std::uint64_t>(kTileOffsets, kLong8, tiles.offsets); | |
| 179 | 3 | ifd.Add<std::uint64_t>(kTileByteCounts, kLong8, tiles.counts); | |
| 180 | } else { | ||
| 181 | 303 | ifd.Add<std::uint32_t>( | |
| 182 | kTileOffsets, kLong, | ||
| 183 | 909 | std::vector<std::uint32_t>(tiles.offsets.begin(), tiles.offsets.end())); | |
| 184 | 303 | ifd.Add<std::uint32_t>( | |
| 185 | kTileByteCounts, kLong, | ||
| 186 | 909 | std::vector<std::uint32_t>(tiles.counts.begin(), tiles.counts.end())); | |
| 187 | } | ||
| 188 | 306 | if (bands > 1) { | |
| 189 | 14 | ifd.Add<std::uint16_t>(kExtraSamples, kShort, std::vector<std::uint16_t>(bands - 1u, 0)); | |
| 190 | } | ||
| 191 | 306 | ifd.Add<std::uint16_t>(kSampleFormat, kShort, std::vector<std::uint16_t>(bands, 3)); | |
| 192 | 306 | } | |
| 193 | |||
| 194 | 28 | std::string NoDataText(float value) { | |
| 195 | 28 | std::array<char, 40> buffer{}; | |
| 196 | 84 | std::snprintf(buffer.data(), buffer.size(), "%.9g", static_cast<double>(value)); | |
| 197 | 56 | return buffer.data(); | |
| 198 | } | ||
| 199 | |||
| 200 | /// Adds the georeferencing, projection keys and no-data value. | ||
| 201 | 306 | void AddGeoTags(IfdWriter& ifd, const Raster<float>& raster, const GeoKeys& geoKeys) { | |
| 202 | 306 | const GeoTransform& transform = raster.Transform(); | |
| 203 | 306 | if (!geoKeys.Empty() || !(transform == GeoTransform{})) { | |
| 204 | const GeoTransform stored = | ||
| 205 | 305 | geoKeys.PixelIsPoint() ? ShiftHalfPixel(transform, 1.0) : transform; | |
| 206 | 305 | if (stored.rowRotation == 0.0 && stored.colRotation == 0.0) { | |
| 207 | 301 | ifd.Add<double>(kModelPixelScale, kDouble, | |
| 208 | 301 | {stored.pixelWidth, -stored.pixelHeight, 0.0}); | |
| 209 | 602 | ifd.Add<double>(kModelTiepoint, kDouble, | |
| 210 | 301 | {0.0, 0.0, 0.0, stored.originX, stored.originY, 0.0}); | |
| 211 | } else { | ||
| 212 | 8 | ifd.Add<double>(kModelTransformation, kDouble, | |
| 213 | 4 | {stored.pixelWidth, stored.rowRotation, 0.0, stored.originX, | |
| 214 | 4 | stored.colRotation, stored.pixelHeight, 0.0, stored.originY, 0.0, 0.0, | |
| 215 | 0.0, 0.0, 0.0, 0.0, 0.0, 1.0}); | ||
| 216 | } | ||
| 217 | } | ||
| 218 | 306 | if (!geoKeys.Empty()) { | |
| 219 | 92 | ifd.Add<std::uint16_t>(kGeoKeyDirectory, kShort, geoKeys.directory); | |
| 220 | 92 | if (!geoKeys.doubleParams.empty()) { | |
| 221 | 59 | ifd.Add<double>(kGeoDoubleParams, kDouble, geoKeys.doubleParams); | |
| 222 | } | ||
| 223 | 92 | if (!geoKeys.asciiParams.empty()) { | |
| 224 | 8 | ifd.AddAscii(kGeoAsciiParams, geoKeys.asciiParams); | |
| 225 | } | ||
| 226 | } | ||
| 227 | 306 | if (raster.NoData()) { | |
| 228 | 28 | ifd.AddAscii(kGdalNoData, NoDataText(*raster.NoData())); | |
| 229 | } | ||
| 230 | 306 | } | |
| 231 | |||
| 232 | } // namespace | ||
| 233 | |||
| 234 | 281 | void WriteGeoTiff(const std::filesystem::path& path, const Raster<float>& raster, | |
| 235 | const GeoKeys& geoKeys, const GeoTiffWriteOptions& options) { | ||
| 236 | 567 | WriteGeoTiff(path, std::vector<const Raster<float>*>{&raster}, geoKeys, options); | |
| 237 | 276 | } | |
| 238 | |||
| 239 | 30 | void WriteGeoTiff(const std::filesystem::path& path, const std::vector<Raster<float>>& bands, | |
| 240 | const GeoKeys& geoKeys, const GeoTiffWriteOptions& options) { | ||
| 241 | 30 | std::vector<const Raster<float>*> pointers; | |
| 242 | 106 | for (const Raster<float>& band : bands) { | |
| 243 | 76 | pointers.push_back(&band); | |
| 244 | } | ||
| 245 | 30 | WriteGeoTiff(path, pointers, geoKeys, options); | |
| 246 | 30 | } | |
| 247 | |||
| 248 | 314 | void WriteGeoTiff(const std::filesystem::path& path, const std::vector<const Raster<float>*>& bands, | |
| 249 | const GeoKeys& geoKeys, const GeoTiffWriteOptions& options) { | ||
| 250 | 314 | const Stopwatch stopwatch; | |
| 251 | 314 | CheckBands(path, bands, options); | |
| 252 | 307 | const Raster<float>& raster = *bands.front(); | |
| 253 | 307 | const std::size_t tile = options.tileSize; | |
| 254 | const std::size_t tiles = | ||
| 255 | 307 | ((raster.Width() + tile - 1) / tile) * ((raster.Height() + tile - 1) / tile) * bands.size(); | |
| 256 | 307 | const bool bigTiff = options.bigTiff || NeedsBigTiff(tiles, tile); | |
| 257 | 307 | const IfdFormat format = IfdFormat::For(bigTiff); | |
| 258 | |||
| 259 | 307 | std::ofstream file(path, std::ios::binary); | |
| 260 | 307 | if (!file) { | |
| 261 | 1 | throw GeoTiffError("Cannot write " + path.string()); | |
| 262 | } | ||
| 263 | 612 | std::vector<std::uint8_t> header = {'I', 'I', static_cast<std::uint8_t>(bigTiff ? 43 : 42), 0}; | |
| 264 | 306 | if (bigTiff) { | |
| 265 | 3 | header.insert(header.end(), {8, 0, 0, 0}); | |
| 266 | } | ||
| 267 | 306 | header.resize(header.size() + format.wordSize, 0); | |
| 268 | 306 | file.write(reinterpret_cast<const char*>(header.data()), | |
| 269 | 306 | static_cast<std::streamsize>(header.size())); | |
| 270 | |||
| 271 | 306 | const TileIndex index = WriteTiles(file, bands, tile, options.compression, header.size()); | |
| 272 | 306 | if (!bigTiff && index.end >= (std::uint64_t{1} << 32) - (1u << 20)) { | |
| 273 | ✗ | throw GeoTiffError("Data too large for classic TIFF; set bigTiff"); | |
| 274 | } | ||
| 275 | 306 | IfdWriter ifd(bigTiff); | |
| 276 | 306 | AddImageTags(ifd, raster, bands.size(), options, index, bigTiff); | |
| 277 | 306 | AddGeoTags(ifd, raster, geoKeys); | |
| 278 | 306 | const std::vector<std::uint8_t> table = ifd.Serialize(index.end); | |
| 279 | 306 | file.write(reinterpret_cast<const char*>(table.data()), | |
| 280 | 306 | static_cast<std::streamsize>(table.size())); | |
| 281 | 306 | file.seekp(static_cast<std::streamoff>(header.size() - format.wordSize)); | |
| 282 | 306 | file.write(reinterpret_cast<const char*>(&index.end), | |
| 283 | 306 | static_cast<std::streamsize>(format.wordSize)); | |
| 284 | 306 | if (!file) { | |
| 285 | ✗ | throw GeoTiffError("Cannot write " + path.string()); | |
| 286 | } | ||
| 287 | 306 | LogInfo("raster", "Wrote ", raster.Width(), "x", raster.Height(), "x", bands.size(), | |
| 288 | 612 | " GeoTIFF ", path.string(), " in ", FormatDuration(stopwatch.Milliseconds())); | |
| 289 | 307 | } | |
| 290 | |||
| 291 | } // namespace terraingen | ||
| 292 |