Signed-off-by: Frank Villaro-Dixon <frank@villaro-dixon.eu>
This commit is contained in:
Frank Villaro-Dixon 2024-04-11 14:17:39 +02:00
parent 817fb10cd7
commit 2a6982afb0
2 changed files with 55 additions and 54 deletions

53
src/dem.rs Normal file
View file

@ -0,0 +1,53 @@
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))
}
pub 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]
}

View file

@ -1,60 +1,8 @@
use std::fmt::Error;
use gdal::errors::GdalError;
use gdal::Dataset;
use gdal;
mod dem;
use crate::dem::elevation_from_coordinates;
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.;