Last-Mile Delivery Zone Optimizer

Build a Dockerized GIS pipeline to cluster Warsaw delivery zones with PostGIS.

Introduction

30 Second Summary

Every time you order food or a package for same-day delivery, a logistics system decides which driver covers your neighborhood. Companies split cities into delivery zones to minimize driving time and balance workload across their fleet.

In this project, you will build a Python application that downloads real commercial locations from OpenStreetMap, stores them in a spatial database, clusters them into optimized delivery zones, and generates an interactive map of Warsaw's delivery network.

What You'll Build

Picture opening an HTML file and seeing an interactive map of Warsaw with five color-coded zones, each filled with clickable delivery points and a legend showing logistics metrics for every zone.

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

  • A PostGIS spatial database running in Docker, loaded with real Warsaw shop and restaurant locations fetched from the Overpass API.
  • A K-Means clustering pipeline that partitions delivery points into optimized zones and calculates distance metrics for each one.
  • An interactive Folium map with zone boundaries, color-coded delivery markers, and an analytics legend you can open in any browser.
  • Secret Mission: Add a depot placement optimizer that calculates the ideal warehouse location based on demand density and visualizes coverage radius on the map.

Are there any prerequisites?

This project assumes you have Python 3.12, Docker Desktop, and VS Code installed on your machine.

You should also have basic familiarity with Python and SQL. No prior experience with geospatial tools is needed.

Before We Start

Before we dive in, let's lock in what you're about to build and why it matters.

Set Up Docker and Python Environment

A delivery zone optimizer needs two foundations: a spatial database that understands geography, and a Python environment loaded with geospatial libraries.

In this step, you'll spin up a PostGIS database inside a Docker container, then build a Python environment with all the packages the project needs. By the end, your Python code will talk directly to PostGIS and confirm the connection is live.

In this step, get ready to:
  • Create a Docker Compose file and start the PostGIS container.
  • Set up a Python virtual environment and install all dependencies.
  • Create a configuration file and verify database connectivity from Python.
Create the project folder and Docker Compose file

Docker Compose lets you define a database server as a simple YAML file. When you run docker compose up, Docker reads that file and launches the container for you.

  • Make sure Docker Desktop is running in the background.
  • In VS Code, open the integrated terminal by pressing Ctrl+` (Windows) or Cmd+` (macOS).
  • Move to your Desktop by running this command:
cd ~/Desktop
  • Create your project folder and move into it by running these commands:
mkdir delivery-zone-optimizer
cd delivery-zone-optimizer
  • Confirm the folder exists by running:
pwd

You should see a path ending in delivery-zone-optimizer.

  • Create the Docker Compose file by running:
touch docker-compose.yml
  • Open docker-compose.yml in VS Code and paste the following configuration:
services:
  postgis:
    image: postgis/postgis:17-3.6-alpine
    environment:
      POSTGRES_DB: logistics_gis
      POSTGRES_USER: gisuser
      POSTGRES_PASSWORD: gispassword
    ports:
      - "5432:5432"
    volumes:
      - pgdata:/var/lib/postgresql/data

volumes:
  pgdata:

What does this configuration do?

This file defines a single service called postgis that runs the postgis/postgis:17-3.6-alpine image. That image bundles PostgreSQL 17 with PostGIS 3.6 on a lightweight Alpine Linux base.

  • The environment block sets the database name, username, and password that PostgreSQL creates on first launch.
  • The ports mapping exposes PostgreSQL on localhost:5432 so your Python scripts can reach it.
  • The pgdata volume persists your data across container restarts.
  • Save docker-compose.yml.
  • Start the PostGIS container by running:
docker compose up -d

The -d flag runs the container in the background so your terminal stays free.

  • Verify the container is running:
docker compose ps

You should see one container with a status of running and port 5432 mapped.

Container not starting?

  • If you see a port conflict error, another process is using port 5432. Stop any local PostgreSQL service, then re-run docker compose up -d.
  • If Docker reports "Cannot connect to the Docker daemon", make sure Docker Desktop is running.
  • help me troubleshoot my Docker Compose PostGIS container not starting
  • Confirm PostGIS is active inside the container by running:
docker compose exec postgis psql -U gisuser -d logistics_gis -c "SELECT PostGIS_Version();"

You should see a version number like 3.6 printed. Your spatial database is live.

✔️ Awesome, I've got everything!

Great. Make sure you've saved docker-compose.yml and the container is running.

ⓧ I'd like to double check the full code

services:
  postgis:
    image: postgis/postgis:17-3.6-alpine
    environment:
      POSTGRES_DB: logistics_gis
      POSTGRES_USER: gisuser
      POSTGRES_PASSWORD: gispassword
    ports:
      - "5432:5432"
    volumes:
      - pgdata:/var/lib/postgresql/data

volumes:
  pgdata:
Set up the Python virtual environment

A virtual environment keeps your project's packages isolated from other Python projects on your machine. This prevents version conflicts between different projects.

  • Create a virtual environment by running:
python3 -m venv venv
  • Activate the virtual environment:

macOS

source venv/bin/activate

Windows

venv\Scripts\activate

You should see (venv) appear at the start of your terminal prompt. This confirms the virtual environment is active.

Don't see (venv) in your prompt?

  • On Windows, you may need to allow script execution first. Run Set-ExecutionPolicy -Scope CurrentUser RemoteSigned and then retry the activate command.
  • Make sure you're in the delivery-zone-optimizer folder where you created the venv folder.
  • help me activate my Python virtual environment
  • Create the requirements file by running:
touch requirements.txt
  • Open requirements.txt in VS Code and paste these dependencies:
osmnx==2.1.0
geopandas==1.1.4
folium==0.20.0
sqlalchemy==2.0.51
geoalchemy2==0.20.0
psycopg2-binary==2.9.12
scikit-learn==1.9.0
shapely==2.1.2

What are these packages?

  • osmnx downloads street networks and points of interest from OpenStreetMap.
  • geopandas extends pandas DataFrames with geometry columns for spatial data.
  • folium generates interactive Leaflet maps as HTML files.
  • sqlalchemy and geoalchemy2 let Python talk to PostGIS using an ORM pattern.
  • scikit-learn provides the K-Means clustering algorithm for zone optimization.
  • shapely handles geometric operations like computing convex hulls and centroids.
  • Save requirements.txt.
  • Install all dependencies by running:
