connectivity_core/
affine.rs

1use std::collections::HashMap;
2use std::convert::From;
3
4/// Affine transform from (row, col) pixel indices to world (x, y) coordinates.
5///
6/// The fields follow the rasterio/`affine` tuple order `(a, b, c, d, e, f)`, i.e.
7/// `tuple(transform)[:6]` in Python, not GDAL's `GetGeoTransform` order:
8/// `x = a * col + b * row + c` and `y = d * col + e * row + f`.
9/// Projected coordinates are in metres; geographic coordinates are longitude (x) and
10/// latitude (y) in degrees.
11#[derive(Debug, Clone, Copy)]
12pub struct Affine {
13    /// Pixel width (`a`).
14    pub x_scale: f64,
15    /// Rotation/skew term (`b`), usually 0.
16    pub x_skew: f64,
17    /// Top-left corner X (`c`).
18    pub x_origin: f64,
19    /// Rotation/skew term (`d`), usually 0.
20    pub y_skew: f64,
21    /// Pixel height (`e`), negative if north-up.
22    pub y_scale: f64,
23    /// Top-left corner Y (`f`).
24    pub y_origin: f64,
25}
26
27impl Affine {
28    /// Convert (row, col) to (x, y) at pixel centre by default.
29    pub fn xy(&self, row: i32, col: i32) -> (f64, f64) {
30        let col_f = (col as f64) + 0.5;
31        let row_f = (row as f64) + 0.5;
32
33        let x = self.x_scale * col_f + self.x_skew * row_f + self.x_origin;
34        let y = self.y_skew * col_f + self.y_scale * row_f + self.y_origin;
35        (x, y)
36    }
37
38    /// Convert world (x, y) to continuous (col, row) pixel coordinates.
39    pub fn col_row(&self, x: f64, y: f64) -> (f64, f64) {
40        let det = self.x_scale * self.y_scale - self.x_skew * self.y_skew;
41        let dx = x - self.x_origin;
42        let dy = y - self.y_origin;
43        let col = (self.y_scale * dx - self.x_skew * dy) / det;
44        let row = (-self.y_skew * dx + self.x_scale * dy) / det;
45        (col, row)
46    }
47}
48
49impl From<(f64, f64, f64, f64, f64, f64)> for Affine {
50    /// From an `(a, b, c, d, e, f)` tuple in rasterio/`affine` order.
51    fn from(t: (f64, f64, f64, f64, f64, f64)) -> Self {
52        Self {
53            x_scale: t.0,
54            x_skew: t.1,
55            x_origin: t.2,
56            y_skew: t.3,
57            y_scale: t.4,
58            y_origin: t.5,
59        }
60    }
61}
62
63/// Level transforms of a tile, keyed by level, anchored like `rastio.read_raster`:
64/// `base * Affine.scale(level) * Affine.translation(col0 // level, row0 // level)`,
65/// with the `affine` package's operation order.
66///
67/// `base` is the full raster's base-level transform and `offsets` the tile's `(row0, col0)`
68/// in base pixels, so level `L` maps the tile's overview grid, whose cell (0, 0) is global
69/// overview cell `(row0 / L, col0 / L)`, and level 1 maps the tile itself.
70pub fn level_transforms(
71    base: impl Into<Affine>,
72    levels: &[usize],
73    offsets: (usize, usize),
74) -> HashMap<i32, Affine> {
75    let Affine { x_scale: a, x_skew: b, x_origin: c, y_skew: d, y_scale: e, y_origin: f } = base.into();
76    levels
77        .iter()
78        .map(|&level| {
79            let l = level as f64;
80            let (sa, sb, sd, se) = (a * l + b * 0.0, a * 0.0 + b * l, d * l + e * 0.0, d * 0.0 + e * l);
81            let (tc, tr) = ((offsets.1 / level) as f64, (offsets.0 / level) as f64);
82            let t = (sa, sb, sa * tc + sb * tr + c, sd, se, sd * tc + se * tr + f);
83            (level as i32, Affine::from(t))
84        })
85        .collect()
86}
connectivity_core/affine.rs

connectivity_core/
affine.rs
