Thanks to a few packages, it is super easy to analyze and visualize geospatial data in Python.

Table of Contents

1. GeoPandas
2. OSMnx
3. Create Figure
4. Examples

GeoPandas

The backbone holds everything together is GeoPandas. It extends the pandas DataFrame with a geometry column that contains a shapely object and turns it into a GeoDataFrame. These have a number of useful functions when it comes to working with geographic data.

OSMnx

The data has to come from somewhere and that is where OpenStreetMap comes in. The data there is maintained by the community and available to everyone. In Python, the OSMnx package makes it especially easy deal with the API and retrieve any desired data.

import numpy as np

# latitude and longitude of Heidelberg
center = (49.41165, 8.68805) 
# size in km
width, height = 4, 4

# calculate the coordinates on a sphere with radius=6371
lat,lon = center
dlat = height / (2 * 6371 * np.pi / 360) / 2
dlon = width / height * dlat / float(np.cos(np.radians(lat)))

# west, south, east, north
bbox = lon-dlon, lat-dlat, lon+dlon ,lat+dlat   

tags = {
    "building": True,
    "highway": True,
    "railway": True,
    "leisure": True,
    "landuse": True,
    "natural": True,
    "waterway": True,
}

We can then pass this to osmnx and download the data we want. Before saving it, we make some small adjustments.

import geopandas as gpd 
import osmnx as ox
from pathlib import Path

gdf = ox.features_from_bbox(bbox, tags)
gdf = gdf.rename_axis(index={"id": "osmid"})
# some objects like rivers extend beyond the original box
gdf = gpd.clip(gdf,bbox)

# to hide things that are under the ground
gdf["visible"] = ~(
    (gdf.get("tunnel") == "yes") |
    (gdf.get("location") == "underground") |
    (gdf.get("layer").astype("float", errors="ignore") < 0) |
    (gdf.get("covered") == "yes")
)

# we only save the columns we later need
gdf = gdf.reindex(columns=list(tags.keys())+["visible","geometry"])
filename = Path("data","OSM_Heidelberg2026.gpkg")
gdf.to_file(filename,index=True)

Create Figure

We can plot the data with gdf.plot(). In order to assign the different components individual colors, we plot them step by step in different layers. This was inspired by prettymaps.

import json
import matplotlib.pyplot as plt 
from shapely.geometry import box

# the layers are defined in this json file
with open(Path("styles",f"historic.json")) as f:
    layout = json.load(f)

fig,ax=plt.subplots(figsize=(10,12))

# draw background and frame
background = gpd.GeoDataFrame(geometry=[box(*gdf.total_bounds)],crs=gdf.crs)
ax = background.plot(ax=ax,color=layout["background"]["style"]["color"],zorder=-1)

for layer in layout:
    tags  = layout[layer]['tags']
    style = layout[layer]['style'].copy()

    # the objects used in this layer are selected with tags
    expr = []
    for key, value in tags.items():
        if value is True: 
            expr.append(f'{key}.notna()')
        elif value is False:
            expr.append(f'{key}.isna()')
        elif isinstance(value,str):
            expr.append(f'{key} == "{value}"')
        elif isinstance(value,list):
            expr.append(f'{key} in {value}')
        else:
            raise ValueError(f'the value for {key} must be bool|str|list')
    subsample = gdf.query('|'.join(expr))

    # buildings can use a random color, drawn from a palette
    if 'palette' in style:
        palette = style.pop('palette')
        style['fc'] = np.random.choice(palette, size=len(subsample))

    # skip empty GeoDataFrames
    if len(subsample) > 0:
        ax = subsample.plot(ax=ax,**style,aspect="auto")  
    else:
        print(f'no elements for {layer}')

xmin, ymin, xmax, ymax = gdf.total_bounds
ax.set(xlim=[xmin,xmax],ylim=[ymin,ymax])
ax.axis('off')

ax.text(0.995, 0.005, '© OpenStreetMap', transform=ax.transAxes, 
        ha='right', va='bottom',fontsize='small',color="black",
        bbox={'fc':'white','ec':'none', 'alpha':0.9, 'pad':1})

plt.show()

Examples

‹ ›