Disallow non-centred tiles and remove --scale argument

README.md

@@ -31,7 +31,7 @@ The primary use case for this project is to find the longest line of sight on th
 Any file format supported by `gdal` should be supported. It must follow this requirements:
 * It must be perfectly square.
 * The width must be divisible by 48.
-* It should be in an AEQD or similar metric projection whose anchor is the centre of the tile.
+* It must be in an AEQD or similar metric projection whose anchor is the centre of the tile.
 * All points must be the same metric distance apart.
 
 The best source of elevation data I have found is here: https://www.viewfinderpanoramas.org/Coverage%20map%20viewfinderpanoramas_org3.htm It's mostly in the `.hgt` format, but can easily be converted to `.bt` using `gdal_translate`.
@@ -57,9 +57,10 @@ Once that's run, you will have access to what's called "ring data", stored by de
   Usage: tvs <COMMAND>
 
   Commands:
-    compute   Run main computations
-    viewshed  Reconstruct a viewshed
-    help      Print this message or the help of the given subcommand(s)
+    compute       Run main computations
+    viewshed      Reconstruct a viewshed
+    post-process  Post-processing tasks
+    help          Print this message or the help of the given subcommand(s)
 
   Options:
     -h, --help
@@ -98,7 +99,7 @@ Once that's run, you will have access to what's called "ring data", stored by de
             - cpu:        Optimised cache-efficient CPU kernel
             - cuda:       TBC
 
-            [default: vulkan]
+            [default: cpu]
 
         --output-dir <Directory to save output to>
             Directory to save results in
@@ -134,30 +135,65 @@ Once that's run, you will have access to what's called "ring data", stored by de
 
             [default: 0.13]
 
-        --thread-count <thread count>
+        --thread-count <Thread count>
             Thread count used for CPU parallelism
 
             [default: 8]
 
+        --disable-image-render
+            Controls line of sight and total viewshed image generation
+
+        --viewsheds-db-path <viewshed storage path>
+            Where to store the viewshed data. Requires build with `--features=ring_data`
+
+            [default: ./viewsheds.db]
+
+        --aoi-point <Lon/lat coords of region>
+            Lon/lat coordinates for a polygon that represents an Area of Interest within the DEM. Points outside this polygon will be ignored
+
+        --database-per-thread
+            Creates a database per thread, managed by an additional worker thread
+
     -h, --help
             Print help (see a summary with '-h')
 
 
 * Viewshed
   Reconstruct a viewshed
 
-  Usage: tvs viewshed <Path to existing output directory> [COORDINATES]...
+  Usage: tvs viewshed <Path to existing database> <Viewshed output directory> [COORDINATES]...
 
   Arguments:
-    <Path to existing output directory>
-            Directory where compute output was saved
+    <Path to existing database>
+            Path or directory where viewshed DB was saved
+
+    <Viewshed output directory>
+            Directory to save viewshed in
 
     [COORDINATES]...
             Coordinates to reconstruct viewsheds for
 
   Options:
     -h, --help
             Print help
