h3ron_ndarray/transform.rs
1use std::ops::Mul;
2
3use geo_types::{Coord, Rect};
4
5use crate::error::Error;
6
7/// Affine Geotransform
8///
9/// Ported from [affine library](https://github.com/sgillies/affine/blob/master/affine/__init__.py) (used by rasterio).
10///
11/// `a`, `b`, `c`, `d`, `e` and `f` are typed as `f64` and are coefficients of an augmented affine
12/// transformation matrix:
13///
14/// ```text
15/// | x' | | a b c | | x |
16/// | y' | = | d e f | | y |
17/// | 1 | | 0 0 1 | | 1 |
18/// ```
19///
20/// `a`, `b`, and `c` are the elements of the first row of the matrix. `d`, `e`, and `f` are the elements of the second row.
21///
22/// Other sources:
23/// * [GDAL geotransform](https://gdal.org/tutorials/geotransforms_tut.html)
24/// * [rasterio 1.0+ vs. GDAL](https://rasterio.readthedocs.io/en/latest/topics/migrating-to-v1.html#affine-affine-vs-gdal-style-geotransforms)
25///
26#[derive(Clone, PartialEq, Debug)]
27pub struct Transform {
28 a: f64,
29 b: f64,
30 c: f64,
31 d: f64,
32 e: f64,
33 f: f64,
34}
35
36impl Transform {
37 #![allow(clippy::many_single_char_names)]
38 pub const fn new(a: f64, b: f64, c: f64, d: f64, e: f64, f: f64) -> Self {
39 Self { a, b, c, d, e, f }
40 }
41
42 /// create from an f64 slice in the ordering used by rasterio
43 pub const fn from_rasterio(transform: &[f64; 6]) -> Self {
44 Self::new(
45 transform[0],
46 transform[1],
47 transform[2],
48 transform[3],
49 transform[4],
50 transform[5],
51 )
52 }
53
54 /// create from an f64 slice in the ordering used by gdal
55 pub const fn from_gdal(transform: &[f64; 6]) -> Self {
56 Self::new(
57 transform[1],
58 transform[2],
59 transform[0],
60 transform[4],
61 transform[5],
62 transform[3],
63 )
64 }
65
66 /// The determinant of the transform matrix
67 pub fn determinant(&self) -> f64 {
68 self.a * self.e - self.b * self.d
69 }
70
71 /// True if this transform is degenerate.
72 ///
73 /// Which means that it will collapse a shape to an effective area
74 /// of zero. Degenerate transforms cannot be inverted.
75 pub fn is_degenerate(&self) -> bool {
76 self.determinant() == 0.0
77 }
78
79 pub fn invert(&self) -> Result<Self, Error> {
80 if self.is_degenerate() {
81 return Err(Error::TransformNotInvertible);
82 }
83
84 let idet = 1.0 / self.determinant();
85 let ra = self.e * idet;
86 let rb = -self.b * idet;
87 let rd = -self.d * idet;
88 let re = self.a * idet;
89 Ok(Self::new(
90 ra,
91 rb,
92 -self.c * ra - self.f * rb,
93 rd,
94 re,
95 -self.c * rd - self.f * re,
96 ))
97 }
98
99 /// Apply the transformation to a coordinate
100 pub fn transform_coordinate(&self, coordinate: &Coord<f64>) -> Coord<f64> {
101 Coord {
102 x: coordinate.x.mul_add(self.a, coordinate.y * self.b) + self.c,
103 y: coordinate.x.mul_add(self.d, coordinate.y * self.e) + self.f,
104 }
105 }
106}
107
108/// apply the transformation to a coordinate
109impl Mul<&Coord<f64>> for &Transform {
110 type Output = Coord<f64>;
111
112 fn mul(self, rhs: &Coord<f64>) -> Self::Output {
113 self.transform_coordinate(rhs)
114 }
115}
116
117/// apply the transformation to a coordinate
118impl Mul<Coord<f64>> for &Transform {
119 type Output = Coord<f64>;
120
121 fn mul(self, rhs: Coord<f64>) -> Self::Output {
122 self.transform_coordinate(&rhs)
123 }
124}
125
126/// apply the transformation to a rect
127impl Mul<&Rect<f64>> for &Transform {
128 type Output = Rect<f64>;
129
130 fn mul(self, rhs: &Rect<f64>) -> Self::Output {
131 Rect::new(self * rhs.min(), self * rhs.max())
132 }
133}
134
135#[cfg(test)]
136mod tests {
137 /*
138 $ gdalinfo data/r.tiff
139 Driver: GTiff/GeoTIFF
140 Files: data/r.tiff
141 Size is 2000, 2000
142 Coordinate System is:
143 GEOGCRS["WGS 84",
144 DATUM["World Geodetic System 1984",
145 ELLIPSOID["WGS 84",6378137,298.257223563,
146 LENGTHUNIT["metre",1]]],
147 PRIMEM["Greenwich",0,
148 ANGLEUNIT["degree",0.0174532925199433]],
149 CS[ellipsoidal,2],
150 AXIS["geodetic latitude (Lat)",north,
151 ORDER[1],
152 ANGLEUNIT["degree",0.0174532925199433]],
153 AXIS["geodetic longitude (Lon)",east,
154 ORDER[2],
155 ANGLEUNIT["degree",0.0174532925199433]],
156 ID["EPSG",4326]]
157 Data axis to CRS axis mapping: 2,1
158 Origin = (8.113770000000001,49.407919999999997)
159 Pixel Size = (0.001196505000000,-0.001215135000000)
160 Metadata:
161 AREA_OR_POINT=Area
162 Image Structure Metadata:
163 COMPRESSION=LZW
164 INTERLEAVE=BAND
165 Corner Coordinates:
166 Upper Left ( 8.1137700, 49.4079200) ( 8d 6'49.57"E, 49d24'28.51"N)
167 Lower Left ( 8.1137700, 46.9776500) ( 8d 6'49.57"E, 46d58'39.54"N)
168 Upper Right ( 10.5067800, 49.4079200) ( 10d30'24.41"E, 49d24'28.51"N)
169 Lower Right ( 10.5067800, 46.9776500) ( 10d30'24.41"E, 46d58'39.54"N)
170 Center ( 9.3102750, 48.1927850) ( 9d18'36.99"E, 48d11'34.03"N)
171 Band 1 Block=2000x4 Type=Byte, ColorInterp=Gray
172 NoData Value=0
173 */
174
175 use geo_types::Coord;
176
177 use crate::transform::Transform;
178
179 fn r_tiff_test_helper(gt: &Transform) {
180 // upper left pixel
181 let px_ul = Coord { x: 0., y: 0. };
182
183 let coord_ul = gt * px_ul;
184 assert_relative_eq!(coord_ul.x, 8.11377);
185 assert_relative_eq!(coord_ul.y, 49.40792);
186
187 let gt_inv = gt.invert().unwrap();
188 let px_ul_back = >_inv * coord_ul;
189 assert_relative_eq!(px_ul_back.x, 0.0);
190 assert_relative_eq!(px_ul_back.y, 0.0);
191 }
192
193 #[test]
194 fn test_r_tiff_from_gdal() {
195 /*
196 Python 3.8.5 (default, Jul 28 2020, 12:59:40)
197 [GCC 9.3.0] on linux
198 >>> from osgeo import gdal
199 >>> ds = gdal.Open("data/r.tiff")
200 >>> ds.GetGeoTransform()
201 (8.11377, 0.0011965049999999992, 0.0, 49.40792, 0.0, -0.001215135)
202 */
203 let gt = Transform::from_gdal(&[
204 8.11377,
205 0.0011965049999999992,
206 0.0,
207 49.40792,
208 0.0,
209 -0.001215135,
210 ]);
211 r_tiff_test_helper(>);
212 }
213
214 #[test]
215 fn test_r_tiff_from_rasterio() {
216 /*
217 Python 3.8.5 (default, Jul 28 2020, 12:59:40)
218 [GCC 9.3.0] on linux
219 >>> import rasterio
220 >>> ds = rasterio.open("data/r.tiff")
221 >>> ds.transform
222 Affine(0.0011965049999999992, 0.0, 8.11377,
223 0.0, -0.001215135, 49.40792)
224 */
225 let gt = Transform::from_rasterio(&[
226 0.0011965049999999992,
227 0.0,
228 8.11377,
229 0.0,
230 -0.001215135,
231 49.40792,
232 ]);
233 r_tiff_test_helper(>);
234 }
235}