🚧 Under construction — VoTechWorks is being built. A free tech program — complimentary, as our commitment to get America working.
voTechWorksGrade 13 · a vocational technology learning center, by urbanicity
newSchoolWorks.com › VoTechWorks › Technology & Design › Data Science & Analytics › GIS and Spatial Data Analysis › Lesson 6
TEC-DATA85-09 · Lesson 6 · 16 clock hours · lab

Python for Spatial Data: GeoPandas, Shapely, and Folium

Objective. Students will read, manipulate, and spatially analyze vector data using GeoPandas and Shapely in a Jupyter Notebook, and render an interactive web map using Folium.

⚠ Pending professional review

This lesson, its worked examples and its practice questions has not yet been verified by an industry-certified professional or by a certification or approval body. It is a study draft and may contain errors. Do not use it as instructions for work on a job site or on live equipment. Always follow your instructor, your supervising professional, and current official references and regulations.

VoTechWorks program approval and accreditation are still in progress. Credentials are issued only by the certifying bodies named in each program, after their own exams.

Python for Spatial Data: GeoPandas, Shapely, and Folium

Geographic Information Systems (GIS) analysis no longer requires expensive desktop software. Python's open-source ecosystem provides everything you need to read, manipulate, analyze, and visualize spatial vector data inside a Jupyter Notebook. In this lab you will use three libraries that work together seamlessly: GeoPandas, Shapely, and Folium.

What Are Vector Data?

Vector data represent geographic features as points, lines, or polygons stored with attribute information. A library location is a point. A road network is a collection of lines. A census tract boundary is a polygon. The GeoJSON format is a common, human-readable text format for vector data, and it is the format you will use in this lab. US Census TIGER/Line shapefiles and open-data portals both distribute GeoJSON files freely.

GeoPandas: The Foundation

GeoPandas extends the familiar Pandas DataFrame by adding a geometry column. The result is called a GeoDataFrame. Every row in a GeoDataFrame holds one feature, its attributes, and its geometry object. You interact with attributes exactly as you would in Pandas — boolean indexing, `groupby`, `merge` — while GeoPandas handles the spatial side.

To read a GeoJSON file:

```python

import geopandas as gpd

tracts = gpd.read_file("census_tracts.geojson")

print(tracts.crs) # always check the coordinate reference system

print(tracts.head())

```

The `crs` attribute tells you the Coordinate Reference System (CRS) — the mathematical model that maps longitude/latitude or projected coordinates to positions on Earth. A common geographic CRS is EPSG:4326 (WGS 84), which stores coordinates in decimal degrees. For area or distance calculations you must reproject to a projected CRS that uses meters, for example EPSG:3857 (Web Mercator) or a state-appropriate EPSG code.

To reproject and compute area:

```python

tracts_m = tracts.to_crs(epsg=3857) # reproject to meters

tracts_m["area_m2"] = tracts_m.geometry.area # area in square meters

```

Never compute area in EPSG:4326; the result will be in squared degrees and will be meaningless.

Shapely: Geometry Under the Hood

Shapely provides the geometry objects — `Point`, `LineString`, `Polygon`, `MultiPolygon` — that populate the geometry column. In most workflows you do not import Shapely directly; GeoPandas calls it automatically. However, understanding Shapely is important when you construct custom geometries or run predicates such as `contains`, `intersects`, or `within`.

```python

from shapely.geometry import Point

pt = Point(-87.6298, 41.8781) # longitude, latitude — note the order

```

Shapely predicates are the engine behind spatial joins.

Spatial Joins with `gpd.sjoin`

A spatial join combines two GeoDataFrames based on a geometric relationship rather than a shared key column. The function `gpd.sjoin` is the primary tool.

```python

libs = gpd.read_file("libraries.geojson").to_crs(epsg=4326)

tracts = gpd.read_file("census_tracts.geojson").to_crs(epsg=4326)

joined = gpd.sjoin(libs, tracts, how="left", predicate="within")

```

