John Snow & Spatial Analytics

🎯 Introduction & Objectives

In his book The Ghost Map: The Story of London's Most Terrifying Epidemic - and How it Changed Science, Cities and the Modern World, author Steven Berlin Johnson describes how Dr. John Snow, with the assistance of Reverend Henry Whitehead utilized spatial analysis to support their theory that cholera is a waterborne disease. Their analysis resulted in a famous map of the London SOHO neighborhood identifying public water pump locations overlaid with cholera deaths from the August 1854 cholera outbreak.

The intuitive presentation of their spatial analysis convinced local officials and the larger medical community that the Broad Street Pump spread cholera and, more importantly, via contaminated water.

Objective

In this lab, you will recreate one of the first recorded uses of spatial analysis utilizing Geographic Information Systems (GIS). For this lab, you will do the following:

• Recreate the Ghost Map.
• Spatially analyze the 1854 London cholera epidemic.
• Create the maps and a summary of deliverables.

---

🗺️ Option 1: QGIS Analysis

In 1854, John Snow was tasked with performing and presenting his spatial analysis using a pen and paper. Today, we can use QGIS to recreate and analyze the original Ghost Map.

Step 1: Georeferencing the Original Map

The Ghost_Map.png must be georeferenced before overlaying supporting GIS data. Its native coordinates are in "model space".

1. Open QGIS and start a blank project.
   
2. Before dealing with the map, add your tabular datasets. Go to Layer > Add Layer > Add Vector Layer...
   
3. Browse for PUMPS.dbf and add it. Repeat for DEATHS.dbf. Close the dialog.
   

   

4. You'll notice nothing is on the map yet because the .dbf files only contain data tables, not geometry yet. However, PUMPS.dbf has coordinates!
   

   

5. Let's geocode the pumps. Open the Processing Toolbox (Gear icon or Processing > Toolbox). Search for and run Create points from table.
   

   

6. Input table: PUMPS. X field: X_BNG. Y field: Y_BNG. Target CRS: EPSG:27700 - OSGB36 / British National Grid. Click Run and close. Rename the output layer to PUMPS Events.
   

   

   

7. Now, select Layer > Georeferencer... from the top menu.
   

   

8. In the Georeferencer window, click File > Open Raster and select Ghost_Map.png.
   

   

9. Click on the Add Point tool. Find a recognizable pump on the PNG (like the Oxford Market Pump near Market Place/Castle Street East) and click it. Tip: You can turn on labels for your PUMPS Events layer in the main QGIS window to help identify points.
   

   

10. When prompted for coordinates, click the button "From Map Canvas". Your QGIS main map window will appear. Make sure PUMPS Events is visible and click the corresponding pump point.
    

    

11. Repeat this process for at least 5 different pumps spread across all four corners of the map. Remember to check the Residual (pixels) column to ensure your points are geographically accurate.
    

    

12. Once you have control points, click the yellow gear icon (Transformation Settings). Choose Polynomial 1 (or Linear) as transformation type, Nearest Neighbour for resampling, and set the Target SRS to EPSG:3857 or EPSG:27700. Check the option to Create world file only if you prefer not to create a new modified TIF, though producing a new, modified TIF is also an option. Save the output raster in your folder. Click OK.
    

    

13. Click the green Play button (Start Georeferencing). It's a good idea to save your GCPs (Ground Control Points) in case you need to troubleshoot later. Close the georeferencer and add the new georeferenced map to your main QGIS map.
    

    

14. Fixing the map projection: Notice how the map might look squished? This is because of the default projection (EPSG:4326). If you want the map to look better, you need to set the project coordinate system to the British National Grid (EPSG:27700) using the button in the bottom right corner.
    

    

    

15. Finally, you can rename layers and rearrange them (by dragging them in the Layers panel) so that the pumps and deaths show up correctly on top of the map.
    

    

Step 2: Tabular Data Analysis

1. We need to bring the deaths onto the map. Run Create points from table again.
2. Input: DEATHS. X field: X_BNG. Y field: Y_BNG. CRS: EPSG:27700. Run and rename to DEATHS Events.
3. Open the attribute table for DEATHS Events (Right-click > Open Attribute Table). Note the COUNT field (number of fatalities) and STREET field.
4. Open the Processing Toolbox and run Statistics by categories.
5. Input layer: DEATHS Events. Field to calculate: COUNT. Field with categories: STREET. Save the output table.
6. Open the resulting table and sort by Sum. Which streets have the most fatalities? Is this conclusive evidence on its own? Probably not, but it supports a broader spatial analysis.

Step 3: Proportional Symbols

1. Double-click the DEATHS Events layer to open Properties, go to Symbology.
2. Change the dropdown at the very top of the window that says "Single Symbol" to Graduated.
   

   

3. Set the Value to COUNT.
4. Note that if you use the default Method of Color, it won't do the job effectively. Instead, set the Method: change it from Color to Size.
   

   

   

5. Classes: 3 or 4. Click the Classify button. Play with the settings (such as the minimum and maximum sizes) until the map reads well.
   

   

6. Change the base marker to a red circle (you can click directly on the symbol icon to do this). Click OK.
7. View the map. Does displaying fatalities using proportional symbols reinforce or discredit John Snow's hypothesis? Are there fatalities clustered near other pumps?
   

   

Step 4: Spatial Statistics

