Vector Analysis

In the previous section, we learned how to ask spatial questions to filter our data (e.g., finding all primary homes inside an evacuation zone).

Now, we will take the final step in vector analysis: Data Fusion. You will learn how to permanently combine the geometries and attributes of entirely different datasets to engineer new, multi-dimensional analytical tables.

1. The Limitation of Pandas Merges¶

In standard data science, when you want to combine two tables, you use pd.merge(). However, a standard table merge requires a shared column (like a unique ID or a matching name) to act as the “key” between the two datasets.

But what happens when your datasets have no matching columns?

Imagine you have a dataset of Alpine huts containing only names and GPS coordinates, and a dataset of Swiss cantons containing population statistics. If you want to know exactly which canton each hut belongs to, a standard Pandas merge is completely useless. The datasets share no text or ID columns.

To solve this, we must use geography itself as the joining key!

2. Spatial Joins (sjoin)¶

A Spatial Join (.sjoin()) transfers the attributes of one layer to another based on their spatial relationship.

If a point falls inside a polygon, GeoPandas will grab the attributes of that polygon and permanently attach them to the point’s row in the attribute table.

Let us answer the question: Which Swiss canton has the highest number of Alpine huts? We will spatially join the cantons to the huts, and then count the results.

The Anatomy of a Spatial Join¶

If you look at the .sjoin() code above, you will notice two very specific arguments: predicate="within" and how="inner". These two parameters give you total control over how the spatial data is combined, and changing them can drastically alter your results.

Let us use the diagrams below to understand exactly how these parameters work.

1. The Spatial Predicate (predicate) The predicate defines the strict geometric rule for the join. It is crucial to pick a rule that makes logical sense (e.g., checking if points are within a polygon makes sense, but checking if points contain a polygon does not).

• within (See Result Table i): This requires the point to be strictly inside the boundary. Look at Point 6: because it sits exactly on the border of Polygon C, it is not considered “inside”, so it is entirely dropped from the results in table i.

• intersects (See Result Table ii): This is more forgiving. If geometries share even a single pixel of space, it is a match. Because Point 6 touches the border of Polygon C, it successfully joins and appears in table ii.

2. The Join Type (how) The join type determines what happens to geometries that completely fail the spatial predicate rule (like Point 3, which is floating out in the middle of nowhere).

• Inner Join (how="inner"): The default, strict option. As seen in tables i and ii, Point 3 is completely deleted from the final dataset because it could not find a matching polygon. If an Alpine hut is located just across the border in Italy, it is dropped.

• Left Outer Join (how="left"): The safe option. Look at result table iii. The “Left” dataset (Layer 1: Points) is prioritized. Every single point is kept. Point 3 survives, but because it has no overlapping polygon, its new columns simply say No data (or NaN in Pandas).

• Right Outer Join (how="right"): Look at result table iv. The “Right” dataset (Layer 2: Polygons) is prioritized. Every single polygon is kept. Polygon D has zero points inside it, but it is kept in the final table with No data in the point columns.

Concept Check: The Disappearing Data¶

Scenario: You have a Point dataset of all 300 hospitals in Switzerland and a Polygon dataset of high-risk avalanche zones. You want to create a safety report that lists every single hospital in the country, with an extra column indicating if it is in an avalanche zone or not.

Which join type must you use to ensure hospitals in safe zones are not deleted from your final table?

A) how="inner"
B) how="left"
C) how="right"

Aggregating the Fused Data¶

Now that every individual hut has a canton name securely attached to it, we can use a standard Pandas .groupby() to count them up!

By fusing the attributes based on location, we seamlessly transition from spatial data back into standard statistical analysis.

3. Dissolving Geometries (dissolve)¶

If .sjoin() is the spatial equivalent of a Pandas merge, then .dissolve() is the spatial equivalent of a Pandas groupby.

Dissolving takes multiple adjacent geometries that share the same attribute and physically melts their borders together into a single, massive geometry. Simultaneously, it allows you to aggregate (sum, average, etc.) the numeric attributes of those melted shapes.

Let us test this. We will load the Swiss Municipalities dataset. First, we will use a spatial join to figure out which canton each municipality belongs to. Then, we will dissolve the thousands of tiny municipalities back into large cantonal boundaries, summing up their population (einwohnerzahl) and surface area along the way!

Because we just engineered entirely new polygons and summed up their data from scratch, it is good practice to verify our math. Let us compare our cantons_rebuilt data against the official statistics stored in the original cantons_gdf:

The numbers match perfectly! Best quality made by SwissTopo. By dissolving the municipalities, we successfully recreated the cantonal dataset from the ground up.

We can visually prove that this operation physically melted the geometries together by plotting our original municipalities next to our newly dissolved cantons_rebuilt layer:

4. Spatial Overlays (overlay)¶

Spatial Joins (sjoin) only transfer attributes; they never alter the actual shapes of your geometries. But what if you need to physically cut a geometry into smaller pieces?

A Spatial Overlay (.overlay()) mathematically intersects two layers, physically slicing the geometries exactly where they cross boundaries to create entirely new shapes.

When performing an overlay, you use the how parameter to control exactly which parts of the sliced geometries are kept in the final dataset. The five main operations are:

• intersection: Keeps only the areas where the two layers overlap.

• union: Keeps all geometries from both layers, splitting them into new distinct pieces wherever their borders cross.

