64 lines
1.9 KiB
Rust
64 lines
1.9 KiB
Rust
|
use std::fmt::Error;
|
||
|
|
|
||
|
|
use gdal::errors::GdalError;
|
||
|
|
use gdal::Dataset;
|
||
|
|
use gdal;
|
||
|
|
|
||
|
|
|
||
|
|
|
||
|
|
|
||
|
|
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))
|
||
|
|
}
|
||
|
|
|
||
|
|
fn elevation_from_coordinates(lat: f64, lon: f64) -> f64 {
|
||
|
|
let file = get_filename_from_latlon(lat, lon);
|
||
|
|
let full_filename = format!("data/{file}");
|
||
|
|
|
||
|
|
println!("file for {lat} {lon} is {full_filename}");
|
||
|
|
|
||
|
|
let ds = Dataset::open(full_filename).unwrap();
|
||
|
|
println!("This {} is in '{}' and has {} bands.", ds.driver().long_name(), ds.spatial_ref().unwrap().name().unwrap(), ds.raster_count());
|
||
|
|
println!("PRojection: {}", ds.projection());
|
||
|
|
|
||
|
|
|
||
|
|
let (px, py) = geo_to_pixel(&ds, lat, lon).unwrap();
|
||
|
|
|
||
|
|
let raster_band = ds.rasterband(1).unwrap();
|
||
|
|
let raster_value = raster_band.read_as::<f64>((px, py), (1, 1), (1, 1), None).unwrap();
|
||
|
|
raster_value.data[0]
|
||
|
|
}
|
||
|
|
|
||
|
|
fn main() {
|
||
|
|
let lat = 46.6;
|
||
|
|
let lon = 6.;
|
||
|
|
println!("ELE {lat} {lon}: {}", elevation_from_coordinates(lat, lon))
|
||
|
|
|
||
|
|
}
|