Spatial aggregation: from HDB transactions to planning areas
Question: How does summarising individual sales by geographical areas change the patterns we describe?
Objectives: check the coordinate system, spatially join HDB transactions to planning subzones, aggregate transaction records, and interpret differences between transaction and area-level statistics.
Related Part 1: Spatial representation and statistical comparisons. Source: HDB resale teaching extract and URA planning subzones in the data catalogue.
This notebook is a guided investigation. Its companion solution is openly accessible. You need pandas, geopandas, shapely and matplotlib. For remote execution, datasets are fetched from a pinned Git revision. No course Google Drive is necessary.
Companion notebook: View the worked solution.
from pathlib import Path
import pandas as pd
import geopandas as gpd
from shapely import wkt
import matplotlib.pyplot as plt
folder = Path("../datasets/singapore_hdb_ura")
if not folder.exists():
folder = Path("src/content/datasets/singapore_hdb_ura")
root = "https://raw.githubusercontent.com/wcchin/geo-stat-vis/8cf3ddb0a7eff7d94018940e91f32d7482f32aaa/src/content/datasets/singapore_hdb_ura/"
def source(filename):
path = folder / filename
return path if path.exists() else root + filename
hdb = pd.read_csv(source("hdb_data_2023.csv.xz"), index_col=0)
subzones = gpd.read_file(source("master-plan-2019-subzone_landuse_area.zip")).to_crs("EPSG:3414")
assert {"geometry","resale_price"}.issubset(hdb.columns)
assert {"PLN_AREA_N"}.issubset(subzones.columns)
print("Transactions", len(hdb), "Subzones", len(subzones))1. Inspect geography and locate transactions¶
The HDB extract stores geometries as WKT. Before joining, ensure that points and polygons have the same projected CRS. The dataset’s point geometry is supplied in Singapore’s projected coordinate system; confirm this from the source notes before comparing distances or areas.
points = gpd.GeoDataFrame(hdb.copy(), geometry=hdb["geometry"].map(wkt.loads), crs="EPSG:3414")
print(points.crs, subzones.crs)
print("Missing geometry:", points.geometry.isna().sum())
print("Valid prices:", points["resale_price"].notna().sum())2. Spatial join¶
A transaction may fall outside a mapped subzone due to coverage or coordinate issues. Count unmatched observations before computing area summaries.
matched = gpd.sjoin(points, subzones[["PLN_AREA_N","geometry"]], how="left", predicate="within")
print("Transactions without planning-area match:", matched["PLN_AREA_N"].isna().sum())
print("Matches:", matched["PLN_AREA_N"].value_counts().head())3. TODO — transaction-level versus area-level values¶
Calculate the median resale price and transaction count for each planning area. Explain why the average of area medians does not equal the median of all individual transactions.
# TODO: group matched transactions by PLN_AREA_N; compute median resale_price and count.
4. TODO — representation and interpretation¶
Draw a bar chart comparing the areas with the most transactions. Should areas with small numbers of records carry the same weight when interpreting average prices?
# TODO: plot transaction counts by planning area for the top 10 areas.
Reflection¶
What information is discarded when transactions are aggregated? Is this an areal rate, a transaction count or an area-level summary of individual sale prices? How might differences in flat type and age change the comparison?
Check your understanding¶
Try the exercises and write down your interpretation before comparing your work with the complete answers.