Do-It-Yourself

In this session, we will practice your skills in mapping with R and Python. Create a new quarto or jupyter notebook document you can edit interactively, and let’s do this!

library(dplyr)
library(ggplot2)
library(sf)
library(viridis)

import geopandas as gpd
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib.colors as mcolors
from matplotlib import cm
from matplotlib.ticker import FuncFormatter

Data preparation

Polygons

For this section, you will have to push yourself out of the comfort zone when it comes to sourcing the data. As nice as it is to be able to pull a dataset directly from the web at the stroke of a url address, most real-world cases are not that straight forward. Instead, you usually have to download a dataset manually and store it locally on your computer before you can get to work.

We are going to use data from the Geographic Data Service (GeoDS) about Liverpool, in particular an extract from the Census. GeoDS is part of Smart Data Research UK, based at UCL, the University of Liverpool, the University of Oxford and the University of Edinburgh — it took over hosting this data after the Consumer Data Research Centre (CDRC) portal it used to live on was retired. You can download a copy of the data at:

Important

Go to the Geographic Data Service’s teaching materials page:

https://data.geods.ac.uk/dataset/teaching-materials

You’ll need to log in (registration is free) to download. Find the “Census Residential Data Pack 2011 for Liverpool (E08000012)” resource on that page and download it.

Once you have the .zip file on your computer, right-click and “Extract all”. The resulting folder will contain all you need. Create a folder called Liverpool in data folder you created in the first Lab.

library(sf)
lsoas <- read_sf("data/Liverpool/Census_Residential_Data_Pack_2011/Local_Authority_Districts/E08000012/shapefiles/E08000012.shp")

lsoas = gpd.read_file("data/Liverpool/Census_Residential_Data_Pack_2011/Local_Authority_Districts/E08000012/shapefiles/E08000012.shp")

Lines

For a line layer, we’ll use roads_osm.gpkg — an OSM-derived roads extract for Liverpool.

roads <- read_sf("data/Liverpool/lines/roads_osm.gpkg")

ggplot() + 
  geom_sf(data = roads, color = 'darkorange') + theme_minimal()

roads = gpd.read_file("data/Liverpool/lines/roads_osm.gpkg")

fig, ax = plt.subplots(1, 1, figsize=(15, 10))
ax.set_title("Roads in Liverpool")
ax.set_axis_off() # we don't need the ticks function
roads.plot(ax=ax, color = 'darkorange', lw = 0.5) #linewidth can be also passed as lw 
# Display the plot
plt.show()

Points

For points, we’ll use points_merseyside.gpkg — points of interest (POIs) in Merseyside. Since this file covers the wider Merseyside region rather than just Liverpool, we clip it down to the lsoas boundary first.

pois_merseyside <- read_sf("data/Liverpool/points/points_merseyside.gpkg")
pois_merseyside <- st_transform(pois_merseyside, st_crs(lsoas)) # match lsoas' CRS before filtering
pois <- st_filter(pois_merseyside, st_as_sfc(st_bbox(lsoas)))

ggplot() + 
  geom_sf(data = pois) + theme_minimal()

pois_merseyside = gpd.read_file("data/Liverpool/points/points_merseyside.gpkg")
pois_merseyside = pois_merseyside.to_crs(lsoas.crs)  # match lsoas' CRS before clipping
pois = pois_merseyside.clip(lsoas.total_bounds)

pois.plot()
plt.show()

Tasks

Task I: Tweak your map

With those three layers, try to complete the following tasks:

Make a map of the Liverpool neighborhoods that includes the following characteristics:
1. Features a title
2. Does not include axes frame
3. Polygons are all in color #525252 and 50% transparent
4. Boundary lines (“edges”) have a width of 0.3 and are of color #B9EBE3
5. Includes a basemap different from the one used in class

Note

Not all of the requirements above are not equally hard to achieve. If you can get some but not all of them, that’s also great! The point is you learn something every time you try.

Task II: Non-spatial manipulations

For this one we will combine some of the ideas we learnt in the previous block with this one.

Focus on the LSOA liverpool layer and use it to do the following:

Calculate the area of each neighbourhood
Find the five smallest areas in the table. Create a new object (e.g. smallest with them only)
Create a multi-layer map of Liverpool where the five smallest areas are coloured in red, and the rest appear in grey.

Task III: Average price per district

districts <- read_sf("data/London/Polygons/districts.shp")
housesales <- read.csv("data/London/Tables/housesales.csv") # import housesales data from csv

housesales_clean = housesales %>%
  filter(price < 500000) %>% 
  st_as_sf(coords = c("greastings", "grnorthing")) %>% 
  st_set_crs(27700)