• difference: Keeps only the areas of the first layer that do not overlap with the second layer.

• symmetric_difference: Keeps the areas from both layers that do not overlap (the exact opposite of intersection).

• identity: Computes a geometric intersection, but keeps all features of the primary input layer, splitting them where they overlap the second layer.

Intersection

Union

Difference

Symmetric Difference

Identity

Applying an Intersection Overlay¶

To demonstrate why this is vital, let us calculate the total length of highways in each canton. A highway is typically a single, continuous line that spans across multiple cantons. If we just did a spatial join, the entire length of the highway would be assigned to just one canton. We must use .overlay() to physically cut the highway lines at the cantonal borders first!

Because .overlay() computationally slices thousands of vertices, it is one of the heaviest operations in GeoPandas. However, it is absolutely essential when dealing with long networks like rivers, roads, or animal migration paths that cross administrative boundaries.

Concept Check: The Overlay Trap¶

Scenario: You have a Polygon representing a large farm. In its attribute table, the area_sqkm column says 10. You use .overlay(how="intersection") to physically cut the farm exactly in half using a municipal border.

If you immediately look at the attribute table of your two newly cut farm pieces, what will their area_sqkm columns say?

A) 5
B) 10
C) NaN (No data)

5. Exercise: Aggregating by Area¶

It is time to put your fusion skills to the test. In the event of a severe emergency, helicopter transport is crucial. You have been asked to determine which Swiss canton has the highest number of registered hospital landing sites.

You will need to use a Spatial Join to combine the points and polygons, and then aggregate the results.

Tasks:

1. Load Data: Load hospital_landing_sites.geojson and swissBoundaries3D_cantons.gpkg. Ensure both are in the metric EPSG:2056 projection.

2. Spatial Join: Perform an sjoin keeping the landing sites as the primary (left) geometry, and joining the canton attributes based on the within predicate.

3. Aggregate: Use .groupby() on the newly attached canton name column and count the number of landing sites per canton.

4. Identify: Sort the results in descending order and print the top 5 cantons.

Sample solution

import geopandas as gpd

# 1. Load the data and project to metric grid
hospitals_gdf = gpd.read_file("hospital_landing_sites.geojson").to_crs(epsg=2056)
cantons_gdf = gpd.read_file("swissBoundaries3D_cantons.gpkg").to_crs(epsg=2056)

# 2. Perform the spatial join (Points within Polygons)
hospitals_with_cantons = gpd.sjoin(
    hospitals_gdf, 
    cantons_gdf, 
    how="inner", 
    predicate="within"
)

# 3. Group by the new canton name column ('name_right') and count
landing_site_counts = hospitals_with_cantons.groupby("name_right").size()

# 4. Sort and display the top 5
print("Top 5 Cantons by Hospital Landing Sites:")
display(landing_site_counts.sort_values(ascending=False).head(5))

# --- BONUS: Visualizing the Aggregation ---

# A. Convert our grouped counts (Series) back into a DataFrame
counts_df = landing_site_counts.reset_index(name="hospital_count")

# B. Merge those counts back onto the original cantons polygons
cantons_with_counts = cantons_gdf.merge(
    counts_df, 
    left_on="name", 
    right_on="name_right", 
    how="left"
)

# Fill any cantons that had zero hospitals with a 0 instead of NaN
cantons_with_counts["hospital_count"] = cantons_with_counts["hospital_count"].fillna(0)

# C. Draw the base map colored by the hospital_count column
ax = cantons_with_counts.plot(
    column="hospital_count", 
    cmap="Blues", 
    edgecolor="darkgrey", 
    linewidth=0.5, 
    figsize=(10, 6),
    legend=True,
    legend_kwds={"label": "Number of Landing Sites"}
)

# D. Overlay the exact hospital point locations in red
hospitals_gdf.plot(ax=ax, color="red", marker="+", markersize=50)

ax.set_title("Hospital Landing Sites per Canton")
ax.set_axis_off()

name_right	count
Bern	17
Graubünden	15
Zürich	13
Vaud	12
St. Gallen	11

6. Summary: Fusing Data¶

Congratulations! You have mastered the core toolkit of vector spatial analysis. You are no longer just looking at maps; you are actively engineering new spatial insights by combining datasets that otherwise could not communicate.

Key takeaways¶

• The Power of Geography: When datasets lack matching ID columns, geography acts as the universal key to connect them.

• Spatial Joins (.sjoin()): Use this to transfer attribute data from one layer to another based on spatial relationships (like points inside polygons) without altering the physical geometries.

• Dissolving (.dissolve()): The spatial equivalent of groupby. Use this to physically melt adjacent boundaries together while aggregating their numeric attributes.

• Overlays (.overlay()): Use this to physically cut and slice geometries at boundaries. Always remember to recalculate your areas and lengths after performing an overlay, as the physical shapes have permanently changed!

What comes next?¶

You have successfully engineered new spatial datasets and aggregated complex statistics across regions. But how do we effectively communicate these numbers to an audience?

If you look back at the hospital map we just created, we used a continuous color gradient. While it looks nice, the human brain actually struggles to interpret the exact difference between two similar shades of blue.

To make our maps truly readable, we need to group our raw numbers into logical “bins” (e.g., “Low”, “Medium”, and “High” density). In the next section, Data Classification, you will learn how to bridge the gap between raw data and clear visual communication using rule-based logic and statistical grouping methods!