Megacity Air Quality Feature Pipeline

Build a spatial SQL pipeline engineering air quality features for 10 megacities.

Introduction

30 Second Summary

Air quality can vary dramatically from one street to the next. A monitor next to a highway breathes in a completely different cocktail of pollutants than one surrounded by parks.

In this project, you will build a geospatial feature engineering pipeline that pulls real air quality station data from 10 of the world's largest megacities and computes spatial features around each monitor using SQL.

What You'll Build

A Jupyter Notebook pipeline that ingests air quality station locations from OpenAQ, downloads geospatial features from OpenStreetMap using OSMnx, loads everything into DuckDB with its spatial extension, and writes spatial SQL queries to compute features like road density, green space ratio, and distance to highways around each monitor.

By the end of this project, you'll have:

  • A DuckDB database with spatial tables containing air quality stations, road networks, green spaces, and industrial zones for 10 global megacities.
  • A set of spatial SQL queries using functions like ST_Buffer, ST_Distance, and ST_Intersects to compute meaningful environmental features.
  • An ML-ready CSV dataset with computed spatial features for every air quality monitor across all 10 cities, ready to feed into a predictive model.
  • Secret Mission: Add temporal air quality patterns by computing seasonal pollution features per city.

Are there any prerequisites?

You should be comfortable with basic Python (variables, functions, loops), basic SQL (SELECT, WHERE, JOIN), and have some familiarity with pandas DataFrames. No prior GIS or geospatial experience is needed.

Before We Start

Predicting air quality requires more than just pollution readings. It requires understanding the physical environment around each monitoring station.

In this step, you will install the Python packages needed for geospatial analysis, verify that DuckDB's spatial extension works, and get a free API key from OpenAQ to access global air quality data.

In this step, get ready to:
  • Install Python packages for geospatial data engineering.
  • Verify DuckDB's spatial extension can create and query geometry objects.
  • Sign up for a free OpenAQ API key to access air quality station data.
Install Python packages

