Most navigation services will consider several variables when looking at how safe a particular road is to ride on. Is the road a highway? How many lanes does the road have? What’s the speed limit? These are all important variables to consider, but one that’s often overlooked is the geometry of the road itself. Sharper turns and corners are statistically less safe. We’ll dive in to our accident data to find out just how much curvature affects the safety of a road.

What is Curvature?

The first question is, what is curvature and how do you measure it? That may sound like a silly question, but the answer is not so obvious. Imagine you’re riding in a car at constant speed with your eyes closed. How would you tell if you’re going around a curve? Well, you’d feel it! We’re all familiar with the acceleration we feel going around an off-ramp, and that’s exactly how we define curvature. For the more math inclined this means all we need to do is parameterize our curve by arc length and then take the second derivative.

Most GIS systems represent roads as a series of point locations, or what geometers would call a polygonal approximation. How do we go from that to a smooth curve we can take derivatives of? Bezier curves! Bezier curves are widely used in computer science to interpolate a smooth curve from a set of points. Our curve needs to meet two very important properties. First, our interpolated curve should pass exactly through the points specified in our polygonal approximation (sometimes called deBoor points). And secondly, our curve should be C², that is the second derivative needs to exist everywhere and it should be continuous. To achieve this we’ll use an interpolation method known as B-splines.

So we have a method for measuring curvature at a point, but we’d like to extend this to measuring how curvy a road is overall. This is a concept known as total curvature. In essence we’d like to “add up” the curvature at each point along the curve. How do we do that? We simply take the integral of the curvature with respect to arc length.

Putting it into Practice

Now that we have a good handle on how to measure curvature, let’s see how it’s done with some real data. For this part we’ll be showing you step-by-step how to measure curvature using OSM data and R. We’ll be using a Debian based system throughout.

We’ll be working with OpenStreetMap data, in particular this extract of metro Atlanta. To load the data we’ll use osm2pgsql. To install it, simply run:

apt-get install osm2pgsql

And you should have osm2pgsql installed as well as Postgres and PostGIS. Next create a spatially enabled database:

postgres=# CREATE USER youruser;
postgres=# CREATE DATABSE osm;
postgres=# \connect osm
postgres=# CREATE EXTENSION POSTGIS;

Now that we have a spatially enabled database we can import our data with osm2pgsql:

osm2pgsql --slim -U youruser -d osm -C 16000 --number-processes 16 atlanta.osm.pbf

Make sure you adjust the -C parameter to match the amount of RAM your system has, likewise make sure to adjust --number-processes to match the number of cores you have. Now go grab a cup of coffee, it’s going to be a couple of hours (or more). Once it’s loaded, fire up R and go ahead and install the packages we’ll need:

install.packages(c("pracma", "sf", "sp"))

This first step is to load all the road geometries we just imported into Postgres, luckily sf will handle this for us:

lines <- st_read("PG:dbname=osm", "planet_osm_roads")

Next we’ll go ahead and define our curvature function:

dot <- function (x, y) as.numeric(x%*%y)
norm <- function(x) sqrt(sum(x^2))
curvature <- function(t, fx, fy) {
    # t is an argument to spline functions fx and fy.
    xp <- fx(t,1); yp <- fy(t,1)
    xpp <- fx(t,2); ypp <- fy(t,2)
    v <- c(xp, yp)
    a <- c(xpp, ypp)
    ((norm(v)^2 * norm(a)^2) - dot(a,v)^2) / norm(v)^3
}

Here, fx and fy will be our parameterized splines, v and a are the first and second derivatives (velocity and acceleration respectively), and k is the actual curvature calculation. So we’ve got a nice formula, but how can we verify that it’s correct? One of the more amazing facts in differential geometry is that the total curvature of any closed curve will always be a multiple of 2π, in particular a circle has total curvature 2π. Let’s test this:

int <- seq(0, 2*pi, length.out = 1000)
fx <- splinefun(int, cos(int), method="natural")
fy <- splinefun(int, sin(int), method="natural")
circle <- function(t) curvature(t, fx, fy)
integrate(Vectorize(circle), 0 , 2*pi)

This gives us a result of 6.279552 with absolute error < 0.00072, exactly what we’d expect. Now we’ll make a handy convenience function to compute curvature from our sfc objects created with st_read:

k <- function(geom) {
    coords <- st_coordinates(geom$way)
    n <- dim(coords)[1]
    fx <- splinefun(1:n, coords[,1], method="natural")
    fy <- splinefun(1:n, coords[,2], method="natural")
    kappa <- function(t) abs(curvature(t, fx, fy))
    result <- tryCatch(integrate(Vectorize(kappa), 0, n, subdivisions = 1e4L)$value, error=function(err) 0)
    c(geom$osm_id, result)
}

You’ll notice we’re taking the absolute value of curvature, this is because we want to be able differentiate the curviness of any closed paths, otherwise they’d all have similar curvature. Now we’ll apply this to all of the roads we imported and grab the road with the greatest absolute curvature:

computed <- t(apply(lines[,], 1, function(x) k(x)))
colnames(computed) <- c("osm_id", "curvature")
computed <- as.data.frame(computed)
computed[which.max(computed$curvature),]

Go grab another cup of coffee and maybe watch a movie, when you get back you’ll see that our maximum curvature is 37.0988, which corresponds to osm_id: 396809075. Take a look below, it is indeed pretty curvy. In our next post we’ll use the techniques we discussed today to explore our accident data and see just how much curvature affects safety.