raster/geo_transform.h
| Line | Branch | Exec | Source |
|---|---|---|---|
| 1 | /// @file | ||
| 2 | /// @brief The affine pixel-to-world transform of a georeferenced raster. | ||
| 3 | #pragma once | ||
| 4 | |||
| 5 | #include "terraingen/util/geometry.h" | ||
| 6 | |||
| 7 | #include <array> | ||
| 8 | #include <stdexcept> | ||
| 9 | |||
| 10 | namespace terraingen { | ||
| 11 | |||
| 12 | /// Affine mapping from pixel to world coordinates, as GDAL's six-coefficient geotransform. | ||
| 13 | /// | ||
| 14 | /// Pixel (0, 0) is the outer corner of the first pixel, so centres sit at half-integer | ||
| 15 | /// coordinates. World coordinates are in the CRS's units (metres when projected). North-up rasters | ||
| 16 | /// have a negative `pixelHeight` and no rotation. The default maps pixels to world one to one, as | ||
| 17 | /// GDAL does for ungeoreferenced rasters. | ||
| 18 | struct GeoTransform { | ||
| 19 | 888 | double originX = 0.0; ///< World x of pixel coordinate (0, 0). | |
| 20 | 683 | double pixelWidth = 1.0; ///< World x change per column. | |
| 21 | 675 | double rowRotation = 0.0; ///< World x change per row. | |
| 22 | 674 | double originY = 0.0; ///< World y of pixel coordinate (0, 0). | |
| 23 | 670 | double colRotation = 0.0; ///< World y change per column. | |
| 24 | 669 | double pixelHeight = 1.0; ///< World y change per row; negative for north-up. | |
| 25 | |||
| 26 | /// A north-up transform with square pixels. | ||
| 27 | /// @param originX World x of the raster's west edge. | ||
| 28 | /// @param originY World y of the raster's north edge. | ||
| 29 | /// @param pixelSize Pixel edge length in world units; positive. | ||
| 30 | 2691 | static GeoTransform NorthUp(double originX, double originY, double pixelSize) { | |
| 31 | 2691 | return GeoTransform{originX, pixelSize, 0.0, originY, 0.0, -pixelSize}; | |
| 32 | } | ||
| 33 | |||
| 34 | /// A transform from coefficients in the order `GDALDataset::GetGeoTransform` gives them. | ||
| 35 | 272 | static GeoTransform FromGdal(const std::array<double, 6>& c) { | |
| 36 | 272 | return GeoTransform{c[0], c[1], c[2], c[3], c[4], c[5]}; | |
| 37 | } | ||
| 38 | |||
| 39 | /// The coefficients in GDAL's order, as passed to `GDALDataset::SetGeoTransform`. | ||
| 40 | 4 | std::array<double, 6> ToGdal() const { | |
| 41 | 4 | return {originX, pixelWidth, rowRotation, originY, colRotation, pixelHeight}; | |
| 42 | } | ||
| 43 | |||
| 44 | /// Maps pixel coordinates, 0 at the first pixel's outer edge, to world coordinates. | ||
| 45 | 134254117 | Vec2d PixelToWorld(double col, double row) const { | |
| 46 | 134254117 | return {originX + col * pixelWidth + row * rowRotation, | |
| 47 | 134254117 | originY + col * colRotation + row * pixelHeight}; | |
| 48 | } | ||
| 49 | |||
| 50 | /// Maps world coordinates to pixel coordinates (column, row). | ||
| 51 | /// @throws std::domain_error if the transform is not invertible. | ||
| 52 | 1026516 | Vec2d WorldToPixel(double x, double y) const { | |
| 53 | 1026516 | const double det = pixelWidth * pixelHeight - rowRotation * colRotation; | |
| 54 | 1026516 | if (det == 0.0) { | |
| 55 | 2 | throw std::domain_error("GeoTransform is not invertible"); | |
| 56 | } | ||
| 57 | 1026514 | const double dx = x - originX; | |
| 58 | 1026514 | const double dy = y - originY; | |
| 59 | 1026514 | return {(dx * pixelHeight - dy * rowRotation) / det, | |
| 60 | 1026514 | (dy * pixelWidth - dx * colRotation) / det}; | |
| 61 | } | ||
| 62 | |||
| 63 | 888 | bool operator==(const GeoTransform&) const = default; | |
| 64 | }; | ||
| 65 | |||
| 66 | } // namespace terraingen | ||
| 67 |