pip install -r requirements.txt

This will download and install all eight packages plus their dependencies. It may take a minute or two.

Installation failed?

  • If you see errors about missing C libraries (especially for shapely or geopandas), make sure your (venv) prefix is visible. Installing without the venv active may cause permission errors.
  • On macOS, if psycopg2-binary fails, ensure Xcode Command Line Tools are installed by running xcode-select --install.
  • help me fix my pip install errors for the delivery zone optimizer project

✔️ Awesome, I've got everything!

Great. Make sure you've saved requirements.txt and all packages installed without errors.

ⓧ I'd like to double check the full code

osmnx==2.1.0
geopandas==1.1.4
folium==0.20.0
sqlalchemy==2.0.51
geoalchemy2==0.20.0
psycopg2-binary==2.9.12
scikit-learn==1.9.0
shapely==2.1.2
Create the configuration file and verify connectivity

A central configuration file keeps your database connection string and project constants in one place. Every module in the project will import from this file instead of hardcoding values.

  • Create the config file by running:
touch config.py
  • Open config.py in VS Code and paste this configuration:
DATABASE_URL = "postgresql://gisuser:gispassword@localhost:5432/logistics_gis"
CITY = "Warsaw, Poland"
NUM_ZONES = 5

What do these constants control?

  • DATABASE_URL is the connection string matching the credentials you set in docker-compose.yml.
  • CITY tells OSMnx which city to download data for.
  • NUM_ZONES sets how many delivery zones the K-Means algorithm will create.
  • Save config.py.

✔️ Awesome, I've got everything!

Great. Make sure you've saved config.py with all three constants.

ⓧ I'd like to double check the full code

DATABASE_URL = "postgresql://gisuser:gispassword@localhost:5432/logistics_gis"
CITY = "Warsaw, Poland"
NUM_ZONES = 5

Now let's prove the full chain works: Python connects through SQLAlchemy to the PostGIS database running in Docker.

Before you run the next command, what do you expect to see if everything is wired up correctly?

  • Test the database connection from Python by running:
python -c "from sqlalchemy import create_engine, text; from config import DATABASE_URL; engine = create_engine(DATABASE_URL); conn = engine.connect(); print(conn.execute(text('SELECT PostGIS_Version()')).scalar()); conn.close()"

You should see a PostGIS version number (like 3.6) printed in your terminal. This confirms Python can reach the PostGIS database through SQLAlchemy using your config.

Seeing a connection refused error?

  • Run docker compose ps to confirm the container is still running. If it stopped, run docker compose up -d again.
  • Make sure your (venv) is active. Without the virtual environment, Python won't find the sqlalchemy package.
  • Verify that config.py uses the exact same credentials as docker-compose.yml (database logistics_gis, user gisuser, password gispassword).
  • help me fix my Python connection to PostGIS

Your infrastructure is ready. You have a spatial database running in Docker and a Python environment with every package the project needs. Next up, you'll pull real delivery location data from OpenStreetMap into this database.

Fetch Warsaw Data from OpenStreetMap

Your PostGIS database is running and Python can connect to it. Now you need real delivery data to work with.

In this step, you will pull real-world commercial locations from OpenStreetMap using the OSMnx library. These locations (shops, restaurants, cafes, pharmacies) will serve as delivery destination points for your zone optimizer.

In this step, get ready to:
  • Download Warsaw's street network graph from OpenStreetMap.
  • Fetch commercial points of interest as delivery destinations.
  • Verify the downloaded data contains real Warsaw locations with coordinates.
Download Warsaw's street network

OSMnx talks to the Overpass API behind the scenes to download geographic data from OpenStreetMap. The first thing you need is Warsaw's street network, which represents all the roads a delivery driver could travel on.

  • Create a new file called fetch_data.py in your project folder.
  • Add the imports and the street network function by pasting this code:
import osmnx as ox
import geopandas as gpd
from config import CITY


def fetch_street_network():
    """Download Warsaw street network from OSM via Overpass API."""
    print(f"Downloading street network for {CITY}...")
    G = ox.graph.graph_from_place(CITY, network_type="drive")
    print(f"Network downloaded: {len(G.nodes)} nodes, {len(G.edges)} edges")
    return G

What does this code do?

import osmnx as ox loads the OSMnx library, which wraps the OpenStreetMap Overpass API into simple Python calls.

ox.graph.graph_from_place(CITY, network_type="drive") geocodes "Warsaw, Poland" to find its boundary polygon, then downloads every driveable road within that boundary as a network graph. The network_type="drive" filter means you only get roads that cars can use, which is exactly what a delivery fleet needs.

  • Save fetch_data.py.
  • Confirm the file was created by running:
ls fetch_data.py

You should see fetch_data.py listed in the output.

Fetch commercial POIs as delivery points

A delivery zone optimizer needs delivery destinations. You will use OSMnx to download real commercial locations tagged in OpenStreetMap: shops, restaurants, cafes, and pharmacies.

  • Add the fetch_delivery_points() function below fetch_street_network() in fetch_data.py:
def fetch_delivery_points():
    """Download commercial POIs as delivery destinations from OSM."""
    print(f"Downloading commercial POIs for {CITY}...")
    tags = {"shop": True, "amenity": ["restaurant", "cafe", "pharmacy"]}
    pois = ox.features.features_from_place(CITY, tags=tags)

    # Filter to point geometries only
    pois = pois[pois.geometry.type == "Point"].copy()
    pois = pois[["name", "shop", "amenity", "geometry"]].copy()
    pois["category"] = pois["shop"].fillna(pois["amenity"])
    pois = pois[["name", "category", "geometry"]].reset_index(drop=True)

    print(f"Found {len(pois)} delivery points")
    return pois

What does this code do?

  • tags = {"shop": True, "amenity": ["restaurant", "cafe", "pharmacy"]} tells OSMnx to fetch ALL shops (any type) plus specifically restaurants, cafes, and pharmacies.
  • ox.features.features_from_place(CITY, tags=tags) downloads matching features within Warsaw's boundary as a GeoDataFrame.
  • The filter pois[pois.geometry.type == "Point"] keeps only point geometries, discarding polygons (building outlines) that would not work as delivery coordinates.
  • The category column merges the shop and amenity columns into one label so every point has a single category string.
  • Save fetch_data.py.
