Skip to content

Commit

Permalink
[geoprof] allow f32
Browse files Browse the repository at this point in the history
  • Loading branch information
@JayKickliter
JayKickliter committed Oct 4, 2023
1 parent b1759ea commit 84cc992
Show file tree
Hide file tree
Showing 5 changed files with 175 additions and 58 deletions.
17 changes: 14 additions & 3 deletions Cargo.lock
View file

Some generated files are not rendered by default. Learn more about how customized files appear on GitHub.

1 change: 1 addition & 0 deletions Cargo.toml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -15,6 +15,7 @@ dashmap = "5.5.3"
env_logger = "0.10.0"
fast-math = "0.1.1"
geo = "0.26.0"
itertools = "0.11.0"
log = "0.4.20"
memmap2 = "0.7.1"
num-traits = "0.2.16"
Expand Down
2 changes: 2 additions & 0 deletions geopath/Cargo.toml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8,6 +8,8 @@ anyhow = "1"
clap = { workspace = true }
env_logger = { workspace = true }
geo = { workspace = true }
itertools = { workspace = true }
num-traits = { workspace = true }
rfprop = { git = "https://github.com/JayKickliter/Signal-Server", branch = "master" }
serde = { workspace = true }
serde_json = { workspace = true }
Expand Down
210 changes: 155 additions & 55 deletions geopath/src/main.rs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,11 +2,14 @@ mod options;

