Spatial Relationships

In the previous sections, we learned how to load shapes, project them onto a metric grid, and calculate their areas and buffers. However, we only manipulated one dataset at a time.

In this section, we transition to answering complex spatial questions by comparing two entirely different datasets. You will learn how to evaluate the physical relationships between geometries and use those relationships to filter data and find nearest neighbors.

1. Asking Spatial Questions¶

When we look at a map, our brains instantly recognize spatial relationships. We can easily see if a dot is inside a square or if two lines cross each other. Computers, however, only see lists of coordinates.

To teach the computer how to recognize these visual relationships, spatial data science uses Topological Spatial Relations. These are mathematical rules that describe how the interior, boundary, and exterior of two geometric shapes interact.

When you ask GeoPandas a spatial question (e.g., “Are these two shapes touching?”), it evaluates these mathematical rules under the hood and returns a simple Boolean answer: True or False. We refer to these evaluation methods as spatial predicates.

While there are many possible mathematical combinations, GeoPandas provides 10 core methods to test these relationships.

The Core Three¶

You will use these three spatial predicates for the vast majority of your spatial queries:

• .intersects(): Returns True if the geometries share at least one point in common. This is the broadest and most commonly used query. (If a shape touches, crosses, or is within another shape, it also intersects it).

• .within(): Returns True if geometry A is completely inside the boundary of geometry B.

• .contains(): Returns True if geometry A completely encloses geometry B. (This is the exact inverse of .within()).

Other Useful Predicates¶

Depending on the specific rules of your analysis, you might need more precise constraints:

• .disjoint(): Returns True if the geometries are completely separated. (The exact opposite of .intersects()).

• .touches(): Returns True if the geometries share a border/edge, but their interiors do not overlap at all.

• .overlaps(): Returns True if the geometries’ interiors intersect partially, but neither is completely inside the other.

• .equals(): Returns True if both geometries are topologically perfectly identical.

• .crosses(): Returns True if the geometries intersect, but their intersection creates a shape of a lower dimension (e.g., a line crossing a polygon creates a line, not a polygon).

• .covers() / .covered_by(): Very similar to contains and within, but slightly more lenient because they allow the shapes to share exact boundaries.

2. Intersects, Contains, and Within¶

Let us test these core relationship methods using two distinct datasets: the polygonal borders of the Swiss cantons, and the point coordinates of the MeteoSwiss pollen monitoring stations.

First, we will load both datasets, ensure they share the exact same Coordinate Reference System (EPSG:2056), and extract a single canton and a single station to compare.

The Geometry Methods¶

Because canton_zh is a Shapely Polygon and station_zh is a Shapely Point, they possess built-in methods to test their relationships against one another.

To understand how these methods work, we must test them in both directions. Does the order in which we ask the question matter? Let us ask the computer six specific questions:

Output:

Canton contains station: True
Station contains canton: False

Station within canton: True
Canton within station: False

Station intersects canton: True
Canton intersects station: True

This output reveals a critical rule of spatial querying: Order matters immensely for most predicates! * Asymmetric Methods: .contains() and .within() are directional. A massive polygon can contain a tiny coordinate point, but a point can never contain a polygon. Notice, however, that they are perfect inverse methods of each other. Use whichever makes your code read more naturally.

• Symmetric Methods: .intersects() is bi-directional. If shape A touches shape B, then shape B must be touching shape A. The order does not matter.

Concept Check: The Rhine River¶

You have a LineString geometry representing the entire Rhine River (flowing from the Swiss Alps all the way to the Netherlands) and a Polygon geometry representing the country of Switzerland (swiss_poly).

Based on the mathematical rules of spatial predicates, which of the following queries will return True?

A) rhine_line.within(swiss_poly)
B) swiss_poly.within(rhine_line)
C) swiss_poly.contains(rhine_line)
D) rhine_line.intersects(swiss_poly)

Scaling Up: The Spatial Loop¶

Checking one station against one canton is easy. But what if we want to calculate a brand new regional metric for our dataset?

For this example, we will calculate the average altitude of the pollen stations inside each individual canton. (You can imagine altitude here as a tangible substitute for any real-world scientific measurement you might want to aggregate, such as daily pollen concentrations, temperature, or rainfall!)

