Introduction to GeoPandas and its Python ecosystem

A talk from the OpenGeoHub Summer School 2022.

The ecosystem of packages for spatial data handling and analysis in Python is extensive and covers both vector and raster analytics from small to large distributed data. This talk covers only a small part, focusing on vector data processing with GeoPandas at its core. First, it covers what GeoPandas is and how it relates to other packages and combines them into a user-friendly API. Then it looks at what GeoPandas enables with a light introduction to PySAL Python Spatial Analysis Library for spatial statistics and modelling. The final part touches on the scalability issue and introduces how to handle parallel computation and big data using GeoPandas and Dask. All of that is illustrated by hands-on tasks on real-world data.

Fleischmann, Martin: Introduction to GeoPandas and its Python ecosystem. OpenGeoHub Summer School 2022 – KISTE project workshop, OpenGeoHub Foundation, 2022. https://doi.org/10.5446/59414

How to create a vector-based web map hosted on GitHub

This is the map we have created for the Urban Grammar AI project. It is created using open source software stack and hosted on GitHub, for free.

This post will walk you through the whole process of generation of the map, step by step, so you can create your own. It is a bit longer than usual, so a quick outline for better orientation:

  1. Vector tiles
    1. Zoom levels
    2. Tile format and compression
    3. Generate tiles
  2. Leaflet.js map
    1. Load leaflet.js
    2. Create the map
    3. Add a base map
    4. Add vector tiles
    5. Style the tiles
  3. Github Pages
  4. Further styling and interactivity
    1. Labels on top of the map
    2. Popup information
    3. Conditional formatting based on the zoom level
    4. Legend and details

By the end of this tutorial, you will be able to take your vector data and turn them into a fully-fledged vector map hosted on GitHub with no cost involved, for everyone to enjoy.

Vector tiles

Web maps are (usually) based on tiles. That means that every dataset is turned into a large set of tiny squares representing small portions of the map. And if you want to zoom in and out, you need those tiles for every zoom level. So the first step in our process will be the creation of vector tiles.

The easiest option is to use tippecanoe, a tool created by Mapbox that allows you to create vector tiles from any GeoJSON input. You should export only the columns you need and make sure that the geometry is in EPSG:4326 projection (lat/long) – that is what tippecanoe expects. Something like this will do if you have your data in geopandas.GeoDataFrame:

gdf[["signature_type", "geometry"]].to_crs(4326).to_file(
    "signatures_4326.geojson"
)

Now you have a signatures_4326.geojson with geometry and a column encoding the signature type. It means that it is time to use tippecanoe (check the installation instructions).

tippecanoe is a command-line tool, so fire up your favourite terminal, or use it within a Jupyter notebook (as we did). There are some options to set but let’s talk about the most important ones (you can check the rest in the documentation).

Zoom levels

Every map has a range of zoom levels at which the shown data is meaningful. If you’re comparing states, you don’t need to zoom to a level of a street, and if you’re showing building-level data, you don’t need a zoom to see the whole continent.

Zoom level 0 means you are zoomed out to see the whole world. It is a single tile. Zoom level 1 is a bit closer and is split into 4 tiles. Zoom level 2 has 16 tiles, level 3 has 64 and so on. Zoom level 16, which shows a street-level detail, consists of 4 294 967 296 tiles if you want to create them for the whole world. That is a lot of data, so you need to be mindful of which detail you actually need because the difference between zooming up to level 15 and up to level 16 is huge. OpenStreetMap has a nice guide on zoom levels if you are interested in more details.

For the signature geometry that represents neighbourhood-level classification, level 15 feels like the reasonable maximum, so you have the first option – -z15.

Tile format and compression

In later steps, you will use the Leaflet.VectorGrid plugin to load the tiles. To make that work, you need them as uncompressed protobuf files (not Mapbox vector tiles), which gives you another two options – -no-tile-compression' to control the compression and –output-to-directorytellingtippecanoeto use directories, instead ofmbtiles`.

We have used a few more options, but you may not need those.

Generate tiles

You have all you need to create the tiles:

tippecanoe -z15 \
           --no-tile-compression \
           --output-to-directory=tiles/ \
           signatures_4326.geojson

The command above sets the maximum zoom to 15 so you don’t create overly detailed tiles (z15), disables tile compression (--no-tile-compression) and saves the result to the tiles/ directory (--output-to-directory=tiles/), based on your GeoJSON file (signatures_4326.geojson).

This step may take a bit of time, but once it is done, you will see thousands of files in the tiles/ folder. It is time to show them on a map.

Leaflet.js map

One of the best tools for interactive maps is a leaflet.js, an open-source JavaScript library. With the right plugins, it is an ideal tool for almost any web map.

The whole map and a JavaScript code will be served from a single HTML file. Let’s start with some basic structure, with comments as a placeholder you will fill later:

<!DOCTYPE html>
<html>

<head>
    <title>Spatial Signatures in Great Britain</title>
    <meta charset="utf-8" />
    <meta name="viewport" content="width=device-width, initial-scale=1.0">
        <!-- load leaflet.js -->
