diff --git a/Cargo.lock b/Cargo.lock index 48192ba..d0dc4a3 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -468,9 +468,9 @@ dependencies = [ [[package]] name = "geo" -version = "0.30.0" +version = "0.31.0" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "4416397671d8997e9a3e7ad99714f4f00a22e9eaa9b966a5985d2194fc9e02e1" +checksum = "2fc1a1678e54befc9b4bcab6cd43b8e7f834ae8ea121118b0fd8c42747675b4a" dependencies = [ "earcutr", "float_next_after", @@ -663,24 +663,24 @@ checksum = "2304e00983f87ffb38b55b444b5e3b60a884b5d30c0fca7d82fe33449bbe55ea" [[package]] name = "i_float" -version = "1.7.0" +version = "1.15.0" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "85df3a416829bb955fdc2416c7b73680c8dcea8d731f2c7aa23e1042fe1b8343" +checksum = "010025c2c532c8d82e42d0b8bb5184afa449fa6f06c709ea9adcb16c49ae405b" dependencies = [ - "serde", + "libm", ] [[package]] name = "i_key_sort" -version = "0.2.0" +version = "0.6.0" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "347c253b4748a1a28baf94c9ce133b6b166f08573157e05afe718812bc599fcd" +checksum = "9190f86706ca38ac8add223b2aed8b1330002b5cdbbce28fb58b10914d38fc27" [[package]] name = "i_overlay" -version = "2.0.5" +version = "4.0.7" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "0542dfef184afdd42174a03dcc0625b6147fb73e1b974b1a08a2a42ac35cee49" +checksum = "413183068e6e0289e18d7d0a1f661b81546e6918d5453a44570b9ab30cbed1b3" dependencies = [ "i_float", "i_key_sort", @@ -691,19 +691,18 @@ dependencies = [ [[package]] name = "i_shape" -version = "1.7.0" +version = "1.14.0" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "0a38f5a42678726718ff924f6d4a0e79b129776aeed298f71de4ceedbd091bce" +checksum = "1ea154b742f7d43dae2897fcd5ead86bc7b5eefcedd305a7ebf9f69d44d61082" dependencies = [ "i_float", - "serde", ] [[package]] name = "i_tree" -version = "0.8.3" +version = "0.16.0" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "155181bc97d770181cf9477da51218a19ee92a8e5be642e796661aee2b601139" +checksum = "35e6d558e6d4c7b82bc51d9c771e7a927862a161a7d87bf2b0541450e0e20915" [[package]] name = "image" @@ -1358,14 +1357,13 @@ dependencies = [ "image", "itertools 0.14.0", "memchr", - "proj4rs", "rayon", "rusqlite", - "serde", "serde_json", "tempfile", "tracing", "tracing-subscriber", + "tvs-lib", ] [[package]] @@ -1439,6 +1437,17 @@ dependencies = [ "tracing-log", ] +[[package]] +name = "tvs-lib" +version = "0.1.0" +dependencies = [ + "color-eyre", + "geo", + "proj4rs", + "rstar 0.12.2", + "serde", +] + [[package]] name = "typenum" version = "1.19.0" diff --git a/Cargo.toml b/Cargo.toml index e951514..553e4cd 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -2,8 +2,15 @@ resolver = "2" members = [ "crates/total-viewsheds", + "crates/tvs-lib", ] +[workspace.dependencies] +color-eyre = "0.6.5" +geo = "0.31.0" +serde = "1.0.219" +serde_json = "1.0.143" + # Canonical lints for whole crate # # Official docs: diff --git a/benchmarks/run.sh b/benchmarks/run.sh index e27a085..3102024 100755 --- a/benchmarks/run.sh +++ b/benchmarks/run.sh @@ -30,10 +30,10 @@ rm $sqlite_db_path || true # Do the calculations time cargo run --features ring_data --release -- \ compute "$PROJECT_ROOT/benchmarks/cardiff.tiff" \ - --rings-per-km 3 \ --process all \ --viewsheds-db-path $sqlite_db_path \ - --thread-count 1 + --thread-count 1 \ + --angle-subdivisions 1 viewshed_file="output/viewsheds/-3.122999906539917-51.48979949951172.json" rm $viewshed_file || true diff --git a/crates/total-viewsheds/Cargo.toml b/crates/total-viewsheds/Cargo.toml index 9a34a34..a3b8c35 100644 --- a/crates/total-viewsheds/Cargo.toml +++ b/crates/total-viewsheds/Cargo.toml @@ -8,20 +8,19 @@ rust-version = "1.89" [dependencies] bytemuck = "1.22.0" clap = { version = "4.5.42", features = ["derive"] } -color-eyre = "0.6.5" +color-eyre.workspace = true gdal = "0.19.0" -geo = "0.30.0" +geo.workspace = true geojson = "0.24.2" image = { version = "0.25.6", default-features = false, features = ["png"] } itertools = "0.14.0" -proj4rs = { version = "0.1.8", features = ["aeqd"] } -serde = "1.0.219" -serde_json = "1.0.143" -tracing = { version = "0.1.41" } -tracing-subscriber = { version = "0.3.19", features = ["env-filter"] } -rusqlite = { version = "0.38.0", features = ["bundled"] } memchr = "2.8.0" +serde_json.workspace = true +rusqlite = { version = "0.38.0", features = ["bundled"] } rayon = "1.12.0" +tracing = { version = "0.1.41" } +tracing-subscriber = { version = "0.3.19", features = ["env-filter"] } +tvs-lib = { path = "../tvs-lib", version = "0.1.0" } [dev-dependencies] tempfile = "3.20.0" diff --git a/crates/total-viewsheds/src/compute/area_of_interest.rs b/crates/total-viewsheds/src/compute/area_of_interest.rs index 0b6bd6a..9bbcdf4 100644 --- a/crates/total-viewsheds/src/compute/area_of_interest.rs +++ b/crates/total-viewsheds/src/compute/area_of_interest.rs @@ -21,15 +21,15 @@ impl Pruner { /// Convert user-provided lon/lat coordinates to a polygon. pub fn lonlat_coords_to_polygon( points: Vec<(f32, f32)>, - metadata: &crate::storage::metadata::MetaData, + metadata: &tvs_lib::metadata::MetaData, ) -> Result { let mut vertices = Vec::new(); for intrest_point in points { - let lonlat = crate::projection::LonLatCoord( + let lonlat = tvs_lib::projector::LonLatCoord( geo::coord!(x: intrest_point.0.into(), y: intrest_point.1.into()), ); - let dem_coord = crate::projection::Converter::lonlat_to_dem_coord(metadata, lonlat)?; + let dem_coord = tvs_lib::projector::Convert::lonlat_to_dem_coord(metadata, lonlat)?; let dem_point = geo::point!(x: dem_coord.0.x, y: dem_coord.0.y); vertices.push(dem_point); } @@ -52,10 +52,10 @@ impl Pruner { reason = "We're only dealing with a max of the DEM's width" )] /// Convert a DEM 1D index to a 2D coordinate. - pub fn convert_dem_id_to_coord(&self, dem_id: i64) -> crate::dem::Coordinate { + pub fn convert_dem_id_to_coord(&self, dem_id: i64) -> tvs_lib::dem::Coordinate { let x = dem_id.rem_euclid(self.width.into()) as f64; let y = dem_id.div_euclid(self.width.into()) as f64; - crate::dem::Coordinate(geo::coord! {x: x, y: y}) + tvs_lib::dem::Coordinate(geo::coord! {x: x, y: y}) } } @@ -65,10 +65,10 @@ mod tests { fn make_pruner() -> Pruner { let width = 300; - let metadata = crate::storage::metadata::MetaData { + let metadata = tvs_lib::metadata::MetaData { width, scale: 100.0, - centre: crate::projection::LonLatCoord((-3.1791, 51.4816).into()), + centre: tvs_lib::projector::LonLatCoord((-3.1791, 51.4816).into()), ..Default::default() }; let cardiff_10km = vec![ diff --git a/crates/total-viewsheds/src/compute/kernel.rs b/crates/total-viewsheds/src/compute/kernel.rs index 9e1350a..3522e43 100644 --- a/crates/total-viewsheds/src/compute/kernel.rs +++ b/crates/total-viewsheds/src/compute/kernel.rs @@ -209,11 +209,11 @@ mod test { let dem = &crate::tests::fixtures::bigger_dem(); let db_worker = crate::storage::worker::Worker::new_noop(); let width = 4; - let metadata = crate::storage::metadata::MetaData { + let metadata = tvs_lib::metadata::MetaData { width, scale: 1.0, max_line_of_sight: width, - centre: crate::projection::LonLatCoord((0.0, 0.0).into()), + centre: tvs_lib::projector::LonLatCoord((0.0, 0.0).into()), neighbourhood_size: 0, angle_subdivisions: 1, }; diff --git a/crates/total-viewsheds/src/compute/unrolled_los.rs b/crates/total-viewsheds/src/compute/unrolled_los.rs index fb78eaf..b77ddc7 100644 --- a/crates/total-viewsheds/src/compute/unrolled_los.rs +++ b/crates/total-viewsheds/src/compute/unrolled_los.rs @@ -17,7 +17,7 @@ fn generate_distances(max_los: usize, refraction: f32, scale: f32) -> (Vec, .map(|step| { let distance = (step as f32) * scale; let adjustment = - (distance * distance * adjusted_refraction) / crate::projection::EARTH_DIAMETER; + (distance * distance * adjusted_refraction) / tvs_lib::projector::EARTH_DIAMETER; (distance, adjustment) }) diff --git a/crates/total-viewsheds/src/main.rs b/crates/total-viewsheds/src/main.rs index 1f38d7e..e39f937 100644 --- a/crates/total-viewsheds/src/main.rs +++ b/crates/total-viewsheds/src/main.rs @@ -30,7 +30,6 @@ use std::mem; use tracing_subscriber::{Layer as _, layer::SubscriberExt as _, util::SubscriberInitExt as _}; mod config; -mod dem; mod dump_usage; mod los_pack; mod post_process; @@ -65,7 +64,6 @@ mod compute { mod storage { pub mod db; pub mod engine; - pub mod metadata; pub mod segments; pub mod worker; } @@ -76,10 +74,14 @@ mod output { pub mod bresenham; pub mod png; pub mod tiff; - pub mod viewshed; -} -mod projection; + /// Load, parse and reconstruct euclidean polygon viewsheds from their raw polar segments. + pub mod viewsheds { + pub mod join; + pub mod viewshed; + pub mod visible_polygon; + } +} fn main() -> Result<()> { color_eyre::install()?; @@ -91,15 +93,15 @@ fn main() -> Result<()> { config::Commands::Compute(compute_config) => compute(compute_config)?, config::Commands::Viewshed(viewshed_config) => { for coordinate in &viewshed_config.coordinates { - let geo_coord = projection::LonLatCoord( + let geo_coord = tvs_lib::projector::LonLatCoord( geo::coord! {x: f64::from(coordinate.0), y: f64::from(coordinate.1)}, ); - let (_, viewshed) = crate::output::viewshed::Viewshed::reconstruct( + let (_, viewshed) = crate::output::viewsheds::viewshed::Viewshed::reconstruct( viewshed_config.db_path.clone(), geo_coord, )?; - crate::output::viewshed::Reconstructor::save( + crate::output::viewsheds::viewshed::Viewshed::save( viewshed, &viewshed_config.output_dir, geo_coord, @@ -145,7 +147,7 @@ fn compute(config: &config::Compute) -> Result<()> { let max_line_of_sight_as_points = tile.width.div_euclid(3); - let mut dem = crate::dem::DEM::new( + let mut dem = tvs_lib::dem::DEM::new( tile.centre, tile.width, tile.scale, @@ -159,7 +161,7 @@ fn compute(config: &config::Compute) -> Result<()> { // Free up RAM drop(tile); - let dem_metadata = crate::storage::metadata::MetaData { + let dem_metadata = tvs_lib::metadata::MetaData { width: dem.width, scale: dem.scale, max_line_of_sight: max_line_of_sight_as_points, diff --git a/crates/total-viewsheds/src/output/ascii.rs b/crates/total-viewsheds/src/output/ascii.rs index 97d2adb..dc2d92d 100644 --- a/crates/total-viewsheds/src/output/ascii.rs +++ b/crates/total-viewsheds/src/output/ascii.rs @@ -3,7 +3,7 @@ #![cfg(test)] #![expect(clippy::indexing_slicing, reason = "This code is mostly for tests")] -use crate::output::viewshed::Viewshed; +use crate::output::viewsheds::viewshed::Viewshed; pub fn make_viewshed( elevations: &[i16], @@ -16,14 +16,14 @@ pub fn make_viewshed( x: viewshed_pov.x - dem_half_width, y: -(viewshed_pov.y - dem_half_width), }; - let coord_lonlat = crate::projection::Converter { base: dem.centre } + let coord_lonlat = tvs_lib::projector::Convert::new(dem.centre) .to_degrees(viewshed_pov_metric) .unwrap(); crate::run::test::compute(&mut dem, config.clone()); let (pov_coord, mut viewshed) = Viewshed::reconstruct(config.viewsheds_db_path, coord_lonlat).unwrap(); - let viewsheder = crate::output::viewshed::Viewshed { + let viewsheder = crate::output::viewsheds::viewshed::Viewshed { dem: &dem, pov_coord, }; @@ -43,9 +43,9 @@ pub fn make_viewshed( let mut maybe_from = None; for coordinate in line.coords_mut() { let projected = viewsheder - .convert_viewshed_coord_to_dem_coord(crate::output::viewshed::Coordinate( - *coordinate, - )) + .convert_viewshed_coord_to_dem_coord( + crate::output::viewsheds::viewshed::Coordinate(*coordinate), + ) .unwrap(); if maybe_from.is_none() { maybe_from = Some(projected); diff --git a/crates/total-viewsheds/src/output/tiff.rs b/crates/total-viewsheds/src/output/tiff.rs index 8bc3ce4..1de965d 100644 --- a/crates/total-viewsheds/src/output/tiff.rs +++ b/crates/total-viewsheds/src/output/tiff.rs @@ -3,7 +3,7 @@ use color_eyre::Result; /// Save an array of `f32`s (total surfaces, longest lines of sight) to a `.tiff` file. -pub fn save(dem: &crate::dem::DEM, data: &[f32], path: &std::path::PathBuf) -> Result<()> { +pub fn save(dem: &tvs_lib::dem::DEM, data: &[f32], path: &std::path::PathBuf) -> Result<()> { let driver = gdal::DriverManager::get_driver_by_name("GTiff")?; let mut dataset = driver.create_with_band_type::( diff --git a/crates/total-viewsheds/src/output/viewsheds/join.rs b/crates/total-viewsheds/src/output/viewsheds/join.rs new file mode 100644 index 0000000..be47313 --- /dev/null +++ b/crates/total-viewsheds/src/output/viewsheds/join.rs @@ -0,0 +1,32 @@ +//! Join the visible segments from each computed angle into the viewshed. + +impl super::visible_polygon::VisiblePolygon { + /// Convert polar segments to euclidean polygons. + pub fn parse_polar_segments( + data: &[Vec], + scale: f32, + ) -> geo::MultiPolygon { + let angle_count = data.len(); + let mut polygons = Vec::new(); + for (anglish, segments) in data.iter().enumerate() { + #[expect( + clippy::as_conversions, + clippy::cast_precision_loss, + reason = "The angle count should never strain the f32 mantissa" + )] + let (anglish_f32, angle_count_f32) = { (anglish as f32, angle_count as f32) }; + let polygoner = Self { + scale, + current_angle: (anglish_f32 / angle_count_f32) * 360.0, + }; + for segment in segments { + let opening = u32::from(segment.start()); + let closing = u32::from(segment.start() + segment.distance()); + let polygon = polygoner.make_visible_polygon(opening, closing); + polygons.push(polygon); + } + } + + geo::unary_union(polygons.iter()) + } +} diff --git a/crates/total-viewsheds/src/output/viewshed.rs b/crates/total-viewsheds/src/output/viewsheds/viewshed.rs similarity index 60% rename from crates/total-viewsheds/src/output/viewshed.rs rename to crates/total-viewsheds/src/output/viewsheds/viewshed.rs index 184cb8d..a085a82 100644 --- a/crates/total-viewsheds/src/output/viewshed.rs +++ b/crates/total-viewsheds/src/output/viewsheds/viewshed.rs @@ -1,7 +1,6 @@ -//! Reconstruct _individual_ viewsheds, not total viewsheds. +//! A viewshed is a (multi)polygon representing all the visibile terrain visible from a given point. use color_eyre::eyre::Result; -use geo::BooleanOps as _; /// A viewshed-based coordinate is projected to a metric system where the anchor is the viewshed's /// point of view. The other option would be a metric projection with an anchor in the DEM centre, @@ -13,22 +12,22 @@ pub struct Coordinate(pub geo::Coord); /// `Viewshed` pub struct Viewshed<'viewshed> { /// The DEM used to compute the final data. - pub dem: &'viewshed crate::dem::DEM, + pub dem: &'viewshed tvs_lib::dem::DEM, /// Coordinate of the observer for the viewshed we want to reconstruct. - pub pov_coord: crate::dem::Coordinate, + pub pov_coord: tvs_lib::dem::Coordinate, } impl Viewshed<'_> { /// Reconstruct a viewshed. pub fn reconstruct( db_path: std::path::PathBuf, - clicked_lonlat: crate::projection::LonLatCoord, - ) -> Result<(crate::dem::Coordinate, geo::MultiPolygon)> { + requested_lonlat: tvs_lib::projector::LonLatCoord, + ) -> Result<(tvs_lib::dem::Coordinate, geo::MultiPolygon)> { let db = crate::storage::db::DB::new(db_path)?; let metadata = db.load_metadata()?; tracing::debug!("Using metadata: {:?}", metadata); - let dem = crate::dem::DEM::new( + let dem = tvs_lib::dem::DEM::new( metadata.centre, metadata.width, metadata.scale, @@ -36,7 +35,7 @@ impl Viewshed<'_> { )?; let dem_coord = - crate::projection::Converter::lonlat_to_dem_coord(&metadata, clicked_lonlat)?; + tvs_lib::projector::Convert::lonlat_to_dem_coord(&metadata, requested_lonlat)?; let dem_id = dem.dem_coord_to_id(dem_coord); let (segments, pov_dem_coord) = match metadata.neighbourhood_size { @@ -60,7 +59,7 @@ impl Viewshed<'_> { clippy::cast_precision_loss, reason = "I assume that we never hit the 52 bit mantissa limit" )] - let coordinate = crate::dem::Coordinate( + let coordinate = tvs_lib::dem::Coordinate( geo::coord! { x: x as f64, y: y as f64, @@ -75,169 +74,37 @@ impl Viewshed<'_> { pov_dem_coord ); + let start = std::time::Instant::now(); let viewshed = Viewshed { dem: &dem, pov_coord: pov_dem_coord, }; - let mut reconstructor = Reconstructor::new(&viewshed, 0.0)?; - let polygon = reconstructor.parse_polar_segments(&segments); - Ok((pov_dem_coord, polygon)) - } - - /// Convert from the viewshed projection to DEM coordinates. - #[cfg(test)] - pub fn convert_viewshed_coord_to_dem_coord( - &self, - viewshed_coord: Coordinate, - ) -> Result { - let scale = f64::from(self.dem.scale); - let origin = crate::projection::Converter { - base: self.dem.centre, - } - .to_degrees((self.pov_coord.0.x, self.pov_coord.0.y).into())?; - let flipped = Coordinate(geo::Coord { - x: viewshed_coord.0.x, - y: -viewshed_coord.0.y, - }); - let projected_coord = crate::projection::Converter::change_metric_origin( - origin, - // The path back to (0,0) is exactly the opposite of the viewshed's point of view. - -self.pov_coord.0 * scale, - flipped.0 * scale, - )?; - Ok(projected_coord / scale) - } -} - -/// `Reconstructor` -// TODO: Find a way to make this part of [`Viewshed`]. -pub struct Reconstructor<'viewshed> { - /// Data for the entire viewshed. - viewshed: &'viewshed Viewshed<'viewshed>, - /// The current sector angle - current_angle: f32, -} - -impl<'viewshed> Reconstructor<'viewshed> { - /// Instantiate the reconstructor for a single angle. - /// - /// The reason we only reconstruct for a single angle is that the sector data (for the angle) - /// is most useful exposed as an iterator. And it's not so easy with lifetimes to keep around - /// the raw data and an iterator for each angle. - fn new(viewshed: &'viewshed Viewshed<'viewshed>, angle: f32) -> Result { let pov_id = viewshed.dem.dem_coord_to_id(viewshed.pov_coord); - let reconstructor = Self { - viewshed, - current_angle: angle, - }; - if !viewshed.dem.is_point_computable(pov_id) { - color_eyre::eyre::bail!( - "Point of view ({:?}) is not calculable", - reconstructor.viewshed.pov_coord - ); + color_eyre::eyre::bail!("Point of view ({:?}) is not calculable", viewshed.pov_coord); } - Ok(reconstructor) - } - - /// Convert polar segements to `GeoJson`. - pub fn parse_polar_segments( - &mut self, - data: &[Vec], - ) -> geo::MultiPolygon { - let mut viewshed_so_far = geo::MultiPolygon::empty(); - let angle_count = data.len(); - for (anglish, segments) in data.iter().enumerate() { - #[expect( - clippy::as_conversions, - clippy::cast_precision_loss, - reason = "The angle count should never strain the f32 mantissa" - )] - let (anglish_f32, angle_count_f32) = { (anglish as f32, angle_count as f32) }; - self.current_angle = (anglish_f32 / angle_count_f32) * 360.0; - for segment in segments { - let opening = u32::from(segment.start()); - let closing = u32::from(segment.start() + segment.distance()); - let polygon = self.make_visible_polygon(opening, closing); - viewshed_so_far = viewshed_so_far.union(&polygon); - } - } - - viewshed_so_far - } - - /// Convert an index along a line of sight into a coordinate. - fn index_to_coordinate(&self, index: u32) -> Coordinate { - let radians = self.current_angle.to_radians(); - let distance = f64::from(index); - - Coordinate(geo::coord! { - x: distance * f64::from(radians.cos()), - y: distance * f64::from(radians.sin()) - }) - } + let multi_polygon = super::visible_polygon::VisiblePolygon::parse_polar_segments( + &segments, + viewshed.dem.scale, + ); - /// Rotate a point about the centre of the viewshed. - #[expect( - clippy::suboptimal_flops, - reason = "I think readability is more important?" - )] - fn rotate_by(point: Coordinate, angle: f64) -> geo::Coord { - let dx = point.0.x; - let dy = point.0.y; - let cos = angle.to_radians().cos(); - let sin = angle.to_radians().sin(); - geo::coord! { - x: dx * cos - dy * sin, - y: dx * sin + dy * cos - } - } + tracing::debug!("Viewshed reconstructed in {:?}.", start.elapsed()); - /// Make a single polygon representing a visible region of the planet. - fn make_visible_polygon(&self, opening_index: u32, closing_index: u32) -> geo::Polygon { - let opening_coord = self.index_to_coordinate(opening_index); - let closing_coord = self.index_to_coordinate(closing_index); - - let spread = 0.5001f64; - let bottom_left = Self::rotate_by(opening_coord, spread); - let bottom_right = Self::rotate_by(opening_coord, -spread); - let top_left = Self::rotate_by(closing_coord, spread); - let top_right = Self::rotate_by(closing_coord, -spread); - - let scale = f64::from(self.viewshed.dem.scale); - - geo::Polygon::new( - geo::LineString(vec![ - bottom_left * scale, - bottom_right * scale, - top_right * scale, - top_left * scale, - bottom_left * scale, - ]), - vec![], - ) + Ok((pov_dem_coord, multi_polygon)) } - /// Save the viewshed to disk. #[expect( - clippy::panic_in_result_fn, clippy::panic, - reason = "The closures expects () so I don't think there's any other way?" + reason = "The closures expect () so I don't think there's any other way?" )] - pub fn save( + /// Convert the local metric coordinates of the viewshed to WGS84 lon/lat coordinates. + fn convert_viewshed_coords_to_lonlat( mut viewshed: geo::MultiPolygon, - output_directory: &std::path::Path, - viewshed_latlon: crate::projection::LonLatCoord, - ) -> Result<()> { - let filename = format!("{}-{}.json", viewshed_latlon.0.x, viewshed_latlon.0.y); - let directory = output_directory.join("viewsheds"); - std::fs::create_dir_all(&directory)?; - let path = directory.join(filename); - let projector = crate::projection::Converter { - base: viewshed_latlon, - }; + viewshed_latlon: tvs_lib::projector::LonLatCoord, + ) -> geo::MultiPolygon { + let projector = tvs_lib::projector::Convert::new(viewshed_latlon); for point in viewshed.iter_mut() { point.exterior_mut(|line| { @@ -274,147 +141,56 @@ impl<'viewshed> Reconstructor<'viewshed> { } }); } - let json = geojson::GeoJson::from(&viewshed).to_string(); - std::fs::write(path, json)?; - Ok(()) - } -} - -#[cfg(test)] -mod test { - use super::*; - use crate::output::ascii::assert_viewshed; - - fn builder<'viewshed>(viewshed: &'viewshed Viewshed, angle: f32) -> Reconstructor<'viewshed> { - Reconstructor::new(viewshed, angle).unwrap() + viewshed } - #[derive(Debug)] - struct VisiblePolygonFor { - pov: geo::Coord, - angle: f32, - opening_index: u32, - closing_index: u32, - } - - fn make_visible_polygon_for(setup: &VisiblePolygonFor) -> Vec { - let dem = crate::run::test::make_dem(&crate::tests::fixtures::single_peak_dem()); - let viewshed = Viewshed { - dem: &dem, - pov_coord: crate::dem::Coordinate(setup.pov), - }; - let viewsheder = builder(&viewshed, setup.angle); - let polygon = viewsheder.make_visible_polygon(setup.opening_index, setup.closing_index); - - let mut polygon_as_dem_coords = Vec::new(); - for coord in &polygon.exterior().0 { - let converted_coord = viewsheder - .viewshed - .convert_viewshed_coord_to_dem_coord(Coordinate(*coord)) - .unwrap(); - polygon_as_dem_coords.push(round_coordinate(converted_coord)); - } - polygon_as_dem_coords - } + /// Save the viewshed to disk. + pub fn save( + viewshed: geo::MultiPolygon, + output_directory: &std::path::Path, + viewshed_latlon: tvs_lib::projector::LonLatCoord, + ) -> Result<()> { + let filename = format!("{}-{}.json", viewshed_latlon.0.x, viewshed_latlon.0.y); + let directory = output_directory.join("viewsheds"); + std::fs::create_dir_all(&directory)?; + let path = directory.join(filename); - fn round(float: f64) -> f64 { - let factor = 10f64.powi(7); - (float * factor).round() / factor - } + let viewshed_lonlat = Self::convert_viewshed_coords_to_lonlat(viewshed, viewshed_latlon); + let json = geojson::GeoJson::from(&viewshed_lonlat).to_string(); + std::fs::write(path, json)?; - fn round_coordinate(coordinate: geo::Coord) -> geo::Coord { - geo::coord! { - x: round(coordinate.x), - y: round(coordinate.y), - } + Ok(()) } +} - // Guide for the following tests: - // - // 0 1 2 3 4 5 6 7 8 - // 0 . . . . . . . . . - // 1 . . . . . .d . . . - // 2 . . . . .a . ) . . - // 3 . . . . . ( . c. . - // 4 . . . . o . b. . . - // 5 . . . . . . . . . - // 6 . . . . . . . . . - // 7 . . . . . . . . . - // 8 . . . . . . . . . - // - mod from_centre_to_top_right { - use super::*; - - const POV: geo::Coord = geo::coord! {x: 4.0, y: 4.0}; - const ANGLE: f32 = 45.0; - - // The polygon we're making is `abcd` from the above guide. - #[test] - fn making_a_visible_polygon() { - assert_eq!( - make_visible_polygon_for(&VisiblePolygonFor { - pov: POV, - angle: ANGLE, - opening_index: 1, - closing_index: 2, - }), - vec![ - (4.7009067, 3.2867503), - (4.7132503, 3.2990939), - (5.4265022, 2.5981862), - (5.401815, 2.5734989), - (4.7009067, 3.2867503) - ] - .into_iter() - .map(Into::into) - .collect::>() - ); - } +#[cfg(test)] +impl Viewshed<'_> { + /// Convert from the viewshed projection to DEM coordinates. + pub fn convert_viewshed_coord_to_dem_coord( + &self, + viewshed_coord: Coordinate, + ) -> Result { + let scale = f64::from(self.dem.scale); + let origin = tvs_lib::projector::Convert::new(self.dem.centre) + .to_degrees((self.pov_coord.0.x, self.pov_coord.0.y).into())?; + let flipped = Coordinate(geo::Coord { + x: viewshed_coord.0.x, + y: -viewshed_coord.0.y, + }); + let projected_coord = tvs_lib::projector::Convert::change_metric_origin( + origin, + // The path back to (0,0) is exactly the opposite of the viewshed's point of view. + -self.pov_coord.0 * scale, + flipped.0 * scale, + )?; + Ok(projected_coord / scale) } +} - // Guide for the following tests: - // - // 0 1 2 3 4 5 6 7 8 - // 0 . . . . . . . . . - // 1 . . . . . . . . . - // 2 . . . . . . . . . - // 3 . . . . . . . . . - // 4 . . . . . . . . . - // 5 . . . o .a . . . . - // 6 . . . . ( .d . . . - // 7 . . . . b. ) . . . - // 8 . . . . . c. . . . - // - mod from_bottom_left_to_bottom_right { - use super::*; - - const POV: geo::Coord = geo::coord! {x: 3.0, y: 5.0}; - const ANGLE: f32 = 135.0 + 180.0; - - // The polygon we're making is `abcd` from the above guide. - #[test] - fn making_a_visible_polygon() { - assert_eq!( - make_visible_polygon_for(&VisiblePolygonFor { - pov: POV, - angle: ANGLE, - opening_index: 1, - closing_index: 2, - }), - vec![ - (3.7132498, 5.7009093), - (3.7009061, 5.7132529), - (4.4018138, 6.4265049), - (4.4265011, 6.4018176), - (3.7132498, 5.7009093) - ] - .into_iter() - .map(Into::into) - .collect::>() - ); - } - } +#[cfg(test)] +mod test { + use crate::output::ascii::assert_viewshed; const SUMMIT_VIEWSHED: [&str; 12] = [ "████████████████████████", @@ -526,7 +302,7 @@ mod test { let temp_db_for_viewsheds = tempfile::NamedTempFile::new().unwrap(); let config = crate::run::Config { - dem_metadata: crate::storage::metadata::MetaData { + dem_metadata: tvs_lib::metadata::MetaData { neighbourhood_size: neighbourhood_size.try_into().unwrap(), ..crate::run::test::big_dem_metadata() }, diff --git a/crates/total-viewsheds/src/output/viewsheds/visible_polygon.rs b/crates/total-viewsheds/src/output/viewsheds/visible_polygon.rs new file mode 100644 index 0000000..3887360 --- /dev/null +++ b/crates/total-viewsheds/src/output/viewsheds/visible_polygon.rs @@ -0,0 +1,203 @@ +//! Create euclidean polygons from polar segments. + +/// `VisiblePolygon` +pub struct VisiblePolygon { + /// Scale of DEM data. + pub scale: f32, + /// The current sector angle. + pub current_angle: f32, +} + +impl VisiblePolygon { + /// Convert an index along a line of sight into a coordinate. + fn index_to_coordinate(&self, index: u32) -> super::viewshed::Coordinate { + let radians = self.current_angle.to_radians(); + let distance = f64::from(index); + + super::viewshed::Coordinate(geo::coord! { + x: distance * f64::from(radians.cos()), + y: distance * f64::from(radians.sin()) + }) + } + + /// Rotate a point about the centre of the viewshed. + #[expect( + clippy::suboptimal_flops, + reason = "I think readability is more important?" + )] + fn rotate_by(point: super::viewshed::Coordinate, angle: f64) -> geo::Coord { + let dx = point.0.x; + let dy = point.0.y; + let cos = angle.to_radians().cos(); + let sin = angle.to_radians().sin(); + geo::coord! { + x: dx * cos - dy * sin, + y: dx * sin + dy * cos + } + } + + /// Make a single polygon representing a visible region of the planet. + pub fn make_visible_polygon(&self, opening_index: u32, closing_index: u32) -> geo::Polygon { + let opening_coord = self.index_to_coordinate(opening_index); + let closing_coord = self.index_to_coordinate(closing_index); + + let spread = 0.5001f64; + let bottom_left = Self::rotate_by(opening_coord, spread); + let bottom_right = Self::rotate_by(opening_coord, -spread); + let top_left = Self::rotate_by(closing_coord, spread); + let top_right = Self::rotate_by(closing_coord, -spread); + + let scale = f64::from(self.scale); + + geo::Polygon::new( + geo::LineString(vec![ + bottom_left * scale, + bottom_right * scale, + top_right * scale, + top_left * scale, + bottom_left * scale, + ]), + vec![], + ) + } +} + +#[cfg(test)] +mod test { + fn builder( + viewshed: &crate::output::viewsheds::viewshed::Viewshed, + angle: f32, + ) -> crate::output::viewsheds::visible_polygon::VisiblePolygon { + crate::output::viewsheds::visible_polygon::VisiblePolygon { + scale: viewshed.dem.scale, + current_angle: angle, + } + } + + #[derive(Debug)] + struct VisiblePolygonFor { + pov: geo::Coord, + angle: f32, + opening_index: u32, + closing_index: u32, + } + + fn make_visible_polygon_for(setup: &VisiblePolygonFor) -> Vec { + let dem = crate::run::test::make_dem(&crate::tests::fixtures::single_peak_dem()); + let viewshed = crate::output::viewsheds::viewshed::Viewshed { + dem: &dem, + pov_coord: tvs_lib::dem::Coordinate(setup.pov), + }; + let viewsheder = builder(&viewshed, setup.angle); + let polygon = viewsheder.make_visible_polygon(setup.opening_index, setup.closing_index); + + let mut polygon_as_dem_coords = Vec::new(); + for coord in &polygon.exterior().0 { + let converted_coord = viewshed + .convert_viewshed_coord_to_dem_coord( + crate::output::viewsheds::viewshed::Coordinate(*coord), + ) + .unwrap(); + polygon_as_dem_coords.push(round_coordinate(converted_coord)); + } + polygon_as_dem_coords + } + + fn round(float: f64) -> f64 { + let factor = 10f64.powi(7); + (float * factor).round() / factor + } + + fn round_coordinate(coordinate: geo::Coord) -> geo::Coord { + geo::coord! { + x: round(coordinate.x), + y: round(coordinate.y), + } + } + + // Guide for the following tests: + // + // 0 1 2 3 4 5 6 7 8 + // 0 . . . . . . . . . + // 1 . . . . . .d . . . + // 2 . . . . .a . ) . . + // 3 . . . . . ( . c. . + // 4 . . . . o . b. . . + // 5 . . . . . . . . . + // 6 . . . . . . . . . + // 7 . . . . . . . . . + // 8 . . . . . . . . . + // + mod from_centre_to_top_right { + use super::*; + + const POV: geo::Coord = geo::coord! {x: 4.0, y: 4.0}; + const ANGLE: f32 = 45.0; + + // The polygon we're making is `abcd` from the above guide. + #[test] + fn making_a_visible_polygon() { + assert_eq!( + make_visible_polygon_for(&VisiblePolygonFor { + pov: POV, + angle: ANGLE, + opening_index: 1, + closing_index: 2, + }), + vec![ + (4.7009067, 3.2867503), + (4.7132503, 3.2990939), + (5.4265022, 2.5981862), + (5.401815, 2.5734989), + (4.7009067, 3.2867503) + ] + .into_iter() + .map(Into::into) + .collect::>() + ); + } + } + + // Guide for the following tests: + // + // 0 1 2 3 4 5 6 7 8 + // 0 . . . . . . . . . + // 1 . . . . . . . . . + // 2 . . . . . . . . . + // 3 . . . . . . . . . + // 4 . . . . . . . . . + // 5 . . . o .a . . . . + // 6 . . . . ( .d . . . + // 7 . . . . b. ) . . . + // 8 . . . . . c. . . . + // + mod from_bottom_left_to_bottom_right { + use super::*; + + const POV: geo::Coord = geo::coord! {x: 3.0, y: 5.0}; + const ANGLE: f32 = 135.0 + 180.0; + + // The polygon we're making is `abcd` from the above guide. + #[test] + fn making_a_visible_polygon() { + assert_eq!( + make_visible_polygon_for(&VisiblePolygonFor { + pov: POV, + angle: ANGLE, + opening_index: 1, + closing_index: 2, + }), + vec![ + (3.7132498, 5.7009093), + (3.7009061, 5.7132529), + (4.4018138, 6.4265049), + (4.4265011, 6.4018176), + (3.7132498, 5.7009093) + ] + .into_iter() + .map(Into::into) + .collect::>() + ); + } + } +} diff --git a/crates/total-viewsheds/src/run.rs b/crates/total-viewsheds/src/run.rs index 017251b..3b3919d 100644 --- a/crates/total-viewsheds/src/run.rs +++ b/crates/total-viewsheds/src/run.rs @@ -11,7 +11,7 @@ where /// User configuration. pub config: Config, /// The Digital Elevation Model that we're computing. - pub dem: &'compute mut crate::dem::DEM, + pub dem: &'compute mut tvs_lib::dem::DEM, /// Keeps track of the cumulative surfaces from every angle. pub total_surfaces: Vec, /// Keeps track of the longest lines of sight. @@ -38,7 +38,7 @@ pub struct Config { /// Where to store the viewshed pub viewsheds_db_path: PathBuf, /// Metadata about the DEM and compute run. - pub dem_metadata: crate::storage::metadata::MetaData, + pub dem_metadata: tvs_lib::metadata::MetaData, /// Polygon for Area of Interest pub area_of_interest: geo::Polygon, /// Should a database aggregate per thread? @@ -51,7 +51,7 @@ pub struct Config { impl<'compute> Compute<'compute> { /// Instantiate. - pub fn new(config: Config, dem: &'compute mut crate::dem::DEM) -> Result { + pub fn new(config: Config, dem: &'compute mut tvs_lib::dem::DEM) -> Result { if Self::is_process_viewsheds(&config.process) && !cfg!(any(test, feature = "ring_data")) { color_eyre::eyre::bail!( "Viewshed storage is only supported with the ring_data feature, \ @@ -160,10 +160,10 @@ pub mod test { use super::*; use googletest::prelude::*; - pub fn make_dem(elevations: &[i16]) -> crate::dem::DEM { + pub fn make_dem(elevations: &[i16]) -> tvs_lib::dem::DEM { let width = elevations.len().isqrt() as u32; - let mut dem = crate::dem::DEM::new( - crate::projection::LonLatCoord((33.33, 33.33).into()), + let mut dem = tvs_lib::dem::DEM::new( + tvs_lib::projector::LonLatCoord((33.33, 33.33).into()), width, 1.0, width / 3, @@ -173,18 +173,18 @@ pub mod test { dem } - pub fn default_metadata() -> crate::storage::metadata::MetaData { - crate::storage::metadata::MetaData { + pub fn default_metadata() -> tvs_lib::metadata::MetaData { + tvs_lib::metadata::MetaData { scale: 1.0, ..Default::default() } } - pub fn big_dem_metadata() -> crate::storage::metadata::MetaData { - crate::storage::metadata::MetaData { + pub fn big_dem_metadata() -> tvs_lib::metadata::MetaData { + tvs_lib::metadata::MetaData { width: 12, max_line_of_sight: 4, - centre: crate::projection::LonLatCoord((33.33, 33.33).into()), + centre: tvs_lib::projector::LonLatCoord((33.33, 33.33).into()), ..crate::run::test::default_metadata() } } @@ -207,7 +207,7 @@ pub mod test { } } - pub fn compute(dem: &mut crate::dem::DEM, config: Config) -> Compute<'_> { + pub fn compute(dem: &mut tvs_lib::dem::DEM, config: Config) -> Compute<'_> { let mut compute = Compute::new(config, dem).unwrap(); compute.run().unwrap(); compute @@ -301,7 +301,7 @@ pub mod test { let compute_very_refraction = compute( &mut dem_for_very_refraction, super::Config { - refraction: -crate::projection::EARTH_DIAMETER, + refraction: -tvs_lib::projector::EARTH_DIAMETER, dem_metadata: crate::run::test::big_dem_metadata(), ..default_config(&temp_db_for_very_refraction) }, @@ -325,7 +325,7 @@ pub mod test { let compute_small_scale = compute( &mut dem_for_small_scale, super::Config { - dem_metadata: crate::storage::metadata::MetaData { + dem_metadata: tvs_lib::metadata::MetaData { scale: 0.01, ..crate::run::test::big_dem_metadata() }, @@ -348,7 +348,7 @@ pub mod test { let compute_big_scale = compute( &mut dem_for_big_scale, super::Config { - dem_metadata: crate::storage::metadata::MetaData { + dem_metadata: tvs_lib::metadata::MetaData { scale: 100.0, ..crate::run::test::big_dem_metadata() }, diff --git a/crates/total-viewsheds/src/storage/db.rs b/crates/total-viewsheds/src/storage/db.rs index bd29467..b440e67 100644 --- a/crates/total-viewsheds/src/storage/db.rs +++ b/crates/total-viewsheds/src/storage/db.rs @@ -26,7 +26,7 @@ impl DB { } /// Save the metadata. - pub fn save_metadata(&self, metadata: &super::metadata::MetaData) -> Result<()> { + pub fn save_metadata(&self, metadata: &tvs_lib::metadata::MetaData) -> Result<()> { tracing::debug!("Saving metadata to {:?}...", self.connection.path()); self.connection.execute( @@ -49,13 +49,13 @@ impl DB { } /// Load the metadata. - pub fn load_metadata(&self) -> Result { + pub fn load_metadata(&self) -> Result { tracing::debug!("Loading metadata from {:?}...", self.connection.path()); let metadata_string: String = self.connection .query_row("SELECT json FROM metadata", [], |row| row.get(0))?; - let metadata: super::metadata::MetaData = serde_json::from_str(&metadata_string)?; + let metadata: tvs_lib::metadata::MetaData = serde_json::from_str(&metadata_string)?; tracing::info!("Loaded metadata: {metadata:?}"); Ok(metadata) diff --git a/crates/total-viewsheds/src/tile.rs b/crates/total-viewsheds/src/tile.rs index 964214b..e38ab3d 100644 --- a/crates/total-viewsheds/src/tile.rs +++ b/crates/total-viewsheds/src/tile.rs @@ -15,7 +15,7 @@ pub struct Tile { /// All the elevation data. pub data: Vec, /// The lat/lon coordinates of the centre of the tile. - pub centre: crate::projection::LonLatCoord, + pub centre: tvs_lib::projector::LonLatCoord, } impl Tile { @@ -126,7 +126,9 @@ impl Tile { } /// Get the lat/lon centre of the tile by querying the projection's definition. - fn get_centre_by_projection(dataset: &gdal::Dataset) -> Result { + fn get_centre_by_projection( + dataset: &gdal::Dataset, + ) -> Result { const ALLOWED_DEVIATION: f64 = 0.001; // 1 milimetre let geometric_centre = Self::get_metric_centre(dataset)?; if geometric_centre.x.abs() > ALLOWED_DEVIATION @@ -155,7 +157,7 @@ impl Tile { .and_then(|value| value.parse::().ok()) .context("Couldn't find `lon_0` in projection defintion")?; - Ok(crate::projection::LonLatCoord( + Ok(tvs_lib::projector::LonLatCoord( geo::coord! { x: lon_0, y: lat_0 }, )) } @@ -174,7 +176,7 @@ mod test { let tile = Tile::load(&config).unwrap(); assert_eq!( tile.centre, - crate::projection::LonLatCoord(geo::Coord { + tvs_lib::projector::LonLatCoord(geo::Coord { x: -0.1278, y: 51.5074 }) diff --git a/crates/total-viewsheds/src/workers.rs b/crates/total-viewsheds/src/workers.rs index e1f4c55..7f84de9 100644 --- a/crates/total-viewsheds/src/workers.rs +++ b/crates/total-viewsheds/src/workers.rs @@ -77,7 +77,7 @@ fn compilation_worker( /// a noop worker pub fn init_worker>( path: P, - meta_data: &crate::storage::metadata::MetaData, + meta_data: &tvs_lib::metadata::MetaData, is_db_worker: bool, ) -> Result { if !is_db_worker { diff --git a/crates/tvs-lib/Cargo.toml b/crates/tvs-lib/Cargo.toml new file mode 100644 index 0000000..a6d005d --- /dev/null +++ b/crates/tvs-lib/Cargo.toml @@ -0,0 +1,14 @@ +[package] +name = "tvs-lib" +version = "0.1.0" +edition = "2024" + +[dependencies] +color-eyre.workspace = true +geo.workspace = true +proj4rs = { version = "0.1.9", features = ["aeqd"] } +rstar = "0.12.2" +serde.workspace = true + +[lints] +workspace = true diff --git a/crates/total-viewsheds/src/dem.rs b/crates/tvs-lib/src/dem.rs similarity index 92% rename from crates/total-viewsheds/src/dem.rs rename to crates/tvs-lib/src/dem.rs index b326564..e769d93 100644 --- a/crates/total-viewsheds/src/dem.rs +++ b/crates/tvs-lib/src/dem.rs @@ -12,10 +12,12 @@ use color_eyre::Result; /// * Each integer maps exactly to a point in the DEM data. /// /// These coordinates are used by the kernel and tests. +#[expect(clippy::exhaustive_structs, reason = "This should never change")] #[derive(Debug, Clone, Copy, PartialEq)] pub struct Coordinate(pub geo::Coord); /// `DEM` +#[non_exhaustive] pub struct DEM { /// All the elevation data. pub elevations: Vec, @@ -30,7 +32,7 @@ pub struct DEM { /// The size of each point in meters. pub scale: f32, /// The geographic location of the centre of the DEM tile. - pub centre: crate::projection::LonLatCoord, + pub centre: crate::projector::LonLatCoord, /// The maximum distance in terms of points to search. pub max_los_as_points: u32, /// The total number of points that can have full viewsheds calculated for them. @@ -39,8 +41,12 @@ pub struct DEM { impl DEM { /// `Instantiate` + /// + /// # Errors + /// If width is not divisible by 3 + #[inline] pub fn new( - centre_latlon: crate::projection::LonLatCoord, + centre_latlon: crate::projector::LonLatCoord, width: u32, scale: f32, max_line_of_sight_as_points: u32, @@ -77,6 +83,8 @@ impl DEM { } /// Is the DEM ID within the DEM such that a valid viewshed can be made for it? + #[inline] + #[must_use] pub fn is_point_computable(&self, dem_id: u32) -> bool { let coord = self.convert_dem_id_to_coord(dem_id).0; let lower = f64::from(self.max_los_as_points); @@ -85,6 +93,8 @@ impl DEM { } /// Convert a DEM 1D index to a 2D coordinate. + #[inline] + #[must_use] pub fn convert_dem_id_to_coord(&self, dem_id: u32) -> Coordinate { let x = f64::from(dem_id.rem_euclid(self.width)); let y = f64::from(dem_id.div_euclid(self.width)); @@ -92,6 +102,8 @@ impl DEM { } /// Convert a DEM coordinate to a DEM ID. + #[inline] + #[must_use] pub fn dem_coord_to_id(&self, coord: Coordinate) -> u32 { let x = coord.0.x.round(); let yish = coord.0.y.round() * f64::from(self.width); @@ -113,6 +125,7 @@ impl DEM { reason = "We don't want to output GBs of data!" )] impl std::fmt::Debug for DEM { + #[inline] fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result { f.debug_struct("DEM") .field("width", &self.width) diff --git a/crates/tvs-lib/src/lib.rs b/crates/tvs-lib/src/lib.rs new file mode 100644 index 0000000..cde6056 --- /dev/null +++ b/crates/tvs-lib/src/lib.rs @@ -0,0 +1,6 @@ +//! Code we want to share amongst projects. + +pub mod dem; +pub mod metadata; +pub mod projector; +pub mod utils; diff --git a/crates/total-viewsheds/src/storage/metadata.rs b/crates/tvs-lib/src/metadata.rs similarity index 68% rename from crates/total-viewsheds/src/storage/metadata.rs rename to crates/tvs-lib/src/metadata.rs index c1cdd97..f245ddd 100644 --- a/crates/total-viewsheds/src/storage/metadata.rs +++ b/crates/tvs-lib/src/metadata.rs @@ -1,6 +1,10 @@ //! Struct for storing essential data about the underlying DEM for which viewsheds are created. /// Metadata about the viewshed data. +#[expect( + clippy::exhaustive_structs, + reason = "It's just nice to contsruct using struct syntax" +)] #[derive(serde::Serialize, serde::Deserialize, Default, Debug, Clone)] pub struct MetaData { /// The width of the 2D grid of elevation data. The algorithm requires that the grid be square, @@ -8,14 +12,16 @@ pub struct MetaData { pub width: u32, /// The diameter in meters each point of the data covers. pub scale: f32, - /// The maximum line of sight, in points). + /// The maximum line of sight (in meters) that was used to calculate the ring data. It is needed + /// to instantiate the `DEM` struct and therefore reconstruct the bands of sight used to create + /// the ring data. pub max_line_of_sight: u32, /// The lat/lon coordinates for the centre of the 2D DEM grid. Used for accurately converting /// between degree and metric coordinate systems. - pub centre: crate::projection::LonLatCoord, + pub centre: crate::projector::LonLatCoord, /// The size of the region (in raster points) within which we will find the viewsheds with the /// largest surface area. Used for reducing the final size of viewshed data saved to disk. pub neighbourhood_size: u32, - /// The number of angle subdivisions used + /// The number of angle subdivisions used. pub angle_subdivisions: u32, } diff --git a/crates/total-viewsheds/src/projection.rs b/crates/tvs-lib/src/projector.rs similarity index 62% rename from crates/total-viewsheds/src/projection.rs rename to crates/tvs-lib/src/projector.rs index 65700d5..b9be800 100644 --- a/crates/total-viewsheds/src/projection.rs +++ b/crates/tvs-lib/src/projector.rs @@ -1,39 +1,104 @@ //! Project coordinates between different systems. -use color_eyre::Result; +use color_eyre::{Result, eyre::ContextCompat as _}; -/// Diameter of the Earth in meters. So that some points are not visible simply -/// by virtue of the earth's spherical shape. -pub const EARTH_DIAMETER: f32 = 12_742_000.0; +/// The radius of the planet in kilometers. +pub const EARTH_RADIUS: f32 = 6378.137; -// TODO: Rename to `LonLatCoord`. -/// A latitude/longtitude coordinate. +/// The diameter of the planet in meters. +pub const EARTH_DIAMETER: f32 = 12_756_274.0; + +/// A longtitude/latitude coordinate. +#[expect(clippy::exhaustive_structs, reason = "This should never change")] #[derive(Debug, Clone, Copy, serde::Serialize, serde::Deserialize, PartialEq, Default)] pub struct LonLatCoord(pub geo::Coord); +impl LonLatCoord { + /// Parse coordinates from a string. + /// + /// # Errors + /// On parsing errors. + #[inline] + pub fn parse(coordinates: &str) -> Result { + let mut parts = coordinates.split(','); + let lon_str = parts.next().context("missing longtitude")?.trim(); + let lat_str = parts.next().context("missing latitude")?.trim(); + + Ok(Self(geo::Coord { + x: lon_str.parse()?, + y: lat_str.parse()?, + })) + } +} + +impl rstar::Point for LonLatCoord { + type Scalar = f64; + const DIMENSIONS: usize = 2; + + #[inline] + fn generate(mut generator: impl FnMut(usize) -> Self::Scalar) -> Self { + Self(geo::coord! { + x: generator(0), + y: generator(1), + }) + } + + #[inline] + fn nth(&self, index: usize) -> Self::Scalar { + match index { + 0 => self.0.x, + 1 => self.0.y, + #[expect(clippy::unreachable, reason = "This is just a 2D coordinate.")] + _ => unreachable!(), + } + } + + #[inline] + fn nth_mut(&mut self, index: usize) -> &mut Self::Scalar { + match index { + 0 => &mut self.0.x, + 1 => &mut self.0.y, + #[expect(clippy::unreachable, reason = "This is just a 2D coordinate.")] + _ => unreachable!(), + } + } +} + /// Convert between different coordinate system. -pub struct Converter { - /// The lat/lon base coordinates for the AEQD mercator projected coordinates. +#[non_exhaustive] +pub struct Convert { + /// The lon/lat base coordinates for the AEQD mercator projected coordinates. pub base: LonLatCoord, } -impl Converter { +impl Convert { + /// Instantiate. + #[inline] + #[must_use] + pub const fn new(base: LonLatCoord) -> Self { + Self { base } + } + /// The projection description for lat/lon. - pub fn degrees_projection() -> Result { - let string = "+proj=latlong +datum=WGS84"; + fn degrees_projection() -> Result { + let string = "+proj=latlong +datum=WGS84 +over"; Ok(proj4rs::Proj::from_proj_string(string)?) } /// The projection description for the AEQD metric projection. fn meters_projection(&self) -> Result { let string = format!( - "+proj=aeqd +lat_0={} +lon_0={} +datum=WGS84", + "+proj=aeqd +lat_0={} +lon_0={} +datum=WGS84 +over", self.base.0.y, self.base.0.x ); Ok(proj4rs::Proj::from_proj_string(&string)?) } + #[inline] /// Convert from degrees to the AEQD metric projection. + /// + /// # Errors + /// When projection conversion fails. pub fn to_meters(&self, source: LonLatCoord) -> Result { let mut converted = (source.0.x.to_radians(), source.0.y.to_radians(), 0.0f64); proj4rs::transform::transform( @@ -45,7 +110,11 @@ impl Converter { Ok(geo::coord! { x: converted.0, y: converted.1 }) } + #[inline] /// Convert from the AEQD metric projection to degrees. + /// + /// # Errors + /// When projection conversion fails. pub fn to_degrees(&self, source: geo::Coord) -> Result { let mut converted = (source.x, source.y, 0.0f64); proj4rs::transform::transform( @@ -59,10 +128,30 @@ impl Converter { )) } + #[inline] + #[must_use] + /// Calculate the width of a degree in meters at the given latitude. + pub fn meters_per_degree(latitude: f64) -> f32 { + #[expect( + clippy::cast_possible_truncation, + clippy::as_conversions, + reason = "It's the only way" + )] + (EARTH_RADIUS + * 1000.0 + * latitude.to_radians().cos() as f32 + * (std::f32::consts::PI / 180.0)) + .abs() + } + /// Convert a lat/lon to a DEM coordinate. + /// + /// # Errors + /// If projection fails. + #[inline] pub fn lonlat_to_dem_coord( - metadata: &crate::storage::metadata::MetaData, - latlon: crate::projection::LonLatCoord, + metadata: &crate::metadata::MetaData, + latlon: LonLatCoord, ) -> Result { let width = f64::from(metadata.width - 1); let scale = f64::from(metadata.scale); @@ -82,9 +171,12 @@ impl Converter { Ok(dem_coord) } - #[cfg(test)] - /// Chante the anchor of the AEQD projection. This just gives slightly more accuracy when + /// Change the anchor of the AEQD projection. This just gives slightly more accuracy when /// reconstructing larger viewsheds on larger DEMs. + /// + /// # Errors + /// If projection fails. + #[inline] pub fn change_metric_origin( // The lat/lon of the DEM's top-left corner source_degrees_anchor: LonLatCoord, @@ -130,7 +222,7 @@ mod test { x: -2.5879, y: 51.4545, }); - let converter = Converter { base }; + let converter = Convert { base }; assert_eq!( converter.to_meters(base).unwrap(), geo::Coord { x: 0.0, y: 0.0 } @@ -143,7 +235,7 @@ mod test { x: -2.5879, y: 51.4545, }); - let converter = Converter { base }; + let converter = Convert { base }; assert_eq!( converter .to_meters(LonLatCoord(geo::Coord { @@ -164,7 +256,7 @@ mod test { x: -2.5879, y: 51.4545, }); - let converter = Converter { base }; + let converter = Convert { base }; assert_eq!( converter.to_degrees(geo::Coord { x: 0.0, y: 0.0 }).unwrap(), LonLatCoord(geo::Coord { @@ -180,7 +272,7 @@ mod test { x: -2.5879, y: 51.4545, }); - let converter = Converter { base }; + let converter = Convert { base }; assert_eq!( converter .to_degrees(geo::Coord { @@ -194,22 +286,4 @@ mod test { }) ); } - - #[test] - fn latlon_to_dem_coord() { - let centre = crate::projection::LonLatCoord((-33.33f64, 12.34f64).into()); - let metadata = crate::storage::metadata::MetaData { - width: 102, - scale: 5.0, - max_line_of_sight: 250, - centre, - neighbourhood_size: 0, - angle_subdivisions: 1, - }; - - assert_eq!( - Converter::lonlat_to_dem_coord(&metadata, centre).unwrap(), - crate::dem::Coordinate((50.5f64, 50.5f64).into()) - ); - } } diff --git a/crates/tvs-lib/src/utils.rs b/crates/tvs-lib/src/utils.rs new file mode 100644 index 0000000..0dd6e15 --- /dev/null +++ b/crates/tvs-lib/src/utils.rs @@ -0,0 +1,77 @@ +//! Helper code. + +use color_eyre::Result; + +#[inline] +/// Convert a lon/lat to a DEM ID. +/// +/// # Errors +/// +/// If projection errors. +pub fn lonlat_to_dem_id( + metadata: &crate::metadata::MetaData, + latlon: crate::projector::LonLatCoord, +) -> Result { + let width_f64 = f64::from(metadata.width + 1); + let scale = f64::from(metadata.scale); + let coord_metric = crate::projector::Convert { + base: metadata.centre, + } + .to_meters(latlon)?; + let offset = (width_f64 * scale) / 2.0f64; + #[expect( + clippy::as_conversions, + clippy::cast_possible_truncation, + reason = "The coordinates should always fit in `u64`" + )] + let (x, y) = { + ( + ((coord_metric.x + offset) / scale) as i64, + ((-coord_metric.y + offset) / scale) as i64, + ) + }; + let dem_id = (y * i64::from(metadata.width)) + x; + Ok(dem_id) +} + +#[inline] +/// Convert a DEM ID to a lon/lat. +/// +/// # Errors +/// +/// If projection errors. +pub fn dem_id_to_lonlat( + metadata: &crate::metadata::MetaData, + dem_id: i64, +) -> Result { + let width_f64 = f64::from(metadata.width + 1); + let scale = f64::from(metadata.scale); + let offset = (width_f64 * scale) / 2.0f64; + + #[expect( + clippy::as_conversions, + clippy::cast_precision_loss, + reason = "These are just lon/lat coordinates" + )] + let (x_raster, y_raster) = { + ( + dem_id.rem_euclid(metadata.width.into()) as f64, + dem_id.div_euclid(metadata.width.into()) as f64, + ) + }; + #[expect( + clippy::suboptimal_flops, + reason = "We don't need the perfomance and this reads better" + )] + let coord_metric = geo::coord! { + x: (x_raster * scale) - offset, + y: (-y_raster * scale) + offset + }; + + let coord_lonlat = crate::projector::Convert { + base: metadata.centre, + } + .to_degrees(coord_metric)?; + + Ok(coord_lonlat) +}