Creating Beautiful Isochrone Maps with Python

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 Process

The map generation process involves several key steps:

  1. Loading and processing elevation data
  2. Calculating travel times
  3. Combining travel times with elevation data
  4. Adding city labels and street networks

Let’s go through each step in detail.

1. Elevation Data Processing

First, we load and process the Digital Elevation Model (DEM) data. This involves several transformations:

Original DEM 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:

Transformed DEM 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:

Masked DEM 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:

Hillshaded DEM 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

2. Travel Time Calculation

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.

Travel Times 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:

Dissolved Travel Times 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")

3. Combining Data

The most interesting part is combining the travel time data with the elevation data to create a beautiful visualization:

Combined Shape and DEM 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

4. Final Map

The final map includes:

Final Map 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')

Technical Implementation

The implementation uses several Python libraries:

The code is structured in a modular way, with separate classes for:

Key technical aspects:

  1. Coordinate Systems: The code handles multiple coordinate systems:
  2. Data Processing:
  3. Visualization:

Conclusion

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.

References