Run the script and verify the data

Now you will add a main block that calls both functions and prints sample output so you can confirm everything works.

  • Add this main block at the bottom of fetch_data.py:
if __name__ == "__main__":
    G = fetch_street_network()
    pois = fetch_delivery_points()
    print(pois.head())

What does this code do?

The if __name__ == "__main__" guard means this code only runs when you execute the file directly. It calls both fetch functions in sequence, then prints the first 5 rows of the delivery points GeoDataFrame so you can inspect the data.

  • Save fetch_data.py.

Before you run this, what do you expect to see? The script downloads real geographic data from the internet, so think about what kind of output a street network and a list of commercial locations would produce.

  • Run the data fetch script:
python fetch_data.py

The download may take 30-60 seconds depending on your internet connection. You should see output showing the street network statistics (thousands of nodes and edges) followed by a count of delivery points and a table of the first 5 rows with columns for name, category, and geometry.

Why are my numbers different?

OpenStreetMap is a live, community-edited map. The exact node/edge counts and POI totals will vary as contributors add and update data. As long as you see thousands of nodes/edges and hundreds (or thousands) of delivery points, your script is working correctly.

Seeing a timeout or connection error?

  • The Overpass API can be slow during peak hours. Wait a minute and try running the script again.
  • Check your internet connection is active.
  • If you see a geocoding error, make sure CITY = "Warsaw, Poland" in your config.py is spelled exactly as shown.

Help me debug a connection or timeout error when running fetch_data.py with OSMnx.

✔️ Awesome, I've got everything!

Great. Double check you have saved fetch_data.py and your terminal showed the street network stats and delivery points table.

ⓧ I'd like to double check the full code

Your fetch_data.py has three sections built in order. Scroll up and confirm each is present:

  • The imports block (import osmnx as ox, import geopandas as gpd, from config import CITY) at the top of the file.
  • The fetch_street_network() function that calls ox.graph.graph_from_place() with network_type="drive".
  • The fetch_delivery_points() function that uses ox.features.features_from_place() and filters to point geometries.
  • The if __name__ == "__main__" block at the bottom calling both functions and printing pois.head().

You now have real Warsaw delivery data flowing into your pipeline. Next up, you will store these points in your PostGIS database so you can run spatial queries against them.

Store Data in PostGIS

In the previous step, you downloaded real delivery locations from OpenStreetMap into a Python GeoDataFrame. That's great for quick exploration, but a GeoDataFrame lives only in memory.

To persist your data, share it between scripts, and run powerful spatial queries like "find all points within 2km," you need a spatial database. PostGIS extends PostgreSQL with geometry types and spatial functions. In this step, you'll define your database schema, load the delivery points into PostGIS, and run your first spatial query.

In this step, get ready to:
  • Define SQLAlchemy models with PostGIS geometry columns.
  • Load delivery point data from your GeoDataFrame into PostGIS.
  • Run a spatial proximity query using ST_DWithin.
Create database models with geometry columns

Your database needs two tables: one for individual delivery points (with a POINT geometry) and one for delivery zones (with a POLYGON geometry). SQLAlchemy provides the ORM layer, and GeoAlchemy2 adds the spatial column types.

  • Create a new file called database.py in your project folder.
  • Add the imports, engine, session factory, and the DeliveryPoint model by pasting this code:
from sqlalchemy import Column, Integer, String, Float, create_engine
from sqlalchemy.orm import DeclarativeBase, sessionmaker
from geoalchemy2 import Geometry
from config import DATABASE_URL

# Create database engine and session factory
engine = create_engine(DATABASE_URL)
Session = sessionmaker(bind=engine)


class Base(DeclarativeBase):
    pass


# Delivery point: a single commercial location
class DeliveryPoint(Base):
    __tablename__ = "delivery_points"

    id = Column(Integer, primary_key=True)
    name = Column(String, nullable=True)
    category = Column(String, nullable=True)
    geom = Column(Geometry("POINT", srid=4326))
    zone = Column(Integer, nullable=True)

What does this code do?

The engine connects to your PostGIS container using the connection string from config.py. Session creates database sessions for running queries.

DeclarativeBase is SQLAlchemy 2.0's modern base class for ORM models. Each model class maps to a database table.

The geom column uses GeoAlchemy2's Geometry("POINT", srid=4326) type. SRID 4326 is the standard coordinate system for GPS latitude/longitude.

  • Now add the DeliveryZone model and the init_db() function below the code you just pasted:
# Delivery zone: a cluster of delivery points
class DeliveryZone(Base):
    __tablename__ = "delivery_zones"

    id = Column(Integer, primary_key=True)
    zone_id = Column(Integer, unique=True)
    centroid = Column(Geometry("POINT", srid=4326))
    area = Column(Geometry("POLYGON", srid=4326))
    num_points = Column(Integer)
    avg_distance_km = Column(Float)


def init_db():
    Base.metadata.create_all(engine)
    print("Database tables created successfully.")


if __name__ == "__main__":
    init_db()

What does the zone model store?

DeliveryZone stores each zone's boundary as a POLYGON, its center as a POINT, and summary statistics (number of points, average distance). You'll populate this table in the next step when you run clustering.

init_db() calls Base.metadata.create_all(engine) which tells SQLAlchemy to create all tables that don't already exist in the database.

Before you run this, what do you think the output will be?

  • Save database.py and run it to create the tables:
python database.py

You should see:

Database tables created successfully.
  • Verify the tables exist directly in PostGIS by running:
docker exec -it $(docker ps -q) psql -U gisuser -d logistics_gis -c "\dt"

You should see delivery_points and delivery_zones listed as tables in the output.

Don't see your tables?

  • Make sure your Docker container is still running. Run docker ps to confirm.
  • Check that config.py has the correct DATABASE_URL matching your container credentials.
  • If you see a connection refused error, wait a few seconds for the container to finish starting up and try again.

Why can't I connect to my PostGIS database from Python?

✔️ Awesome, I've got everything!

Great. Make sure you've saved database.py and confirmed both tables appear in the \dt output.

