1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153

use std::fs::File;
use std::io::{BufRead, BufReader, BufWriter, Write};
use std::process::Command;

use anyhow::Result;

use abstutil::{must_run_cmd, Timer};
use geom::{Distance, PolyLine};
use map_model::raw::{OriginalRoad, RawMap};

pub fn add_data(map: &mut RawMap, timer: &mut Timer) -> Result<()> {
    timer.start("add elevation data");

    timer.start("generate input");
    let ids = generate_input(map)?;
    timer.stop("generate input");

    timer.start("run elevation_lookups");
    std::fs::create_dir_all("elevation_output")?;
    std::fs::create_dir_all("data/input/shared/elevation")?;
    let pwd = std::env::current_dir()?.display().to_string();
    must_run_cmd(
        // Because elevation_lookups has so many dependencies, just depend on Docker.
        Command::new("docker")
            .arg("run")
            // Bind the input directory to the temporary place we just created
            .arg("--mount")
            .arg(format!(
                "type=bind,source={}/elevation_input,target=/elevation/input,readonly",
                pwd
            ))
            // We want to cache the elevation data sources in A/B Street's S3 bucket, so bind to
            // our data/input/shared directory.
            .arg("--mount")
            .arg(format!(
                "type=bind,source={}/data/input/shared/elevation,target=/elevation/data",
                pwd
            ))
            .arg("--mount")
            .arg(format!(
                "type=bind,source={}/elevation_output,target=/elevation/output",
                pwd
            ))
            .arg("-t")
            // TODO Upload this to Docker Hub, so it's easier to distribute
            .arg("elevation_lookups_lidar")
            .arg("python3")
            .arg("main.py")
            .arg("query"),
    );
    timer.stop("run elevation_lookups");

    timer.start("grab output");
    scrape_output(map, ids)?;
    timer.stop("grab output");

    // Clean up temporary files
    std::fs::remove_file("elevation_input/query")?;
    std::fs::remove_dir("elevation_input")?;
    std::fs::remove_file("elevation_output/query")?;
    std::fs::remove_dir("elevation_output")?;

    timer.stop("add elevation data");
    Ok(())
}

fn generate_input(map: &RawMap) -> Result<Vec<OriginalRoad>> {
    std::fs::create_dir_all("elevation_input")?;
    let mut f = BufWriter::new(File::create("elevation_input/query")?);
    let mut ids = Vec::new();
    for (id, r) in &map.roads {
        // TODO Handle cul-de-sacs
        if let Ok(pl) = PolyLine::new(r.center_points.clone()) {
            ids.push(id.clone());
            // Sample points every meter along the road
            let mut pts = Vec::new();
            let mut dist = Distance::ZERO;
            while dist <= pl.length() {
                let (pt, _) = pl.dist_along(dist).unwrap();
                pts.push(pt);
                // Smaller gives more detail, but is slower.
                dist += Distance::meters(5.0);
            }
            // Always ask for the intersection
            if *pts.last().unwrap() != pl.last_pt() {
                pts.push(pl.last_pt());
            }
            for (idx, gps) in map.gps_bounds.convert_back(&pts).into_iter().enumerate() {
                write!(f, "{},{}", gps.x(), gps.y())?;
                if idx != pts.len() - 1 {
                    write!(f, " ")?;
                }
            }
            writeln!(f)?;
        }
    }
    Ok(ids)
}

fn scrape_output(map: &mut RawMap, ids: Vec<OriginalRoad>) -> Result<()> {
    let num_ids = ids.len();
    let mut cnt = 0;
    for (line, id) in BufReader::new(File::open("elevation_output/query")?)
        .lines()
        .zip(ids)
    {
        cnt += 1;
        let line = line?;
        let mut values = Vec::new();
        for x in line.split('\t') {
            let x = x.parse::<f64>()?;
            if !x.is_finite() {
                // TODO Warn
                continue;
            }
            values.push(Distance::meters(x));
        }
        if values.len() != 4 {
            error!("Elevation output line \"{}\" doesn't have 4 numbers", line);
            continue;
        }
        // TODO Also put total_climb and total_descent on the roads
        map.intersections.get_mut(&id.i1).unwrap().elevation = values[0];
        map.intersections.get_mut(&id.i2).unwrap().elevation = values[1];
    }
    if cnt != num_ids {
        bail!("Output had {} lines, but we made {} queries", cnt, num_ids);
    }

    // Calculate the incline for each road here, before the road gets trimmed for intersection
    // geometry. If we did this after trimming, we'd miss some of the horizontal distance.
    for (id, road) in &mut map.roads {
        let rise = map.intersections[&id.i2].elevation - map.intersections[&id.i1].elevation;
        let run = PolyLine::unchecked_new(road.center_points.clone()).length();
        if !(rise / run).is_finite() {
            // TODO Warn?
            continue;
        }
        road.percent_incline = rise / run;
        // Per https://wiki.openstreetmap.org/wiki/Key:incline#Common_.26_extreme_inclines, we
        // shouldn't often see values outside a certain range. Adjust this when we import
        // somewhere exceeding this...
        if road.percent_incline.abs() > 0.3 {
            error!(
                "{} is unexpectedly steep! Incline is {}%",
                id,
                road.percent_incline * 100.0
            );
        }
    }

    Ok(())
}