Creating isochrones: what is the optimal way?
Go to the Streamlit app.
Accessibility
A common concept in Transportation and Urban Planning is the concept of Accessibility to services and amenities; not to be confused with Universal Accessibility which is more related to inclusive design. Examples of accessibility can be:
- What percentage of households are within a 5-minute drive to a pharmacy?
- A Transportation Master Plan objective aims to have 100% of households at a maximum of 500 m distance walk from a bus stop by redesigning their bus network in the next 2 years.
- How many people live 20 minutes, 30 minutes, 45 minutes, 1 hour, and more than an hour away from the university by cycling or public transportation?
These questions often translate to Accessibility or Reach Map analyses that aim to create a single — or a series of — isochrones (lines/areas of equal travel time or distance); similar to contour lines on topo maps (lines of equal elevation). The term isochrone is used interchangeably to refer to both time and distance. An advantage of isochrones is that they transform a routing problem into a spatial join problem. Instead of having to calculate routes between a source and a multitude of destinations (POIs) which can be time-consuming, you can instead intersect the isochrone (polygon/area of equal distance or travel time) with a POI layer represented by points or building polygons (spatial join) and know with good accuracy how far each POI lies from the start point for a given mode of travel.
Isochrones
However, An aspect of reach maps that is often overlooked is how the accessibility surface is actually computed or constructed. The final result is a polygon but the underlying calculations to determine the accessible area(s) are graph/network calculations (commonly ego graph or ego network). These calculations yield a subset of nodes (intersections) and links/edges (roads) that are within the required distance or travel time.
An additional step is needed to transform or approximate a shape (line or polygon) from those nodes and links. The correct choice of methodology or algorithm to approximate a shape can depend on the:
- Use case
— Are you calculating a single isochrone or do you need nested contour lines?
— Is the area of the isochrone important in the calculation? Maybe you have a KPI that is designed around an area metric and not a number or percentage of destinations reached - Location and network characteristics
— Urban vs rural areas
— Link and node density, connectivity, circuity, etc.
— Average width of streets, average perpendicular dept at which amenity/building access points are, etc.
There are two main groups of methodologies to approximate accessibility polygons:
- Computational Geometry methods (vector calculations)
- Spatial/Grid interpolation methods (raster calculations)
About Ego Graphs
An Ego Graph is a subgraph made up of nodes and their corresponding edges that are within a given “distance” from a starting node. It is purely a Graph Theory concept and the distance can be any type of weight. However, applied to a transportation or routing context, the ego graph can be limiting. Graphs are well suited to model a road network with streets and intersections but in real life, most trips do not start or end at intersections, they begin somewhere along the road; a driveway, a pedestrian access leading to a building, etc. Modeling every access point using nodes and edges would be ideal, but also difficult to maintain, and would probably cause performance issues at some point. Most public road network datasets will not have that level of detail and routing engines have built-in functionality to be able to merge new nodes temporarily to calculate routes without polluting the base network.
In the context of isochrones, the ego graph can be limiting in areas where we have big blocks with long stretches of road between intersections. For example, when computing the 500 m ego graph from a node, if the last accessible node is 400 m away and the node after that is 600 meters away, then the ego graph will end at the 400 m node and we will lose 100 m of accessibility (see image below).
To counter this problem, we can use interpolation to approximate where the 500 m would be and insert an artificial node and capture the real accessible distance on the edge. Thankfully, GeoNetworkX has a method to accomplish just that called the extended_ego_graph. The image below illustrates the difference between the ego graph and the extended ego graph in the vicinity of the Beaubien metro Station in Montreal where long rectangular blocks can be seen.
Isochrone Calculator Streamlit App
The Isochrone Calculator Streamlit app can calculate isochrones for anywhere in the world that is mapped in OpenStreetMap. The objective of the app is to allow the user to experiment with different methods of approximating the accessibility surface, understand the differences between them, tweak the parameters of each method, and determine which method works best for a specific type of location.
It uses OSMnx to download OSM data as networks and NetworkX, and GeoNetworkX for graph calculations. Geospatial data manipulation is done using GeoPandas and Rasterio (for vector and raster data respectively) and Pillow is used for image format manipulation. SciPy is used for the spatial/grid interpolation (to avoid using GDAL which would have been easier but installing GDAL is really a hit or miss each time, and I wanted to keep the package management simple for the sake of hosting on Streamlit). User interaction and visualization on web maps are achieved using Folium and PyDeck.
This article is meant to be a tutorial for the Streamlit app. It introduces the concepts used in the computation of isochrones and explains how different parameters can affect the results. As mentioned above, apart from algorithms and their parameters, the location itself can affect the result depending on the methodology used. I have selected 9 locations with different urban fabrics that will help illustrate how each method can affect the results and which isochrone shape approximation methods work best for each setting.
The Convex Hull
The Convex Hull is a polygon that encloses all the points in a set and has no concave corners. It can be thought of as a rubber band around the points. It is the simplest shape to calculate among the shapes considered in this article and app. Furthermore, all GIS software and geospatial programming packages will come with built-in functionality to create a convex hull given a set of points.
Advantages
- Easy to calculate
- Most if not all GIS software commercial or open-source will have built-in functionality to calculate a convex hull given a set of points
- Can be used to create isochrones at intervals (concentric lines/contour lines) instead of grid/spatial interpolation (covered later in the article) since the shapes are simple, smooth, and don’t have jagged edges or parts sticking out.
Disadvantages
- The shape will never have holes in it and consequently might not always represent the most correct area but this could also be an advantage depending on the situation or use case;
- Is greatly affected by the shape of the network; might miss some nodes and links (or parts of them) that were actually accessible or it might include some nodes and links that were not accessible. Works better with regular grid networks and might underperform in irregular pattern network settings.
The Concave Hull or Alpha Shape
The Concave Hull or the Alpha Shape is more difficult to compute because it requires user input or judgment. As opposed to the convex hull, the concave hull can have concave corners and theoretically can yield a more accurate shape along nodes, sort of like placing a flexible bag over the points and sucking out the air. However, not all sets of points will have the same density or spatial distribution, and applying the same parameters to every situation can yield different and undesirable results.
The alpha parameter determines how tight or concave the shape will be but another common method to calculate the concave hull is to compute the Delaunay Triangulation of the set of points and filter out the triangles (simplexes) based on their circumcircle radius. The Streamlit app borrows the function already available in GeoNetworkX which filters the triangles based on a percentile of the triangles’ radii.
As can be seen in the image above, using the 85th percentile for all the locations does not yield a good result all around. In fact, all of the sites have edges left out of the isochrone shape using the 85th percentile and some sites even have multiple disconnected polygons. It is always possible to artificially add nodes along each accessible edge to have a denser starting point set but that makes for another topic not covered in this article and app. This method can potentially result in shapes that have holes, but plugging the holes is also a question of judgment. The Streamlit app has a toggle for removing holes in the shapes.
One way to obtain a better concave shape for each case without having to inspect the results visually is to use recursion. The recursion and the conditions to break it would be:
If the shape is a multi-polygon (more than 1 polygon) or not all edges are within the resulting shape, then increase the percentile value by 1 and recompute the concave shape. Continue until one or both conditions are met or the percentile reaches 100.
Using the 100th percentile effectively yields the equivalent of the convex shape. The Streamlit app allows the user to impose the single polygon and the edges containment conditions at once or just one of them. The animation below gives an idea of what the Delaunay Triangulation looks like and how triangles are added or filtered out depending on the percentile value of two of the locations.
The image below shows the results of imposing both conditions on the Beaubien Metro and Manila sites. In the first case, for both conditions to be respected, we had to include a 100% of the triangles which is effectively the convex hull. In the second case, the conditions are met at the 95th percentile yielding a slightly more accurate shape compared to the convex hull.
Advantages
- Can yield more accurate shapes than the convex hull
- Shapes can have holes. If there is a vacant/undeveloped parcel of land, or a water body within the isochrone, having them represented by holes can be advantageous
Disadvantages
- Requires user input or judgment to select the best parameters
- Sometimes the best concave shape is the convex shape which results in a waste of time and computing resources
The Link Offset
The Link Offset method is probably the most precise method of the three but it’s also the most limited. It is more precise because it uses the accessible links and not the accessible nodes to build the shape, making it independent of network node density.
For example, looking at the Sherbrooke (rural area) case, some links will always be partially included or partially excluded from the isochrone when using methods based on nodes/points unless a large buffer is applied at the end of the calculation; whereas with the link offset method, it will always be included.
The link offset method is essentially applying a buffer to the accessible edges of the network. Like the concave hull, user input is required to determine the optimal offset/buffer distance. As opposed to the point methods which yield a relatively full shape with no or few holes in them, the link offset method will yield a polygon with potentially too many holes in it depending on the buffer distance. Judgment is required to determine if the holes should be filled or not but for the sake of this tutorial, we will not be plugging the holes to better illustrate the strengths and weaknesses of the method. In the Streamlit app, however, the user has the choice to fill in the holes.
Since the link offset depends on buffering lines and polylines, it will be affected by the join and endcap methods used. If your POI layer is represented by buildings, then you are likely to intersect the buildings when using an arbitrary buffer distance but if your POI layer is represented by points, then you need to make sure the buffer distance is large enough to contain the points.
If you use a high value for the buffer, depending on the endcap method, the shape might include more nodes or links that are actually inaccessible. Most GIS software or packages use round joins and endcaps by default; sometimes customization is available and sometimes it’s not.
The image below shows the link offset method applied using a 25 m buffer with round endcaps. At 25 meters, most of the buildings are intersected except for the courtyard design cases (Paris and Barcelona). The rural single-family home neighborhood of Sherbrooke also has a few cases where the buildings are barely intersected at 25 m. What if we do not have access to a buildings layer (as is the case with Mexico City where only half the buildings in the specific area are mapped in OSM)?
The image below focuses on Paris, Barcelona, and Sherbrooke. To simulate the absence of a buildings layer, we calculate the centroid of each building to create a fictional POI layer. Some of the centroids (blue dots) are left out of the isochrone shape.
Advantages
- Yields the most accurate shape because it uses the links and not the nodes
Disadvantages
- The correct buffering distance requires knowledge of the area or else some destinations might be left out
- The shape is likely to have many holes (but can be easily filled)
- The endcap buffering method will affect the results (high buffer value using round endcaps will include inaccessible areas at the extremities)
- Results will be drastically different depending on if the extended ego graph is used or the simple ego graph.
- Not well suited for computing contour lines
The Grid/Spatial Interpolation
The Spatial Interpolation method is mostly used for creating isochrones with long travel times or for creating contour lines (isochrones at regular intervals). This method is usually applied using a major trip generator like a university or a hospital as its source and it is not well suited to calculate accessibility areas for households for example (thousands of calculations).
There are many spatial interpolation algorithms but two of the more common ones are Inverse Distance Weighting (IDW) and Triangular Irregular Network (TIN) also known as Linear Interpolation. There is an abundance of information online about each of these methods so I will not go into the details but the important things to note are:
- The TIN algorithm is closely related to Delaunay Triangulation and to the Convex Hull and it cannot extrapolate/estimate data beyond the convex hull. The resulting raster image is rather simple and therefore is easy to extract pretty contours from.
- The IDW algorithm requires a decay parameter (how fast does the effect of a known sample point decrease with distance) which is set by the user, already making the IDW more demanding compared to the TIN. It yields much more elaborate surfaces compared to the TIN which makes it more accurate but also makes for more difficult contour extraction often requiring an additional step to filter, clean, and smooth out the contours. At this moment, I have not made the decay parameter tweakable in the Streamlit app.
As with any interpolation involving a regular grid, a grid cell size must be determined. In other words how wide do we want each pixel to be which in turn means how many meters on the ground will each pixel represent? The pixel size is usually a tradeoff between what we are trying to measure and how crisp of an image we want to get.
What are we measuring? Since we are measuring distances in a transportation/travel time context, data at 1 m accuracy is not necessary, we can probably say 20 m accuracy is good enough, but it might make for a fairly pixelated image.
How crisp of an image do we need? 20 m does not make a clear enough image and 1 m accuracy is not necessary, so we can opt for 5 meters resolution for example.
The user can experiment with the cell size in the Streamlit app and see how it affects the results.
Advantages
- The best choice for creating isochrones at intervals (contour lines) normally involving longer travel times or distances
Disadvantages
- The interpolation algorithm chosen can affect the result. Some algorithms work well with sparse or irregular data while others do not
What about different modes of travel?
All the examples in this article and the Streamlit app use the walking network to illustrate the concepts. For each mode of travel, we need to use its respective network.
For walking, the network would consist of all links where walking is permitted. Sometimes, the sidewalks would be represented as links of their own, parallel to the links representing roads, while other times an attribute of the road (link) would indicate if walking is permitted or not. For cycling it would be all links where cycling is permitted and not just bike lanes. In both these cases, the extended ego graph can be used.
In the case of car travel, the network would consist of links where driving is permitted. In this case, the extended ego graph should be used with caution because in the case of freeways, for example, trips cannot realistically end between two junctions (entrances or exits).
A similar constraint exists for public transportation. The network is the transit network and trips can only start and end at transit stops. Generally, walking to and from stops is included in the calculation, in which case we have a hybrid approach:
- First, we calculate the ego graph using the transit network which will yield a set of stops (nodes) that are within the target distance.
- Then, we calculate walking isochrones (or cycling, or another mode) from each accessible stop where the target distance for each mini isochrone is a function of the remaining total distance.
Conclusion
In conclusion, isochrones are an efficient way to measure accessibility and the procedure is similar for all modes of transportation as long as the correct network data is available for a given mode. For bi-modal or multi-modal accessibility, the procedure is not as straightforward but does share many similarities.
Any method covered in this article will yield a good enough result but understanding why and how different methods work can come in handy for edge cases. In addition, having the ability to calculate the extended ego graph (either through expert knowledge or through access to third-party software packages that integrate well with your existing ecosystems) minimizes the risk involved in choosing one method over another because the lack of detail in a network is no longer a problem.
In this tutorial and in the corresponding Streamlit app, we used the distance as cost or weight. It is also common to use travel time as weight and calculate isochrones in terms of time instead of distance. In this case, an average travel speed per mode can be applied to calculate the time required to travel each link in the network/graph. If congestion data or elevation data is available, it can be factored into the cost.
Lastly, judgment was mentioned a lot to determine the best step forward but it would also be possible to train a model that would predict the best method to use. Input to the model could be the nodes and links and their spatial properties (what was used in this article), network properties of the ego graph, and external data such as buildings, POIs, and their properties.
