115 lines
3.2 KiB
Rust
115 lines
3.2 KiB
Rust
use std::sync::Arc;
|
|
|
|
use gdal::errors::GdalError;
|
|
use gdal::Dataset;
|
|
|
|
use tracing::{debug, debug_span, error, info};
|
|
|
|
use moka::future::Cache;
|
|
|
|
pub struct MyDataset {
|
|
pub ds: Dataset,
|
|
}
|
|
unsafe impl Send for MyDataset {}
|
|
unsafe impl Sync for MyDataset {}
|
|
|
|
pub struct DatasetRepository {
|
|
cache: Cache<String, Arc<MyDataset>>,
|
|
basedir: String,
|
|
}
|
|
|
|
unsafe impl Send for DatasetRepository {}
|
|
unsafe impl Sync for DatasetRepository {}
|
|
impl DatasetRepository {
|
|
pub fn new(basedir: String) -> Self {
|
|
let c = Cache::builder()
|
|
.max_capacity(100)
|
|
// Create the cache.
|
|
.build();
|
|
|
|
DatasetRepository { cache: c, basedir }
|
|
}
|
|
|
|
async fn get(&self, filename: String) -> Option<Arc<MyDataset>> {
|
|
let full_filename = format!("{}/{filename}", self.basedir);
|
|
|
|
if !self.cache.contains_key(&full_filename) {
|
|
info!("Will open {full_filename} because not in cache!");
|
|
let ds = Dataset::open(full_filename.clone());
|
|
match ds {
|
|
Err(_) => {
|
|
error!("File not present");
|
|
return None;
|
|
}
|
|
Ok(ds) => {
|
|
let mds = Arc::new(MyDataset { ds });
|
|
self.cache.insert(full_filename.clone(), mds).await;
|
|
}
|
|
}
|
|
}
|
|
|
|
self.cache.get(&full_filename).await
|
|
}
|
|
}
|
|
|
|
impl Clone for DatasetRepository {
|
|
fn clone(&self) -> Self {
|
|
Self {
|
|
basedir: self.basedir.clone(),
|
|
cache: self.cache.clone(),
|
|
}
|
|
}
|
|
}
|
|
|
|
pub async fn elevation_from_coordinates(
|
|
dr: &DatasetRepository,
|
|
lat: f64,
|
|
lon: f64,
|
|
) -> Result<Option<f64>, GdalError> {
|
|
let span = debug_span!("req", lat=%lat, lon=%lon);
|
|
let _guard = span.enter();
|
|
|
|
let filename = get_filename_from_latlon(lat, lon);
|
|
debug!(filename, "filename");
|
|
|
|
let ds = &match dr.get(filename).await {
|
|
Some(x) => x,
|
|
None => return Ok(None),
|
|
}
|
|
.ds;
|
|
|
|
let (px, py) = geo_to_pixel(ds, lat, lon)?;
|
|
|
|
let raster_band = ds.rasterband(1)?;
|
|
let raster_value = raster_band.read_as::<f64>((px, py), (1, 1), (1, 1), None)?;
|
|
Ok(Some(raster_value.data[0]))
|
|
}
|
|
|
|
fn get_filename_from_latlon(lat: f64, lon: f64) -> String {
|
|
// Calculate the rounded values for latitude and longitude
|
|
let rounded_lat = lat.floor();
|
|
let rounded_lon = lon.floor();
|
|
|
|
// Convert the rounded values to strings with leading zeros if necessary
|
|
let lat_deg = format!("{:02}", rounded_lat.abs());
|
|
let lon_deg = format!("{:03}", rounded_lon.abs());
|
|
|
|
// Determine the hemisphere prefixes for latitude and longitude
|
|
let lat_prefix = if rounded_lat >= 0.0 { "N" } else { "S" };
|
|
let lon_prefix = if rounded_lon >= 0.0 { "E" } else { "W" };
|
|
|
|
// Construct the filename
|
|
let filename = format!(
|
|
"Copernicus_DSM_30_{}{}_00_{}{}_00_DEM.tif",
|
|
lat_prefix, lat_deg, lon_prefix, lon_deg
|
|
);
|
|
|
|
filename
|
|
}
|
|
|
|
fn geo_to_pixel(dataset: &Dataset, lat: f64, lon: f64) -> Result<(isize, isize), GdalError> {
|
|
let transform = dataset.geo_transform()?;
|
|
let x_pixel = ((lon - transform[0]) / transform[1]).round() as isize;
|
|
let y_pixel = ((lat - transform[3]) / transform[5]).round() as isize;
|
|
Ok((x_pixel, y_pixel))
|
|
}
|