Chapter 2: Visualizing and modelling spatial data¶

Tomas Beuzen, 2021

Chapter Outline¶

Chapter Learning Objectives
Imports
1. Spatial visualization
2. Spatial modelling

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.

Imports¶

import warnings
import keplergl
import numpy as np
import osmnx as ox
import pandas as pd
import geopandas as gpd
import plotly.express 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
plt.style.use('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
      )
ubc.head()

	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")
plt.title("UBC");

../_images/chapter2_spatial-viz-and-modelling_14_0.png

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)
point_office

../_images/chapter2_spatial-viz-and-modelling_16_0.svg

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

ubc[ubc.contains(point_office)]

	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")
plt.title("UBC");

../_images/chapter2_spatial-viz-and-modelling_22_0.png

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
plt.axis("off")
plt.title("UBC");

../_images/chapter2_spatial-viz-and-modelling_24_0.png

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 plotly.express 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,
                           mapbox_style="open-street-map")
fig.update_layout(margin=dict(l=0, r=0, t=30, b=10))

Python for Geospatial Analysis