To do this, we can loop through every single canton and ask the computer to find which stations fall .within() its borders. If it finds any, we will calculate their average altitude and save it directly into the cantons_gdf attribute table.

3. Filtering by Location¶

In the previous section, we used a for loop to test geometries one by one. While loops are useful, they are very slow for large datasets. The true power of GeoPandas shines when we apply spatial predicates to entire DataFrames at once.

In early chapters, you learned how to filter Pandas DataFrames using logical conditions (e.g., data[data["population"] > 1000]). In GeoPandas, you can use geography itself as the logical condition!

Let us answer a critical safety question: If a nuclear emergency required a 16 km evacuation zone around the Beznau nuclear power plant, how many primary homes in the surrounding communes would be affected?

We will load the Swiss housing inventory (which contains communal polygons and household counts in a column called ZWG_3010) and the Nuclear Power Plant dataset.

Now, instead of writing a loop, we will ask the entire housing_gdf dataset which communal polygons intersect with our evacuation zone. This generates a Boolean mask—a simple list of True and False values—which we can use to instantly filter the rows!

With just two lines of spatial filtering logic, we instantly calculated demographic vulnerability based purely on geographic location!

To prove that our spatial filter worked perfectly, let us visualize the results. We will plot all the communes in light grey, overlay our filtered impacted_communes in red, and draw the 16km evac_zone boundary on top.

4. Nearest Neighbor¶

Sometimes geometries do not touch or overlap at all, but you still need to know their spatial relationship. A common task is finding the absolute closest feature between two different datasets.

Instead of writing complex mathematical loops to calculate the distance between every single point, GeoPandas provides the incredibly fast .sjoin_nearest() method (Spatial Join Nearest). It calculates the distance from every geometry in your primary dataset to the geometries in a secondary dataset, and immediately joins the attribute data of the closest match.

Let us find the closest hospital landing site for every Swiss Alpine Club (SAC) mountain hut in case of a helicopter medical emergency.

Output:

This operation requires a massive amount of background math, but thanks to advanced spatial indexing (KD Trees) working under the hood, GeoPandas can cross-reference and match thousands of points in milliseconds.

Visualizing the Connections

To truly understand what the computer just did, we can draw a map that explicitly connects every mountain hut to its newly assigned hospital.

Do not worry about the plotting syntax right now! We will formally introduce advanced mapping (matplotlib) in the next session. Treat this as an optional “sneak peek” into what you will be able to do next week.

To draw these connecting lines, we use the index_right column (which .sjoin_nearest() automatically creates to remember which hospital it matched with) to retrieve the exact hospital coordinates.

from shapely.geometry import LineString
import matplotlib.pyplot as plt

# 1. Bring the hospital's geometry back into our results table
emergency_plan = emergency_plan.merge(
    hospitals_gdf[["geometry"]], 
    left_on="index_right", 
    right_index=True, 
    suffixes=("", "_hospital")
)

# 2. Create a brand new "flight path" line connecting the hut to the hospital
emergency_plan["flight_path"] = emergency_plan.apply(
    lambda row: LineString([row["geometry"], row["geometry_hospital"]]), axis=1
)

# 3. Tell GeoPandas to use our new flight paths as the active geometry for plotting
flight_paths_gdf = emergency_plan.set_geometry("flight_path")

# 4. Plot the results!
fig, ax = plt.subplots(figsize=(10, 6))

# Plot the flight paths in grey
flight_paths_gdf.plot(ax=ax, color="grey", linewidth=0.5, linestyle="--");

# Plot the huts as blue triangles
huts_gdf.plot(ax=ax, color="blue", marker="^", markersize=15, label="SAC Huts");

# Plot the hospitals as red crosses
hospitals_gdf.plot(ax=ax, color="red", marker="+", markersize=60, label="Hospitals");

ax.set_title("Emergency Helicopter Flight Paths (Nearest Hospital)")
ax.axis("off")
plt.legend()
plt.show()

Advanced Tip: Search Limits¶

Helicopters have limited fuel, and in medical emergencies, time is critical. Look at the Mischabelhütte in the output above—it is almost 20 kilometers away from the nearest hospital! What if a rescue company only wants to find hospitals that are within a strict 15-kilometer flight radius?

