In this post, I’ll walk you through the process of creating beautiful isochrone maps using Python. An isochrone map shows areas that can be reached within a certain travel time from a central point. We’ll combine this with elevation data to create visually stunning maps that show both travel times and topography.
The map generation process involves several key steps:
Let’s go through each step in detail.
First, we load and process the Digital Elevation Model (DEM) data. This involves several transformations:
The original Digital Elevation Model (DEM) data
The DEM processing is handled by the _process_dem
method in our Isochrone
class:
def _process_dem(self):
"""Process DEM data including transformation, masking and hillshading."""
if self.kwargs.use_elevation:
bounds = self.get_shp_bounds()
elevation_file = self.get_elevation_dem(bounds)
else:
elevation_file = self.kwargs.dem_file
with rio.open(elevation_file) as src:
# Transform DEM to target CRS
working_dem = transform_dem(src, self.config["dem"]["dest_crs"])
# Resize if needed
if self.kwargs.reduce_size != 1.0:
working_dem = resize_dem(working_dem, self.kwargs.reduce_size)
# Mask to area of interest
working_dem = mask_dem(working_dem, self.admin_shp["geometry"])
# Apply hillshading
transformed_shaded = hillshade(working_dem.read(1),
azimuth=315,
angle_altitude=45)
return recreate_dataset_reader(transformed_shaded,
working_dem.meta,
1)
The DEM is then transformed to match our target coordinate system:
The DEM after coordinate system transformation
The transformation is handled by the transform_dem
function:
def transform_dem(dem_dataset, target_crs):
"""Transform DEM to target coordinate system."""
if dem_dataset.crs == target_crs:
return dem_dataset
# Calculate new transform and dimensions
transform, width, height = calculate_default_transform(
dem_dataset.crs, target_crs,
dem_dataset.width, dem_dataset.height,
*dem_dataset.bounds
)
# Create new dataset with transformed data
kwargs = dem_dataset.meta.copy()
kwargs.update({
'crs': target_crs,
'transform': transform,
'width': width,
'height': height
})
return recreate_dataset_reader(dem_dataset.read(), kwargs)
Next, we mask the DEM to focus on our area of interest:
The DEM after masking to our target area
The masking process uses the mask_dem
function:
def mask_dem(dem_dataset, geometries):
"""Mask DEM using provided geometries."""
out_image, out_transform = rio_mask(
dem_dataset,
geometries,
crop=True,
all_touched=True
)
out_meta = dem_dataset.meta.copy()
out_meta.update({
"height": out_image.shape[1],
"width": out_image.shape[2],
"transform": out_transform
})
return recreate_dataset_reader(out_image, out_meta)
Finally, we create a hillshade effect to better visualize the terrain:
The DEM with hillshading applied for better terrain visualization
The hillshading is calculated using the hillshade
function:
def hillshade(array, azimuth=315, angle_altitude=45):
"""Calculate hillshade for the DEM."""
x, y = np.gradient(array)
slope = np.pi/2 - np.arctan(np.sqrt(x*x + y*y))
aspect = np.arctan2(-x, y)
azimuthrad = azimuth * np.pi / 180
altituderad = angle_altitude * np.pi / 180
hillshade = (np.sin(altituderad) * np.sin(slope) +
np.cos(altituderad) * np.cos(slope) *
np.cos((azimuthrad - np.pi/2) - aspect))
return hillshade * 255
We calculate travel times from a central point to all areas in our region. This is done using the OSRM (Open Source Routing Machine) service, which provides accurate routing information.
Initial travel time calculations for each area
The travel time calculation is handled by the OSRM
class:
class OSRM:
def __init__(self, base_url, config, kwargs):
self._base_url = base_url
self._origin_coords = (config["origin_city"]["lat"],
config["origin_city"]["lon"])
self._config = config
self._kwargs = kwargs
def calculate_travel_time(self, gdf):
"""Calculate travel times for all areas."""
# Initialize fields
gdf = self._init_gdf_fields(gdf)
# Calculate for each area
for _, row in tqdm(gdf.iterrows(),
total=len(gdf),
desc="Calculating travel times"):
geoid = row["DEBKG_ID"]
# Convert UTM to lat/lon
utm_x, utm_y = row["geometry"].centroid.coords[0]
lat, lon = utm.to_latlon(utm_x, utm_y, 32, 'U')
# Get route information
result = self._get_distance_duration(lat, lon)
if result is None:
continue
distance, duration = result
gdf = self._update_gdf(gdf, geoid, lat, lon,
distance, duration)
return gdf
The travel times are then dissolved into larger areas for better visualization:
Travel times dissolved into larger areas
The dissolution process is handled by the dissolve_gdf
method:
@staticmethod
def dissolve_gdf(gdf, based_on_column):
"""Dissolve GeoDataFrame based on column."""
gdf["grouping"] = gdf[based_on_column]
return gdf.dissolve(by="grouping")
The most interesting part is combining the travel time data with the elevation data to create a beautiful visualization:
The final combined visualization showing both travel times and elevation
The combination process uses the combine_shape_dem
method:
def combine_shape_dem(self, admin_gdf, dem):
"""Combine shape data with DEM for visualization."""
unique_colors = len(admin_gdf["traveltime_h"].unique())
cmap = LinearSegmentedColormap.from_list(
"Custom",
[(255/255, 215/255, 0/255), (0/255, 87/255, 183/255)],
N=unique_colors
)
base_image = np.full((dem.shape[0], dem.shape[1], 3),
np.array([0, 0, 0]))
# Process each area
for _, area in tqdm(admin_gdf.iterrows(),
total=len(admin_gdf),
desc="Processing areas"):
out_image, _ = rio_mask(dem, [area["geometry"]],
crop=False,
all_touched=False)
# Calculate color based on travel time
color_index, ccode = self._get_colorcode_val(
area["traveltime_h"],
area["traveltime_h30"],
unique_colors
)
# Apply color to area
rgba = cmap(ccode)
rgb = np.array(rgba[:3])
base_color = np.full((dem.shape[0], dem.shape[1], 3), rgb)
image_data = np.expand_dims(out_image[0, :, :], axis=2)
tinted_base_color = np.rint(base_color * image_data).astype("int")
base_image = base_image + tinted_base_color
return base_image
The final map includes:
The final isochrone map with all elements combined
The final visualization is created using the plot_map
function:
def plot_map(colored_image, street_graph_edges, city_label_gdf,
eu_gdf, config, legend_used_colors):
"""Create the final map visualization."""
fig, ax = plt.subplots(figsize=config["map"]["size"])
# Plot base image
ax.imshow(colored_image)
# Add street network
if street_graph_edges is not None:
street_graph_edges.plot(ax=ax,
linewidth=0.5,
alpha=0.5)
# Add city labels
for idx, row in city_label_gdf.iterrows():
ax.annotate(row['city'],
(row.geometry.x, row.geometry.y),
fontsize=config["map"]["fontsize"])
# Add legend
plot_legend(ax, legend_used_colors)
plt.savefig(config["map"]["output_file"],
dpi=300,
bbox_inches='tight')
The implementation uses several Python libraries:
geopandas
for handling geographic datarasterio
for processing elevation datamatplotlib
for visualizationosmnx
for street network dataOSRM
for travel time calculationsThe code is structured in a modular way, with separate classes for:
Isochrone
: Main class handling the map generation processOSRM
: Client for travel time calculationsKey technical aspects:
Creating isochrone maps is a powerful way to visualize accessibility and travel times. By combining this with elevation data and proper visualization techniques, we can create informative and beautiful maps that help understand the relationship between travel times and terrain.
The code is available on GitHub, and you can use it to create your own isochrone maps for any region of interest.