</head>

<body style='margin:0'>
        <!-- div containing map -->
        <!-- specification of leaflet map -->
</body>

</html>

Load leaflet.js

To load the leaflet.js, you need to get its source script and CSS styles. Both are usually fetched over the web. You can then replace <!-- load leaflet.js --> with the following snippet:

<!-- load leaflet.js -->
<!-- load CSS -->
<link rel="stylesheet" href="https://unpkg.com/leaflet@1.7.1/dist/leaflet.css"/>

<!-- load JS -->
<script src="https://unpkg.com/leaflet@1.7.1/dist/leaflet.js"></script>

However, the basic leaflet.js is not able to work with vector tiles, so you need to load a Leaflet.VectorGrid extension on top. Just add this to the snippet:

<!-- load VectorGrid extension -->
<script src="https://unpkg.com/leaflet.vectorgrid@1.3.0/dist/Leaflet.VectorGrid.bundled.js"></script> 

At this moment, you have everything you need to create the map in your file.

Create the map

You need two things – a div object where the map will live and a bit of JavaScript specifying what the map should do. The former is easy:

<!-- div containing map -->
<div id="map" style="width: 100vw; height: 100vh; background: #fdfdfd"></div>

This div is empty, but the important bit is id="map", which will link it to the JavaScript. It also has some basic styles ensuring that it covers the whole screen (width: 100vw; height: 100vh;) and has an almost-white background (background: #fdfdfd).

The second part is the actual JavaScript generating the map.

<script>
// create map
var map = L.map('map', {
    center: [53.4107, -2.9704],
    minZoom: 6,
    maxZoom: 15,
    zoomControl: true,
    zoom: 12,
});

// other map specifications
</script>

Let’s unpack this. This creates a variable called map (var map), which is a leaflet map object (L.map). The object is linked to the HTML element (in this case div) with the id map you created above (“map”). It also has a dictionary with some basic specification, that is setting the centre of the map to Liverpool (center: [53.4107, -2.9704] where the numbers are latitude and longitude), limiting zoom levels to 6-15 (minZoom: 6 and maxZoom: 15), turning on zoom control (+- buttons) (zoomControl: true) and setting the default zoom to 12 (zoom: 12). The zoom level 6, selected as a minimum, is approximately the level on which you see the whole of Great Britain. Since spatial signatures cover only GB, there’s no need to zoom further away. Maximum zoom 15 mirrors the maximum zoom of tiles we generated earlier. You would get no tiles if you zoomed closer than that.

The map object is ready! There’s nothing there now, but that will change soon.

Add a base map

It is nice to have some base map in there before adding other layers, just to make sure everything works. You can add a default OpenStreetMap or something more subtle as a CartoDB Positron map.

// add background basemap
var mapBaseLayer = L.tileLayer(
    'https://{s}.basemaps.cartocdn.com/rastertiles/light_all/{z}/{x}/{y}.png', {
        attribution: '(C) OpenStreetMap contributors (C) CARTO'
    }
).addTo(map);

This creates a leaflet tile layer (L.tileLayer) with the URL of CartoDB raster XYZ tiles and attribution of the origin of tiles. Finally, it adds the layer to our map (.addTo(map)). Now you can actually see a real map.

Add vector tiles

Now for the key part – vector tiles. We need to define quite a few things here.

URL

Let’s first set an URL to each tile.

// get vector tiles URL
var mapUrl = "tiles/{z}/{x}/{y}.pbf";

Tiles are in the tiles/ folder where tippecanoe dumped them. They are organised into folders, following the XYZ tiling standard. That means that they are organised according to zoom level ({z}) and tile coordinates within each level ({x} and {y}). You then need to specify the URL to tiles with placeholders that will be automatically filled by leaflet.js. Here you save the URL to the variable mapUrl.

Options

Another variable you may want to specify is a dictionary with options for the tiles.

// define options of vector tiles
var mapVectorTileOptions = {
    rendererFactory: L.canvas.tile,
    interactive: true,
    attribution: '(C) Martin Fleischmann',
    maxNativeZoom: 15,
    minZoom: 6,
};

The snippet creates a variable (mapVectorTileOptions) and specifies that the tiles should be rendered using leaflet (rendererFactory: L.canvas.tile), they should allow interactivity (interactive: true) as we will need it later, it sets the attribution (attribution: '(C) Martin Fleischmann') and min and max zoom levels.

Vector tile layer

Now you can define the vector tile layer itself.

var mapPbfLayer = new L.VectorGrid.Protobuf(
    mapUrl, mapVectorTileOptions
)

This creates a new variable (mapPbfLayer) to store vector grid protobuf layer (new L.VectorGrid.Protobuf) based on the URL and tile options specified above.

Add to map

Finally, you can add the created layer to the map.

// add VectorGrid layer to map
mapPbfLayer.addTo(map);

At this point, you should be able to see your vector tiles on your map with no styles. However, to be able to use it locally, you need to start a server first. An easy way is to do it with Python. In the terminal, navigate to the folder with your html file (something like cd path/to/my_folder) and start the http server:

python -m http.server 8000

If you now open http://localhost:8000 in your browser, you should be able to open your html file (it will open automatically if it is called index.html) and see your map. It should look like this:

Style the tiles

The GeoJSON used to create these tiles has a single attribute column called signature_type. What you want to do now is to use this column to style the map to your liking.

Let’s assume that you know which RGB colours should apply to which signature type code. Let’s save those in a dictionary, so you can use them later.

// mapping of colors to signature types
const cmap = {
    "0_0": "#f2e6c7",
    "1_0": "#8fa37e",
    "3_0": "#d7a59f",
    "4_0": "#d7ded1",
    "5_0": "#c3abaf",
    "6_0": "#e4cbc8",
    "7_0": "#c2d0d9",
    "8_0": "#f0d17d",
    "2_0": "#678ea6",
    "2_1": "#94666e",
    "2_2": "#efc758",
    "9_0": "#3b6e8c",
    "9_1": "#333432",
    "9_2": "#ab888e",
    "9_3": "#c1c1c0",
    "9_4": "#bc5b4f",
    "9_5": "#a7b799",
    "9_6": "#c1c1c0",
    "9_7": "#c1c1c0",
    "9_8": "#c1c1c0"
};

This dictionary is easy left side encodes the signature type code, right side the colour as a HEX code.

Having colours, let’s define a new variable with styles and unpack it a bit.

// define styling of vector tiles
var vectorTileStyling = {
    signatures_4326: function(properties) {
        return ({
            fill: true,
                        fillColor: cmap[properties.signature_type],
            fillOpacity: 0.9,
            weight: 1,
            color: "#ffffff",
            opacity: 1.0,
        });
    }
}

You need to create a function mapping the styles to geometries based on the signature type property. For layer signatures_4326 (that matches the name of input GeoJSON), you shall define a function that can use the layer properties (function(properties)) and specify that we want to fill each polygon with a colour (fill: true), that the fill colour should use the signature type to retrieve the colour from the dictionary defined above (fillColor: cmap[properties.signature_type]) with the fill opacity 0.8 to see a bit of background through (fillOpacity: 0.9`). It also specifies the outline of each polygon as a white line of a weight 1.

Then you just need to edit the vector tile options (the mapVectorTileOptions dictionary from above) and add tile styling.

        // ...
        minZoom: 6,
        vectorTileLayerStyles: vectorTileStyling,  // this line is new
};

At this point, the basic map is done, and you can try to get it online or style it a bit more (scroll below for that).

This is the whole code of the page for reference.

<!DOCTYPE html>
<html>

<head>
    <title>Spatial Signatures in Great Britain</title>
    <meta charset="utf-8" />
    <meta name="viewport" content="width=device-width, initial-scale=1.0">

    <!-- load leaflet.js -->
    <!-- load CSS -->
    <link rel="stylesheet" href="https://unpkg.com/leaflet@1.7.1/dist/leaflet.css" />

    <!-- load JS -->
    <script src="https://unpkg.com/leaflet@1.7.1/dist/leaflet.js"></script>

    <!-- load VectorGrid extension -->
    <script src="https://unpkg.com/leaflet.vectorgrid@1.3.0/dist/Leaflet.VectorGrid.bundled.js"></script>
</head>

<body style='margin:0'>
    <!-- div containing map -->
    <div id="map" style="width: 100vw; height: 100vh; background: #fdfdfd"></div>
    <!-- specification of leaflet map -->
    <script>
        // create map
        var map = L.map('map', {
            center: [53.4107, -2.9704],
            minZoom: 6,
            maxZoom: 15,
            zoomControl: true,
            zoom: 12,
        });

        // add background basemap
        var mapBaseLayer = L.tileLayer(
            'https://{s}.basemaps.cartocdn.com/rastertiles/light_all/{z}/{x}/{y}.png', {
                attribution: '(C) OpenStreetMap contributors (C) CARTO'
            }
        ).addTo(map);

        // get vector tiles URL
        var mapUrl = "tiles/{z}/{x}/{y}.pbf";

        // mapping of colors to signature types
        const cmap = {
            "0_0": "#f2e6c7",
            "1_0": "#8fa37e",
            "3_0": "#d7a59f",
            "4_0": "#d7ded1",
            "5_0": "#c3abaf",
            "6_0": "#e4cbc8",
            "7_0": "#c2d0d9",
            "8_0": "#f0d17d",
            "2_0": "#678ea6",
            "2_1": "#94666e",
            "2_2": "#efc758",
            "9_0": "#3b6e8c",
            "9_1": "#333432",
            "9_2": "#ab888e",
            "9_3": "#c1c1c0",
            "9_4": "#bc5b4f",
            "9_5": "#a7b799",
            "9_6": "#c1c1c0",
            "9_7": "#c1c1c0",
            "9_8": "#c1c1c0"
        };

        // define styling of vector tiles
        var vectorTileStyling = {
            signatures_4326: function(properties) {
                return ({
                    fill: true,
                    fillColor: cmap[properties.signature_type],
                    fillOpacity: 0.9,
                    weight: 1,
                    color: "#ffffff",
                    opacity: 1.0,
                });
            }
        }

        // define options of vector tiles
        var mapVectorTileOptions = {
            rendererFactory: L.canvas.tile,
            interactive: true,
            attribution: '(C) Martin Fleischmann',
            maxNativeZoom: 15,
            minZoom: 6,
            vectorTileLayerStyles: vectorTileStyling,
        };

        // create VectorGrid layer
        var mapPbfLayer = new L.VectorGrid.Protobuf(
            mapUrl, mapVectorTileOptions
        )

        // add VectorGrid layer to map
        mapPbfLayer.addTo(map);
    </script>
</body>

</html>

Github Pages

You can use GitHub Pages to host the whole map online. GitHub has a size limit for individual files of 100MB. Considering that the original GeoJSON has 2.2GB, that may sound like a problem. But remember that the data are now cut to small tiles, where the size of each is a few KB. You want to create a repository on GitHub, upload your tiles/ folder and the HTML file you have just created, called index.html and set up GitHub Pages hosting.

Once your data are uploaded, create an empty file called .nojekyll that tells GitHub Pages to use the repository as an HTML code as it is (what it really says is “do not use Jekyll to compile the page”) and go to Settings > Pages and enable GitHub Pages.

🎉 Your web map is live and accessible by anyone online. GitHub Pages settings show you the exact link. It is as easy as that.

Tip: Once you have everything set on GitHub, use github.dev to make the further changes to the HTML file avoid download and indexing of the whole repository (just hit . when you are on GitHub in your repository).

Our map lives in the urbangrammarai/great-britain repository if you want to check the source code and is available on https://urbangrammarai.xyz/great-britain/ (we use the custom domain on GitHub Pages).

Further styling and interactivity

You may have noticed that our map has a couple of features on top of this. Let’s quickly go through those.

Labels on top of the map

You can use another XYZ base map that has just labels and load it on top of your tiles to see the names of places (or roads or something else). The order of layers matches the order in which you use .addTo(map).

Adding CartoDB Positron labels (this need to be after mapPbfLayer.addTo(map):

// add labels basemap
var mapBaseLayer = L.tileLayer(
    'https://{s}.basemaps.cartocdn.com/rastertiles/light_only_labels/{z}/{x}/{y}.png', {
        attribution: '(C) OpenStreetMap contributors (C) CARTO'
    }
).addTo(map);

You may want to show some additional information when you click on the geometry. Let’s link show the name of the signature type based on its code. First, create a dictionary with the info:

// mapping of names to signature types to be shown in popup on click
const popup_info = {
    "0_0": "Countryside agriculture",
    "1_0": "Accessible suburbia",
    "3_0": "Open sprawl",
    "4_0": "Wild countryside",
    "5_0": "Warehouse/Park land",
    "6_0": "Gridded residential quarters",
    "7_0": "Urban buffer",
    "8_0": "Disconnected suburbia",
    "2_0": "Dense residential neighbourhoods",
    "2_1": "Connected residential neighbourhoods",
    "2_2": "Dense urban neighbourhoods",
    "9_0": "Local urbanity",
    "9_1": "Concentrated urbanity",
    "9_2": "Regional urbanity",
    "9_4": "Metropolitan urbanity",
    "9_5": "Hyper concentrated urbanity",
};

Then you need to link the on click popup element to the VectorGrid object.

// create VectorGrid layer
var mapPbfLayer = new L.VectorGrid.Protobuf(
    mapUrl, mapVectorTileOptions
).on('click', function (e) { // this line and below
    L.popup()
        .setContent(popup_info[e.layer.properties.signature_type])
        .setLatLng(e.latlng)
        .openOn(map);
});

Conditional formatting based on the zoom level

The white outline can be a bit think when you zoom out, but there is a way how you can control that. You can just specify the condition to a variable and pass the variable to the tile styling dictionary.

// define styling of vector tiles
var vectorTileStyling = {
    signatures_4326: function(properties, zoom) {
                // set line weight conditionally based on zoom
        var weight = 0;
        if (zoom > 12) {
            weight = 1.0;
        }
        return ({
            fill: true,
            fillColor: cmap[properties.signature_type],
            fillOpacity: 0.9,
            weight: weight, // pass the weight variable instead of a value
            color: "#ffffff",
            opacity: 1.0,
        });
    }
}

The function signature has a new zoom attribute you can use in the condition.

Legend and details

A legend and details are just custom HTML objects detached from the leaflet map. You can check the source code to see how they’re done.

Understanding the structure of cities through the lens of data

The workshop organised together with James D. Gaboardi during the Spatial Data Science Symposium 2022 is now available online. See the recording below and access the workshop material on Github from which you can even run the code online, in your browser.

Annotation

Martin & James will walk you through the fundamentals of analysis of the structure of cities. You will learn what can be measured, how to do that in Python, and how to use those numbers to capture the patterns that make up cities. Using a real-life example, you will learn every step of the method, from data retrieval to detection of individual patterns.

Introducing Dask-GeoPandas for scalable spatial analysis in Python

Using Python for data science is usually a great experience, but if you’ve ever worked with pandas or GeoPandas, you may have noticed that they use only a single core of your processor. Especially on larger machines, that is a bit of a sad situation.

Developers came up with many solutions to scale pandas, but the one that seems to take the lead is Dask. Dask (specifically dask.dataframe as Dask can do much more) creates a partitioned data frame, where each partition is a single pandas.DataFrame. Each of them can be then processed in parallel and combined when necessary. On top of that, the whole pipeline can be scaled to a cluster of machines and can deal with out-of-core computation, i.e. with datasets that do not fit the memory.

Today, we announce the release of Dask-GeoPandas 0.1.0, a new Python package that extends dask.dataframe in the same way GeoPandas extends pandas, bringing the support for geospatial data to Dask. That means geometry columns and spatial operations but also spatial partitioning, ensuring that geometries that are close in space are within the same partition, necessary for efficient spatial indexing.

The project has been in development for quite some time. The original exploration of bridging Dask and GeoPandas started almost 5 years ago by Matt Rocklin, the author of Dask. Later, in 2020, Julia Signell revised the idea and created the foundations of the current project. Since then, GeoPandas maintainers have taken over and led the recent development.

What is awesome about Dask-GeoPandas? First, you can do your spatial analysis in parallel, making sure all available resources are used (no more sad idle cores!), turning your workflow into faster and more efficient ones. You can also use Dask-GeoPandas to process data that do not fit your machine’s memory as Dask comes with a support of out-of-core computation. Finally, you can distribute the work across many machines in a cluster. And all that with almost the same familiar GeoPandas APIs.

The latest evolution of underlying libraries powering GeoPandas ensures that it is efficient in terms of utilisation of resources but also performant within each partition. For example, unlike GeoPandas, where the use of pygeos, a new vectorised interface to GEOS is optional, Dask-GeoPandas requires it. Similarly, it depends on pyogrio, a vectorised interface to GDAL, to read geospatial file formats.

At this moment, Dask-GeoPandas can do a lot of what GeoPandas can, with some limitations. When your code involves individual geometries, without assessing a relationship between them (like computing a centroid or area), you should be able to use it directly. When you need to work out some relationships, you can try (still a bit limited) sjoin or make use of spatial partitions and spatial indexing.

But not everything is ready. For example, overlapping computation needed for use cases like accessibility or K-nearest neighbour analyses is not yet implemented, PostGIS IO is not done, and some overlay operations are implemented only partially (sjoin) or not at all (overlay, sjoin_nearest). But the 0.1.0 release is just a start.

You can try it yourself, installing via conda (or mamba) or from PyPI (but see the instructions, GeoPandas can be tricky to install using pip).

mamba install dask-geopandas
pip install dask-geopandas 

The best starting point to learn how Dask-GeoPandas works is the documentation, but this is the gist:

import geopandas
import dask_geopandas

df = geopandas.read_file(
    geopandas.datasets.get_path("naturalearth_lowres")
)
dask_df = dask_geopandas.from_geopandas(df, npartitions=8)

dask_df.geometry.area.compute()

The code creates a dask_geopandas.GeoDataFrame with 8 partitions because I have 8 cores and compute each polygon’s area in parallel, giving almost 8x speedup compared to the vanilla GeoPandas.

You can also check my latest post comparing Dask-GeoPandas performance on a large spatial join with PostGIS and cuSpatial (GPU) implementations.

If you want to help, have questions or ideas, you are always welcome. Just head over to Github or Gitter and say hi!

Dask-GeoPandas vs PostGIS vs GPU: Performance and Spatial Joins

Paul Ramsey saw a spatial join done using a GPU and tried to do the same with PostGIS, checking how fast that is compared to the GPU-based RAPIDS.AI solution. Since Paul used parallelisation in PostGIS, I got curious how fast Dask-GeoPandas is on the same task.

So, I gave it a go.

import download
import geopandas
import dask_geopandas
import dask.dataframe
from dask.distributed import Client, LocalCluster

Let’s download the data using Paul’s query, to ensure we work with the same CSV.

curl "https://phl.carto.com/api/v2/sql?filename=parking_violations&format=csv&skipfields=cartodb_id,the_geom,the_geom_webmercator&q=SELECT%20*%20FROM%20parking_violations%20WHERE%20issue_datetime%20%3E=%20%272012-01-01%27%20AND%20issue_datetime%20%3C%20%272017-12-31%27" > phl_parking.csv

And then download and unzip the neighbourhoods shapefile.

download.download(
    "https://github.com/azavea/geo-data/raw/master/Neighborhoods_Philadelphia/Neighborhoods_Philadelphia.zip",
    "Neighborhoods_Philadelphia", 
    kind="zip"
)

Paul used a machine with 8 cores. Since I use a machine with 16 cores, I’ll create a local cluster limited to 8 workers. That should be as close to Paul’s machine as I can get without using some virtual one. Keep in mind that this distorts the benchmark as we use different processors with different performances. But the point here is to get a sense of how fast can Dask-based solution be compared to PostGIS and the original GPU code.

client = Client(
    LocalCluster(
        n_workers=8, 
        threads_per_worker=1
    )
)

With Dask, we create the whole pipeline to create a task graph and then run it all, so we won’t have the timings for individual steps, just the total one.

Read parking data CSV into a partitioned data frame (25MB per partition).

ddf = dask.dataframe.read_csv(
    "phl_parking.csv", 
    blocksize=25e6, 
    assume_missing=True
)

Create point geometry and assign it to the data frame, creating dask_geopandas.GeoDataFrame.

ddf = ddf.set_geometry(
    dask_geopandas.points_from_xy(
        ddf, 
        x="lon", 
        y="lat", 
        crs=4326
    )
)

Read neighbourhood polygons and reproject to EPSG:4326 (same as parking data).

neigh = geopandas.read_file(
    "Neighborhoods_Philadelphia"
).to_crs(4326)

Create the spatial join.

joined = dask_geopandas.sjoin(ddf, neigh, predicate="within")

Finally, let’s compute the result.

%%time
r = joined.compute()

Time on a local cluster with 8 workers and 1 thread per worker to pretend it is an 8-core CPU:

CPU times: user 9.34 s, sys: 2.09 s, total: 11.4 s
Wall time: 21.3 s

The complete pipeline took 21.3 seconds, including sending all data to a single process, in the end, to create a single partition joined GeoDataFrame. Usually, that is unnecessary as you work with the data directly in Dask. It does take a few seconds guessing from the Dask Dashboard.

Let’s compare it to the PostGIS solution:

  • Reading in the 9M records from CSV takes about 29 seconds
  • Making a second copy while creating a geometry column takes about 24 seconds
  • The final query running with 4 workers takes 24 seconds

That gives us a total of 77 seconds compared to 21 seconds using Dask-GeoPandas. It’s still slower than 13 seconds using RAPIDS.AI although that covers only the join itself, not reading and creating geometry, so my sense is that it will be almost equal. One aspect that makes the difference between Dask and PostGIS is that our pipeline is parallelised at every step – reading the CSV, creating points, generating spatial index (that is done under the hood in sjoin), the actual join.

While Paul was using the 8-core machine, PostGIS actually utilised only 4 cores (I am not sure why). Let’s try to run our code limited to 4 workers as well.

CPU times: user 9.53 s, sys: 2 s, total: 11.5 s
Wall time: 28.4 s

28 seconds is a bit slower than before, but still quite fast!

When comparing PostGIS and GPU solutions, Paul says

Basically, it is very hard to beat a bespoke performance solution with a general-purpose tool. Yet, PostgreSQL/PostGIS comes within “good enough” range of a high end GPU solution, so that counts as a “win” to me.

At the moment, Dask-GeoPandas is somewhere between PostGIS and bespoke solutions. It does not offer as many functions as PostGIS, but it is designed as a general-purpose tool. So I would say that we are all winners here.

The notebook is available here.

EDIT (Mar 24, 2022): See also Dewey Dunnington’s follow-up expanding the comparison to R.

Methodological Foundation of a Numerical Taxonomy of Urban Form

The final paper based on my PhD thesis is (finally!) out in the Environment and Planning B: Urban Analytics and City Science. We looked into ways of identifying patterns of urban form and came up with the Methodological foundation of a numerical taxonomy of urban form. You can read it on the journal website (open access).

We use urban morphometrics (i.e. data-driven methods) to derive a classification of urban form in Prague and Amsterdam, and you can check the results in online interactive maps – http://martinfleis.github.io/numerical-taxonomy-maps/ or below. (Check the layers toggle!)

The paper explores the method that can eventually support the creation of a taxonomy of urban types in a similar way you know from biology. We even borrowed the foundations of the method from biology.

We measure many variables based on building footprints and street networks (using the momepy Python package) and use Gaussian Mixture Model clustering to get urban tissue types independently in both cities. Then we apply Ward’s hierarchical clustering to build a taxonomy of these types.

The code is available, and the repo even includes a clean Jupyter notebook with the complete method, so you can apply it to your data if you wish. https://github.com/martinfleis/numerical-taxonomy-paper. If you instead want to play with our data, it is available as well https://doi.org/10.6084/m9.figshare.16897102.

Abstract

Cities are complex products of human culture, characterised by a startling diversity of visible traits. Their form is constantly evolving, reflecting changing human needs and local contingencies, manifested in space by many urban patterns. Urban morphology laid the foundation for understanding many such patterns, largely relying on qualitative research methods to extract distinct spatial identities of urban areas. However, the manual, labour-intensive and subjective nature of such approaches represents an impediment to the development of a scalable, replicable and data-driven urban form characterisation. Recently, advances in geographic data science and the availability of digital mapping products open the opportunity to overcome such limitations. And yet, our current capacity to systematically capture the heterogeneity of spatial patterns remains limited in terms of spatial parameters included in the analysis and hardly scalable due to the highly labour-intensive nature of the task. In this paper, we present a method for numerical taxonomy of urban form derived from biological systematics, which allows the rigorous detection and classification of urban types. Initially, we produce a rich numerical characterisation of urban space from minimal data input, minimising limitations due to inconsistent data quality and availability. These are street network, building footprint and morphological tessellation, a spatial unit derivative of Voronoi tessellation, obtained from building footprints. Hence, we derive homogeneous urban tissue types and, by determining overall morphological similarity between them, generate a hierarchical classification of urban form. After framing and presenting the method, we test it on two cities – Prague and Amsterdam – and discuss potential applications and further developments. The proposed classification method represents a step towards the development of an extensive, scalable numerical taxonomy of urban form and opens the way to more rigorous comparative morphological studies and explorations into the relationship between urban space and phenomena as diverse as environmental performance, health and place attractiveness.

Capturing the Structure of Cities with Data Science

During the Spatial Data Science Conference 2021, I had a chance to deliver a workshop illustrating the application of PySAL and momepy in understanding the structure of cities. The recording is now available for everyone. The materials are available on my GitHub and you can even run the whole notebook in your browser using the MyBinder service.

xyzservices: a unified source of XYZ tile providers in Python

A Python ecosystem offers numerous tools for the visualisation of data on a map. A lot of them depend on XYZ tiles, providing a base map layer, either from OpenStreetMap, satellite or other sources. The issue is that each package that offers XYZ support manages its own list of supported providers.

We have built xyzservices package to support any Python library making use of XYZ tiles. I’ll try to explain the rationale why we did that, without going into the details of the package. If you want those details, check its documentation.

The situation

Let me quickly look at a few popular packages and their approach to tile management – contextily, folium, ipyleaflet and holoviews.

contextily

contextily brings contextual base maps to static geopandas plots. It comes with a dedicated contextily.providers module, which contains a hard-coded list of providers scraped from the list used by leaflet (as of version 1.1.0).

folium

folium is a lightweight interface to a JavaScript leaflet library. It providers built-in support for 6 types of tiles and allows passing any XYZ URL and its attribution to a map. It means that it mostly relies on external sources of tile providers.

ipyleaflet

ipyleaflet brings leaflet support to Jupyter notebooks and comes with a bit more options than folium. It has a very similar approach as contextily does – it has a hard-coded list of about 37 providers in its basemaps module.

holoviews

holoviews provides a Python interface to the Bokeh library and its list of supported base maps is also hard-coded.

A similar situation is in other packages like geemap or leafmap.

Each package has to maintain the list of base maps, ensure that they all work, respond to users requiring more, update links… That is a lot of duplicated maintenance burden. We think it is avoidable.

The vision

All XYZ tile providers have a single lightweight home and a clean API supporting the rest of the ecosystem. All the other packages use the same resource, one which is tested and expanded by a single group of maintainers.

We have designed xyzservices to be exactly that. It is a Python package that has no dependencies and only a single purpose – to collect and process metadata of tile providers.

We envisage a few potential use cases.

The first – packages like contextily and geopandas will directly support xyzservices.TileProvider object when specifying tiles. Nothing else is needed, contextily will fetch the data it needs (final tile URL, an attribution, zoom and extent limits) from the object. In the code form:

import xyzservices.providers as xyz
from contextily import add_basemap

add_basemap(ax, source=xyz.CartoDB.Positron)

The second option is wrapping xyzservices.providers into a custom API providing, for example, an interactive selection of tiles.

The third one is using different parts of a TileProvider individually when passing the information. This option can be currently used, for example, with folium:

import folium
import xyzservices.providers as xyz

tiles = xyz.CartoDB.Positron

folium.Map(
    location=[53.4108, -2.9358],
    tiles=tiles.build_url(),
    attr=tiles.attribution,
)

The last one is the most versatile. The xyzservices comes with a JSON file used as a storage of all the metadata. The JSON is automatically installed to share/xyzservices/providers.json where it is available for any other package without depending on xyzservices directly.

We hope to cooperate with maintainers of other existing packages and move most of the functionality around XYZ tiles that can be reused to xyzservices. We think that it will:

  1. Remove the burden from individual developers. Any package will just implement an interface to Python or JSON API of xyzservices.
  2. Expand the list easy-to-use tiles for users. xyzservices currently has over 200 providers, all of which should be available for users across the ecosystem, without the need to individually hard-code them in every package.

While this discussion started in May 2020 (thanks @darribas!), the initial version of the package is out now and installable from PyPI and conda-forge. We hope to have as many developers as possible on board to allow for the consolidation of the ecosystem in the future.

Evolution of Urban Patterns: Urban Morphology as an Open Reproducible Data Science

We have a new paper published in the Geographical Analysis on the opportunities current developments in geographic data science within the Python ecosystem offer to urban morphology. To sum up – there’s a lot to play with and if you’re interested in the quantification of urban form, there’s no better choice for you at the moment.

Urban morphology (study of urban form) is historically a qualitative discipline that only recently expands into more data science-ish waters. We believe that there’s a lot of potential in this direction and illustrate it on the case study looking into the evolution of urban patterns, i.e. how different aspects of urban form has changed over time.

The paper is open access (yay!) (links to the online version and a PDF), the research is fully reproducible (even in your browser thanks to amazing MyBinder!) with all code on GitHub.

Short summary

We have tried to map all the specialised open-source tools urban morphologists can use these days, which resulted in this nice table. The main conclusion? Most of them are plug-ins for QGIS or ArcGIS, hence depend on pointing and clicking. A tricky thing to reproduce properly.

Table 1 from the paper

We prefer code-based science, so we took the Python ecosystem and put it in the test. We have designed a fully reproducible workflow based on GeoPandas, OSMnx, PySAL and momepy to sample and study 42 places around the world, developed at very different times.

We measured a bunch of morphometric characters (indicators for individual aspects of form) and looked at their change in time. And there is a lot to look at. There are significant differences not only in scale (on the figure below) but in other aspects as well. One interesting observation – it seems that we have forgotten how to make a properly connected, dense grid.

As we said in the paper, “Switching to a code-based analysis may be associated with a steep learning curve. However, not everyone needs to reach the developer level as the data science ecosystem aims to provide a middle ground user level. That is a bit like Lego—the researcher learns how to put pieces together and then find pieces they need to build a house.

We think that moving from QGIS to Python (or R), as daunting as it may seem to some, is worth it. It helps us overcome the reproducibility crisis science is going through, the crisis caused, among other things, by relying on pointing and clicking too much.

The open research paradigm, based on open platforms and transparent community-led governance, has the potential to democratize science and remove unnecessary friction caused by the lack of cooperation between research groups while bringing additional transparency to research methods and outputs.

page 18

Abstract

The recent growth of geographic data science (GDS) fuelled by increasingly available open data and open source tools has influenced urban sciences across a multitude of fields. Yet there is limited application in urban morphology—a science of urban form. Although quantitative approaches to morphological research are finding momentum, existing tools for such analyses have limited scope and are predominantly implemented as plug-ins for standalone geographic information system software. This inherently restricts transparency and reproducibility of research. Simultaneously, the Python ecosystem for GDS is maturing to the point of fully supporting highly specialized morphological analysis. In this paper, we use the open source Python ecosystem in a workflow to illustrate its capabilities in a case study assessing the evolution of urban patterns over six historical periods on a sample of 42 locations. Results show a trajectory of change in the scale and structure of urban form from pre-industrial development to contemporary neighborhoods, with a peak of highest deviation during the post-World War II era of modernism, confirming previous findings. The wholly reproducible method is encapsulated in computational notebooks, illustrating how modern GDS can be applied to urban morphology research to promote open, collaborative, and transparent science, independent of proprietary or otherwise limited software.

Fleischmann, M., Feliciotti, A. and Kerr, W. (2021), Evolution of Urban Patterns: Urban Morphology as an Open Reproducible Data Science. Geogr Anal. https://doi.org/10.1111/gean.12302

Talk at ISUF 2021: Classifying urban form at a national scale

I had a chance to present our ongoing work on the classification of the (built) environment in Great Britain during the International Seminar on Urban Form 2021, which was held virtually in Glasgow. I was presenting the classification of urban form, one component of Spatial Signatures we’re developing as part of the Urban Grammar AI project together with Dani Arribas-Bel. The video of the presentation is attached below, as well as the abstract.

Classifying urban form at a national scale: The case of Great Britain

There is a pressing need to monitor urban form and function in ways that can feed into better planning and management of cities. Both academic and policymaking communities have identified the need for more spatially and temporally detailed, consistent, and scalable evidence on the nature and evolution of urban form. Despite impressive progress, the literature can achieve only two of those characteristics simultaneously. Detailed and consistent studies do not scale well because they tend to rely on small-scale, ad-hoc datasets that offer limited coverage. Until recently, consistent and scalable research has only been possible by using simplified measures that inevitably miss much of the nuance and richness behind the concept of urban form.

This paper outlines the notion of “spatial signatures”, a characterisation of space based on form and function, and will specifically focus on its form component. Whilst spatial signature sits between the purely morphological and purely functional description of the built environment, its form-based component reflects the morphometric definition of urban tissue, the distinct structurally homogenous area of a settlement. The proposed method employs concepts of “enclosures” and “enclosed tessellation” to derive indivisible hierarchical geographies based on physical boundaries (streets, railway, rivers, coastline) and building footprints to delineate such tissues in the built fabric. Each unit is then characterised by a comprehensive set of data-driven morphometric characters feeding into an explicitly spatial contextual layer, which is used as an input of cluster analysis.

The classification based on spatial signatures is applied to the entirety of Great Britain on a fine grain scale of individual tessellation cells and released as a fully reproducible open data product. The results provide a unique input for local authorities to drive planning and decision-making and for the wider research community as data input.

Video