You can make your nearest neighbor search even faster by passing the max_distance parameter:

# Only find hospitals if they are within 15,000 meters
emergency_plan = huts_gdf.sjoin_nearest(
    hospitals_gdf, 
    distance_col="distance_m",
    max_distance=15000 
)

(Note: If you ran this code, the Mischabelhütte would simply be dropped from the final results, as it has no hospitals within that 15km search radius!)

5. Exercise: Point in Polygon¶

It is time to test your spatial querying skills. We have GPS tracking data from wolves across Switzerland, but we are only interested in studying their movements when they are safely inside protected national park boundaries.

Your task is to use a spatial predicate as a Boolean mask to filter the GPS tracks.

Tasks:

1. Load Data: Load the two datasets and ensure both are in EPSG:2056.

   • Real-World Tip: Both of these files actually contain multiple spatial layers! To suppress GeoPandas’ warnings, you must specify exactly which layer to open. Use layer="Schneeschuhwanderwege" for the wolf tracks (apparently our wolves like to follow the snowshoe trails!), and layer="N2026_Revision_Park_20250430" for the parks.

2. Merge Parks: The parks dataset contains dozens of individual polygons. To make our query faster, merge them into a single massive “protected zone” geometry using all_parks = parks_gdf.geometry.union_all().

3. Filter by Location: Create a Boolean mask by checking which geometries in your wolf tracks dataset are .within() the all_parks geometry.

4. Apply Mask: Filter the wolf GeoDataFrame using your mask to create a new variable called protected_wolves.

5. Count: Use len() to see how many GPS tracks fell entirely inside the park boundaries compared to the original dataset.

Sample solution

import geopandas as gpd

# 1. Load Data (Specifying the exact layer prevents warnings!)
wolves_gdf = gpd.read_file("wolf_tracks_switzerland.gpkg", layer="Schneeschuhwanderwege").to_crs(epsg=2056)
parks_gdf = gpd.read_file("swiss_parks.shp.zip", layer="N2026_Revision_Park_20250430").to_crs(epsg=2056)

# 2. Merge Parks into a single massive geometry
all_parks = parks_gdf.geometry.union_all()

# 3. Create the spatial boolean mask
inside_park_mask = wolves_gdf.geometry.within(all_parks)

# 4. Filter the dataframe using the mask
protected_wolves = wolves_gdf[inside_park_mask]

# 5. Print counts
print(f"Original GPS tracks: {len(wolves_gdf)}")
print(f"Tracks strictly inside parks: {len(protected_wolves)}")

# Optional Plotting
ax = parks_gdf.plot(color="lightgreen", edgecolor="darkgreen", figsize=(10, 6));
wolves_gdf.plot(ax=ax, color="lightgrey", markersize=1, alpha=0.5);
protected_wolves.plot(ax=ax, color="red", markersize=5);
ax.set_title("Wolf Tracks Inside Swiss National Parks")
ax.axis("off")

6. Summary: Spatial Querying¶

In this section, you unlocked the ability to link disparate datasets together purely based on where they exist in the real world. From calculating nuclear evacuation zones to routing emergency helicopters and tracking wolves, you have learned the foundational logic of all advanced geographic information systems.

Key takeaways¶

• Topological Rules: Spatial predicates are mathematical evaluations that return True or False depending on how shapes physically interact (or don’t).

• The Core Three: Memorize .intersects(), .contains(), and .within(). They are your primary tools for querying space. Always remember that order matters when using asymmetric predicates like contains and within!

• Geographic Filtering: Just like filtering text or numbers in standard Pandas, you can use a spatial predicate to generate a Boolean mask and filter an entire GeoDataFrame at once—no slow for loops required.

• Nearest Neighbor: When items do not touch, use .sjoin_nearest() to rapidly calculate the closest feature between two datasets and merge their attributes. Remember to use max_distance to set logical real-world limits (like maximum flight ranges) and speed up your code.

What comes next?¶

We have successfully used spatial relationships to filter data (e.g., keeping only the wolves inside the park). But what if we want to permanently merge the columns of two overlapping datasets together?

In the final section of this chapter, Spatial Joins, we will explore how to fuse datasets based on their spatial relationships to create entirely new, multi-dimensional analytical tables!