This attaches each library point to the census tract that contains it. Both GeoDataFrames must share the same CRS before the join — always verify with `assert libs.crs == tracts.crs`. After joining, you count libraries per tract with `groupby` and merge the counts back to the tract GeoDataFrame for mapping.

Choropleth Mapping

A choropleth map shades polygons by a numeric variable — here, the count of libraries per tract. GeoPandas supports a quick static choropleth via Matplotlib:

```python

tracts.plot(column="library_count", cmap="YlOrRd", legend=True)

```

For an interactive map you use Folium.

Folium: Interactive Web Maps

Folium wraps the Leaflet.js JavaScript library, allowing you to build interactive maps entirely in Python and export them as self-contained HTML files viewable in any web browser — no server required.

```python

import folium

m = folium.Map(location=[41.88, -87.63], zoom_start=11)

folium.Choropleth(

geo_data=tracts,

data=tracts,

columns=["GEOID", "library_count"],

key_on="feature.properties.GEOID",

fill_color="YlOrRd",

legend_name="Libraries per Census Tract"

).add_to(m)

m.save("libraries_map.html")

```

Add `GeoJsonTooltip` to display the tract ID and library count when a user hovers over a polygon. Open the exported HTML file in your browser to verify tooltips and pan/zoom behavior at multiple zoom levels.

Lab Safety and Best Practices

This is a software lab; physical hazards are minimal. Follow these professional practices:

Key terms8

GeoDataFrameA GeoPandas data structure that extends a Pandas DataFrame with a geometry column, enabling spatial operations alongside tabular attribute management.
Coordinate Reference System (CRS)A mathematical framework that defines how coordinates map to real positions on the Earth's surface; examples include EPSG:4326 (geographic, degrees) and EPSG:3857 (projected, meters).
Spatial JoinAn operation that combines two GeoDataFrames by matching features based on a geometric relationship — such as 'within' or 'intersects' — rather than a shared attribute key.
ShapelyA Python library that provides and manipulates geometric objects (Point, LineString, Polygon) and supports spatial predicates such as contains, intersects, and within.
Choropleth MapA thematic map that shades or colors polygons according to the value of a numeric variable, making spatial patterns in aggregated data visually apparent.
GeoJSONAn open, text-based format for encoding geographic features and their attributes as JSON; widely used by open-data portals and web mapping tools.
ReprojectionThe process of transforming a GeoDataFrame's coordinates from one CRS to another, required before computing measurements such as area or distance in metric units.
FoliumA Python library that generates interactive Leaflet.js web maps and exports them as standalone HTML files, supporting layers, tooltips, and choropleths.

Worked examples3

1. A GeoDataFrame of census tracts is in EPSG:4326. A tract's polygon has an area of 0.0021 square degrees as reported by Shapely. After reprojecting to EPSG:3857, the same polygon reports an area of 1,850,000 square meters. Convert that area to square kilometers.
  1. Identify the area in square meters: 1,850,000 m².
  2. Recall the conversion: 1 km² = 1,000,000 m².
  3. Divide: 1,850,000 ÷ 1,000,000 = 1.85 km².
Answer: The tract area is 1.85 square kilometers. The value in square degrees (0.0021) is meaningless for area measurement and must never be reported.
2. After a spatial join, you use groupby to count libraries per census tract. The groupby result shows tract '17031010100' has 3 libraries and tract '17031010200' has 0 libraries (it never appeared in the join output). You merge the count Series back to the tract GeoDataFrame with a left merge. What value will tract '17031010200' have in the library_count column, and how should you handle it before mapping?
  1. A left merge keeps all rows from the tract GeoDataFrame, even tracts absent from the groupby result.
  2. Tracts with no matching libraries will receive NaN in the library_count column after the merge.
  3. Fill NaN with 0 using: tracts['library_count'] = tracts['library_count'].fillna(0).
  4. Convert the column to integer type if needed: tracts['library_count'] = tracts['library_count'].astype(int).