ⓧ I'd like to double check the full code

Your database.py should have these sections in order. Scroll up and confirm each is present:

  • Imports (sqlalchemy, geoalchemy2, config) and engine/Session setup.
  • Base(DeclarativeBase) class.
  • DeliveryPoint model with id, name, category, geom, zone columns.
  • DeliveryZone model with id, zone_id, centroid, area, num_points, avg_distance_km columns.
  • init_db() function and if __name__ == "__main__" block.
Load delivery points into PostGIS

Now that your tables exist, you need a function to take the GeoDataFrame from fetch_data.py and insert each point into PostGIS. The key challenge is converting Shapely geometry objects into a format PostGIS can store.

  • Create a new file called analytics.py in your project folder.
  • Add the imports and the load_points_to_db() function by pasting this code:
import numpy as np
from sklearn.cluster import KMeans
from shapely.geometry import Point, MultiPoint
from geoalchemy2.shape import from_shape, to_shape
from sqlalchemy import text
from database import engine, Session, DeliveryPoint, DeliveryZone, init_db
from config import NUM_ZONES


def load_points_to_db(pois_gdf):
    """Load delivery points GeoDataFrame into PostGIS."""
    init_db()
    session = Session()

    # Clear existing points to avoid duplicates
    session.query(DeliveryPoint).delete()
    session.commit()

    for idx, row in pois_gdf.iterrows():
        point = DeliveryPoint(
            name=row.get("name", f"Point_{idx}"),
            category=row.get("category", "unknown"),
            geom=from_shape(row.geometry, srid=4326),
        )
        session.add(point)

    session.commit()
    count = session.query(DeliveryPoint).count()
    print(f"Loaded {count} delivery points into PostGIS")
    session.close()

What does from_shape do?

from_shape() converts a Shapely geometry object (like a Point with lat/lon) into a WKBElement that PostGIS can store. You pass srid=4326 to tag the geometry with its coordinate system.

The function iterates over each row in the GeoDataFrame, wraps the geometry with from_shape(), and adds it to the session. A single session.commit() at the end writes all points to the database in one transaction.

  • Save analytics.py and test the load function by running:
python -c "from fetch_data import fetch_delivery_points; from analytics import load_points_to_db; pois = fetch_delivery_points(); load_points_to_db(pois)"

This fetches the Warsaw POI data and loads it into your PostGIS database. You should see output showing the number of delivery points downloaded, followed by a confirmation line like:

Loaded 4523 delivery points into PostGIS

Your exact number will vary depending on current OpenStreetMap data for Warsaw.

Seeing an import error?

  • If you see a ModuleNotFoundError, make sure your virtual environment is activated before running the command.
  • If you see a connection error, confirm the Docker container is running with docker ps.

Why is my data not loading into PostGIS?

Run a spatial proximity query

The real power of PostGIS is spatial queries. You can ask the database questions like "find all delivery points within 2 kilometers of Warsaw's center" and get results sorted by distance.

  • In analytics.py, add the spatial_query_example() function below load_points_to_db():
def spatial_query_example():
    """Demonstrate PostGIS spatial queries - points near Warsaw center."""
    with engine.connect() as conn:
        result = conn.execute(text("""
            SELECT name, category,
                   ST_Distance(
                       geom::geography,
                       ST_SetSRID(ST_MakePoint(21.0122, 52.2297), 4326)::geography
                   ) / 1000 AS dist_km
            FROM delivery_points
            WHERE ST_DWithin(
                geom::geography,
                ST_SetSRID(ST_MakePoint(21.0122, 52.2297), 4326)::geography,
                2000
            )
            ORDER BY dist_km
            LIMIT 10;
        """))

        print("\n--- Points within 2km of Warsaw Center ---")
        for row in result:
            print(f"  {row.name} ({row.category}): {row.dist_km:.2f} km")

What do these PostGIS functions do?

  • ST_MakePoint(21.0122, 52.2297) creates a point at Warsaw's center (longitude, latitude).
  • ST_DWithin(..., 2000) filters to points within 2000 meters of that center. The ::geography cast tells PostGIS to calculate real-world distances in meters.
  • ST_Distance computes the actual distance in meters, then divides by 1000 for kilometers.

What do you think you'll see when you run this query against your loaded data?

  • Save analytics.py and run the spatial query:
python -c "from analytics import spatial_query_example; spatial_query_example()"

You should see a list of 10 delivery points sorted by distance from Warsaw's center, each showing the location name, category, and distance in kilometers.

--- Points within 2km of Warsaw Center ---
  Biedronka (convenience): 0.12 km
  Starbucks (cafe): 0.34 km
  Apteka (pharmacy): 0.51 km
  ...

Your exact results will differ based on current OpenStreetMap data for Warsaw.

Query returning no results?

  • Make sure you ran the load step first. Run the full pipeline in one command if needed.
  • Check that the coordinates are in the correct order. PostGIS uses ST_MakePoint(longitude, latitude), not latitude first.

Why is my PostGIS spatial query returning empty results?

  • As a final check, verify the data count directly in PostGIS by running:
docker exec -it $(docker ps -q) psql -U gisuser -d logistics_gis -c "SELECT COUNT(*) FROM delivery_points;"

You should see a count matching the number printed during the load step.

✔️ Awesome, I've got everything!

Great. Make sure you've saved both database.py and analytics.py.

ⓧ I'd like to double check the full code

Your database.py should have these sections in order. Scroll up and confirm each is present:

  • Imports (sqlalchemy, geoalchemy2, config) and engine/Session setup.
  • Base(DeclarativeBase) class.
  • DeliveryPoint model with id, name, category, geom, zone columns.
  • DeliveryZone model with id, zone_id, centroid, area, num_points, avg_distance_km columns.
  • init_db() function and if __name__ == "__main__" block.

Your analytics.py should have these sections in order:

  • Imports (numpy, sklearn, shapely, geoalchemy2, sqlalchemy, database, config).
  • load_points_to_db() function that clears existing points and inserts from a GeoDataFrame.
  • spatial_query_example() function with the ST_DWithin query.

Your delivery data is now living in a real spatial database. Next up, you'll use K-Means clustering to divide these points into optimized delivery zones.

Build the Spatial Analytics Module