+
+
+* Post-process
+  Post-processing tasks
+
+  Usage: tvs post-process [OPTIONS] <Path to existing databases>
+
+  Arguments:
+    <Path to existing databases>
+            Directory where viewshed DBs were saved
+
+  Options:
+        --threads <Number of threads to use>
+            Number of threads to use. Defaults to letting Rayon choose the thread count
+
+    -h, --help
+            Print help
+
 ```
 
 ## Building Vulkan shader

benchmarks/run.sh

@@ -32,8 +32,6 @@ rm $sqlite_db_path || true
 # Do the calculations
 time cargo run --features ring_data --release -- \
 	compute "$PROJECT_ROOT/benchmarks/cardiff.tiff" \
-	--scale 100 \
-	--centre-from-projection \
 	--rings-per-km 3 \
 	--backend "$backend" \
 	--process all \

crates/total-viewsheds/src/config.rs

@@ -72,10 +72,6 @@ pub struct Compute {
     )]
     pub output_dir: std::path::PathBuf,
 
-    /// Override the calculated DEM points scale from the DEM file. Units in meters.
-    #[arg(long, value_name = "DEM scale (meters)")]
-    pub scale: Option<f32>,
-
     /// What to compute.
     #[arg(
         long,
@@ -116,15 +112,6 @@ pub struct Compute {
     #[arg(long, value_name = "Render image", default_value = "false")]
     pub disable_image_render: bool,
 
-    /// Derive the tile's centre from the tile's anchored projection. This can be more accurate for
-    /// large metric-projected tiles.
-    #[arg(
-        long,
-        value_name = "Get centre from projection",
-        default_value = "false"
-    )]
-    pub centre_from_projection: bool,
-
     /// Where to store the viewshed data. Requires build with `--features=ring_data`.
     #[arg(
         long,
@@ -142,6 +129,7 @@ pub struct Compute {
         value_name = "Lon/lat coords of region")
     ]
     pub aoi_point: Vec<(f32, f32)>,
+
     /// Creates a database per thread, managed by an additional worker thread
     #[arg(
         long,

crates/total-viewsheds/src/cpu/unrolled_los.rs

@@ -87,7 +87,8 @@ impl<const UNROLL: usize, const VECTOR_WIDTH: usize> UnrolledVectorLos<UNROLL, V
         assert_eq!(
             max_los % VECTOR_WIDTH,
             0,
-            "the maximum line of sight must be divisible by {VECTOR_WIDTH} for vectorization"
+            "the maximum line of sight ({max_los}) \
+            must be divisible by {VECTOR_WIDTH} for vectorization"
         );
 
         let (distances, adjustments) = generate_distances(max_los, refraction, scale);

crates/total-viewsheds/src/dem.rs

@@ -43,22 +43,14 @@ impl DEM {
         centre_latlon: crate::projection::LonLatCoord,
         width: u32,
         scale: f32,
-        max_line_of_sight: u32,
+        max_line_of_sight_as_points: u32,
     ) -> Result<Self> {
         let size: u64 = u64::from(width) * u64::from(width);
-        #[expect(
-            clippy::cast_possible_truncation,
-            clippy::as_conversions,
-            clippy::cast_sign_loss,
-            clippy::cast_precision_loss,
-            reason = "This shouldn't be a problem in most sane cases"
-        )]
-        let max_los_as_points = (max_line_of_sight as f32 / scale) as u32;
         let max_possible_los_as_points = width.div_euclid(2);
 
-        if max_los_as_points > max_possible_los_as_points {
+        if max_line_of_sight_as_points > max_possible_los_as_points {
             color_eyre::eyre::bail!(
-                "The maximum line of sight ({max_line_of_sight}m) is longer than the maximum \
+                "The maximum line of sight ({max_line_of_sight_as_points}) is longer than the maximum \
                 distance that any point can completely see ({}m).",
                 f64::from(max_possible_los_as_points) * f64::from(scale)
             );
@@ -78,23 +70,12 @@ impl DEM {
             size,
             scale,
             centre: centre_latlon,
-            max_los_as_points,
+            max_los_as_points: max_line_of_sight_as_points,
             computable_points_count,
         };
         Ok(dem)
     }
 
-    #[expect(
-        clippy::as_conversions,
-        clippy::cast_possible_truncation,
-        clippy::cast_sign_loss,
-        reason = "The scale is positive and never be beyond the limts of f32"
-    )]
-    /// Cast the scale into an integer
-    pub const fn scale_u32(&self) -> u32 {
-        self.scale as u32
-    }
-
     /// Is the DEM ID within the DEM such that a valid viewshed can be made for it?
     pub fn is_point_computable(&self, dem_id: u32) -> bool {
         let coord = self.convert_dem_id_to_coord(dem_id).0;

crates/total-viewsheds/src/main.rs

@@ -140,7 +140,19 @@ fn compute(config: &config::Compute) -> Result<()> {
     }
 
     let mut tile = tile::Tile::load(config)?;
-    let scale = config.scale.unwrap_or(tile.scale);
+
+    let max_line_of_sight_as_points = tile.width.div_euclid(3);
+
+    let mut dem = crate::dem::DEM::new(
+        tile.centre,
+        tile.width,
+        tile.scale,
+        max_line_of_sight_as_points,
+    )?;
+
+    dem.elevations = mem::take(&mut tile.data);
+
+    tracing::debug!("Created DEM: {dem:?}");
 
     #[expect(
         clippy::as_conversions,
@@ -149,30 +161,18 @@ fn compute(config: &config::Compute) -> Result<()> {
         clippy::cast_precision_loss,
         reason = "Sign loss and truncation aren't relevant"
     )]
-    let max_line_of_sight_metres = (tile.width.div_euclid(3) as f32 * scale) as u32;
-
-    let mut dem = crate::dem::DEM::new(tile.centre, tile.width, scale, max_line_of_sight_metres)?;
-
-    dem.elevations = mem::take(&mut tile.data);
+    let max_line_of_sight_as_metres = (max_line_of_sight_as_points as f32 * tile.scale) as u32;
 
     // Free up RAM
     drop(tile);
 
-    tracing::debug!("Created DEM: {dem:?}");
-
-    let max_line_of_sight = if matches!(config.backend, config::Backend::CPU) {
-        dem.max_los_as_points
-    } else {
-        max_line_of_sight_metres
-    };
-
     let dem_metadata = crate::cpu::storage::metadata::MetaData {
         width: dem.width,
         scale: dem.scale,
-        max_line_of_sight,
+        max_line_of_sight: max_line_of_sight_as_points,
         reserved_ring_size: run::compute::Compute::ring_count_per_band(
             config.rings_per_km,
-            dem.max_los_as_points * dem.scale_u32(),
+            max_line_of_sight_as_metres,
         ),
         centre: dem.centre,
     };

crates/total-viewsheds/src/run/compute.rs

@@ -194,7 +194,7 @@ impl<'compute> Compute<'compute> {
         Ok(crate::cpu::storage::metadata::MetaData {
             width: self.dem.width,
             scale: self.dem.scale,
-            max_line_of_sight: self.dem.max_los_as_points * self.dem.scale_u32(),
+            max_line_of_sight: self.dem.max_los_as_points,
             reserved_ring_size: usize::try_from(self.constants.reserved_rings_per_band)?,
             centre: self.dem.centre,
         })

crates/total-viewsheds/src/tile.rs

@@ -6,8 +6,6 @@
 
 use color_eyre::eyre::{ContextCompat as _, Result};
 
-use crate::projection;
-
 /// The basic raw data needed to compute total viewsheds.
 pub struct Tile {
     /// The width of the tile.
@@ -43,11 +41,7 @@ impl Tile {
             );
         }
 
-        let centre = if config.centre_from_projection {
-            Self::get_centre_by_projection(&dataset)?
-        } else {
-            Self::get_centre_by_raster(&dataset)?
-        };
+        let centre = Self::get_centre_by_projection(&dataset)?;
 
         let lon = centre.0.x;
         let lat = centre.0.y;
@@ -103,10 +97,8 @@ impl Tile {
         }
     }
 
-    /// Get the lat/lon of the centre of the tile simply by dividing the extent in half. This isn't
-    /// ideal as there's no guarantee that it's the same centre that the creator used to generate
-    /// the tile.
-    fn get_centre_by_raster(dataset: &gdal::Dataset) -> Result<crate::projection::LonLatCoord> {
+    /// Get the metric centre of the tile simply by dividing the extent in half.
+    fn get_metric_centre(dataset: &gdal::Dataset) -> Result<geo::Coord> {
         let (width, height) = dataset.raster_size();
 
         #[expect(
@@ -126,24 +118,27 @@ impl Tile {
             x_top_left + (centre_x * pixel_width),
             y_top_left - (centre_y * pixel_height),
         );
-        let crs = &dataset.spatial_ref()?.to_proj4()?;
-
-        let mut converted = (x_world, y_world, 0.0f64);
-        proj4rs::transform::transform(
-            &proj4rs::Proj::from_proj_string(crs)?,
-            &projection::Converter::degrees_projection()?,
-            &mut converted,
-        )?;
-
-        Ok(crate::projection::LonLatCoord(geo::coord! {
-            x: converted.0.to_degrees(),
-            y: converted.1.to_degrees()
-        }))
+
+        Ok(geo::coord! {
+            x: x_world,
+            y: y_world
+        })
     }
 
-    /// Get the lat/lon centre of the tile by querying the projection's definition. This is more
-    /// likely to guarantee that the tile's centre matches the creator's intended centre.
+    /// Get the lat/lon centre of the tile by querying the projection's definition.
     fn get_centre_by_projection(dataset: &gdal::Dataset) -> Result<crate::projection::LonLatCoord> {
+        const ALLOWED_DEVIATION: f64 = 0.001; // 1 milimetre
+        let geometric_centre = Self::get_metric_centre(dataset)?;
+        if geometric_centre.x.abs() > ALLOWED_DEVIATION
+            || geometric_centre.y.abs() > ALLOWED_DEVIATION
+        {
+            color_eyre::eyre::bail!(
+                "Tile centre ({},{}) must be within {ALLOWED_DEVIATION} of 0,0.",
+                geometric_centre.x,
+                geometric_centre.y
+            );
+        }
+
         let projection = &dataset.spatial_ref()?.to_proj4()?;
 
         let lat_0: f64 = projection
@@ -174,7 +169,6 @@ mod test {
     fn get_centre_by_projection() {
         let config = crate::config::Compute {
             input: "../../benchmarks/samples/aeqd_10x10.tiff".into(),
-            centre_from_projection: true,
             ..Default::default()
         };
         let tile = Tile::load(&config).unwrap();
@@ -188,18 +182,18 @@ mod test {
     }
 
     #[test]
-    fn get_centre_by_raster() {
-        let config = crate::config::Compute {
-            input: "../../benchmarks/samples/aeqd_10x10.tiff".into(),
-            ..Default::default()
-        };
-        let tile = Tile::load(&config).unwrap();
-        assert_eq!(
-            tile.centre,
-            crate::projection::LonLatCoord(geo::Coord {
-                x: -0.1278,
-                y: 51.5074
-            })
-        );
+    fn get_metric_centre() {
+        let dataset = gdal::Dataset::open("../../benchmarks/samples/aeqd_10x10.tiff").unwrap();
+        let centre = Tile::get_metric_centre(&dataset).unwrap();
+        assert_eq!(centre, geo::Coord { x: 0.0, y: 0.0 });
+    }
+
+    #[test]
+    fn bad_centre() {
+        let dataset = gdal::Dataset::open("../../benchmarks/samples/aeqd_bad_centre.tiff").unwrap();
+        let error = Tile::get_centre_by_projection(&dataset)
+            .unwrap_err()
+            .to_string();
+        assert!(error.contains("Tile centre"), "Wrong error: '{error}'");
     }
 }