Answer: Tract '17031010200' will have NaN after the merge. Fill it with 0 using fillna(0) so that zero-library tracts display correctly on the choropleth instead of appearing as missing data.
3. You have 45 libraries spread across 12 census tracts. After the spatial join and groupby, the sum of all library_count values equals 43, not 45. What are two likely causes, and how do you investigate each?
  1. Cause 1 — CRS mismatch: if the library points and tract polygons had different CRS values when sjoin ran, some points may have fallen outside all polygons. Investigate by verifying libs.crs == tracts.crs before the join; reproject one GeoDataFrame to match the other if they differ.
  2. Cause 2 — Points on or outside boundaries: two library points may lie exactly on a tract boundary or slightly outside all tracts (data quality issue). Investigate by checking the sjoin result for rows where the tract GEOID is NaN (indicating no match); map those unmatched points to see if they fall outside the tract layer extent.
  3. Fix: after confirming the same CRS, use predicate='within' or try predicate='intersects' for boundary-straddling points, and inspect unmatched features before deciding how to handle them.
Answer: Two missing libraries likely result from either a CRS mismatch causing incorrect alignment, or library points that lie on polygon boundaries or outside the tract layer. Always assert matching CRS before a spatial join and check for unmatched rows in the output.

The math in this lesson1

Learners convert tract areas from square meters to square kilometers in a projected CRS.
N-Q.1
Brush up in the math halls at newSchoolWorks.com:
In class: Code-along: instructor opens a Jupyter Notebook and walks through reading a GeoJSON file into a GeoDataFrame, inspecting CRS, filtering rows with Pandas-style boolean indexing, and computing a new column (e.g., area in square meters using projected CRS); students type and run each cell, adding markdown commentary. · Independent coding exercise: students use GeoPandas to perform a spatial join between a point GeoDataFrame (e.g., library locations) and a polygon GeoDataFrame (census tracts), compute the count of libraries per tract, and merge the result back to the tract GeoDataFrame for choropleth mapping. · Folium mapping task: students build an interactive Folium map displaying their choropleth layer with a tooltip showing tract ID and library count, add a tile layer, and export the map as an HTML file; they open it in a browser and screenshot three zoom levels.
Materials: Python 3.x with GeoPandas, Shapely, Folium, Matplotlib installed (via course Anaconda environment), Jupyter Notebook, GeoPandas documentation (geopandas.org), Folium documentation (python-visualization.github.io/folium), US Census TIGER/Line GeoJSON files, library location GeoJSON from city open-data portal
Check: Students submit their Jupyter Notebook (.ipynb) with all cells executed and the exported Folium HTML file. Instructor checks that the spatial join produces correct feature counts, the choropleth classification is applied correctly, and the Folium map loads and displays tooltips; pass mark: all three deliverables correct, notebook runs top-to-bottom without errors.

Practice check8

Answer each question, then check it. Nothing is recorded; this is your own practice.

Data Manipulation and Analysis
1. Before computing the area of polygon features in a GeoPandas GeoDataFrame, what step is always required?
Data Manipulation and Analysis
2. Which GeoPandas function combines two GeoDataFrames based on a geometric relationship such as 'within' or 'intersects'?
Data Integrity and Quality
3. A student runs gpd.sjoin(libs, tracts, how='left', predicate='within') but finds that many library points have NaN for tract attributes in the output. What is the most likely cause?
Data Manipulation and Analysis
4. After a spatial join and groupby count, a student merges the library count back to the tract GeoDataFrame. Tracts with no libraries show NaN in the count column. What is the correct next step?
Data Visualization
5. What does the Folium library produce as its primary output?
Data Workflow and Reproducibility
6. A student's Jupyter Notebook runs without errors when cells are executed individually but fails when the instructor restarts the kernel and runs all cells top-to-bottom. What is the most likely explanation?
Spatial Data Concepts
7. Which EPSG code represents the WGS 84 geographic coordinate reference system that stores coordinates in decimal degrees of longitude and latitude?
Data Integrity and Ethics
8. In the context of this lab, why is recording the source URL and download date of GeoJSON files in a notebook markdown cell considered a best practice?
Record this lesson on my transcript. Sign in with your key and enrollment email to record lessons (enroll free if you have no key). This page stays open to read and practise.
← Previous lessonProgram planNext lesson →