districts = gpd.read_file("data/London/Polygons/districts.shp")
housesales = pd.read_csv("data/London/Tables/housesales.csv")

housesales_f = housesales[housesales['price'] < 500000]

housesales_gdf = gpd.GeoDataFrame(
    housesales_f, geometry=gpd.points_from_xy(housesales_f.greastings, housesales_f.grnorthing), crs="EPSG:27700"
)

This one is a bit more advanced, so don’t despair if you can’t get it on your first try. It relies on the London data you used in the Lab. Here is the questions for you to answer:

What is the district with the highest housing prices in London?

Answering this questions involve 3 steps:

1. Performing a spatial join (st_join) between the district layer (polygons) and the households (points).

2. Aggregating the data at district level: group_by & summarise()

3. Figure out the district with the highest price

Really try not to open the answer below right away.

R Answer

Spatial overlay between points and polygons. We use left = FALSE so that only districts with at least one matching sale are kept — this mirrors the Python answer below, which uses an inner join, and avoids NA prices sneaking into the average for districts with no sales at all.

housesales_districts <- st_join(districts, housesales_clean, left = FALSE)

Aggregate at district level

housesales_districts_agg <- housesales_districts %>% 
  group_by(DIST_CODE, DIST_NAME) %>% # group at district level
  summarise(count_sales = n(),  # create count
            mean_price = mean(price)) # average price

`summarise()` has grouped output by 'DIST_CODE'. You can override using the
`.groups` argument.

head(housesales_districts_agg)

Simple feature collection with 6 features and 4 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: 515484.9 ymin: 156480.8 xmax: 554503.8 ymax: 198355.2
Projected CRS: OSGB36 / British National Grid
# A tibble: 6 × 5
# Groups:   DIST_CODE [6]
  DIST_CODE DIST_NAME           count_sales mean_price                  geometry
  <chr>     <chr>                     <int>      <dbl>             <POLYGON [m]>
1 00AB      Barking and Dagenh…          38     91802. ((550817 184196, 550814 …
2 00AC      Barnet                       83    169662. ((526830.3 187535.5, 526…
3 00AD      Bexley                       82    119276. ((552373.5 174606.9, 552…
4 00AE      Brent                        49    174498. ((524661.7 184631, 52466…
5 00AF      Bromley                     124    142468. ((533852.2 170129, 53385…
6 00AG      Camden                       37    224889. ((531410.7 181576.1, 531…

Python Answer

# Spatial join (inner: matches the R answer's left = FALSE above)
housesales_districts = housesales_gdf.sjoin(districts, how="inner", predicate='intersects')

# Aggregate at district level
housesales_districts_agg = housesales_districts.groupby(['DIST_CODE', 'DIST_NAME']).agg(
    count_sales=('price', 'size'),  # count number of sales
    mean_price=('price', 'mean')    # calculate average price
).reset_index()

# Merge the aggregated data back with the original GeoDataFrame to retain geometry
housesales_districts_agg = districts[['DIST_CODE', 'geometry']].drop_duplicates().merge(
    housesales_districts_agg, on='DIST_CODE'
)

Once that’s done, create a map of the data

R Answer

Use ggplot if you are working in R and if you’re feeling adventurous the function scale_fill_viridis() to make your map look especially good.

# map housesales by wards
map3 <- ggplot()+
  geom_sf(data = housesales_districts_agg, inherit.aes = FALSE, aes(fill = mean_price)) + # add the district level housing price 
  scale_fill_viridis("Price", direction = -1, labels = scales::dollar_format(prefix = "£"), option = "magma" )+ # change the legend scale to £ and the colour to magma
  xlab("") +
  ylab("") +
  theme_minimal() # choose a nicer theme https://ggplot2.tidyverse.org/reference/ggtheme.html
map3

Python Answer

# Function to format currency with £
def currency(x, pos):
    return f'£{int(x):,}'

# Create a colormap similar to Viridis magma
cmap = plt.colormaps.get_cmap('magma_r')

# Plotting the map
fig, ax = plt.subplots(figsize=(10, 10))

# Plot the district level housing price
housesales_districts_agg.plot(column='mean_price', ax=ax, legend=True, 
                              cmap=cmap, edgecolor='black')

# Modify the legend scale to £
formatter = FuncFormatter(currency)
cbar = ax.get_figure().get_axes()[1]  # Get the color bar
cbar.yaxis.set_major_formatter(formatter)

# Remove x and y axis labels
ax.set_xlabel('')
ax.set_ylabel('')

# Set theme to minimal
ax.set_axis_off()

# Show the plot
plt.show()