Chapter 2: Visualizing and modelling spatial data

Tomas Beuzen, 2021

Chapter Learning Objectives

  • Make informed choices about how to plot your spatial data, e.g., scattered, polygons, 3D, etc..

  • Plot spatial data using libraries such as geopandas, plotly, and keplergl.

  • Interpolate unobserved spatial data using deterministic methods such as nearest-neighbour interpolation.

  • Interpolate data from one set of polygons to a different set of polygons using areal interpolation.


import warnings
import keplergl
import numpy as np
import osmnx as ox
import pandas as pd
import geopandas as gpd
import as px
from skgstat import Variogram
import matplotlib.pyplot as plt
from shapely.geometry import Point
from pykrige.ok import OrdinaryKriging
from scipy.interpolate import NearestNDInterpolator
from tobler.area_weighted import area_interpolate
# Custom functions
from scripts.utils import pixel2poly
# Plotting defaults'ggplot')
px.defaults.height = 400; px.defaults.width = 620
plt.rcParams.update({'font.size': 16, 'axes.labelweight': 'bold', 'figure.figsize': (6, 6), 'axes.grid': False})

1. Spatial visualization

1.1. Geopandas

We saw last chapter how to easily plot geospatial data using the geopandas method .plot(). This workflow is useful for making quick plots, exploring your data, and easily layering geometries. Let’s import some data of UBC buildings using osmnx (our Python API for accessing OpenStreetMap data) and make a quick plot:

ubc = (ox.geometries_from_place("University of British Columbia, Canada", tags={'building':True})
         .loc[:, ["geometry"]]                 # just keep the geometry column for now
         .query("geometry.type == 'Polygon'")  # only what polygons (buidling footprints)
         .assign(Label="Building Footprints")  # assign a label for later use
         .reset_index(drop=True)               # reset to 0 integer indexing
geometry Label
0 POLYGON ((-123.25526 49.26695, -123.25506 49.2... Building Footprints
1 POLYGON ((-123.25328 49.26803, -123.25335 49.2... Building Footprints
2 POLYGON ((-123.25531 49.26859, -123.25493 49.2... Building Footprints
3 POLYGON ((-123.25403 49.26846, -123.25408 49.2... Building Footprints
4 POLYGON ((-123.25455 49.26906, -123.25398 49.2... Building Footprints

Recall that we can make a plot using the .plot() method on a GeoDataFrame:

ax = ubc.plot(figsize=(8, 8), column="Label", legend=True,
              edgecolor="0.2", markersize=200, cmap="rainbow")

Say I know the “point” location of my office but I want to locate the building footprint (a “polygon”). That’s easily done with geopandas!

First, I’ll use shapely (the Python geometry library geopandas is built on) to make my office point, but you could also use the geopandas function gpd.points_from_xy() like we did last chapter:

point_office = Point(-123.2522145, 49.2629555)

Now, I can use the .contains() method to find out which building footprint my office resides in:

geometry Label
48 POLYGON ((-123.25217 49.26345, -123.25196 49.2... Building Footprints

Looks like it’s index 48! I’m going to change the label of that one to “Tom’s Office”:

ubc.loc[48, "Label"] = "Tom's Office"

Now let’s make a plot!

ax = ubc.plot(figsize=(8, 8), column="Label", legend=True,
              edgecolor="0.2", markersize=200, cmap="rainbow")

We can add more detail to this map by including a background map. For this, we need to install the contextily package. Note that most web providers use the Web Mercator projection, “EPSG:3857” (interesting article on that here) so I’ll convert to that before plotting:

import contextily as ctx

ax = (ubc.to_crs("EPSG:3857")
         .plot(figsize=(10, 8), column="Label", legend=True,
               edgecolor="0.2", markersize=200, cmap="rainbow")
ctx.add_basemap(ax, source=ctx.providers.OpenStreetMap.Mapnik)  # I'm using OSM as the source. See all provides with ctx.providers

1.2. Plotly

The above map is nice, but it sure would be helpful in some cases to have interactive functionality for our map don’t you think? Like the ability to zoom and pan? Well there are several packages out there that can help us with that, but plotly is one of the best for this kind of mapping. plotly supports plotting of maps backed by MapBox (a mapping and location data cloud platform).

Here’s an example using the function px.choropleth_mapbox():

fig = px.choropleth_mapbox(ubc, geojson=ubc.geometry, locations=ubc.index, color="Label",
                           center={"lat": 49.261, "lon": -123.246}, zoom=12.5,
fig.update_layout(margin=dict(l=0, r=0, t=30, b=10))