Skip to main content

rusterize/geo/
raster.rs

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