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

use std::fs::File;

use anyhow::Result;
use byteorder::{BigEndian, ReadBytesExt};

use geom::{Distance, LonLat};

// TODO Make sure this actually works, and link to references describing the format.
// - It's VERY wrong! Look at 19th and Boyer, and Delmar.
// TODO Use http://gis.ess.washington.edu/data/raster/tenmeter/byquad/seattle/index.html or
// something else instead?

// Assuming the 3-second arcs
const GRID_DIM: usize = 1201;

pub struct Elevation {
    lon_offset: f64,
    lat_offset: f64,
    data: Vec<i16>,
}

impl Elevation {
    pub fn load(path: &str) -> Result<Elevation> {
        println!("Reading elevation data from {}", path);
        let mut f = File::open(path)?;

        let mut e = Elevation {
            // TODO dont hardcode
            lon_offset: -122.0,
            lat_offset: 47.0,
            data: Vec::with_capacity(GRID_DIM.pow(2)),
        };
        // TODO off by one?
        for _ in 0..GRID_DIM.pow(2) {
            e.data.push(f.read_i16::<BigEndian>()?);
        }
        Ok(e)
    }

    pub fn get(&self, pt: LonLat) -> Distance {
        // TODO assert the (lon, lat) match the offsets
        // TODO not tons of confidence in any of this.
        // TODO interpolate from the 4 matching tiles?
        let x = ((pt.x() - self.lon_offset).abs() * (GRID_DIM as f64)) as usize;
        let y = ((pt.y() - self.lat_offset).abs() * (GRID_DIM as f64)) as usize;
        let i = x + (y * GRID_DIM);
        if i >= self.data.len() {
            panic!("srtm lookup for out-of-bounds {}", pt);
        }
        Distance::meters(f64::from(self.data[i]))
    }
}