Various cleanups after the new elevation code

docs/credits.md

@@ -18,6 +18,8 @@ Special thanks to [Hadrien Salat](https://github.com/HSalat) for the logo and
 [Stuart Lynn](https://github.com/stuartlynn) for help with pmtiles, both from
 [the Alan Turing Institute](https://www.turing.ac.uk).
 
+Thanks to [KDKasonde](https://github.com/KDKasonde) for adding elevation support.
+
 ## References / inspiration:
 
 - [Propensity to Cycle Tool](https://www.pct.bike) / [NPTScot](https://nptscot.github.io)

docs/output.md

@@ -12,7 +12,8 @@ The main mode of the tool outputs a GeoJSON FeatureCollection, with each LineStr
 - `way` is the OSM way ID of the road
 - `node1` and `node2` are the OSM node IDs bounding this road segment. Intermediate nodes of a curvy way (of degree 2, with no other connecting roads) are not used.
 - `count` represents the sum of trips along the segment. This is equal to the number of trips crossing the segment when the uptake model is "Identity", and something weighted for other uptake models.
-- `cost` is the cost for crossing this segment for routing
+- `forward_cost` and `backward_cost` are the costs for crossing this segment in each direction for routing
+- `slope` is the slope as a percent (3% grade encoded as `3.0`) in the forwards direction
 - `lts` is the Level of Traffic Stress for the segment, based on the chosen configuration. `0` means not allowed, `1` is suitable for children, and `4` is high stress.
 - `nearby_amenities` is the number of shops and amenities that're closest to this segment.
 

docs/tutorial_examples.md

@@ -123,12 +123,15 @@ First let's tell od2net to use a custom command for cost. Edit `config.json` and
 },
 ```
 
-Then create a new file called `cost.py`, copying from a different [example cost.py](https://github.com/Urban-Analytics-Technology-Platform/od2net/blob/main/examples/edinburgh/example_cost.py). The program gets a JSON array of dictionaries in STDIN and needs to print a JSON array of integer numbers as a result. Each dictionary input gives you:
+Then create a new file called `cost.py`, copying from a different [example cost.py](https://github.com/Urban-Analytics-Technology-Platform/od2net/blob/main/examples/edinburgh/example_cost.py). The program gets a JSON array of dictionaries in STDIN and needs to print a JSON array of pairs of integer numbers as a result. Each dictionary input gives you:
 
 - `length_meters`
 - `tags`, a JSON dictionary with the raw OSM tags
 - `lts` as a number 0 to 4, with 0 representing "cyclists not allowed here"
 - `nearby_amenities`, the number of shops that're closest to this road
+- Optional `slope`, the percent grade in the forwards direction
+
+The output is a corresponding JSON array with a pair of costs for each edge. The pair is `[forward_cost, backward_cost]`, which may be the same.
 
 The example does something very boring -- if [highway = residential](https://wiki.openstreetmap.org/wiki/Tag:highway%3Dresidential), just return length. Otherwise, multiply by 10, meaning all other roads will be **heavily** penalized. Write something more interesting here!
 

examples/edinburgh/example_cost.py

@@ -7,14 +7,15 @@ def calculate(edge):
     length_meters = edge["length_meters"]
     lts = edge["lts"]
     nearby_amenities = edge["nearby_amenities"]
+    slope = edge["slope"]
 
     # Return None to not use the edge at all
 
     if tags["highway"] == "residential":
-        return round(length_meters)
+        return [round(length_meters), round(length_meters)]
     else:
         # Strongly avoid non-residential roads
-        return round(10 * length_meters)
+        return [round(10 * length_meters), round(10 * length_meters)]
 
 
 # Read an array of JSON dictionaries from STDIN

examples/york/cost.py

@@ -7,14 +7,15 @@ def calculate(edge):
     length_meters = edge["length_meters"]
     lts = edge["lts"]
     nearby_amenities = edge["nearby_amenities"]
+    slope = edge["slope"]
 
     # Return None to not use the edge at all
 
     if tags["highway"] == "residential":
-        return round(length_meters)
+        return [round(length_meters), round(length_meters)]
     else:
         # Strongly avoid non-residential roads
-        return round(10 * length_meters)
+        return [round(10 * length_meters), round(10 * length_meters)]
 
 
 # Read an array of JSON dictionaries from STDIN

od2net/src/config.rs

@@ -16,7 +16,9 @@ pub struct InputConfig {
 
     pub lts: LtsMapping,
 
-    pub dem: Option<String>,
+    /// Path to a GeoTIFF file with elevation data. It must use WGS84 coordinates and have heights
+    /// in units of meters.
+    pub elevation_geotiff: Option<String>,
 }
 
 #[derive(Clone, Serialize, Deserialize)]

od2net/src/main.rs

@@ -81,9 +81,8 @@ fn main() -> Result<()> {
                 };
 
                 println!("That failed ({err}), so generating it from {osm_path}");
-                let dem_geotiff = if let Some(ref dem) = config.dem {
-                    let dem_file_path = format!("{directory}/input/{}", dem);
-                    Some(fs_err::read(&dem_file_path)?)
+                let geotiff_bytes = if let Some(ref filename) = config.elevation_geotiff {
+                    Some(fs_err::read(&format!("{directory}/input/{filename}"))?)
                 } else {
                     None
                 };
@@ -93,7 +92,7 @@ fn main() -> Result<()> {
                     &config.lts,
                     &mut config.cost,
                     &mut timer,
-                    dem_geotiff,
+                    geotiff_bytes,
                 )?;
 
                 timer.start(format!("Saving to {bin_path}"));

od2net/src/network/create_from_osm.rs

@@ -23,7 +23,7 @@ impl Network {
         lts: &LtsMapping,
         cost: &mut CostFunction,
         timer: &mut Timer,
-        dem_input_buffer: Option<Vec<u8>>,
+        geotiff_bytes: Option<Vec<u8>>,
     ) -> Result<Network> {
         timer.start("Make Network from xml or pbf");
         timer.start("Scrape OSM data");
@@ -81,10 +81,14 @@ impl Network {
         }
         timer.stop();
 
-        if let Some(buffer) = dem_input_buffer {
+        if let Some(bytes) = geotiff_bytes {
             timer.start("Calculate elevation for all edges");
-            let mut geo_tiff = GeoTiffElevation::new(Cursor::new(buffer));
-            network.calculate_elevation(&mut geo_tiff)?;
+            let mut geotiff = GeoTiffElevation::new(Cursor::new(bytes));
+            let progress = utils::progress_bar_for_count(network.edges.len());
+            for (_, edge) in &mut network.edges {
+                progress.inc(1);
+                edge.apply_elevation(&mut geotiff);
+            }
             timer.stop();
         }
 
@@ -117,23 +121,6 @@ impl Network {
 
         Ok(())
     }
-
-    pub fn calculate_elevation(
-        &mut self,
-        elevation_data: &mut GeoTiffElevation<Cursor<Vec<u8>>>,
-    ) -> Result<()> {
-        let progress = utils::progress_bar_for_count(self.edges.len());
-        for (_, edge) in &mut self.edges {
-            let elevation_details = edge.apply_elevation(elevation_data);
-            if let Some((slope, slope_factor)) = elevation_details {
-                progress.inc(1);
-                edge.slope = Some(slope);
-                edge.slope_factor = Some(slope_factor);
-            };
-        }
-
-        Ok(())
-    }
 }
 
 struct Way {

od2net/src/network/mod.rs

@@ -118,10 +118,9 @@ pub struct Edge {
     pub way_id: WayID,
     pub tags: Tags,
     geometry: Vec<Position>,
-    // slope as a percentage, for example a 3% slope is represented as 3.0.
+    // A 3% slope is represented as 3.0.
     pub slope: Option<f64>,
-    // slope factor is the value we will multiply the cost by to account for the
-    // slope of a given edge. The factor is given for traversing the edge in both directions.
+    // A factor to multiply cost by in the (forwards, backwards) direction
     pub slope_factor: Option<(f64, f64)>,
     // Storing the derived field is negligible for file size
     pub length_meters: f64,
@@ -137,87 +136,80 @@ pub struct Edge {
 }
 
 impl Edge {
-    pub fn apply_elevation<R: Read + Seek + Send>(
-        &self,
-        elevation_data: &mut GeoTiffElevation<R>,
-    ) -> Option<(f64, (f64, f64))> {
-        let slope = self.get_slope(elevation_data)?;
-
-        let length = self.length_meters;
-
-        let forward_slope_factor = Edge::calculate_slope_factor(slope, length);
-        let backward_slope_factor = Edge::calculate_slope_factor(-slope, length);
-
-        Some((slope, (forward_slope_factor, backward_slope_factor)))
-    }
-
-    /// This function takes in a slope and length and will calculate a slope factor
-    /// an explanation of the logic used can be found here:  https://github.com/U-Shift/Declives-RedeViaria/blob/main/SpeedSlopeFactor/SpeedSlopeFactor.md#speed-slope-factor-1
-    /// instead of using the slope_factor to divide the speed of a rider, we instead use it
-    /// multiplicatively on the cost to augment it before routing
-    fn calculate_slope_factor(slope: f64, length: f64) -> f64 {
-        let g = if 13.0 >= slope && slope > 10.0 && length > 15.0 {
-            4.0
-        } else if slope < 8.0 && slope <= 10.0 && length > 30.0 {
-            4.5
-        } else if slope < 5.0 && slope <= 8.0 && length > 60.0 {
-            5.0
-        } else if slope < 3.0 && slope <= 5.0 && length > 120.0 {
-            6.0
-        } else {
-            7.0
+    /// Sets `slope` and `slope_factor`.
+    pub fn apply_elevation<R: Read + Seek + Send>(&mut self, geotiff: &mut GeoTiffElevation<R>) {
+        let Some(slope) = self.get_slope(geotiff) else {
+            // Silently return if this fails
+            return;
         };
-
-        let slope_factor = if slope < -30.0 {
-            1.5
-        } else if slope < 0.0 && slope >= -30.0 {
-            1.0 + 2.0 * 0.7 * slope / 13.0 + 0.7 * slope * slope / 13.0 / 13.0
-        } else if slope <= 20.0 && slope >= 0.0 {
-            1.0 + slope * slope / g / g
-        } else {
-            10.0
-        };
-
-        slope_factor
+        self.slope = Some(slope);
+        self.slope_factor = Some((
+            calculate_slope_factor(slope, self.length_meters),
+            calculate_slope_factor(-slope, self.length_meters),
+        ));
     }
 
-    fn get_slope<R: Read + Seek + Send>(
-        &self,
-        elevation_data: &mut GeoTiffElevation<R>,
-    ) -> Option<f64> {
-        let first_node = self.geometry[0].to_degrees();
-        let second_node = self.geometry[self.geometry.len() - 1].to_degrees();
-
-        let first_node_height =
-            elevation_data.get_height_for_lon_lat(first_node.0 as f32, first_node.1 as f32)?;
+    fn get_slope<R: Read + Seek + Send>(&self, geotiff: &mut GeoTiffElevation<R>) -> Option<f64> {
+        let (lon1, lat1) = self.geometry[0].to_degrees();
+        let (lon2, lat2) = self.geometry.last().unwrap().to_degrees();
 
-        let second_node_height =
-            elevation_data.get_height_for_lon_lat(second_node.0 as f32, second_node.1 as f32)?;
+        let height1 = geotiff.get_height_for_lon_lat(lon1 as f32, lat1 as f32)?;
+        let height2 = geotiff.get_height_for_lon_lat(lon2 as f32, lat2 as f32)?;
 
-        let slope = (second_node_height - first_node_height) / self.length_meters as f32 * 100.0;
+        let slope = (height2 - height1) / (self.length_meters as f32) * 100.0;
         Some(slope.into())
     }
 }
 
+/// This returns a factor to multiply cost by, to adjust the speed of a cyclist. See
+/// <https://github.com/U-Shift/Declives-RedeViaria/blob/main/SpeedSlopeFactor/SpeedSlopeFactor.md#speed-slope-factor-1>.
+fn calculate_slope_factor(slope: f64, length: f64) -> f64 {
+    // Ported from https://github.com/U-Shift/Declives-RedeViaria/blob/5b5680ba769ab57f0fe061fd16c626cec66a0452/SpeedSlopeFactor/SpeedSlopeFactor.Rmd#L114
+    let g = if slope > 3.0 && slope <= 5.0 && length > 120.0 {
+        6.0
+    } else if slope > 5.0 && slope <= 8.0 && length > 60.0 {
+        5.0
+    } else if slope > 8.0 && slope <= 10.0 && length > 30.0 {
+        4.5
+    } else if slope > 10.0 && slope <= 13.0 && length > 15.0 {
+        4.0
+    } else {
+        7.0
+    };
+
+    // TODO Check this one again
+    let slope_factor = if slope < -30.0 {
+        1.5
+    } else if slope < 0.0 && slope >= -30.0 {
+        1.0 + 2.0 * 0.7 * slope / 13.0 + 0.7 * slope * slope / 13.0 / 13.0
+    } else if slope <= 20.0 && slope >= 0.0 {
+        1.0 + slope * slope / g / g
+    } else {
+        10.0
+    };
+
+    slope_factor
+}
+
 #[cfg(test)]
 mod tests {
-    use super::Edge;
+    use super::*;
 
     #[test]
     fn speed_slope_test() {
         let speed_flat = 15.0;
         let slope = 3.0;
         let length = 50.0;
-        let slope_factor = Edge::calculate_slope_factor(slope, length);
+        let slope_factor = calculate_slope_factor(slope, length);
         let slope_speed = speed_flat / slope_factor;
         let delta = slope_speed - 12.67241;
-        assert!(delta < f64::EPSILON);
+        assert!(delta < 0.00001);
 
         let slope = -8.0;
         let length = 100.0;
-        let slope_factor = Edge::calculate_slope_factor(slope, length);
+        let slope_factor = calculate_slope_factor(slope, length);
         let slope_speed = speed_flat / slope_factor;
         let delta = slope_speed - 37.17009;
-        assert!(delta < f64::EPSILON);
+        assert!(delta < 0.00001);
     }
 }

od2net/src/network/output.rs

@@ -62,29 +62,21 @@ impl Edge {
         properties.insert("way".to_string(), JsonValue::from(self.way_id.0));
         properties.insert("node1".to_string(), JsonValue::from(node1.0));
         properties.insert("node2".to_string(), JsonValue::from(node2.0));
-        if let (Some(forward_cost), Some(backward_cost)) = (self.forward_cost, self.backward_cost) {
-            // Either both costs are defined or none of them are so the expression above should
-            // be enough.
+        if let Some(forward_cost) = self.forward_cost {
             properties.insert(
                 "forward_cost".to_string(),
                 serde_json::to_value(forward_cost).unwrap(),
             );
+        }
+        if let Some(backward_cost) = self.backward_cost {
             properties.insert(
                 "backward_cost".to_string(),
                 serde_json::to_value(backward_cost).unwrap(),
             );
-            properties.insert(
-                "length".to_string(),
-                serde_json::to_value(self.length_meters).unwrap(),
-            );
-
-            // put the slope in here since if the cost is None there is no need to pull in the
-            // slope.
-            if let Some(slope) = self.slope {
-                properties.insert("slope".to_string(), serde_json::to_value(slope).unwrap());
-            };
         }
-
+        if let Some(slope) = self.slope {
+            properties.insert("slope".to_string(), serde_json::to_value(slope).unwrap());
+        };
         properties.insert("lts".to_string(), serde_json::to_value(self.lts).unwrap());
         properties.insert(
             "nearby_amenities".to_string(),