geotiff_core/
transform.rs

1//! Geo-transform: pixel coordinates to/from geographic coordinates.
2
3use crate::crs::RasterType;
4
5/// An affine geo-transform mapping pixel (col, row) to map (x, y).
6///
7/// Follows the GDAL convention:
8/// ```text
9/// x = origin_x + col * pixel_width + row * skew_x
10/// y = origin_y + col * skew_y     + row * pixel_height
11/// ```
12///
13/// For north-up images, `skew_x` and `skew_y` are 0 and `pixel_height` is negative.
14#[derive(Debug, Clone, Copy)]
15pub struct GeoTransform {
16    pub origin_x: f64,
17    pub pixel_width: f64,
18    pub skew_x: f64,
19    pub origin_y: f64,
20    pub skew_y: f64,
21    pub pixel_height: f64,
22}
23
24impl GeoTransform {
25    /// Build from ModelTiepoint (tag 33922) and ModelPixelScale (tag 33550).
26    pub fn from_tiepoint_and_scale(tiepoint: &[f64; 6], pixel_scale: &[f64; 3]) -> Self {
27        Self::from_tiepoint_and_scale_with_raster_type(
28            tiepoint,
29            pixel_scale,
30            RasterType::PixelIsArea,
31        )
32    }
33
34    /// Build from ModelTiepoint and ModelPixelScale using the GeoTIFF raster type.
35    ///
36    /// The returned transform is normalized to a corner-based affine transform so
37    /// bounds and pixel-space math stay consistent for both PixelIsArea and
38    /// PixelIsPoint rasters.
39    pub fn from_tiepoint_and_scale_with_raster_type(
40        tiepoint: &[f64; 6],
41        pixel_scale: &[f64; 3],
42        raster_type: RasterType,
43    ) -> Self {
44        // tiepoint: [I, J, K, X, Y, Z]
45        // pixel_scale: [ScaleX, ScaleY, ScaleZ]
46        let pixel_offset = match raster_type {
47            RasterType::PixelIsPoint => 0.5,
48            RasterType::PixelIsArea | RasterType::Unknown(_) => 0.0,
49        };
50        Self {
51            origin_x: tiepoint[3] - (tiepoint[0] + pixel_offset) * pixel_scale[0],
52            pixel_width: pixel_scale[0],
53            skew_x: 0.0,
54            origin_y: tiepoint[4] + (tiepoint[1] + pixel_offset) * pixel_scale[1],
55            skew_y: 0.0,
56            pixel_height: -pixel_scale[1],
57        }
58    }
59
60    /// Build from a 4x4 ModelTransformation matrix (tag 34264), row-major.
61    pub fn from_transformation_matrix(matrix: &[f64; 16]) -> Self {
62        Self {
63            origin_x: matrix[3],
64            pixel_width: matrix[0],
65            skew_x: matrix[1],
66            origin_y: matrix[7],
67            skew_y: matrix[4],
68            pixel_height: matrix[5],
69        }
70    }
71
72    /// Create from origin + pixel size (north-up, no skew).
73    pub fn from_origin_and_pixel_size(
74        origin_x: f64,
75        origin_y: f64,
76        pixel_width: f64,
77        pixel_height: f64,
78    ) -> Self {
79        Self {
80            origin_x,
81            pixel_width,
82            skew_x: 0.0,
83            origin_y,
84            skew_y: 0.0,
85            pixel_height,
86        }
87    }
88
89    /// Convert pixel coordinates (col, row) to map coordinates (x, y).
90    pub fn pixel_to_geo(&self, col: f64, row: f64) -> (f64, f64) {
91        let x = self.origin_x + col * self.pixel_width + row * self.skew_x;
92        let y = self.origin_y + col * self.skew_y + row * self.pixel_height;
93        (x, y)
94    }
95
96    /// Convert map coordinates (x, y) to pixel coordinates (col, row).
97    ///
98    /// Returns `None` if the transform is degenerate (zero determinant).
99    pub fn geo_to_pixel(&self, x: f64, y: f64) -> Option<(f64, f64)> {
100        let det = self.pixel_width * self.pixel_height - self.skew_x * self.skew_y;
101        if det.abs() < 1e-15 {
102            return None;
103        }
104        let dx = x - self.origin_x;
105        let dy = y - self.origin_y;
106        let col = (self.pixel_height * dx - self.skew_x * dy) / det;
107        let row = (-self.skew_y * dx + self.pixel_width * dy) / det;
108        Some((col, row))
109    }
110
111    /// Returns the geographic bounds (min_x, min_y, max_x, max_y) for an image
112    /// of the given width and height.
113    pub fn bounds(&self, width: u32, height: u32) -> [f64; 4] {
114        let corners = [
115            self.pixel_to_geo(0.0, 0.0),
116            self.pixel_to_geo(width as f64, 0.0),
117            self.pixel_to_geo(0.0, height as f64),
118            self.pixel_to_geo(width as f64, height as f64),
119        ];
120        let min_x = corners.iter().map(|c| c.0).fold(f64::INFINITY, f64::min);
121        let max_x = corners
122            .iter()
123            .map(|c| c.0)
124            .fold(f64::NEG_INFINITY, f64::max);
125        let min_y = corners.iter().map(|c| c.1).fold(f64::INFINITY, f64::min);
126        let max_y = corners
127            .iter()
128            .map(|c| c.1)
129            .fold(f64::NEG_INFINITY, f64::max);
130        [min_x, min_y, max_x, max_y]
131    }
132
133    /// Serialize to a tiepoint + pixel_scale pair (for north-up, no-skew images).
134    /// Returns `None` if there is skew (use `to_transformation_matrix` instead).
135    ///
136    /// This uses PixelIsArea semantics. Use
137    /// [`Self::to_tiepoint_and_scale_with_raster_type`] when writing a GeoTIFF
138    /// with an explicit raster type.
139    pub fn to_tiepoint_and_scale(&self) -> Option<([f64; 6], [f64; 3])> {
140        self.to_tiepoint_and_scale_with_raster_type(RasterType::PixelIsArea)
141    }
142
143    /// Serialize to a tiepoint + pixel_scale pair using the GeoTIFF raster type.
144    ///
145    /// The transform is stored internally as a corner-based affine transform.
146    /// PixelIsPoint tiepoints, however, refer to pixel centers, so the emitted
147    /// tiepoint is shifted by half a pixel to roundtrip through
148    /// [`Self::from_tiepoint_and_scale_with_raster_type`] without changing the
149    /// normalized transform.
150    pub fn to_tiepoint_and_scale_with_raster_type(
151        &self,
152        raster_type: RasterType,
153    ) -> Option<([f64; 6], [f64; 3])> {
154        if self.skew_x.abs() > 1e-15 || self.skew_y.abs() > 1e-15 {
155            return None;
156        }
157        let scale = [self.pixel_width, -self.pixel_height, 0.0];
158        let pixel_offset = match raster_type {
159            RasterType::PixelIsPoint => 0.5,
160            RasterType::PixelIsArea | RasterType::Unknown(_) => 0.0,
161        };
162        let tiepoint = [
163            0.0,
164            0.0,
165            0.0,
166            self.origin_x + pixel_offset * scale[0],
167            self.origin_y - pixel_offset * scale[1],
168            0.0,
169        ];
170        Some((tiepoint, scale))
171    }
172
173    /// Serialize to a 4x4 transformation matrix (row-major).
174    pub fn to_transformation_matrix(&self) -> [f64; 16] {
175        [
176            self.pixel_width,
177            self.skew_x,
178            0.0,
179            self.origin_x,
180            self.skew_y,
181            self.pixel_height,
182            0.0,
183            self.origin_y,
184            0.0,
185            0.0,
186            0.0,
187            0.0,
188            0.0,
189            0.0,
190            0.0,
191            1.0,
192        ]
193    }
194}
195
196#[cfg(test)]
197mod tests {
198    use super::*;
199    use crate::crs::RasterType;
200
201    #[test]
202    fn tiepoint_and_scale_roundtrip() {
203        let tp = [0.0, 0.0, 0.0, -180.0, 90.0, 0.0];
204        let scale = [0.1, 0.1, 0.0];
205        let gt = GeoTransform::from_tiepoint_and_scale(&tp, &scale);
206
207        let (x, y) = gt.pixel_to_geo(0.0, 0.0);
208        assert!((x - (-180.0)).abs() < 1e-10);
209        assert!((y - 90.0).abs() < 1e-10);
210
211        let (x2, y2) = gt.pixel_to_geo(10.0, 10.0);
212        assert!((x2 - (-179.0)).abs() < 1e-10);
213        assert!((y2 - 89.0).abs() < 1e-10);
214
215        let (col, row) = gt.geo_to_pixel(x2, y2).unwrap();
216        assert!((col - 10.0).abs() < 1e-10);
217        assert!((row - 10.0).abs() < 1e-10);
218    }
219
220    #[test]
221    fn bounds_calculation() {
222        let tp = [0.0, 0.0, 0.0, 0.0, 10.0, 0.0];
223        let scale = [1.0, 1.0, 0.0];
224        let gt = GeoTransform::from_tiepoint_and_scale(&tp, &scale);
225        let bounds = gt.bounds(10, 10);
226        assert!((bounds[0] - 0.0).abs() < 1e-10);
227        assert!((bounds[1] - 0.0).abs() < 1e-10);
228        assert!((bounds[2] - 10.0).abs() < 1e-10);
229        assert!((bounds[3] - 10.0).abs() < 1e-10);
230    }
231
232    #[test]
233    fn pixel_is_point_tiepoint_is_normalized_to_outer_bounds() {
234        let tp = [0.0, 0.0, 0.0, 100.0, 200.0, 0.0];
235        let scale = [2.0, 2.0, 0.0];
236        let gt = GeoTransform::from_tiepoint_and_scale_with_raster_type(
237            &tp,
238            &scale,
239            RasterType::PixelIsPoint,
240        );
241
242        let (min_x, max_y) = gt.pixel_to_geo(0.0, 0.0);
243        assert!((min_x - 99.0).abs() < 1e-10);
244        assert!((max_y - 201.0).abs() < 1e-10);
245
246        let (center_x, center_y) = gt.pixel_to_geo(0.5, 0.5);
247        assert!((center_x - 100.0).abs() < 1e-10);
248        assert!((center_y - 200.0).abs() < 1e-10);
249    }
250
251    #[test]
252    fn to_tiepoint_and_scale_roundtrips() {
253        let gt = GeoTransform::from_origin_and_pixel_size(-180.0, 90.0, 0.1, -0.1);
254        let (tp, scale) = gt.to_tiepoint_and_scale().unwrap();
255        let gt2 = GeoTransform::from_tiepoint_and_scale(&tp, &scale);
256        assert!((gt2.origin_x - gt.origin_x).abs() < 1e-10);
257        assert!((gt2.origin_y - gt.origin_y).abs() < 1e-10);
258        assert!((gt2.pixel_width - gt.pixel_width).abs() < 1e-10);
259        assert!((gt2.pixel_height - gt.pixel_height).abs() < 1e-10);
260    }
261
262    #[test]
263    fn pixel_is_point_tiepoint_and_scale_roundtrips_normalized_transform() {
264        let gt = GeoTransform::from_origin_and_pixel_size(99.0, 201.0, 2.0, -2.0);
265        let (tp, scale) = gt
266            .to_tiepoint_and_scale_with_raster_type(RasterType::PixelIsPoint)
267            .unwrap();
268        assert!((tp[3] - 100.0).abs() < 1e-10);
269        assert!((tp[4] - 200.0).abs() < 1e-10);
270
271        let gt2 = GeoTransform::from_tiepoint_and_scale_with_raster_type(
272            &tp,
273            &scale,
274            RasterType::PixelIsPoint,
275        );
276        assert!((gt2.origin_x - gt.origin_x).abs() < 1e-10);
277        assert!((gt2.origin_y - gt.origin_y).abs() < 1e-10);
278        assert!((gt2.pixel_width - gt.pixel_width).abs() < 1e-10);
279        assert!((gt2.pixel_height - gt.pixel_height).abs() < 1e-10);
280    }
281
282    #[test]
283    fn skewed_transform_returns_none_for_tiepoint_scale() {
284        let gt = GeoTransform {
285            origin_x: 0.0,
286            pixel_width: 1.0,
287            skew_x: 0.5,
288            origin_y: 0.0,
289            skew_y: 0.0,
290            pixel_height: -1.0,
291        };
292        assert!(gt.to_tiepoint_and_scale().is_none());
293    }
294
295    #[test]
296    fn transformation_matrix_roundtrips() {
297        let gt = GeoTransform::from_origin_and_pixel_size(100.0, 200.0, 0.5, -0.5);
298        let matrix = gt.to_transformation_matrix();
299        let gt2 = GeoTransform::from_transformation_matrix(&matrix);
300        assert!((gt2.origin_x - gt.origin_x).abs() < 1e-10);
301        assert!((gt2.origin_y - gt.origin_y).abs() < 1e-10);
302        assert!((gt2.pixel_width - gt.pixel_width).abs() < 1e-10);
303        assert!((gt2.pixel_height - gt.pixel_height).abs() < 1e-10);
304    }
305}
geotiff_core/transform.rs

geotiff_core/
transform.rs
