In Part IV of this series, we fetched OpenStreetMap data and displayed it in the map. This time, we will dive into the world (no pun intended) of geometric geodesy. Our objective is to figure out what is the nearest way from our current location.

Photo by Suhash Villuri on Unsplash

Geodesy is the scientific discipline that deals with the measurement and representation of the Earth, and it was kick-started when the ancient Greek philosopher Eratosthenes estimated the circumference of the planet with great precision. Two thousand years later, most people still think the Earth is round!

We’ll be specifically dealing with spherical trigonometry, which is like ordinary trigonometry turned to 11. If you’ve ever had to find the distance between geographical places, you probably had the luxury of letting a GIS-enabled database or a GIS library do the grunt work for you. I don’t have either of that. Well, there are a few relevant Rust crates out there, but none that seem to do everything I need. Which is totally fine, because implementing this sounds like a fun challenge. And yes, I’ve considered contributing to the geo crate, but it needs to mature a bit in terms of coordinate system support. It currently doesn’t quite separate cartesian and polar coordinates, and an implementation of distance between objects will have to know which coordinate system it is dealing with.

At least I have access to the necessary mathematics. Wherever I look, people tend to refer to this great resource. Turns out it actually has all the answers I need.

Let’s start simple and tackle the distance between two points. There are several ways to do it, but the haversine formula seems to yield the best mix of short and long range accuracy as well as performance. Implemented in Rust it looks like this:

const R: f64 = 6371008.8; // mean Earth radiuspub fn distance(c1: &Coord, c2: &Coord) -> f64 {
    let (phi1, phi2) = (c1.phi(), c2.phi());
    let (lambda1, lambda2) = (c1.lambda(), c2.lambda());
    let delta_phi = phi2 - phi1;
    let delta_lambda = lambda2 - lambda1;
    let a = (delta_phi / 2.0).sin().powi(2)
        + phi1.cos() * phi2.cos() 
        * (delta_lambda / 2.0).sin().powi(2);    R * 2.0 * a.sqrt().atan2((1.0 - a).sqrt())
}impl Coord {
    fn phi(self: &Coord) -> f64 {
        self.lat.to_radians()
    }    fn lambda(self: &Coord) -> f64 {
        self.lon.to_radians()
    }
}

So far, so good, but how do we go about calculating the distance between a point and a line? It’s considerably harder to find a formula for that. We could cheat, of course. We could just calculate the distance to each of the nodes that make up our ways, and take the shortest distance. In most cases that would provide a good enough result, but not always. First, it would sometimes simply give the wrong answer. Consider the situation below. We are clearly nearest to the way on the left, but by only considering the nearest nodes we would be matching the way on the right.

Why does the app tell me I’m on that other road?

Second, and perhaps more problematic, after a fork in the road, the nearest node would be the fork itself, and it would be impossible to tell which way we have chosen to follow, as illustrated below.

Did we turn left or right at the fork? Hard to tell if we’re only considering the nearest node.

A smarter way to cheat would be to distribute many virtual nodes along each line, and then find the nearest one of all those. Not very efficient or elegant, but it would work. A reasonable plan B if calculating the nearest way directly proved impossible.

Luckily, it isn’t impossible. Consider the map below. I have drawn a line from Bergen to Trondheim. Now we want to find the distance from Ålesund to this line, and the same from Stavanger. While our mathematical guide doesn’t give us the answer directly, it gives us the tools to get there. One of them is the cross-track distance, which basically tells us how much off track we are. It corresponds with the green lines on the map. The second tool is the related along-track distance. It tells us how far from Bergen in the direction toward Trondheim do we have to travel in order to get to the green line. In the case of Ålesund, that’s roughly halfway to Trondheim. In the case of Stavanger, that’s backwards, which means that the along-track distance is negative.

I’ll spare you all the details, but the function calculating along-track distance looks like this:

pub fn along_track_distance(
    c1: &Coord, c2: &Coord, c3: &Coord) -> f64 {

let x = (angular_distance(c1, c3).sin() 
        * delta_theta(c1, c2, c3).sin()).asin();    if x.cos().abs() > f64::EPSILON {
        R * (angular_distance(c1, c3).cos() / x.cos())
            .acos()
            .copysign(delta_theta(c1, c2, c3).cos())
    } else {
        0.0
    }
}

In order to correctly report a negative distance, I have included a fix found in a Python library.

We can now find the nearest point on the line. If the along-track distance is negative (as is the case with Stavanger), the nearest point must be the beginning of the line (Bergen). Similarly, if the along-track distance is longer than the length of the line, the nearest point must be the end (Trondheim). In all other cases, start at the beginning, follow the bearing of the line for the appropriate distance, and you’ll end up at the nearest point.