This project uses several specialized libraries. GeoPandas handles geographic data in DataFrames. OSMnx downloads map features from OpenStreetMap. DuckDB gives you a fast SQL engine with a spatial extension for geometry operations.

  • Press Cmd+Space (macOS) or the Windows key (Windows) to open your search bar.
  • Type Visual Studio Code and press Enter to open it.
  • Open the integrated terminal by pressing Ctrl+` (macOS and Windows).
  • Move to your Desktop by running this command:
cd ~/Desktop
  • Create your project folder and move into it by running these commands:
mkdir megacity-air-quality
cd megacity-air-quality
  • Confirm the folder exists by running:
ls

You should see the megacity-air-quality folder listed in the output (or an empty listing if you ran ls from inside the new folder).

  • Install all the project dependencies by running this command:
pip install duckdb geopandas osmnx requests folium shapely jupyter

What does this command do?

  • duckdb is a fast, in-process SQL database engine with a spatial extension for geometry operations.
  • geopandas extends pandas DataFrames to support geographic data (points, lines, polygons).
  • osmnx downloads geospatial features (roads, parks, buildings) directly from OpenStreetMap.
  • requests makes HTTP requests to the OpenAQ API.
  • folium creates interactive maps for visual validation.
  • shapely provides geometric operations that GeoPandas relies on under the hood.
  • jupyter enables running Jupyter Notebooks inside VS Code.

You should see a series of Successfully installed messages once the installation completes.

Seeing permission errors?

If you see a permissions error, try running the command with the --user flag: pip install --user duckdb geopandas osmnx requests folium shapely jupyter.

If pip is not found, make sure Python is correctly installed and available on your system PATH.

help me troubleshoot pip install errors on my system

Verify the DuckDB spatial extension

DuckDB's spatial extension adds geometry types and Spatial SQL functions (like ST_Buffer, ST_Distance, ST_Intersects) to DuckDB. You need to confirm it installs and loads correctly before building the full pipeline.

  • Open the Command Palette by pressing Cmd+Shift+P (macOS) or Ctrl+Shift+P (Windows).
  • Type Create: New Jupyter Notebook and select it from the dropdown.
  • Save the notebook immediately using Cmd+S (macOS) or Ctrl+S (Windows).
  • Name the file air_quality_pipeline.ipynb and save it inside your megacity-air-quality folder.

What is a Jupyter Notebook?

A Jupyter Notebook is an interactive document that combines code, output, and notes in a single file. You write code in "cells" and run them one at a time. This makes it ideal for data pipelines where you want to inspect results at each stage.

  • If VS Code prompts you to select a kernel, click Select Kernel in the top-right of the notebook and choose your Python environment.
  • In the first cell of your notebook, paste the following code and run it by clicking the play button to the left of the cell (or pressing Shift+Enter):
import duckdb

# Connect to DuckDB in-memory
con = duckdb.connect()

# Install and load the spatial extension
con.execute("INSTALL spatial;")
con.execute("LOAD spatial;")

# Test with a point geometry (New Delhi coordinates)
result = con.execute("SELECT ST_Point(77.2245, 28.6358) AS test_point;").fetchone()
print(result)

What does this code do?

  • duckdb.connect() creates an in-memory database. No files on disk yet.
  • INSTALL spatial downloads the spatial extension. LOAD spatial activates it for this session.
  • ST_Point(77.2245, 28.6358) creates a geometry point at New Delhi's coordinates (longitude, latitude). If this returns a geometry object, the extension is working.

You should see output containing POINT (77.2245 28.6358) printed below the cell.

Not seeing a geometry point?

If you see an error about the spatial extension not being found, make sure you installed duckdb via pip in the same Python environment that VS Code is using as the notebook kernel.

Check the kernel shown in the top-right of your notebook. If it does not match the environment where you ran pip install, click it and switch to the correct Python interpreter.

help me debug why DuckDB spatial extension is not loading in my Jupyter notebook

Get your OpenAQ API key

OpenAQ provides free access to global air quality data from thousands of monitoring stations. You need an API key to authenticate your requests.

  • Open https://explore.openaq.org in your browser.
  • Click Sign Up and create a free account.
  • After signing in, navigate to your account settings to find your API key.
  • Copy the API key. You will paste it in the next step.

Why do I need an API key?

The OpenAQ API uses the X-API-Key header to authenticate requests. The free tier is generous and more than enough for this project.

  • In your notebook, add a new cell below the first one by clicking the + Code button at the top of the notebook.
  • Store your API key by pasting this code into the new cell:
# Your OpenAQ API key (free tier)
API_KEY = "your-api-key-here"
  • Replace your-api-key-here with the API key you copied from your OpenAQ account.
  • Run the cell by pressing Shift+Enter.

What does this code do?

This stores your API key in a Python variable called API_KEY. Later cells will use this variable to authenticate requests to the OpenAQ API. The cell runs without visible output, which is expected.

  • Save your notebook with Cmd+S (macOS) or Ctrl+S (Windows).

Forgot where to find your API key?

Go back to https://explore.openaq.org, sign in, and check your account settings or profile page. The key should be listed there.

help me find my OpenAQ API key after signing up

✔️ Awesome, I've got everything!

Great. Double check you have saved air_quality_pipeline.ipynb with both cells (the spatial test and the API key variable).

ⓧ I'd like to double check the full code

Here is what your entire notebook should look like at this point (two cells):

Cell 1:

import duckdb

# Connect to DuckDB in-memory
con = duckdb.connect()

# Install and load the spatial extension
con.execute("INSTALL spatial;")
con.execute("LOAD spatial;")

# Test with a point geometry (New Delhi coordinates)
result = con.execute("SELECT ST_Point(77.2245, 28.6358) AS test_point;").fetchone()
print(result)

Cell 2:

# Your OpenAQ API key (free tier)
API_KEY = "your-api-key-here"

Your environment is ready. Next up, you will pull real air quality station data from 10 of the world's largest megacities using the OpenAQ API.

Ingest Air Quality Station Data from OpenAQ

Your pipeline needs real-world air quality monitoring data to work with. The OpenAQ platform provides free, open access to air quality station locations across the globe.

In this step, you will query the OpenAQ API for PM2.5 monitoring stations in the world's 10 largest megacities, load the results into a DuckDB table, and add spatial geometry points so the data is ready for geospatial analysis.

In this step, get ready to:
  • Define 10 target megacities with their geographic coordinates.
  • Fetch PM2.5 station data from the OpenAQ API for each city.
  • Load stations into DuckDB with spatial geometry points and verify the data.
Define target cities and query the OpenAQ API

Each city needs a center coordinate so the API can search for stations within a radius of that point. You will use the 10 largest urban agglomerations by population.

  • In your notebook, create a new cell and define the CITIES dictionary by pasting the following code:
import requests
import pandas as pd
import time

# 10 largest megacities with center coordinates (lat, lon)
CITIES = {
    'Tokyo': {'lat': 35.6762, 'lon': 139.6503},
    'Delhi': {'lat': 28.6139, 'lon': 77.2090},
    'Shanghai': {'lat': 31.2304, 'lon': 121.4737},
    'Dhaka': {'lat': 23.8103, 'lon': 90.4125},
    'Sao Paulo': {'lat': -23.5505, 'lon': -46.6333},
    'Mexico City': {'lat': 19.4326, 'lon': -99.1332},
    'Cairo': {'lat': 30.0444, 'lon': 31.2357},
    'Beijing': {'lat': 39.9042, 'lon': 116.4074},
    'Mumbai': {'lat': 19.0760, 'lon': 72.8777},
    'Osaka': {'lat': 34.6937, 'lon': 135.5023},
}

print(f"Defined {len(CITIES)} target cities")

What does this code do?

The CITIES dictionary maps each megacity name to its geographic center as a latitude/longitude pair.

The OpenAQ API uses these coordinates to search for air quality stations within a given radius. You also import requests for HTTP calls, pandas for tabular data, and time for rate limiting between API calls.

  • Run the cell. You should see Defined 10 target cities printed.

Now you need a function that queries the OpenAQ API for each city. The GET /v3/locations endpoint accepts coordinates (latitude,longitude) and radius (in meters, max 25000) as query parameters.

  • Create a new cell and define the fetch_stations function by pasting this code:
def fetch_stations(cities, api_key):
    """Query OpenAQ v3 API for PM2.5 stations near each city center."""
    base_url = "https://api.openaq.org/v3/locations"
    headers = {"X-API-Key": api_key}
    all_stations = []

    for city_name, coords in cities.items():
        print(f"Fetching stations for {city_name}...")
        params = {
            "coordinates": f"{coords['lat']},{coords['lon']}",
            "radius": 25000,
            "parameters_id": 2,  # PM2.5
            "limit": 100,
        }

        response = requests.get(base_url, headers=headers, params=params)
        response.raise_for_status()
        data = response.json()

        # Parse each location from the results
        for location in data.get("results", []):
            all_stations.append({
                "city": city_name,
                "station_id": location["id"],
                "station_name": location.get("name", "Unknown"),
                "latitude": location["coordinates"]["latitude"],
                "longitude": location["coordinates"]["longitude"],
                "country": location["country"]["name"],
            })

        # Rate limit to avoid overwhelming the API
        time.sleep(1)

    return all_stations

print("fetch_stations function defined")

What does this code do?

The function loops through each city, sending a GET request to the OpenAQ /v3/locations endpoint with the city's coordinates and a 25km radius.

The parameters_id=2 filter ensures you only get stations that measure PM2.5 (fine particulate matter). The function parses each result's id, name, coordinates, and country fields from the JSON response.

A 1-second delay between requests prevents rate limiting from the API.

  • Run the cell. You should see fetch_stations function defined printed.

Seeing a 401 or 403 error?

This means your API key is invalid or missing. Double-check that the API_KEY variable in your earlier cell contains your actual key from the OpenAQ account settings.

Make sure you ran the cell that defines API_KEY before running this cell. Help me debug my OpenAQ API authentication error

Build the stations DataFrame

With the function defined, you can now call it to collect station data from all 10 cities and organize the results into a pandas DataFrame.

  • Create a new cell and run the fetch function by pasting and running this code:
# Fetch stations from all 10 cities
all_stations = fetch_stations(CITIES, API_KEY)

# Build a DataFrame from the collected station data
stations_df = pd.DataFrame(all_stations)
print(f"Total stations collected: {len(stations_df)}")
print(f"\nStations per city:")
print(stations_df['city'].value_counts())
stations_df.head()

What should I see?

The cell will print progress messages as it queries each city. This takes about 10-15 seconds due to the rate limiting.

Once complete, you should see the total number of stations collected (varies, but expect 50-200+ total) and a breakdown showing at least 5-10 stations per city. Some cities like Delhi and Beijing may have more stations than others due to denser monitoring networks.

  • Confirm that stations_df has columns: city, station_id, station_name, latitude, longitude, country.

Getting fewer stations than expected?

Some cities may return fewer stations if their monitoring network is sparse. As long as you have at least a few stations per city, the pipeline will work.

If a city returns 0 stations, the API might not have PM2.5 monitors in that radius. This is normal for some cities. Help me troubleshoot low station counts from OpenAQ

Load stations into DuckDB with geometry

Now you will persist this data into a DuckDB database file and add a spatial geometry column. This converts each station's lat/lon pair into a proper geometry point that DuckDB's spatial functions can operate on.

  • Create a new cell and set up the DuckDB connection and table by pasting this code:
import duckdb

# Connect to a persistent DuckDB database file
con = duckdb.connect('megacity_pipeline.duckdb')

# Install and load the spatial extension
con.execute("INSTALL spatial;")
con.execute("LOAD spatial;")

# Create the air_quality_stations table from the DataFrame
con.execute("CREATE OR REPLACE TABLE air_quality_stations AS SELECT * FROM stations_df")

# Add a geometry column computed from coordinates
con.execute("""
    ALTER TABLE air_quality_stations
    ADD COLUMN geometry GEOMETRY
""")
con.execute("""
    UPDATE air_quality_stations
    SET geometry = ST_Point(longitude, latitude)
""")

print("Table 'air_quality_stations' created with geometry column!")
con.execute("SELECT COUNT(*) FROM air_quality_stations").fetchone()

What does this code do?

The duckdb.connect('megacity_pipeline.duckdb') call creates a persistent database file on disk. Unlike the in-memory connection from Step 1, this file persists between notebook sessions.

After loading the DataFrame as a table, the code adds a GEOMETRY column and populates it using ST_Point(longitude, latitude). This creates a proper spatial point object for each station that DuckDB's spatial SQL functions can use in later steps.

  • Run the cell. You should see the confirmation message and a tuple showing the total row count.

Seeing a catalog error about the spatial extension?

If you get an error about the spatial extension not being found, make sure you ran INSTALL spatial; before LOAD spatial;. The install only needs to happen once per machine, but the load is required each time you open a new connection.

Help me fix a DuckDB spatial extension error

Finally, verify the data distribution across cities to make sure your ingestion captured stations from all 10 megacities.

  • Create a new cell and run the verification query:
# Verify station distribution across cities
result = con.execute("""
    SELECT city, COUNT(*) as station_count
    FROM air_quality_stations
    GROUP BY city
    ORDER BY COUNT(*) DESC
""").fetchdf()

print("Stations per city:")
print(result)
print(f"\nTotal stations: {result['station_count'].sum()}")
print(f"Cities covered: {len(result)}")

# Check for any null geometries
nulls = con.execute("""
    SELECT COUNT(*) FROM air_quality_stations
    WHERE geometry IS NULL
""").fetchone()[0]
print(f"Null geometries: {nulls}")

What should I see?

You should see a table listing all 10 cities with their station counts, ordered from highest to lowest. The total stations should match the earlier count from your DataFrame.

The null geometries check should return 0. This confirms every station has valid coordinates that were successfully converted to spatial points.

  • Confirm that you see all 10 cities in the output with at least a few stations each.
  • Confirm that null geometries shows 0.

✔️ Awesome, I've got everything!

Great. Double check you have saved your notebook.

ⓧ I'd like to double check the full code

Here is the complete code for this step. Make sure your notebook cells match:

import requests
import pandas as pd
import time

# 10 largest megacities with center coordinates (lat, lon)
CITIES = {
    'Tokyo': {'lat': 35.6762, 'lon': 139.6503},
    'Delhi': {'lat': 28.6139, 'lon': 77.2090},
    'Shanghai': {'lat': 31.2304, 'lon': 121.4737},
    'Dhaka': {'lat': 23.8103, 'lon': 90.4125},
    'Sao Paulo': {'lat': -23.5505, 'lon': -46.6333},
    'Mexico City': {'lat': 19.4326, 'lon': -99.1332},
    'Cairo': {'lat': 30.0444, 'lon': 31.2357},
    'Beijing': {'lat': 39.9042, 'lon': 116.4074},
    'Mumbai': {'lat': 19.0760, 'lon': 72.8777},
    'Osaka': {'lat': 34.6937, 'lon': 135.5023},
}

print(f"Defined {len(CITIES)} target cities")
def fetch_stations(cities, api_key):
    """Query OpenAQ v3 API for PM2.5 stations near each city center."""
    base_url = "https://api.openaq.org/v3/locations"
    headers = {"X-API-Key": api_key}
    all_stations = []

    for city_name, coords in cities.items():
        print(f"Fetching stations for {city_name}...")
        params = {
            "coordinates": f"{coords['lat']},{coords['lon']}",
            "radius": 25000,
            "parameters_id": 2,  # PM2.5
            "limit": 100,
        }

        response = requests.get(base_url, headers=headers, params=params)
        response.raise_for_status()
        data = response.json()

        # Parse each location from the results
        for location in data.get("results", []):
            all_stations.append({
                "city": city_name,
                "station_id": location["id"],
                "station_name": location.get("name", "Unknown"),
                "latitude": location["coordinates"]["latitude"],
                "longitude": location["coordinates"]["longitude"],
                "country": location["country"]["name"],
            })

        # Rate limit to avoid overwhelming the API
        time.sleep(1)

    return all_stations

print("fetch_stations function defined")
# Fetch stations from all 10 cities
all_stations = fetch_stations(CITIES, API_KEY)

# Build a DataFrame from the collected station data
stations_df = pd.DataFrame(all_stations)
print(f"Total stations collected: {len(stations_df)}")
print(f"\nStations per city:")
print(stations_df['city'].value_counts())
stations_df.head()
import duckdb

# Connect to a persistent DuckDB database file
con = duckdb.connect('megacity_pipeline.duckdb')

# Install and load the spatial extension
con.execute("INSTALL spatial;")
con.execute("LOAD spatial;")

# Create the air_quality_stations table from the DataFrame
con.execute("CREATE OR REPLACE TABLE air_quality_stations AS SELECT * FROM stations_df")

# Add a geometry column computed from coordinates
con.execute("""
    ALTER TABLE air_quality_stations
    ADD COLUMN geometry GEOMETRY
""")
con.execute("""
    UPDATE air_quality_stations
    SET geometry = ST_Point(longitude, latitude)
""")

print("Table 'air_quality_stations' created with geometry column!")
con.execute("SELECT COUNT(*) FROM air_quality_stations").fetchone()
# Verify station distribution across cities
result = con.execute("""
    SELECT city, COUNT(*) as station_count
    FROM air_quality_stations
    GROUP BY city
    ORDER BY COUNT(*) DESC
""").fetchdf()

print("Stations per city:")
print(result)
print(f"\nTotal stations: {result['station_count'].sum()}")
print(f"Cities covered: {len(result)}")

# Check for any null geometries
nulls = con.execute("""
    SELECT COUNT(*) FROM air_quality_stations
    WHERE geometry IS NULL
""").fetchone()[0]
print(f"Null geometries: {nulls}")

Your air quality station data is now safely stored in DuckDB with proper spatial geometry. Next, you will download geospatial features from OpenStreetMap to enrich each station's context.

Download and Load Geospatial Features from OpenStreetMap

You have air quality stations loaded into DuckDB, but station coordinates alone do not tell the whole story. What surrounds each monitor matters: busy highways generate exhaust, parks filter particulates, and factories emit pollutants.

In this step, you will download road networks, green spaces, and industrial zones from OpenStreetMap for all 10 cities. You will then convert these geometries into a format DuckDB understands and load them as spatial tables.

In this step, get ready to:
  • Download road networks, green spaces, and industrial zones from OpenStreetMap using OSMnx.
  • Convert GeoDataFrame geometries to WKB format and load them into DuckDB spatial tables.
  • Verify geometry types and row counts with spatial SQL queries.
Download road networks from OpenStreetMap

OSMnx is a Python library that simplifies downloading geospatial data from OpenStreetMap. It handles the Overpass API pagination for you and returns clean GeoDataFrames.

  • In a new cell in your notebook, import OSMnx and download road features for all 10 cities by running this code:
import osmnx as ox
import geopandas as gpd

# Download road networks for each city (15km radius from center)
print("Downloading road networks...")
roads_list = []

for city_name, info in CITIES.items():
    center = (info['lat'], info['lon'])
    try:
        # Fetch all highway features within 15km of city center
        roads = ox.features.features_from_point(
            center, tags={'highway': True}, dist=15000
        )
        roads['city'] = city_name
        # Keep city, geometry, and highway type columns
        roads_list.append(roads[['city', 'geometry', 'highway']].reset_index(drop=True))
        print(f"  {city_name}: {len(roads)} road segments")
    except Exception as e:
        print(f"  {city_name}: Error - {e}")

roads_gdf = gpd.GeoDataFrame(pd.concat(roads_list, ignore_index=True))
print(f"\nTotal road segments: {len(roads_gdf)}")

What does this code do?

ox.features.features_from_point() downloads OpenStreetMap features within a bounding box around a center point. The tags={'highway': True} parameter fetches all road types (motorways, residential streets, footpaths).

The dist=15000 parameter creates a 15km bounding box from the center. This captures the urban core where air quality monitors are concentrated without overwhelming the API with massive administrative boundary areas.

This download takes a few minutes

Fetching road data for 10 global cities is network-intensive. Expect the cell to run for 3-5 minutes depending on your connection speed. The progress prints let you see which city is currently downloading.

  • Run the cell and wait for it to complete.
  • You should see a count of road segments for each city printed in the output, followed by a total count (likely in the tens of thousands).

Seeing timeout or connection errors?

The Overpass API (which OSMnx uses behind the scenes) can be busy during peak hours. Wait a minute and re-run the cell. If a specific city keeps failing, that city's data will be skipped thanks to the try/except block.

help me fix OSMnx timeout errors when downloading OpenStreetMap data

Download green spaces and industrial zones

Green spaces absorb pollutants, while industrial zones produce them. These two feature types complete the environmental context around each air quality station.

  • In a new cell, download green spaces and industrial zones for all 10 cities by running this code:
# Download green spaces for each city
print("Downloading green spaces...")
green_list = []

for city_name, info in CITIES.items():
    center = (info['lat'], info['lon'])
    try:
        # Fetch parks, forests, grass, and meadows
        greens = ox.features.features_from_point(
            center,
            tags={'leisure': 'park', 'landuse': ['forest', 'grass', 'meadow']},
            dist=15000
        )
        greens['city'] = city_name
        green_list.append(greens[['city', 'geometry']].reset_index(drop=True))
        print(f"  {city_name}: {len(greens)} green spaces")
    except Exception as e:
        print(f"  {city_name}: Error - {e}")

green_spaces_gdf = gpd.GeoDataFrame(pd.concat(green_list, ignore_index=True))
print(f"\nTotal green spaces: {len(green_spaces_gdf)}")

# Download industrial zones for each city
print("\nDownloading industrial zones...")
industrial_list = []

for city_name, info in CITIES.items():
    center = (info['lat'], info['lon'])
    try:
        # Fetch areas tagged as industrial land use
        industrial = ox.features.features_from_point(
            center,
            tags={'landuse': 'industrial'},
            dist=15000
        )
        industrial['city'] = city_name
        industrial_list.append(industrial[['city', 'geometry']].reset_index(drop=True))
        print(f"  {city_name}: {len(industrial)} industrial zones")
    except Exception as e:
        print(f"  {city_name}: Error - {e}")

industrial_zones_gdf = gpd.GeoDataFrame(pd.concat(industrial_list, ignore_index=True))
print(f"\nTotal industrial zones: {len(industrial_zones_gdf)}")

What do these tags mean?

The tags dictionary tells OSMnx which OpenStreetMap features to fetch. For green spaces, {'leisure': 'park', 'landuse': ['forest', 'grass', 'meadow']} gets the union of all matching features: parks OR forests OR grass OR meadows.

For industrial zones, {'landuse': 'industrial'} fetches areas explicitly tagged as industrial land use in OpenStreetMap. The coverage varies by city depending on how thoroughly local mappers have tagged their area.

  • Run the cell and wait for it to complete.
  • You should see green space and industrial zone counts for each city. Some cities may have fewer industrial zones than others depending on OpenStreetMap coverage.

Some cities show 0 results or errors?

OpenStreetMap coverage varies by region. Some cities (like Tokyo and Osaka) have excellent coverage, while others may have fewer tagged features. As long as most cities return data, your pipeline is working correctly.

help me understand why some cities return no green spaces or industrial zones from OSMnx

Convert geometries to WKB and load into DuckDB

DuckDB cannot directly read GeoPandas geometry objects. You need to convert them to WKB (Well-Known Binary), a standard binary encoding for geometries that DuckDB can deserialize with ST_GeomFromWKB().

  • In a new cell, convert the GeoDataFrames to WKB and load them into DuckDB by running this code:
# Convert geometries to WKB (Well-Known Binary) for DuckDB ingestion
roads_gdf['geometry_wkb'] = roads_gdf['geometry'].apply(lambda g: g.wkb)
green_spaces_gdf['geometry_wkb'] = green_spaces_gdf['geometry'].apply(lambda g: g.wkb)
industrial_zones_gdf['geometry_wkb'] = industrial_zones_gdf['geometry'].apply(lambda g: g.wkb)

# Convert highway column to string (OSMnx may return lists for some features)
roads_gdf['highway_str'] = roads_gdf['highway'].apply(
    lambda x: x[0] if isinstance(x, list) else str(x)
)

# Create spatial tables in DuckDB
con.execute("""
    CREATE OR REPLACE TABLE roads (
        city VARCHAR,
        geometry GEOMETRY,
        highway VARCHAR
    )
""")

con.execute("""
    CREATE OR REPLACE TABLE green_spaces (
        city VARCHAR,
        geometry GEOMETRY
    )
""")

con.execute("""
    CREATE OR REPLACE TABLE industrial_zones (
        city VARCHAR,
        geometry GEOMETRY
    )
""")

# Load roads into DuckDB using ST_GeomFromWKB to convert WKB bytes to geometry
roads_load = roads_gdf[['city', 'geometry_wkb', 'highway_str']].copy()
con.execute("""
    INSERT INTO roads
    SELECT city, ST_GeomFromWKB(geometry_wkb), highway_str
    FROM roads_load
""")

# Load green spaces
green_load = green_spaces_gdf[['city', 'geometry_wkb']].copy()
con.execute("""
    INSERT INTO green_spaces
    SELECT city, ST_GeomFromWKB(geometry_wkb)
    FROM green_load
""")

# Load industrial zones
industrial_load = industrial_zones_gdf[['city', 'geometry_wkb']].copy()
con.execute("""
    INSERT INTO industrial_zones
    SELECT city, ST_GeomFromWKB(geometry_wkb)
    FROM industrial_load
""")

print("All spatial tables loaded into DuckDB.")
print(f"  Roads: {con.execute('SELECT COUNT(*) FROM roads').fetchone()[0]} rows")
print(f"  Green spaces: {con.execute('SELECT COUNT(*) FROM green_spaces').fetchone()[0]} rows")
print(f"  Industrial zones: {con.execute('SELECT COUNT(*) FROM industrial_zones').fetchone()[0]} rows")

What is the WKB transfer pattern?

GeoPandas stores geometries as Shapely objects. DuckDB needs its own native GEOMETRY type. WKB acts as the bridge between them.

  • geometry.wkb converts each Shapely geometry to a binary blob.
  • ST_GeomFromWKB() deserializes that blob back into a DuckDB GEOMETRY during INSERT.
  • DuckDB's replacement scan reads the pandas DataFrame variable directly in SQL, so no explicit registration is needed.
  • Run the cell.
  • You should see confirmation that all three tables are loaded, with row counts printed for each. The roads table should have the highest count.

Seeing a conversion error on ST_GeomFromWKB?

Make sure the spatial extension is still loaded in your DuckDB connection. If your kernel restarted, you need to re-run INSTALL spatial; LOAD spatial; from your earlier setup cell.

help me fix ST_GeomFromWKB conversion errors in DuckDB when loading OSMnx geometries

  • Verify the tables loaded correctly by running this spatial SQL in a new cell:
# Verify road counts per city
print("Road segments per city:")
print(con.execute("SELECT city, COUNT(*) FROM roads GROUP BY city").fetchdf())

# Verify geometry types for green spaces
print("\nGreen space geometry types:")
print(con.execute("SELECT ST_GeometryType(geometry) FROM green_spaces LIMIT 5").fetchdf())

# Verify geometry types for industrial zones
print("\nIndustrial zone geometry types:")
print(con.execute("SELECT ST_GeometryType(geometry) FROM industrial_zones LIMIT 5").fetchdf())

What should I see?

The first query shows a table with each city and its road segment count. Cities like Tokyo and Delhi should have thousands of road segments each.

The geometry type queries should show values like POLYGON or MULTIPOLYGON for green spaces and industrial zones. Roads will typically be LINESTRING or MULTILINESTRING geometries.

  • Run the cell and confirm you see row counts for all 10 cities and valid geometry types.

Missing cities in the road count query?

If a city does not appear, its OSMnx download likely failed earlier. Scroll up to the download cell's output and check for error messages on that city. You can re-run just the failing city's download by adjusting the loop or running a single-city call.

help me troubleshoot why some cities are missing from my DuckDB roads table after OSMnx download

✔️ Awesome, I've got everything!

Great. Make sure you have saved your notebook before moving on.

ⓧ I'd like to double check the full code

Here is all the code you added in this step across your notebook cells:

import osmnx as ox
import geopandas as gpd

# Download road networks for each city (15km radius from center)
print("Downloading road networks...")
roads_list = []

for city_name, info in CITIES.items():
    center = (info['lat'], info['lon'])
    try:
        # Fetch all highway features within 15km of city center
        roads = ox.features.features_from_point(
            center, tags={'highway': True}, dist=15000
        )
        roads['city'] = city_name
        # Keep city, geometry, and highway type columns
        roads_list.append(roads[['city', 'geometry', 'highway']].reset_index(drop=True))
        print(f"  {city_name}: {len(roads)} road segments")
    except Exception as e:
        print(f"  {city_name}: Error - {e}")

roads_gdf = gpd.GeoDataFrame(pd.concat(roads_list, ignore_index=True))
print(f"\nTotal road segments: {len(roads_gdf)}")

# --- Next cell ---

# Download green spaces for each city
print("Downloading green spaces...")
green_list = []

for city_name, info in CITIES.items():
    center = (info['lat'], info['lon'])
    try:
        # Fetch parks, forests, grass, and meadows
        greens = ox.features.features_from_point(
            center,
            tags={'leisure': 'park', 'landuse': ['forest', 'grass', 'meadow']},
            dist=15000
        )
        greens['city'] = city_name
        green_list.append(greens[['city', 'geometry']].reset_index(drop=True))
        print(f"  {city_name}: {len(greens)} green spaces")
    except Exception as e:
        print(f"  {city_name}: Error - {e}")

green_spaces_gdf = gpd.GeoDataFrame(pd.concat(green_list, ignore_index=True))
print(f"\nTotal green spaces: {len(green_spaces_gdf)}")

# Download industrial zones for each city
print("\nDownloading industrial zones...")
industrial_list = []

for city_name, info in CITIES.items():
    center = (info['lat'], info['lon'])
    try:
        # Fetch areas tagged as industrial land use
        industrial = ox.features.features_from_point(
            center,
            tags={'landuse': 'industrial'},
            dist=15000
        )
        industrial['city'] = city_name
        industrial_list.append(industrial[['city', 'geometry']].reset_index(drop=True))
        print(f"  {city_name}: {len(industrial)} industrial zones")
    except Exception as e:
        print(f"  {city_name}: Error - {e}")

industrial_zones_gdf = gpd.GeoDataFrame(pd.concat(industrial_list, ignore_index=True))
print(f"\nTotal industrial zones: {len(industrial_zones_gdf)}")

# --- Next cell ---

# Convert geometries to WKB (Well-Known Binary) for DuckDB ingestion
roads_gdf['geometry_wkb'] = roads_gdf['geometry'].apply(lambda g: g.wkb)
green_spaces_gdf['geometry_wkb'] = green_spaces_gdf['geometry'].apply(lambda g: g.wkb)
industrial_zones_gdf['geometry_wkb'] = industrial_zones_gdf['geometry'].apply(lambda g: g.wkb)

# Convert highway column to string (OSMnx may return lists for some features)
roads_gdf['highway_str'] = roads_gdf['highway'].apply(
    lambda x: x[0] if isinstance(x, list) else str(x)
)

# Create spatial tables in DuckDB
con.execute("""
    CREATE OR REPLACE TABLE roads (
        city VARCHAR,
        geometry GEOMETRY,
        highway VARCHAR
    )
""")

con.execute("""
    CREATE OR REPLACE TABLE green_spaces (
        city VARCHAR,
        geometry GEOMETRY
    )
""")

con.execute("""
    CREATE OR REPLACE TABLE industrial_zones (
        city VARCHAR,
        geometry GEOMETRY
    )
""")

# Load roads into DuckDB using ST_GeomFromWKB to convert WKB bytes to geometry
roads_load = roads_gdf[['city', 'geometry_wkb', 'highway_str']].copy()
con.execute("""
    INSERT INTO roads
    SELECT city, ST_GeomFromWKB(geometry_wkb), highway_str
    FROM roads_load
""")

# Load green spaces
green_load = green_spaces_gdf[['city', 'geometry_wkb']].copy()
con.execute("""
    INSERT INTO green_spaces
    SELECT city, ST_GeomFromWKB(geometry_wkb)
    FROM green_load
""")

# Load industrial zones
industrial_load = industrial_zones_gdf[['city', 'geometry_wkb']].copy()
con.execute("""
    INSERT INTO industrial_zones
    SELECT city, ST_GeomFromWKB(geometry_wkb)
    FROM industrial_load
""")

print("All spatial tables loaded into DuckDB.")
print(f"  Roads: {con.execute('SELECT COUNT(*) FROM roads').fetchone()[0]} rows")
print(f"  Green spaces: {con.execute('SELECT COUNT(*) FROM green_spaces').fetchone()[0]} rows")
print(f"  Industrial zones: {con.execute('SELECT COUNT(*) FROM industrial_zones').fetchone()[0]} rows")

# --- Next cell ---

# Verify road counts per city
print("Road segments per city:")
print(con.execute("SELECT city, COUNT(*) FROM roads GROUP BY city").fetchdf())

# Verify geometry types for green spaces
print("\nGreen space geometry types:")
print(con.execute("SELECT ST_GeometryType(geometry) FROM green_spaces LIMIT 5").fetchdf())

# Verify geometry types for industrial zones
print("\nIndustrial zone geometry types:")
print(con.execute("SELECT ST_GeometryType(geometry) FROM industrial_zones LIMIT 5").fetchdf())

Your DuckDB database now holds four spatial tables: station points, road networks, green spaces, and industrial zones across 10 megacities. Next up, you will write spatial SQL queries to compute meaningful features from these layers around each air quality monitor.

Compute Spatial Features with SQL

With all your spatial data loaded into DuckDB, you can now answer the real question: what does the physical environment around each air quality monitor look like?

Spatial SQL lets you compute features like road density, green space coverage, and industrial proximity for every station. These features capture the local context that influences pollution levels.

In this step, get ready to:
  • Compute total road length within 500m of each station using buffer analysis.
  • Calculate green space ratio and industrial area within buffers around each station.
  • Join all spatial features into a single ML-ready table called station_features.
Compute road density within a 500m buffer

The first spatial feature measures how much road infrastructure surrounds each monitor. More roads typically means more traffic emissions.

How does buffer analysis work?

A buffer creates a circular zone around a point at a given radius. You then find all features (like roads) that intersect that zone, clip them to the boundary, and measure what remains inside.

Since our coordinates are in EPSG:4326 (degrees), ST_Buffer takes its distance in degrees too. We approximate 500m as 0.0045 degrees (roughly 500 / 111320). For accurate length and distance results in meters, we use the spheroid-aware variants like ST_Length_Spheroid and ST_Distance_Spheroid.

  • Create a new cell in your notebook and add the road density query by pasting the following code:
# Compute total road length within 500m of each air quality station
road_density_query = """
CREATE OR REPLACE TABLE road_density AS
SELECT
    s.station_id,
    s.city,
    -- Create ~500m buffer (in degrees), find intersecting roads,
    -- clip them to the buffer, then measure total length in meters
    COALESCE(
        SUM(ST_Length_Spheroid(
            ST_Intersection(r.geometry, ST_Buffer(s.geometry, 0.0045))
        )),
        0
    ) AS road_length_500m
FROM air_quality_stations s
LEFT JOIN roads r
    ON s.city = r.city
    AND ST_Intersects(r.geometry, ST_Buffer(s.geometry, 0.0045))
GROUP BY s.station_id, s.city;
"""

con.execute(road_density_query)

# Verify road density results per city
result = con.execute("""
    SELECT city,
           COUNT(*) AS stations,
           ROUND(AVG(road_length_500m), 1) AS avg_road_m,
           ROUND(MAX(road_length_500m), 1) AS max_road_m
    FROM road_density
    GROUP BY city
    ORDER BY avg_road_m DESC
""").fetchdf()

print(result)

What does this query do?

  • ST_Buffer(s.geometry, 0.0045) creates a circular zone of approximately 500 meters around each station point.
  • ST_Intersects filters to only roads that overlap with that buffer zone.
  • ST_Intersection clips each road segment to the buffer boundary so you only measure the portion inside.
  • ST_Length_Spheroid measures the clipped road length in meters using an ellipsoidal Earth model.
  • The LEFT JOIN ensures stations with no nearby roads still appear with a value of 0.
  • Run the cell.
  • You should see a table showing average and maximum road lengths per city, with values in the thousands of meters for dense urban areas.

Query running too slowly?

Buffer analysis with spatial joins can take 30-60 seconds across 10 cities. If it takes longer than 2 minutes, check that your roads table has the city column populated correctly. The ON s.city = r.city clause prevents cross-city comparisons that would explode the join size.

help me debug why my spatial join is slow

Compute green space ratio and industrial area

The next features capture environmental context: how much green space surrounds each station (a ratio from 0 to 1) and how much industrial land sits within 500m.

  • Add a new cell with the green space ratio query:
# Compute green space ratio within 1km of each station
green_space_query = """
CREATE OR REPLACE TABLE green_features AS
SELECT
    s.station_id,
    s.city,
    -- Ratio of green area to total buffer area within ~1km
    COALESCE(
        SUM(ST_Area_Spheroid(
            ST_Intersection(g.geometry, ST_Buffer(s.geometry, 0.009))
        )),
        0
    ) / ST_Area_Spheroid(ST_Buffer(s.geometry, 0.009)) AS green_area_ratio_1km
FROM air_quality_stations s
LEFT JOIN green_spaces g
    ON s.city = g.city
    AND ST_Intersects(g.geometry, ST_Buffer(s.geometry, 0.009))
GROUP BY s.station_id, s.city, s.geometry;
"""

con.execute(green_space_query)

# Verify green ratio values are between 0 and 1
green_check = con.execute("""
    SELECT city,
           ROUND(AVG(green_area_ratio_1km), 4) AS avg_green_ratio,
           ROUND(MAX(green_area_ratio_1km), 4) AS max_green_ratio
    FROM green_features
    GROUP BY city
    ORDER BY avg_green_ratio DESC
""").fetchdf()

print(green_check)

What does this query do?

  • ST_Buffer(s.geometry, 0.009) creates an approximately 1km buffer around each station (1000 / 111320 degrees).
  • ST_Area_Spheroid computes area in square meters on an ellipsoidal Earth model.
  • Dividing the total green area by the buffer area gives a ratio between 0 (no green space) and 1 (entirely green).
  • Run the cell.
  • You should see average green ratios per city, typically between 0.01 and 0.3.
  • Add another cell for the industrial area computation:
# Compute total industrial area within 500m of each station
industrial_query = """
CREATE OR REPLACE TABLE industrial_features AS
SELECT
    s.station_id,
    s.city,
    -- Total industrial zone area within ~500m buffer, in square meters
    COALESCE(
        SUM(ST_Area_Spheroid(
            ST_Intersection(i.geometry, ST_Buffer(s.geometry, 0.0045))
        )),
        0
    ) AS industrial_area_500m
FROM air_quality_stations s
LEFT JOIN industrial_zones i
    ON s.city = i.city
    AND ST_Intersects(i.geometry, ST_Buffer(s.geometry, 0.0045))
GROUP BY s.station_id, s.city, s.geometry;
"""

con.execute(industrial_query)

# Verify industrial area results
ind_check = con.execute("""
    SELECT city,
           ROUND(AVG(industrial_area_500m), 1) AS avg_industrial_m2,
           COUNT(*) AS stations
    FROM industrial_features
    GROUP BY city
    ORDER BY avg_industrial_m2 DESC
""").fetchdf()

print(ind_check)

What does this query do?

This follows the same buffer-intersect-measure pattern as road density, but measures area instead of length. Stations near factories or industrial parks will have high industrial_area_500m values (in square meters), while residential-area stations will show zero.

  • Run the cell.
  • Cities with heavy manufacturing (like Delhi, Beijing, Shanghai) should show higher average industrial area values.

Getting all zeros for industrial area?

Some cities may legitimately have very few industrial zones near air quality monitors. If ALL cities show zero, verify your industrial_zones table has rows by running SELECT city, COUNT(*) FROM industrial_zones GROUP BY city.

help me debug why my industrial area query returns all zeros

Compute highway distance and join all features

The final spatial feature measures how far each station is from major highways. Stations close to motorways experience higher pollution from heavy traffic.

  • Add a new cell with the highway distance query:
# Compute distance to nearest major highway for each station
highway_distance_query = """
CREATE OR REPLACE TABLE highway_distance AS
SELECT
    s.station_id,
    s.city,
    -- Filter to major road types and find the closest one
    COALESCE(
        MIN(ST_Distance_Spheroid(s.geometry, r.geometry)),
        15000.0
    ) AS dist_to_nearest_highway_m
FROM air_quality_stations s
LEFT JOIN roads r
    ON s.city = r.city
    AND r.highway IN ('motorway', 'trunk', 'primary')
GROUP BY s.station_id, s.city;
"""

con.execute(highway_distance_query)

# Verify highway distance results
hwy_check = con.execute("""
    SELECT city,
           ROUND(AVG(dist_to_nearest_highway_m), 1) AS avg_dist_m,
           ROUND(MIN(dist_to_nearest_highway_m), 1) AS min_dist_m,
           ROUND(MAX(dist_to_nearest_highway_m), 1) AS max_dist_m
    FROM highway_distance
    GROUP BY city
    ORDER BY avg_dist_m
""").fetchdf()

print(hwy_check)

What does this query do?

  • The JOIN ... AND r.highway IN (...) condition filters roads to only major highways: motorways, trunk roads, and primary roads.
  • ST_Distance_Spheroid returns the geodesic distance in meters between each station and each qualifying road segment.
  • MIN(...) picks the closest highway. The COALESCE(..., 15000.0) assigns 15km as a default if no major highway exists nearby.
  • Run the cell.
  • You should see distance values ranging from a few hundred meters to several thousand meters across cities.

Now join all four feature tables into the final station_features table.

  • Add one more cell to create the joined table and validate it:
# Join all spatial features into the final station_features table
join_query = """
CREATE OR REPLACE TABLE station_features AS
SELECT
    s.station_id,
    s.city,
    s.latitude,
    s.longitude,
    rd.road_length_500m,
    gf.green_area_ratio_1km,
    hd.dist_to_nearest_highway_m,
    ind.industrial_area_500m
FROM air_quality_stations s
JOIN road_density rd ON s.station_id = rd.station_id
JOIN green_features gf ON s.station_id = gf.station_id
JOIN highway_distance hd ON s.station_id = hd.station_id
JOIN industrial_features ind ON s.station_id = ind.station_id;
"""

con.execute(join_query)

# Validate: check for nulls and value ranges
validation = con.execute("""
    SELECT
        COUNT(*) AS total_stations,
        COUNT(*) - COUNT(road_length_500m) AS null_road,
        COUNT(*) - COUNT(green_area_ratio_1km) AS null_green,
        COUNT(*) - COUNT(dist_to_nearest_highway_m) AS null_highway,
        COUNT(*) - COUNT(industrial_area_500m) AS null_industrial,
        ROUND(MIN(road_length_500m), 1) AS min_road,
        ROUND(MAX(road_length_500m), 1) AS max_road,
        ROUND(MIN(green_area_ratio_1km), 4) AS min_green,
        ROUND(MAX(green_area_ratio_1km), 4) AS max_green,
        ROUND(MIN(dist_to_nearest_highway_m), 1) AS min_highway,
        ROUND(MAX(dist_to_nearest_highway_m), 1) AS max_highway
    FROM station_features
""").fetchdf()

print("=== Station Features Validation ===")
print(validation)

# Show sample rows
sample = con.execute("""
    SELECT * FROM station_features LIMIT 5
""").fetchdf()

print("\n=== Sample Rows ===")
print(sample)

What should I see?

The validation query confirms your data quality. All null counts should be 0. Road length values should be positive numbers. Green ratio should be between 0 and 1. Highway distance should range from 0 to 15000 meters.

The sample rows show one row per station with all four spatial feature columns populated.

  • Run the cell.
  • Confirm that total_stations matches the number of stations you loaded earlier, and all null counts show 0.

Seeing fewer rows than expected?

The JOIN (inner join) drops any station missing from one of the feature tables. Check each intermediate table's row count with SELECT COUNT(*) FROM road_density and compare to SELECT COUNT(*) FROM air_quality_stations. If counts differ, a station was lost during grouping.

help me debug why my station_features table has fewer rows than air_quality_stations

✔️ Awesome, I've got everything!

Great. Double check you have saved your notebook.

ⓧ I'd like to double check the full code

Here is the complete code for this step's cells:

# --- Cell 1: Road Density ---

# Compute total road length within 500m of each air quality station
road_density_query = """
CREATE OR REPLACE TABLE road_density AS
SELECT
    s.station_id,
    s.city,
    -- Create ~500m buffer (in degrees), find intersecting roads,
    -- clip them to the buffer, then measure total length in meters
    COALESCE(
        SUM(ST_Length_Spheroid(
            ST_Intersection(r.geometry, ST_Buffer(s.geometry, 0.0045))
        )),
        0
    ) AS road_length_500m
FROM air_quality_stations s
LEFT JOIN roads r
    ON s.city = r.city
    AND ST_Intersects(r.geometry, ST_Buffer(s.geometry, 0.0045))
GROUP BY s.station_id, s.city;
"""

con.execute(road_density_query)

# Verify road density results per city
result = con.execute("""
    SELECT city,
           COUNT(*) AS stations,
           ROUND(AVG(road_length_500m), 1) AS avg_road_m,
           ROUND(MAX(road_length_500m), 1) AS max_road_m
    FROM road_density
    GROUP BY city
    ORDER BY avg_road_m DESC
""").fetchdf()

print(result)


# --- Cell 2: Green Space Ratio ---

# Compute green space ratio within 1km of each station
green_space_query = """
CREATE OR REPLACE TABLE green_features AS
SELECT
    s.station_id,
    s.city,
    -- Ratio of green area to total buffer area within ~1km
    COALESCE(
        SUM(ST_Area_Spheroid(
            ST_Intersection(g.geometry, ST_Buffer(s.geometry, 0.009))
        )),
        0
    ) / ST_Area_Spheroid(ST_Buffer(s.geometry, 0.009)) AS green_area_ratio_1km
FROM air_quality_stations s
LEFT JOIN green_spaces g
    ON s.city = g.city
    AND ST_Intersects(g.geometry, ST_Buffer(s.geometry, 0.009))
GROUP BY s.station_id, s.city, s.geometry;
"""

con.execute(green_space_query)

# Verify green ratio values are between 0 and 1
green_check = con.execute("""
    SELECT city,
           ROUND(AVG(green_area_ratio_1km), 4) AS avg_green_ratio,
           ROUND(MAX(green_area_ratio_1km), 4) AS max_green_ratio
    FROM green_features
    GROUP BY city
    ORDER BY avg_green_ratio DESC
""").fetchdf()

print(green_check)


# --- Cell 3: Industrial Area ---

# Compute total industrial area within 500m of each station
industrial_query = """
CREATE OR REPLACE TABLE industrial_features AS
SELECT
    s.station_id,
    s.city,
    -- Total industrial zone area within ~500m buffer, in square meters
    COALESCE(
        SUM(ST_Area_Spheroid(
            ST_Intersection(i.geometry, ST_Buffer(s.geometry, 0.0045))
        )),
        0
    ) AS industrial_area_500m
FROM air_quality_stations s
LEFT JOIN industrial_zones i
    ON s.city = i.city
    AND ST_Intersects(i.geometry, ST_Buffer(s.geometry, 0.0045))
GROUP BY s.station_id, s.city, s.geometry;
"""

con.execute(industrial_query)

# Verify industrial area results
ind_check = con.execute("""
    SELECT city,
           ROUND(AVG(industrial_area_500m), 1) AS avg_industrial_m2,
           COUNT(*) AS stations
    FROM industrial_features
    GROUP BY city
    ORDER BY avg_industrial_m2 DESC
""").fetchdf()

print(ind_check)


# --- Cell 4: Highway Distance ---

# Compute distance to nearest major highway for each station
highway_distance_query = """
CREATE OR REPLACE TABLE highway_distance AS
SELECT
    s.station_id,
    s.city,
    -- Filter to major road types and find the closest one
    COALESCE(
        MIN(ST_Distance_Spheroid(s.geometry, r.geometry)),
        15000.0
    ) AS dist_to_nearest_highway_m
FROM air_quality_stations s
LEFT JOIN roads r
    ON s.city = r.city
    AND r.highway IN ('motorway', 'trunk', 'primary')
GROUP BY s.station_id, s.city;
"""

con.execute(highway_distance_query)

# Verify highway distance results
hwy_check = con.execute("""
    SELECT city,
           ROUND(AVG(dist_to_nearest_highway_m), 1) AS avg_dist_m,
           ROUND(MIN(dist_to_nearest_highway_m), 1) AS min_dist_m,
           ROUND(MAX(dist_to_nearest_highway_m), 1) AS max_dist_m
    FROM highway_distance
    GROUP BY city
    ORDER BY avg_dist_m
""").fetchdf()

print(hwy_check)


# --- Cell 5: Join All Features ---

# Join all spatial features into the final station_features table
join_query = """
CREATE OR REPLACE TABLE station_features AS
SELECT
    s.station_id,
    s.city,
    s.latitude,
    s.longitude,
    rd.road_length_500m,
    gf.green_area_ratio_1km,
    hd.dist_to_nearest_highway_m,
    ind.industrial_area_500m
FROM air_quality_stations s
JOIN road_density rd ON s.station_id = rd.station_id
JOIN green_features gf ON s.station_id = gf.station_id
JOIN highway_distance hd ON s.station_id = hd.station_id
JOIN industrial_features ind ON s.station_id = ind.station_id;
"""

con.execute(join_query)

# Validate: check for nulls and value ranges
validation = con.execute("""
    SELECT
        COUNT(*) AS total_stations,
        COUNT(*) - COUNT(road_length_500m) AS null_road,
        COUNT(*) - COUNT(green_area_ratio_1km) AS null_green,
        COUNT(*) - COUNT(dist_to_nearest_highway_m) AS null_highway,
        COUNT(*) - COUNT(industrial_area_500m) AS null_industrial,
        ROUND(MIN(road_length_500m), 1) AS min_road,
        ROUND(MAX(road_length_500m), 1) AS max_road,
        ROUND(MIN(green_area_ratio_1km), 4) AS min_green,
        ROUND(MAX(green_area_ratio_1km), 4) AS max_green,
        ROUND(MIN(dist_to_nearest_highway_m), 1) AS min_highway,
        ROUND(MAX(dist_to_nearest_highway_m), 1) AS max_highway
    FROM station_features
""").fetchdf()

print("=== Station Features Validation ===")
print(validation)

# Show sample rows
sample = con.execute("""
    SELECT * FROM station_features LIMIT 5
""").fetchdf()

print("\n=== Sample Rows ===")
print(sample)

Your station_features table now holds four spatial features for every air quality monitor across 10 megacities. Next up, you will export this dataset to CSV and create a visual map to validate everything looks correct.

Export and Validate the ML-Ready Dataset

In the previous step, you computed spatial features for every air quality station across all 10 megacities. Your station_features table now holds the full picture: road density, green space ratio, highway distance, and industrial proximity for each monitor.

Now it is time to validate those features and export the final dataset. You will run summary statistics to confirm the spatial features make geographic sense, export the table as a CSV ready for machine learning, and create an interactive map to visually confirm the pipeline worked end-to-end.

In this step, get ready to:
  • Run summary statistics to validate spatial feature distributions across cities.
  • Export the feature-engineered dataset to a CSV file.
  • Create an interactive Folium map with color-coded station markers and feature tooltips.
Run summary statistics per city

Before exporting, you want confidence that the spatial features reflect real-world geography. Cities with more roads should show higher road density, and cities with more parks should show higher green space ratios.

  • In a new cell in air_quality_pipeline.ipynb, run a summary statistics query by pasting this code:
# Compute average spatial features per city for validation
summary_query = """
SELECT
    city,
    AVG(road_length_500m) AS avg_road_length,
    AVG(green_area_ratio_1km) AS avg_green_ratio,
    AVG(dist_to_nearest_highway_m) AS avg_highway_dist
FROM station_features
GROUP BY city
ORDER BY city
"""

summary_df = con.execute(summary_query).fetchdf()
print(summary_df)

What does this query do?

This query groups all stations by city and computes the average value for each spatial feature. It gives you a quick comparison across cities.

You should see variation that matches real-world expectations. For example, Delhi and Beijing should have relatively high road density values, while cities with more parks should show higher green area ratios.

  • Run the cell.
  • You should see a DataFrame with one row per city and three numeric columns showing the averaged spatial features.

Seeing all zeros or NaN values?

If all values are 0 or NaN, the spatial joins in the previous step may not have matched any geometries to stations. Go back to Step 4 and verify that the station_features table has non-null values by running SELECT * FROM station_features LIMIT 5.

Help me debug why my summary statistics return zeros or NaN

Export the dataset to CSV

With the features validated, you can now export the complete dataset as a CSV file. DuckDB's COPY command writes a table directly to disk without needing to pull data through Python first.

  • Export the station_features table to CSV by running this code in a new cell:
# Export the ML-ready dataset to CSV
con.execute("""
COPY station_features TO 'megacity_air_quality_features.csv' (HEADER, DELIMITER ',')
""")

print("CSV exported: megacity_air_quality_features.csv")

What does COPY TO do?

HEADER includes column names as the first row of the CSV. DELIMITER ',' separates values with commas.

DuckDB writes the file directly from its internal storage, so this is much faster than converting to a pandas DataFrame and calling to_csv(). The DuckDB database file megacity_pipeline.duckdb is already saved to disk automatically since you connected to a file-based database.

  • Run the cell.
  • Confirm the file was created by running this in a new cell:
# Verify the CSV file exists and check the first few rows
import os
print(f"CSV file size: {os.path.getsize('megacity_air_quality_features.csv') / 1024:.1f} KB")

# Quick preview of the exported data
import pandas as pd
preview = pd.read_csv('megacity_air_quality_features.csv')
print(f"Rows: {len(preview)}, Columns: {list(preview.columns)}")
print(preview.head())
  • Run the cell. You should see the file size, row count, and columns matching station_id, city, latitude, longitude, road_length_500m, green_area_ratio_1km, dist_to_nearest_highway_m, industrial_area_500m.

File not found error?

The CSV is written to the same directory where your Jupyter notebook is running. If you changed your working directory at any point, the file may have been saved elsewhere. Run os.getcwd() to check your current directory.

Help me find my exported CSV file

Create an interactive map with Folium

A table of numbers is hard to interpret at global scale. An interactive map lets you visually verify that stations are in the right locations and that feature values make geographic sense.

  • Create a world map with color-coded station markers by running this code in a new cell:
import folium

# Define a color for each city
city_colors = {
    'Tokyo': 'red',
    'Delhi': 'blue',
    'Shanghai': 'green',
    'Dhaka': 'purple',
    'Sao Paulo': 'orange',
    'Mexico City': 'darkred',
    'Cairo': 'cadetblue',
    'Beijing': 'darkgreen',
    'Mumbai': 'darkpurple',
    'Osaka': 'pink'
}

# Fetch all station features from DuckDB
stations_df = con.execute("SELECT * FROM station_features").fetchdf()

# Create a world map centered to show all 10 cities
m = folium.Map(location=[20, 0], zoom_start=2)

# Add a circle marker for each station
for _, row in stations_df.iterrows():
    popup_text = (
        f"Station: {row['station_id']}<br>"
        f"City: {row['city']}<br>"
        f"Road length (500m): {row['road_length_500m']:.1f}<br>"
        f"Green ratio (1km): {row['green_area_ratio_1km']:.3f}<br>"
        f"Dist to highway (m): {row['dist_to_nearest_highway_m']:.1f}<br>"
        f"Industrial area (500m): {row['industrial_area_500m']:.1f}"
    )

    folium.CircleMarker(
        location=[row['latitude'], row['longitude']],
        radius=6,
        color=city_colors.get(row['city'], 'gray'),
        fill=True,
        fill_opacity=0.7,
        popup=folium.Popup(popup_text, max_width=300)
    ).add_to(m)

# Save the map as an HTML file
m.save('megacity_stations_map.html')
print("Map saved as megacity_stations_map.html")

What does this code do?

Folium is a Python library that generates interactive Leaflet maps as HTML files. Each station gets a CircleMarker colored by its city.

The popup for each marker displays all four spatial features you computed. Clicking any marker on the map reveals its station ID, city name, and computed feature values.

  • Run the cell.
  • Open megacity_stations_map.html in your browser. You should see circle markers clustered in 10 distinct city regions across Asia, South America, Africa, and the Middle East.
  • Click a marker to verify the popup shows the station's spatial features.

Map is blank or markers are missing?

Make sure you opened the .html file in a browser, not in VS Code's text editor. Right-click the file in the Explorer sidebar and select Reveal in File Explorer (Windows) or Reveal in Finder (macOS), then double-click to open in your browser.

If markers appear but are all in one spot, check that your latitude and longitude columns were not swapped when creating the air_quality_stations table in Step 2.

Help me debug my Folium map

✔️ Awesome, I've got everything!

Great. Double check you have saved your notebook.

ⓧ I'd like to double check the full code

Here is the complete code added in this step. Make sure your cells match:

# --- Cell 1: Summary statistics ---
summary_query = """
SELECT
    city,
    AVG(road_length_500m) AS avg_road_length,
    AVG(green_area_ratio_1km) AS avg_green_ratio,
    AVG(dist_to_nearest_highway_m) AS avg_highway_dist
FROM station_features
GROUP BY city
ORDER BY city
"""

summary_df = con.execute(summary_query).fetchdf()
print(summary_df)

# --- Cell 2: Export to CSV ---
con.execute("""
COPY station_features TO 'megacity_air_quality_features.csv' (HEADER, DELIMITER ',')
""")

print("CSV exported: megacity_air_quality_features.csv")

# --- Cell 3: Verify CSV ---
import os
print(f"CSV file size: {os.path.getsize('megacity_air_quality_features.csv') / 1024:.1f} KB")

import pandas as pd
preview = pd.read_csv('megacity_air_quality_features.csv')
print(f"Rows: {len(preview)}, Columns: {list(preview.columns)}")
print(preview.head())

# --- Cell 4: Folium map ---
import folium

city_colors = {
    'Tokyo': 'red',
    'Delhi': 'blue',
    'Shanghai': 'green',
    'Dhaka': 'purple',
    'Sao Paulo': 'orange',
    'Mexico City': 'darkred',
    'Cairo': 'cadetblue',
    'Beijing': 'darkgreen',
    'Mumbai': 'darkpurple',
    'Osaka': 'pink'
}

stations_df = con.execute("SELECT * FROM station_features").fetchdf()

m = folium.Map(location=[20, 0], zoom_start=2)

for _, row in stations_df.iterrows():
    popup_text = (
        f"Station: {row['station_id']}<br>"
        f"City: {row['city']}<br>"
        f"Road length (500m): {row['road_length_500m']:.1f}<br>"
        f"Green ratio (1km): {row['green_area_ratio_1km']:.3f}<br>"
        f"Dist to highway (m): {row['dist_to_nearest_highway_m']:.1f}<br>"
        f"Industrial area (500m): {row['industrial_area_500m']:.1f}"
    )

    folium.CircleMarker(
        location=[row['latitude'], row['longitude']],
        radius=6,
        color=city_colors.get(row['city'], 'gray'),
        fill=True,
        fill_opacity=0.7,
        popup=folium.Popup(popup_text, max_width=300)
    ).add_to(m)

m.save('megacity_stations_map.html')
print("Map saved as megacity_stations_map.html")

Your pipeline is complete. You now have a validated, ML-ready CSV dataset and an interactive map proving that every station and feature is in the right place.

Secret mission

Add Temporal Air Quality Features

Your spatial features capture the environment around each station. But pollution is seasonal. Delhi's winter smog is legendary, and Tokyo's summers bring photochemical haze. Can you pull historical monthly PM2.5 data and compute temporal features that reveal seasonal pollution patterns across all 10 cities?

Clean Up Your Resources

Clean Up Your Resources

Decide whether to keep your resources running, pause them to come back later, or delete them entirely.

This project runs entirely locally, so there are no ongoing costs. All files are stored on your machine and your OpenAQ API key is on the free tier.

Resources you used:

  • Jupyter notebook (air_quality_pipeline.ipynb).
  • DuckDB database file (megacity_pipeline.duckdb).
  • CSV exports (megacity_air_quality_features.csv and megacity_air_quality_features_enriched.csv).
  • Folium HTML map (megacity_stations_map.html).
  • OpenAQ API key (free tier, no charges).
  • Python packages installed via pip (duckdb, geopandas, osmnx, requests, folium, shapely, jupyter).

Keep everything running

No action needed. Choose this if you're still actively building or plan to use the DuckDB file and CSVs in a follow-up project.

  • All files remain on your machine at no cost.
  • The megacity_pipeline.duckdb file persists for future reuse.
  • Your OpenAQ API key remains active for future projects.
  • The installed Python packages stay available in your environment.

Pause - I'll come back to this later

Nothing to pause. All resources are local files with no running processes or ongoing costs.

  • Close the Jupyter notebook in VS Code to free up memory.
  • All files stay on disk until you're ready to return.

Delete - I don't want to use this again

Remove all project files from your machine. Your OpenAQ API key remains active and free. There is no need to delete it.

  • Close the Jupyter notebook and any running kernels in VS Code.
  • Delete the project files by running this command in your terminal from the project directory:

macOS

rm air_quality_pipeline.ipynb megacity_pipeline.duckdb megacity_air_quality_features.csv megacity_air_quality_features_enriched.csv megacity_stations_map.html

Windows

del air_quality_pipeline.ipynb, megacity_pipeline.duckdb, megacity_air_quality_features.csv, megacity_air_quality_features_enriched.csv, megacity_stations_map.html
  • Uninstall the Python packages if you no longer need them:
pip uninstall duckdb geopandas osmnx requests folium shapely jupyter -y
  • Confirm the files are gone by running ls (macOS) or dir (Windows) in your project directory. You should no longer see the notebook, database, CSV, or HTML files.

Nice Work!

Nice Work!

You just built a production-quality geospatial feature engineering pipeline that pulls air quality data from 10 of the world's largest megacities and computes spatial features using real SQL.

You've learned how to:

  • Built a multi-city data pipeline that ingests air quality station data from the OpenAQ API and downloads geospatial features (roads, green spaces, industrial zones) from OpenStreetMap using OSMnx.
  • Wrote real spatial SQL queries in DuckDB using functions like ST_Buffer, ST_Distance, ST_Intersects, ST_Area, and ST_Length to compute meaningful features around each air quality monitor.
  • Exported a clean, ML-ready dataset with spatial features like road density, green space ratio, highway distance, and industrial proximity for every station across all 10 cities.
  • Secret Mission: Added temporal air quality patterns by pulling monthly PM2.5 averages and computing seasonal pollution features to create an enriched spatiotemporal dataset.

Ready to quiz yourself?