1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
|
//! Windowed reads from tiled EPSG:4326 GeoTIFFs (EFAS, GHS-POP, Copernicus DEM) without GDAL:
//! only the tiles covering the requested bounds are decoded.
use anyhow::{bail, Context, Result};
use std::{fs::File, io::BufReader, path::Path};
use tiff::decoder::{Decoder, DecodingResult, Limits};
use tiff::tags::Tag;
/// A north-up lon/lat raster held in memory: lon = x0 + col * dx, lat = y0 - row * dy (pixel corners).
pub struct Grid {
pub data: Vec<f32>,
pub width: usize,
pub height: usize,
pub x0: f64,
pub y0: f64,
pub dx: f64,
pub dy: f64,
}
impl Grid {
/// Same georeferencing, new values (e.g. slope computed from a DEM)
pub fn with_data(&self, data: Vec<f32>) -> Grid {
Grid { data, width: self.width, height: self.height, x0: self.x0, y0: self.y0, dx: self.dx, dy: self.dy }
}
pub fn get(&self, row: usize, col: usize) -> f32 {
self.data[row * self.width + col]
}
/// Same lookup as `Area.at` in risk_model.py: floor to a pixel, clipped to the edge.
pub fn at(&self, lon: f64, lat: f64) -> f32 {
let col = ((lon - self.x0) / self.dx).floor().clamp(0.0, (self.width - 1) as f64) as usize;
let row = ((self.y0 - lat) / self.dy).floor().clamp(0.0, (self.height - 1) as f64) as usize;
self.get(row, col)
}
pub fn contains(&self, lon: f64, lat: f64) -> bool {
let col = (lon - self.x0) / self.dx;
let row = (self.y0 - lat) / self.dy;
col >= 0.0 && row >= 0.0 && col < self.width as f64 && row < self.height as f64
}
}
/// Read the pixels of `path` covering lon w..e, lat s..n (clipped to the file).
pub fn read_window(path: &Path, w: f64, s: f64, e: f64, n: f64) -> Result<Grid> {
let file = File::open(path).with_context(|| format!("open {}", path.display()))?;
let mut dec = Decoder::new(BufReader::new(file))?.with_limits(Limits::unlimited());
let (width, height) = dec.dimensions()?;
let scale = dec.get_tag_f64_vec(Tag::ModelPixelScaleTag)?;
let tie = dec.get_tag_f64_vec(Tag::ModelTiepointTag)?;
let (dx, dy) = (scale[0], scale[1]);
let (mut x0, mut y0) = (tie[3] - tie[0] * dx, tie[4] + tie[1] * dy);
// GTRasterTypeGeoKey (1025) = 2 means the tiepoint is a pixel centre (Copernicus DEM); GDAL shifts by half a pixel
if let Ok(keys) = dec.get_tag_u32_vec(Tag::GeoKeyDirectoryTag) {
if keys.chunks(4).skip(1).any(|k| k[0] == 1025 && k[3] == 2) {
x0 -= dx / 2.0;
y0 += dy / 2.0;
}
}
let c0 = ((w - x0) / dx).floor().max(0.0) as usize;
let c1 = (((e - x0) / dx).floor() as i64).min(width as i64 - 1);
let r0 = ((y0 - n) / dy).floor().max(0.0) as usize;
let r1 = (((y0 - s) / dy).floor() as i64).min(height as i64 - 1);
if c1 < c0 as i64 || r1 < r0 as i64 {
bail!("{}: bounds outside raster", path.display());
}
let (c1, r1) = (c1 as usize, r1 as usize);
let (ww, wh) = (c1 - c0 + 1, r1 - r0 + 1);
let (tw, th) = dec.chunk_dimensions();
let (tw, th) = (tw as usize, th as usize);
let across = (width as usize).div_ceil(tw);
let mut data = vec![f32::NAN; ww * wh];
for tr in r0 / th..=r1 / th {
for tc in c0 / tw..=c1 / tw {
let idx = (tr * across + tc) as u32;
let (dw, _) = dec.chunk_data_dimensions(idx);
let vals: Vec<f32> = match dec.read_chunk(idx)? {
DecodingResult::F32(v) => v,
DecodingResult::F64(v) => v.into_iter().map(|x| x as f32).collect(),
_ => bail!("{}: unsupported sample type", path.display()),
};
// edge tiles may come back cropped; infer the row stride from the buffer
let stride = if vals.len() == tw * th { tw } else { dw as usize };
let rows = vals.len() / stride;
for r in 0..rows {
let gr = tr * th + r;
if gr < r0 || gr > r1 {
continue;
}
for c in 0..stride {
let gc = tc * tw + c;
if gc >= c0 && gc <= c1 {
data[(gr - r0) * ww + (gc - c0)] = vals[r * stride + c];
}
}
}
}
}
Ok(Grid { data, width: ww, height: wh, x0: x0 + c0 as f64 * dx, y0: y0 - r0 as f64 * dy, dx, dy })
}
/// scipy.ndimage.minimum_filter(a, size) with the default 'reflect' edges; separable, so rows then columns.
pub fn minimum_filter(g: &Grid, size: usize) -> Vec<f32> {
let h = (size / 2) as i64;
let reflect = |i: i64, n: i64| if i < 0 { -i - 1 } else if i >= n { 2 * n - i - 1 } else { i };
let (w, ht) = (g.width as i64, g.height as i64);
let fmin = |a: f32, b: f32| if a.is_nan() || b < a { b } else { a };
let mut tmp = vec![0f32; g.data.len()];
for r in 0..ht {
for c in 0..w {
let mut m = f32::NAN;
for k in -h..=h {
m = fmin(m, g.data[(r * w + reflect(c + k, w)) as usize]);
}
tmp[(r * w + c) as usize] = m;
}
}
let mut out = vec![0f32; g.data.len()];
for r in 0..ht {
for c in 0..w {
let mut m = f32::NAN;
for k in -h..=h {
m = fmin(m, tmp[(reflect(r + k, ht) * w + c) as usize]);
}
out[(r * w + c) as usize] = m;
}
}
out
}
|