pub fn nearest_point(c1: &Coord, c2: &Coord, c3: &Coord) -> Coord {
    let along_track_distance = along_track_distance(c1, c2, c3);if along_track_distance < 0.0 {
        c1.clone()
    } else if along_track_distance > distance(c1, c2) {
        c2.clone()
    } else {
        destination(c1, bearing(c1, c2), along_track_distance)
    }
}

Why the clone() calls? Most of the time in Rust, we are passing references around. Functions shouldn’t insist on taking ownership of values unless they have to. Well, this function returns a Coord rather than a reference to a Coord because destination() does as well. And destination() does it because it is constructing a brand new Coord. There simply isn’t a value anywhere that it can reference. Both destination() and nearest_point() need to pass on ownership of this new Coord, otherwise it would go out of scope and get destructed. But this means that nearest_point() cannot simply return c1 and c2. That would mean passing ownership of them, but we’re only borrowing them. We have to clone them so that we can return something we actually own.

It’s easy to unit test what we have so far, such as like this:

#[test]
fn test_nearest_point_aalesund() {
    let destination = destination(&BERGEN, 35.93, 212561.3);
    let nearest_point = nearest_point(
        &BERGEN, &TRONDHEIM, &AALESUND);    assert_approx_eq!(nearest_point.lat, destination.lat, 0.001);
    assert_approx_eq!(nearest_point.lon, destination.lon, 0.001);
}#[test]
fn test_nearest_point_stavanger() {
    let nearest_point = nearest_point(
        &BERGEN, &TRONDHEIM, &STAVANGER);    assert_approx_eq!(nearest_point.lat, BERGEN.lat, 0.001);
    assert_approx_eq!(nearest_point.lon, BERGEN.lon, 0.001);
}

With all the tools in place, it’s time to apply them. I would like to draw the shortest line from our current position to each of the ways on the map, and I’d also like to put a separate color on the nearest way. We’ll modify the rendering code:

pub fn render_topology(model: &Model) {
    if let Some(map) = &model.map {
        for way in model.osm.ways.iter() {
            Polyline::new_with_options(
                way.points(&model.osm) ...,
                &JsValue::from_serde(&PolylineOptions {
                    color: if model.nearest_way().id == way.id {
                        "blue"
                    } else {
                        "green"
                    }
                    .into(),
                    weight: 2,
                })
                .expect("Unable to serialize polyline options"),
            )
            .addTo(&map);
        }        if let Some(pos) = &model.position {
            for (destination, _, _) in   
                model.nearest_point_on_each_way().iter() {                Polyline::new_with_options(
                    vec![pos, destination] ...,
                    &JsValue::from_serde(&PolylineOptions {
                        color: "red".into(),
                        weight: 1,
                    })
                    .expect("Unable to serialize polyline options"),
                )
                .addTo(&map);
            }
        }
    }
}

This code relies on model.nearest_way(), as well as model.nearest_point_on_each_way(), a function that returns a list of triples, but we’re only interested in the first part of those, which is the Coord. It is implemented like this:

fn nearest_point_on_each_way(&self) -> Vec<(Coord, f64, &OsmWay)> {
    match &self.position {
        None => vec![],
        Some(pos) => self
            .osm
            .ways
            .iter()
            .map(|way| {
                way.points(&self.osm)
                    .windows(2)
                    .map(|line_segment| {
                        let a = line_segment[0];
                        let b = line_segment[1];
                        let destination =  
                            geo::nearest_point(
                                &a.into(), &b.into(), pos);

let distance = geo::distance(
                            pos, &destination);                        (destination, distance, way)
                    })
                    .min_by(|(_, x, _), (_, y, _)| {
                        x.partial_cmp(y).expect(
                            "Could not compare distances")
                    })
                    .expect("Could not find a nearest distance")
            })
            .collect(),
    }
}

So for each way, we use windows(2) on the node collection to get all line segments from one node to the next. For each line segment, we find our nearest point, and we find the distance to that point. This lets us construct a triple containing the nearest point, the distance, and the way it belongs to. For each way we take just the triple with the minimum distance.

The other function is easy to implement based on that:

fn nearest_way(&self) -> &OsmWay {
    let nearest_points = self.nearest_point_on_each_way();    let (_, _, way) = nearest_points
        .iter()
        .min_by(|(_, x, _), (_, y, _)|    
            x.partial_cmp(y).expect("Could not compare distances"))
        .expect("Could not find a nearest distance");    way
}

We simply find the triple with the minimum distance overall, and return its way.

The app now outputs the following nice little map, which lets us visually validate that everything seems to be in order. There is a red line running to the nearest point on each way. Sometimes the nearest point is a node, and sometimes it is within a line segment. And the nearest way overall is colored blue.

And with that, we’ve achieved our objective for this installment. Next time, we should probably start looking at movement, and updating the app view based on where we are.