use anyhow::Error as AnyError;
use clap::Parser;
use itertools::Itertools;
use num_traits::{AsPrimitive, Float, FromPrimitive};
use options::{Cli, Command as CliCmd};
use rfprop::TerrainProfile;
use serde::Serialize;
use std::{io::Write, path::Path};
use terrain::{
geo::{point, Point},
geo::{coord, point, CoordFloat, Point},
Profile, TileMode, Tiles,
};
use textplots::{Chart, Plot, Shape};
Expand All @@ -16,6 +19,7 @@ fn main() -> Result<(), AnyError> {
let Cli {
tile_dir,
rfprop,
use_f32,
max_step,
earth_curve,
normalize,
Expand All @@ -26,38 +30,88 @@ fn main() -> Result<(), AnyError> {

env_logger::init();

let terrain_profile: CommonProfile = if rfprop {
rfprop::init(Path::new(&tile_dir), false)?;
rfprop::terrain_profile(
cli.start.0.y,
cli.start.0.x,
cli.start.1,
cli.dest.0.y,
cli.dest.0.x,
cli.dest.1,
900e6,
cli.normalize,
)
.into()
if use_f32 {
type C = f32;

let terrain_profile: CommonProfile<C> = if rfprop {
rfprop::init(Path::new(&tile_dir), false)?;
rfprop::terrain_profile(
cli.start.0.y,
cli.start.0.x,
cli.start.1,
cli.dest.0.y,
cli.dest.0.x,
cli.dest.1,
900e6,
cli.normalize,
)
.into()
} else {
let start_point = coord!(x: start.0.x as C, y: start.0.y as C);
let start_alt = start.1 as C;
let max_step = max_step as C;
let dest_point = coord!(x: dest.0.x as C, y: dest.0.y as C);
let dest_alt = dest.1 as C;

eprintln!(
"start_point: {start_point:?}, start_alt: {start_alt}, max_step: {max_step}, dest_point: {dest_point:?}, dest_alt: {dest_alt}"
);

let tile_src = Tiles::new(tile_dir, TileMode::MemMap)?;
Profile::<C>::builder()
.start(start_point)
.start_alt(start_alt)
.max_step(max_step)
.earth_curve(earth_curve)
.normalize(normalize)
.end(dest_point)
.end_alt(dest_alt)
.build(&tile_src)?
.into()
};

match cmd {
CliCmd::Csv => print_csv(terrain_profile),
CliCmd::Plot => plot_ascii(terrain_profile),
CliCmd::Json => print_json(terrain_profile),
CliCmd::Tia => print_tia(terrain_profile),
}
} else {
let tile_src = Tiles::new(tile_dir, TileMode::MemMap)?;
Profile::builder()
.start(start.0)
.start_alt(start.1)
.max_step(max_step)
.earth_curve(earth_curve)
.normalize(normalize)
.end(dest.0)
.end_alt(dest.1)
.build(&tile_src)?
type C = f64;

let terrain_profile: CommonProfile<C> = if rfprop {
rfprop::init(Path::new(&tile_dir), false)?;
rfprop::terrain_profile(
cli.start.0.y,
cli.start.0.x,
cli.start.1,
cli.dest.0.y,
cli.dest.0.x,
cli.dest.1,
900e6,
cli.normalize,
)
.into()
};
} else {
let tile_src = Tiles::new(tile_dir, TileMode::MemMap)?;
Profile::<C>::builder()
.start(coord!(x: start.0.x, y: start.0.y))
.start_alt(start.1)
.max_step(max_step)
.earth_curve(earth_curve)
.normalize(normalize)
.end(coord!(x: dest.0.x, y: dest.0.y))
.end_alt(dest.1)
.build(&tile_src)?
.into()
};

match cmd {
CliCmd::Csv => print_csv(terrain_profile),
CliCmd::Plot => plot_ascii(terrain_profile),
CliCmd::Json => print_json(terrain_profile),
CliCmd::Tia => print_tia(terrain_profile),
match cmd {
CliCmd::Csv => print_csv(terrain_profile),
CliCmd::Plot => plot_ascii(terrain_profile),
CliCmd::Json => print_json(terrain_profile),
CliCmd::Tia => print_tia(terrain_profile),
}
}
}

Expand All @@ -66,7 +120,7 @@ fn main() -> Result<(), AnyError> {
/// ```sh
/// cargo run -- --srtm-dir=data/nasadem/3arcsecond/ --max-step=90 --earth-curve --normalize --start=0,0,100 --dest=0,1,0 csv | tr ',' ' ' > ~/.tmp/plot && gnuplot -p -e "plot for [col=4:5] '~/.tmp/plot' using 1:col with lines"
/// ```
fn print_csv(profile: CommonProfile) -> Result<(), AnyError> {
fn print_csv<T: CoordFloat + std::fmt::Display>(profile: CommonProfile<T>) -> Result<(), AnyError> {
let mut stdout = std::io::stdout().lock();
writeln!(stdout, "Distance,Longitude,Latitude,Los,Elevation")?;
for (((elevation, point), los), distance) in profile
Expand All @@ -86,27 +140,33 @@ fn print_csv(profile: CommonProfile) -> Result<(), AnyError> {
Ok(())
}

fn plot_ascii(profile: CommonProfile) -> Result<(), AnyError> {
fn plot_ascii<T>(profile: CommonProfile<T>) -> Result<(), AnyError>
where
T: CoordFloat + AsPrimitive<f32>,
{
let plot_data: Vec<(f32, f32)> = profile
.terrain_elev_m
.iter()
.enumerate()
.map(|(idx, elev)| (f32::from(idx as u16), *elev as f32))
.map(|(idx, elev)| (f32::from(idx as u16), elev.as_()))
.collect();
Chart::new(300, 150, 0.0, plot_data.len() as f32)
.lineplot(&Shape::Lines(&plot_data))
.display();
Ok(())
}

fn print_json(profile: CommonProfile) -> Result<(), AnyError> {
fn print_json<T>(profile: CommonProfile<T>) -> Result<(), AnyError>
where
T: CoordFloat + Serialize,
{
#[derive(Serialize)]
struct JsonEntry {
location: [f64; 2],
elevation: f64,
struct JsonEntry<T> {
location: [T; 2],
elevation: T,
}

let reshaped: Vec<JsonEntry> = profile
let reshaped: Vec<JsonEntry<T>> = profile
.great_circle
.iter()
.zip(profile.terrain_elev_m.iter())
Expand All @@ -120,9 +180,12 @@ fn print_json(profile: CommonProfile) -> Result<(), AnyError> {
Ok(())
}

fn print_tia(profile: CommonProfile) -> Result<(), AnyError> {
fn print_tia<T>(profile: CommonProfile<T>) -> Result<(), AnyError>
where
T: CoordFloat + FromPrimitive + std::fmt::Display + std::iter::Sum,
{
let tia = terrain_intersection_area(
profile.distances_m.last().unwrap() / profile.distances_m.len() as f64,
&profile.distances_m,
&profile.los_elev_m,
&profile.terrain_elev_m,
);
Expand All @@ -132,31 +195,44 @@ fn print_tia(profile: CommonProfile) -> Result<(), AnyError> {

/// Calculate the positive area of intersection, in m², between the
/// profile (terrain) and the line of sight.
fn terrain_intersection_area(step_size_m: f64, los_vec: &[f64], profile: &[f64]) -> f64 {
los_vec
fn terrain_intersection_area<T>(distances_m: &[T], los_elev_m: &[T], terrain_elev_m: &[T]) -> T
where
T: Float + FromPrimitive + std::iter::Sum,
{
let tia_m2 = los_elev_m
.iter()
.zip(profile.iter())
.map(|(los, prof)| (prof - los).max(0.0) * step_size_m)
.sum::<f64>()
.zip(terrain_elev_m.iter())
.map(|(los, prof)| (*prof - *los).max(T::zero()))
.tuple_windows::<(T, T)>()
.zip(distances_m.iter().tuple_windows::<(&T, &T)>())
.map(|((h_n1, h_n), (d_n1, d_n))| {
let dx = (*d_n - *d_n1).abs();
dx * (h_n + h_n1) / T::from(2).unwrap()
})
.sum::<T>();

// Convert distance from `m^2` to `m*km` to stay compatible with
// DB assumptions.
tia_m2
}

/// A common represention of both native and rfprop profiles.
struct CommonProfile {
great_circle: Vec<Point<f64>>,
distances_m: Vec<f64>,
los_elev_m: Vec<f64>,
terrain_elev_m: Vec<f64>,
struct CommonProfile<T: CoordFloat> {
great_circle: Vec<Point<T>>,
distances_m: Vec<T>,
los_elev_m: Vec<T>,
terrain_elev_m: Vec<T>,
}

impl From<Profile<f64>> for CommonProfile {
impl<T: CoordFloat> From<Profile<T>> for CommonProfile<T> {
fn from(
Profile {
distances_m,
great_circle,
los_elev_m,
terrain_elev_m,
..
}: Profile<f64>,
}: Profile<T>,
) -> Self {
Self {
distances_m,
Expand All @@ -167,14 +243,38 @@ impl From<Profile<f64>> for CommonProfile {
}
}

impl From<rfprop::TerrainProfile> for CommonProfile {
impl From<TerrainProfile> for CommonProfile<f32> {
fn from(
TerrainProfile {
distance,
los,
terrain,
..
}: TerrainProfile,
) -> Self {
let distances_m = distance.iter().map(|val| *val as f32 * 1000.0).collect();
let los_elev_m = los.iter().map(|val| *val as f32).collect();
let terrain_elev_m = los.iter().map(|val| *val as f32).collect();
let great_circle = std::iter::repeat(point!(x: 0.0, y:0.0))
.take(terrain.len())
.collect();
Self {
distances_m,
great_circle,
los_elev_m,
terrain_elev_m,
}
}
}

impl From<TerrainProfile> for CommonProfile<f64> {
fn from(
rfprop::TerrainProfile {
TerrainProfile {
mut distance,
los,
terrain,
..
}: rfprop::TerrainProfile,
}: TerrainProfile,
) -> Self {
distance.iter_mut().for_each(|val| *val *= 1000.0);
let great_circle = std::iter::repeat(point!(x: 0.0, y:0.0))
Expand Down
3 changes: 3 additions & 0 deletions geopath/src/options.rs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -13,6 +13,9 @@ pub struct Cli {
#[arg(long, default_value_t = false)]
pub rfprop: bool,

#[arg(long = "f32", default_value_t = false)]
pub use_f32: bool,

/// Maximum path incremental step size, in meters.
#[arg(short, long, default_value_t = 90.0)]
pub max_step: f64,
Expand Down

0 comments on commit 84cc992

Please sign in to comment.