You have thousands of delivery points stored in PostGIS, but a flat list of locations is not useful for logistics planning. A delivery company needs to partition its service area into zones so that each driver covers a manageable, geographically compact region.

In this step, you will build the analytics engine that transforms raw coordinates into actionable delivery zones. You will use K-Means clustering to group nearby delivery points into 5 zones, compute zone boundaries, and run PostGIS spatial queries to prove your database can answer location-aware questions.

In this step, get ready to:
  • Cluster delivery points into 5 optimized zones using K-Means.
  • Compute convex hull boundaries and zone statistics for each zone.
  • Run PostGIS spatial queries to find delivery points near Warsaw center.
Create the analytics module with data loading

Before running clustering, you need a function to load your fetched GeoDataFrame into PostGIS in a clean, repeatable way. This ensures your analytics pipeline always starts from fresh data.

  • Create a new file called analytics.py in your project folder.
  • Add the imports and load_points_to_db() function by pasting this code:
import numpy as np
from sklearn.cluster import KMeans
from shapely.geometry import Point, MultiPoint
from geoalchemy2.shape import from_shape, to_shape
from sqlalchemy import text
from database import engine, Session, DeliveryPoint, DeliveryZone, init_db
from config import NUM_ZONES


def load_points_to_db(pois_gdf):
    """Load delivery points GeoDataFrame into PostGIS."""
    init_db()
    session = Session()

    # Clear existing points for a clean reload
    session.query(DeliveryPoint).delete()
    session.commit()

    for idx, row in pois_gdf.iterrows():
        point = DeliveryPoint(
            name=row.get("name", f"Point_{idx}"),
            category=row.get("category", "unknown"),
            geom=from_shape(row.geometry, srid=4326),
        )
        session.add(point)

    session.commit()
    count = session.query(DeliveryPoint).count()
    print(f"Loaded {count} delivery points into PostGIS")
    session.close()

What does this code do?

The imports pull in all the tools you need: numpy for array math, KMeans for clustering, Shapely for geometry objects, and GeoAlchemy2 helpers for converting between Python shapes and PostGIS geometries.

The load_points_to_db() function clears any existing data first (so re-runs are safe), then loops through the GeoDataFrame and inserts each row as a DeliveryPoint record with its geometry converted via from_shape().

  • Verify the import works by running:
python -c "from analytics import load_points_to_db; print('Analytics module loaded successfully')"

You should see Analytics module loaded successfully printed with no import errors.

Seeing an ImportError?

  • Make sure your virtual environment is activated. The scikit-learn and geoalchemy2 packages must be installed.
  • Confirm database.py and config.py are in the same folder as analytics.py.

Help me debug an ImportError in my analytics.py file.

Perform K-Means clustering on delivery points

K-Means clustering partitions your delivery points into groups where each point belongs to the nearest cluster center. The algorithm iterates until the centers stabilize, producing compact geographic zones.

  • Add the first half of the perform_clustering() function below load_points_to_db() in analytics.py:
def perform_clustering():
    """Assign delivery points to zones using K-Means clustering."""
    session = Session()
    points = session.query(DeliveryPoint).all()

    # Extract lon/lat coordinates from PostGIS geometries
    coords = np.array([[to_shape(p.geom).x, to_shape(p.geom).y] for p in points])

    # Run K-Means with 5 zones, fixed seed for reproducibility
    kmeans = KMeans(n_clusters=NUM_ZONES, random_state=42, n_init=10)
    labels = kmeans.fit_predict(coords)

    # Assign each delivery point its zone ID
    for point, label in zip(points, labels):
        point.zone = int(label)
    session.commit()

What does this code do?

The function loads all delivery points from PostGIS and converts their geometries back to a coordinate array using to_shape().

Setting n_init=10 runs the algorithm 10 times with different starting points and keeps the best result. The random_state=42 ensures you get the same zones every time you run it. After fitting, each point's zone column gets updated in the database.

  • Complete the perform_clustering() function by adding this code directly below the session.commit() line (still inside the function):
    # Build zone polygons and compute statistics
    session.query(DeliveryZone).delete()
    for zone_id in range(NUM_ZONES):
        zone_points = coords[labels == zone_id]
        centroid = kmeans.cluster_centers_[zone_id]
        # Convex hull wraps a polygon around all zone points
        multi_point = MultiPoint([Point(c) for c in zone_points])
        zone_area = multi_point.convex_hull
        # Approximate distance in km (1 degree ~ 111 km)
        distances = np.sqrt(np.sum((zone_points - centroid) ** 2, axis=1))
        avg_dist = float(np.mean(distances)) * 111

        zone = DeliveryZone(
            zone_id=zone_id,
            centroid=from_shape(Point(centroid), srid=4326),
            area=from_shape(zone_area, srid=4326),
            num_points=len(zone_points),
            avg_distance_km=round(avg_dist, 2),
        )
        session.add(zone)

    session.commit()
    print(f"Clustering complete: {NUM_ZONES} zones created")

    zones = session.query(DeliveryZone).all()
    print("\n--- Zone Summary ---")
    for z in zones:
        print(f"Zone {z.zone_id}: {z.num_points} points, avg distance: {z.avg_distance_km} km")

    session.close()
    return kmeans

What does the zone statistics code do?

For each of the 5 zones, this code computes a convex hull (the smallest polygon enclosing all points in that zone) and calculates the average distance from the zone center to its member points.

The * 111 converts degrees to approximate kilometers. At Warsaw's latitude, 1 degree of latitude is roughly 111 km. This is a useful approximation for city-scale distances.

  • Save analytics.py.
  • Before you run the clustering, what do you think the zone summary will show? Will the 5 zones have roughly equal numbers of points, or will some zones be much larger?
  • Test the clustering by running:
python -c "from analytics import perform_clustering; perform_clustering()"

You should see Clustering complete: 5 zones created followed by a zone summary table listing each zone's point count and average distance.

Why might zones be uneven?

