diff --git a/src/dem.rs b/src/dem.rs new file mode 100644 index 0000000..719a304 --- /dev/null +++ b/src/dem.rs @@ -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::((px, py), (1, 1), (1, 1), None).unwrap(); + raster_value.data[0] +} \ No newline at end of file diff --git a/src/main.rs b/src/main.rs index 9a97e70..18a43c0 100644 --- a/src/main.rs +++ b/src/main.rs @@ -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::((px, py), (1, 1), (1, 1), None).unwrap(); - raster_value.data[0] -} - fn main() { let lat = 46.6; let lon = 6.;