rusterize/geo/
raster.rs

1use std::cmp::Ordering;
2
3use crate::error::{RusterizeError, RusterizeResult};
4use geo::{BoundingRect, Geometry, Rect, coord};
5use ndarray::Array3;
6use num_traits::Num;
7
8/// Contains the spatial information associated with the burned [`geo::Geometry`].
9#[derive(Clone)]
10pub struct RasterInfo {
11    pub ncols: usize,
12    pub nrows: usize,
13    pub xmin: f64,
14    pub xmax: f64,
15    pub ymin: f64,
16    pub ymax: f64,
17    pub xres: f64,
18    pub yres: f64,
19    pub epsg: Option<u16>,
20}
21
22impl RasterInfo {
23    pub(crate) fn build_raster<N>(&self, bands: usize, background: N) -> Array3<N>
24    where
25        N: Num + Copy,
26    {
27        Array3::from_elem((bands, self.nrows, self.ncols), background)
28    }
29}
30
31/// Builder for a [`RasterInfo`] instance.
32/// If extent is not provided, it can be inferred from the [`geo::Geometry`] when building it.
33/// In this case, a half-pixel buffer is applied to avoid missing points on the border.
34/// The logics dictating the final spatial properties of the rasterized geometries follow those of GDAL.
35#[derive(Default)]
36pub struct RasterInfoBuilder {
37    shape: Option<[usize; 2]>,
38    extent: Option<[f64; 4]>,
39    resolution: Option<[f64; 2]>,
40    tap: bool,
41    epsg: Option<u16>,
42}
43
44impl RasterInfoBuilder {
45    pub fn new() -> Self {
46        RasterInfoBuilder::default()
47    }
48
49    /// Build into a [`RasterInfo`] with user-defined extent.
50    pub fn build(self) -> RusterizeResult<RasterInfo> {
51        match self.extent {
52            Some(extent) => {
53                let is_unspecified_extent = extent.iter().all(|x| matches!(x.total_cmp(&0.0), Ordering::Equal));
54                if is_unspecified_extent {
55                    return Err(RusterizeError::ValueError("Unspecified extent (all zeros)."));
56                }
57                self.finalize(extent, false)
58            }
59            None => Err(RusterizeError::RuntimeError(
60                "Extent must be provided for construction. \
61                Use `build_with()` to infer extent from geometries.",
62            )),
63        }
64    }
65
66    /// Same as `build`, but infer extent from the geometry.
67    pub fn build_with(self, geoms: &[Geometry<f64>]) -> RusterizeResult<RasterInfo> {
68        if self.extent.is_some() {
69            return Err(RusterizeError::RuntimeError(
70                "Extent must be inferred from geometries for construction. \
71                Use `build()` to provide a custom extent.",
72            ));
73        }
74
75        let bounds = geoms.iter().fold(None, |acc, geom| {
76            let bounds = geom.bounding_rect();
77
78            match (acc, bounds) {
79                (None, None) => None,
80                (None, Some(r)) | (Some(r), None) => Some(r),
81                (Some(r1), Some(r2)) => Some(Rect::new(
82                    coord! { x: r1.min().x.min(r2.min().x), y: r1.min().y.min(r2.min().y) },
83                    coord! { x: r1.max().x.max(r2.max().x), y: r1.max().y.max(r2.max().y) },
84                )),
85            }
86        });
87
88        if let Some(b) = bounds {
89            self.finalize([b.min().x, b.min().y, b.max().x, b.max().y], true)
90        } else {
91            Err(RusterizeError::RuntimeError("Cannot infer bounding box from geometry."))
92        }
93    }
94
95    fn finalize(
96        self,
97        [mut xmin, mut ymin, mut xmax, mut ymax]: [f64; 4],
98        inferred: bool,
99    ) -> RusterizeResult<RasterInfo> {
100        if self.shape.is_none() && self.resolution.is_none() {
101            return Err(RusterizeError::ValueError(
102                "Must set at least one of `shape` or `resolution`",
103            ));
104        }
105        if self.shape.is_some() && self.resolution.is_some() {
106            return Err(RusterizeError::ValueError(
107                "Shape and resolution are mutually exclusive; provide only one",
108            ));
109        }
110        let has_shape = self.shape.is_some();
111        let has_res = self.resolution.is_some();
112        let [mut nrows, mut ncols] = self.shape.unwrap_or_default();
113        let [mut xres, mut yres] = self.resolution.unwrap_or_default();
114
115        if has_shape && (nrows == 0 || ncols == 0) {
116            return Err(RusterizeError::ValueError("Shape values must be > 0."));
117        }
118
119        if has_res && (xres <= 0.0 || yres <= 0.0) {
120            return Err(RusterizeError::ValueError("Resolution values must be > 0."));
121        }
122
123        if inferred && !self.tap && has_res {
124            xmin -= xres / 2.0;
125            xmax += xres / 2.0;
126            ymin -= yres / 2.0;
127            ymax += yres / 2.0;
128        }
129
130        if !has_res {
131            xres = (xmax - xmin) / ncols as f64;
132            yres = (ymax - ymin) / nrows as f64;
133        } else if self.tap {
134            xmin = (xmin / xres).floor() * xres;
135            xmax = (xmax / xres).ceil() * xres;
136            ymin = (ymin / yres).floor() * yres;
137            ymax = (ymax / yres).ceil() * yres;
138        }
139
140        if !has_shape {
141            nrows = (0.5 + (ymax - ymin) / yres) as usize;
142            ncols = (0.5 + (xmax - xmin) / xres) as usize;
143        }
144
145        Ok(RasterInfo {
146            ncols,
147            nrows,
148            xmin,
149            xmax,
150            ymin,
151            ymax,
152            xres,
153            yres,
154            epsg: self.epsg,
155        })
156    }
157
158    pub fn shape(mut self, nrows: usize, ncols: usize) -> Self {
159        self.shape = Some([nrows, ncols]);
160        self
161    }
162
163    pub fn extent(mut self, xmin: f64, ymin: f64, xmax: f64, ymax: f64) -> Self {
164        self.extent = Some([xmin, ymin, xmax, ymax]);
165        self
166    }
167
168    pub fn resolution(mut self, xres: f64, yres: f64) -> Self {
169        self.resolution = Some([xres, yres]);
170        self
171    }
172
173    pub fn with_target_aligned_pixels(mut self) -> Self {
174        self.tap = true;
175        self
176    }
177
178    pub fn epsg(mut self, epsg: u16) -> Self {
179        self.epsg = Some(epsg);
180        self
181    }
182}
rusterize/geo/raster.rs

rusterize/geo/
raster.rs