1. In the Processing Toolbox, search for Mean coordinate(s) (usually found under Vector analysis). Ensure the Input layer is DEATHS Events and Weight field is set to COUNT. Run it.
2. How far is the Broad Street Pump from this mean center? (Use the Measure Line tool in QGIS to check).
3. In the Processing Toolbox, search for Standard deviational ellipse (sometimes known as Directional Distribution). Input layer: DEATHS Events. Weight matrix: COUNT. Standard Deviation: 1. Run it.
4. Are the majority of fatalities clustered around the Broad Street Pump within this 1 standard deviation ellipse?

Step 5: Vector Analysis (Thiessen Polygons & Spatial Joins)

John Snow needed more than visual reinforcement; he needed straight-line least-distance service areas for each water pump.

1. In the Processing Toolbox, search for Voronoi polygons (the alternative name for Thiessen Polygons).
2. Input layer: PUMPS Events. Run it.
3. Overlay the DEATHS Events points on top of these polygons. Which polygon visually contains the most fatalities?
4. To quantify this, we'll perform a Spatial Join. In the Processing Toolbox, search for Join attributes by location (summary) (or just use basic Join attributes by location).
5. Target layer: Voronoi polygons. Join layer: DEATHS Events. Summaries to calculate: sum (specifically for the COUNT field). Run it.
6. Look at the attribute table of the newly generated layer. Which pump really definitively had the most deaths in its immediate service area?

Step 6: Raster Analysis (Kernel Density)

1. Turn off the Thiessen polygons and Mean Center layers. Leave the DEATHS Events points.
2. Search the Processing Toolbox for Heatmap (Kernel Density Estimation).
3. Point layer: DEATHS Events. Radius: 100 (in meters). Weight from field: COUNT. Pixel size: 2. Run it.
4. Double-click the heatmap layer, go to Symbology, choose Singleband pseudocolor, pick a color ramp (like YlOrRd), and apply it.
5. The Kernel Density analysis creates a smooth surface. Does the resulting surface reinforce or counter Snow's hypothesis?

---

💻 Option 2: VS Code & Python Scripting (Geopandas & Copilot)

Modern spatial analysis can automate these workflows. Use VS Code, Python, and Copilot to replicate Snow's analysis.

Step 1: Environment Setup

1. Open VS Code in your lab folder containing the .dbf files.
2. Install required libraries via terminal:

   pip install geopandas libpysal scipy folium simpledbf matplotlib

Step 2: Scripting with Copilot

Create ghost_map.py. We will use simpledbf to read the legacy ArcGIS dbfs, and geopandas for analysis.

Prompt Copilot (or write) step-by-step:

"""
Goal: Analyze 1854 Cholera Data in Python
1. Read PUMPS.dbf and DEATHS.dbf into pandas DataFrames using Dbf5 from simpledbf.
2. Convert them into GeoDataFrames using X_WGS84 and Y_WGS84. Assign CRS EPSG:4326.
3. Tabular analysis: Group Deaths by 'STREET' and sum 'COUNT'. Print the top 5 streets.
4. Export the data to GeoJSON files or render them onto a Folium map.
"""

Step 3: Reference Complete Script

Here is a working script combining Data Conversion, Voronoi Analysis, and Web Visualization:

import geopandas as gpd
import pandas as pd
from simpledbf import Dbf5
import folium
from folium.plugins import HeatMap

# 1. Load Data
dbf_deaths = Dbf5('DEATHS.dbf').to_dataframe()
dbf_pumps = Dbf5('PUMPS.dbf').to_dataframe()

# 2. To GeoDataFrames
gdf_deaths = gpd.GeoDataFrame(
    dbf_deaths, 
    geometry=gpd.points_from_xy(dbf_deaths.X_WGS84, dbf_deaths.Y_WGS84),
    crs="EPSG:4326"
)
gdf_pumps = gpd.GeoDataFrame(
    dbf_pumps, 
    geometry=gpd.points_from_xy(dbf_pumps.X_WGS84, dbf_pumps.Y_WGS84),
    crs="EPSG:4326"
)

# 3. Tabular Analysis
top_streets = gdf_deaths.groupby('STREET')['COUNT'].sum().sort_values(ascending=False).head()
print("Streets with most fatalities:\n", top_streets)

# 4. Map Setup
m = folium.Map(location=[51.5133, -0.1366], zoom_start=16, tiles='CartoDB Positron')

# Add Heatmap (Kernel Density equivalent)
heat_data = [[row.geometry.y, row.geometry.x, row['COUNT']] for _, row in gdf_deaths.iterrows()]
HeatMap(heat_data, radius=15).add_to(m)

# Add Proportional Symbols
for _, row in gdf_deaths.iterrows():
    folium.CircleMarker(
        location=[row.geometry.y, row.geometry.x],
        radius=row['COUNT'] * 1.5,
        color='red', fill=True, fillOpacity=0.7,
        popup=f"{row['COUNT']} deaths"
    ).add_to(m)

# Add Pumps
for _, row in gdf_pumps.iterrows():
    folium.Marker(
        location=[row.geometry.y, row.geometry.x],
        popup=row['NAME'],
        icon=folium.Icon(color='blue', icon='tint')
    ).add_to(m)

# Save Map
m.save('ghost_map_python.html')
print("\nSuccess! Open ghost_map_python.html in your browser.")

---