K-Means optimizes for geometric compactness, not equal point counts. Areas with dense commercial activity (like Warsaw's city center) produce zones with more points, while suburban zones cover more area but have fewer points. This reflects real-world delivery load imbalance.

Seeing a connection error?

  • Confirm your PostGIS container is still running with docker compose ps.
  • If it stopped, run docker compose up -d to restart it.
  • Make sure the delivery_points table has data from the previous step. If it is empty, re-run your data loading script first.

Help me fix a connection error when running clustering in my analytics module.

Add PostGIS spatial queries and run the full pipeline

Your clustering proves Python can do spatial analytics, but the real power of PostGIS is answering spatial questions directly in the database. A spatial query like "find all delivery points within 2km of a location" runs faster in PostGIS than filtering in Python.

  • Add the spatial_query_example() function below perform_clustering() in analytics.py:
def spatial_query_example():
    """Demonstrate PostGIS spatial queries - points near Warsaw center."""
    with engine.connect() as conn:
        result = conn.execute(text("""
            SELECT name, category,
                   ST_Distance(
                       geom::geography,
                       ST_SetSRID(ST_MakePoint(21.0122, 52.2297), 4326)::geography
                   ) / 1000 AS dist_km
            FROM delivery_points
            WHERE ST_DWithin(
                geom::geography,
                ST_SetSRID(ST_MakePoint(21.0122, 52.2297), 4326)::geography,
                2000
            )
            ORDER BY dist_km
            LIMIT 10;
        """))

        print("\n--- Points within 2km of Warsaw Center ---")
        for row in result:
            print(f"  {row.name} ({row.category}): {row.dist_km:.2f} km")

What do ST_DWithin and ST_Distance do?

ST_DWithin filters rows where the geometry is within a specified distance (2000 meters here) from a reference point. It uses a spatial index for fast lookups.

ST_Distance calculates the exact distance in meters between two geographies. Dividing by 1000 gives you kilometers. The ::geography cast tells PostGIS to use real-world meters instead of raw degree differences.

  • Add the __main__ block at the very bottom of analytics.py:
if __name__ == "__main__":
    from fetch_data import fetch_delivery_points

    pois = fetch_delivery_points()
    load_points_to_db(pois)
    perform_clustering()
    spatial_query_example()
  • Save analytics.py.
  • What do you expect the spatial query to return? Think about what types of businesses sit near Warsaw's city center.
  • Run the full analytics pipeline by running:
python analytics.py

You should see the data download, a confirmation of points loaded, the zone summary table, and a list of the 10 closest delivery points to Warsaw's center with their distances in kilometers.

Script hangs on downloading data?

  • The Overpass API download can take 30-60 seconds for a large city like Warsaw. Wait for it to complete.
  • If it times out, check your internet connection and try again. The OSM Overpass API is free but can be slow during peak hours.

Help me troubleshoot a timeout when running analytics.py in my delivery zone optimizer project.

  • Verify zone data is stored in PostGIS by running:
docker exec -it $(docker compose ps -q postgis) psql -U gisuser -d logistics_gis -c "SELECT zone_id, num_points, avg_distance_km FROM delivery_zones ORDER BY zone_id;"

You should see a table with 5 rows (zone IDs 0 through 4), each showing the number of points and average distance.

What should I see?

The exact numbers depend on the current OSM data for Warsaw, but you should see 5 zones with varying point counts (some zones with hundreds of points, others with fewer) and average distances ranging from roughly 1 to 4 km. Uneven distribution is normal and reflects real commercial density patterns.

✔️ Awesome, I've got everything!

Great. Double check you have saved analytics.py and that your zone summary printed 5 zones.

ⓧ I'd like to double check the full code

Your analytics.py file is too long to display in one block. Scroll up through your step and confirm each section is present in this order:

  • Imports (numpy, KMeans, Point, MultiPoint, from_shape, to_shape, text, database models, and NUM_ZONES).
  • load_points_to_db() function that clears existing data and loads a GeoDataFrame into PostGIS.
  • perform_clustering() function that runs K-Means, assigns zones, computes convex hulls, stores zone statistics, and prints the zone summary.
  • spatial_query_example() function that runs ST_DWithin and ST_Distance queries against Warsaw center.
  • __main__ block that calls fetch_delivery_points(), load_points_to_db(), perform_clustering(), and spatial_query_example() in sequence.

Your analytics module is fully operational. You have delivery points partitioned into zones, zone boundaries stored as polygons, and spatial queries proving PostGIS can answer proximity questions in milliseconds. Next up, you will visualize all of this on an interactive map.

Generate the Interactive Map

Your delivery zones are clustered and stored in PostGIS, but raw database rows don't communicate much to a logistics team.

In this step, you'll turn the spatial data into an interactive HTML map using Folium. Folium generates standalone Leaflet.js maps from Python, so anyone with a browser can pan, zoom, and click through your delivery zone analysis without running any code.

In this step, get ready to:
  • Create visualize.py with zone polygon overlays, centroid markers, and delivery point circles.
  • Add a color-coded analytics legend to the map.
  • Create main.py to orchestrate the full pipeline end-to-end.
Build the map with zone polygon overlays

Folium produces a self-contained HTML file you can open in any browser. You'll start by creating the map canvas and overlaying the convex hull polygons from your delivery_zones table.

  • Create a new file called visualize.py in your project folder.
  • Add the imports, color palette, and the beginning of the create_delivery_map() function by pasting this code:
import folium
from geoalchemy2.shape import to_shape
from database import Session, DeliveryPoint, DeliveryZone

# Color palette for up to 8 zones
COLORS = ["red", "blue", "green", "purple", "orange", "darkred", "darkblue", "darkgreen"]


def create_delivery_map():
    """Generate interactive Folium map with delivery zones and analytics."""
    session = Session()

    # Load all points and zones from PostGIS
    points = session.query(DeliveryPoint).all()
    zones = session.query(DeliveryZone).all()

    # Create base map centered on Warsaw with a clean light basemap
    m = folium.Map(location=[52.2297, 21.0122], zoom_start=12, tiles="cartodbpositron")

What does this code do?

  • COLORS defines a palette so each zone gets a distinct color on the map.
  • The function queries all delivery points and zones from PostGIS using the models you defined in database.py.
  • folium.Map() creates the base map at Warsaw's coordinates (52.2297 N, 21.0122 E) with the cartodbpositron tileset. This provides a minimal light background ideal for data overlays.
  • Add the zone overlay loop directly below the map creation line. This draws each zone's polygon boundary and places a star marker at the centroid:
    # Draw each zone as a semi-transparent polygon overlay
    for zone in zones:
        zone_shape = to_shape(zone.area)
        zone_color = COLORS[zone.zone_id % len(COLORS)]

        folium.GeoJson(
            zone_shape.__geo_interface__,
            style_function=lambda x, color=zone_color: {
                "fillColor": color,
                "color": color,
                "weight": 2,
                "fillOpacity": 0.15,
            },
            tooltip=f"Zone {zone.zone_id}: {zone.num_points} points, avg dist: {zone.avg_distance_km} km",
        ).add_to(m)

        # Star marker at zone centroid with popup analytics
        centroid = to_shape(zone.centroid)
        folium.Marker(
            location=[centroid.y, centroid.x],
            popup=(
                f"Zone {zone.zone_id} Center<br>"
                f"Points: {zone.num_points}<br>"
                f"Avg Distance: {zone.avg_distance_km} km"
            ),
            icon=folium.Icon(color=zone_color, icon="star"),
        ).add_to(m)

What does this code do?

  • to_shape() converts the PostGIS geometry back into a Shapely object so Folium can read it.
  • folium.GeoJson() renders the zone's convex hull polygon with semi-transparent fill (fillOpacity: 0.15). Hovering shows a tooltip with zone stats.
  • The lambda captures zone_color as a default argument. Without this trick, all zones would share the last loop iteration's color.
  • Each centroid gets a star-icon Marker with a clickable popup showing point count and average distance.
Add delivery point markers and the analytics legend

The zone boundaries show the big picture. Now you'll plot every individual delivery point as a small circle, color-matched to its zone.

  • Below the zone loop in visualize.py, add the delivery point plotting loop:
    # Plot each delivery point as a small colored circle
    for point in points:
        pt = to_shape(point.geom)
        color = COLORS[point.zone % len(COLORS)] if point.zone is not None else "gray"

        folium.CircleMarker(
            location=[pt.y, pt.x],
            radius=3,
            color=color,
            fill=True,
            fillOpacity=0.7,
            tooltip=f"{point.name} ({point.category}) - Zone {point.zone}",
        ).add_to(m)

What does this code do?

  • Each delivery point becomes a CircleMarker with radius 3 pixels. Small enough to avoid overlap, large enough to hover and click.
  • The color matches the point's assigned zone, so you can visually confirm the clustering looks sensible.
  • Points with no zone assignment fall back to gray as a safety check.
  • Add the HTML legend, save logic, and standalone entry point below the delivery point loop:
    # Build a fixed-position HTML legend showing per-zone stats
    legend_html = """
    <div style="position: fixed; bottom: 30px; left: 30px; z-index: 1000;
                background-color: white; padding: 15px; border-radius: 5px;
                border: 2px solid grey; font-size: 14px;">
        <b>Delivery Zone Analytics</b><br>
    """
    for zone in zones:
        color = COLORS[zone.zone_id % len(COLORS)]
        legend_html += (
            f'<i style="background:{color};width:12px;height:12px;'
            f'display:inline-block;margin-right:5px;"></i>'
            f"Zone {zone.zone_id}: {zone.num_points} pts, {zone.avg_distance_km} km avg<br>"
        )
    legend_html += "</div>"
    m.get_root().html.add_child(folium.Element(legend_html))

    # Save the map as a standalone HTML file
    output_path = "output/delivery_zones_warsaw.html"
    m.save(output_path)
    print(f"\nMap saved to {output_path}")
    print("Open this file in your browser to explore the interactive map!")

    session.close()
    return m


if __name__ == "__main__":
    create_delivery_map()

What does this code do?

  • The legend is raw HTML injected into the map page via m.get_root().html.add_child(). It floats in the bottom-left corner of the browser window.
  • m.save() writes the entire map (tiles, markers, polygons, legend) as a self-contained HTML file that opens in any browser without Python.
  • The if __name__ block lets you run visualize.py directly to generate the map without the full pipeline.
  • Save visualize.py.
  • Create the output directory and run the visualization by executing these commands:
mkdir -p output
python visualize.py
  • You should see output confirming the map was saved:
Map saved to output/delivery_zones_warsaw.html
Open this file in your browser to explore the interactive map!
  • Open output/delivery_zones_warsaw.html in your browser by double-clicking the file in Finder (macOS) or File Explorer (Windows).
  • You should see an interactive map of Warsaw with 5 colored polygon overlays, star markers at each zone centroid, and hundreds of small colored circles representing delivery points.

Map file is empty or not generated?

  • Make sure your Docker container is still running (docker compose ps). The script queries PostGIS for the data.
  • Confirm you ran the analytics step first. The delivery_zones table must have data or the map will be blank.
  • Check that the output directory exists. The mkdir -p output command creates it if missing.

help me debug why my Folium map is blank or not generating

✔️ Awesome, I've got everything!

Great. Double check you've saved visualize.py and the map renders in your browser.

ⓧ I'd like to double check the full code

The full visualize.py is over 30 lines, so scroll up and confirm each section is present in order:

  • Imports and the COLORS list.
  • create_delivery_map() function opening: session queries and folium.Map() creation.
  • Zone polygon overlay loop with GeoJson and centroid Marker calls.
  • Delivery point CircleMarker loop.
  • HTML legend construction, m.save(), session close, and return.
  • if __name__ == "__main__" block at the bottom.
Orchestrate the full pipeline with main.py

Each module handles one responsibility. Now you'll wire them together into a single entry point that runs the entire workflow from data fetch to map generation.

  • Create a new file called main.py in your project folder.
  • Add the imports and the first two pipeline stages:
import os
from fetch_data import fetch_street_network, fetch_delivery_points
from analytics import load_points_to_db, perform_clustering, spatial_query_example
from visualize import create_delivery_map


def main():
    print("=" * 60)
    print("  Last-Mile Delivery Zone Optimizer - Warsaw")
    print("=" * 60)

    # Fetch real data from OpenStreetMap
    print("\n[1/4] Fetching data from OpenStreetMap...")
    G = fetch_street_network()
    pois = fetch_delivery_points()

    # Load into PostGIS spatial database
    print("\n[2/4] Loading data into PostGIS...")
    load_points_to_db(pois)

What does this code do?

  • main() orchestrates the four pipeline stages in order. The banner and progress messages ([1/4], [2/4], etc.) let you track which stage is running.
  • Stage 1 fetches the OSM street network and delivery POIs. Stage 2 loads them into PostGIS.
  • Add the remaining pipeline stages and the entry point directly below:
    # Cluster and run spatial analysis
    print("\n[3/4] Running spatial analytics...")
    perform_clustering()
    spatial_query_example()

    # Generate the interactive map
    print("\n[4/4] Generating interactive map...")
    os.makedirs("output", exist_ok=True)
    create_delivery_map()

    print("\n" + "=" * 60)
    print("  Analysis complete! Check output/delivery_zones_warsaw.html")
    print("=" * 60)


if __name__ == "__main__":
    main()

What does this code do?

  • Stage 3 runs K-Means clustering and the spatial proximity query. Stage 4 creates the output directory and generates the Folium map.
  • os.makedirs("output", exist_ok=True) ensures the output directory exists without failing if it was already created.
  • The if __name__ block means running python main.py triggers the full pipeline.
  • Save main.py.

Before you run the full pipeline, what do you expect to see in the terminal output? Think about how many stages will print progress messages and what the final line will say.

  • Run the full pipeline by executing:
python main.py

What should I see?

The script prints progress through all four stages. The OSM download (stage 1) may take 1-2 minutes depending on your internet connection. You'll see network statistics, delivery point counts, clustering results, the spatial query output, and finally the map-saved confirmation.

The last lines should read:

============================================================
  Analysis complete! Check output/delivery_zones_warsaw.html
============================================================
  • Open output/delivery_zones_warsaw.html in your browser.
  • You should see an interactive map of Warsaw with 5 color-coded zone polygons, star markers at zone centroids with clickable popups, hundreds of delivery point circles, and a legend in the bottom-left corner showing per-zone statistics.

Pipeline fails partway through?

  • If it fails at stage 1 (fetching OSM data), check your internet connection. The Overpass API needs network access.
  • If it fails at stage 2 (loading into PostGIS), confirm your Docker container is running with docker compose ps.
  • If you see an ImportError, make sure your virtual environment is activated and all packages from requirements.txt are installed.

help me debug why my main.py pipeline is failing

✔️ Awesome, I've got everything!

Great. Confirm you've saved main.py and the pipeline runs end-to-end without errors.

ⓧ I'd like to double check the full code

The full main.py is over 30 lines, so scroll up and confirm each section is present in order:

  • Imports: os, fetch_street_network, fetch_delivery_points, load_points_to_db, perform_clustering, spatial_query_example, create_delivery_map.
  • main() function with the banner print and [1/4] through [4/4] stage labels.
  • Stage 1 calls fetch_street_network() and fetch_delivery_points().
  • Stage 2 calls load_points_to_db(pois).
  • Stage 3 calls perform_clustering() and spatial_query_example().
  • Stage 4 calls os.makedirs("output", exist_ok=True) then create_delivery_map().
  • Closing banner and if __name__ == "__main__" block at the bottom.

You've built a complete geospatial data pipeline from raw OpenStreetMap data all the way to an interactive browser-based visualization. Next up, the Secret Mission challenges you to compute the optimal depot location for your delivery network.

Secret mission

Add a Depot Placement Optimizer

Your 5 delivery zones are optimized, but where should the warehouse go? In this bonus challenge, you'll compute the ideal single-depot location based on demand density, calculate coverage radius statistics, and visualize the depot with distance rings on your interactive map.

Clean Up Your Resources

Clean Up Your Resources

This project runs entirely locally, so there are no ongoing cloud costs. Decide whether to keep your resources running, pause them to come back later, or delete them entirely.

Resources you used:

  • Docker container running PostGIS (uses minimal disk and RAM while running).
  • Docker volume pgdata storing the logistics_gis database.
  • Python virtual environment with installed packages (local disk only).
  • Project files on your Desktop (Python scripts, Docker Compose file, output/ folder with the HTML map).

Keep everything running

No action needed. Choose this if you're still actively building or experimenting with the analytics pipeline.

  • Your PostGIS container stays running and ready for queries.
  • You can re-run main.py anytime to regenerate the map with different parameters.
  • Edit config.py to change NUM_ZONES or CITY and re-run for different results.

Pause - I'll come back to this later

Shut down the Docker container to free up memory. Your data persists in the volume so you can resume later.

  • In your terminal, navigate to the project folder containing docker-compose.yml.
  • Stop the PostGIS container by running:
docker compose stop
  • Your database data stays in the pgdata volume. When you're ready to resume, run docker compose up -d from the same folder.
  • Deactivate your Python virtual environment if it is still active:
deactivate

Delete - I don't want to use this again

Remove all project resources and start fresh. This deletes the container, the database volume, your virtual environment, and all project files.

  • In your terminal, navigate to the project folder containing docker-compose.yml.
  • Remove the PostGIS container and its database volume by running:
docker compose down -v
  • Deactivate your virtual environment if it is still active:
deactivate
  • Move up one level and delete the project folder. On macOS/Linux:
cd ~/Desktop
rm -rf delivery-zone-optimizer
  • On Windows (PowerShell):
cd ~\Desktop
Remove-Item -Recurse -Force delivery-zone-optimizer
  • Verify the folder is gone by running ls (macOS/Linux) or dir (Windows) on your Desktop. You should no longer see the project folder.

That's a Wrap!

That's a Wrap!

Nice work. You just built a full spatial data pipeline that pulls real-world delivery locations from OpenStreetMap, clusters them into optimized logistics zones, and visualizes everything on an interactive map.

You've learned how to:

  • Build a Dockerized PostGIS spatial database and design geometry-aware models with SQLAlchemy + GeoAlchemy2 to store and query real geospatial data.
  • Fetch live location data from OpenStreetMap using OSMnx, then apply K-Means clustering to partition delivery points into optimized zones with PostGIS spatial queries for proximity analysis.
  • Generate an interactive Folium map with color-coded zone overlays, delivery markers, centroid analytics, and an HTML legend.
  • Secret Mission: Computed optimal depot placement using demand-weighted centroids and visualized 2km/5km coverage radius circles on the map.

Ready to quiz yourself?