Perform a spatial join

Use a GeometryOps predicate in a FlexiJoins join to match rows by their geometries.

FlexiJoins is a separate package. If you installed only GeoDataFrames, install and load it before following this guide:

import Pkg; Pkg.add("FlexiJoins")
using FlexiJoins

Join points to polygons

For a point-in-polygon join, put points on the left and use GeometryOps.within against polygons on the right:

using GeoDataFrames
using FlexiJoins
using NaturalEarth

map_units = select(
    DataFrame(naturalearth("admin_0_map_units", 10)),
    :ADMIN => :country,
    :NAME => :map_unit,
    :geometry,
)
cities = select(
    DataFrame(naturalearth("populated_places", 50)),
    :NAME => :city,
    :geometry,
)

zones = subset(
    map_units,
    :country => ByRow(
        name -> name in ["Belgium", "Luxembourg", "Netherlands"],
    ),
    :map_unit => ByRow(!=("Caribbean Netherlands")),
)
points = subset(
    cities,
    :city => ByRow(
        name -> name in ["Amsterdam", "Berlin", "Brussels", "Luxembourg", "Paris"],
    ),
)

joined = innerjoin((points, zones), by_pred(:geometry, GeometryOps.within, :geometry))
joined

3×5 DataFrame

Row	city	geometry	country	map_unit	geometry_1
	String	Point…	String	String	Union…
1	Luxembourg	2D Point	Luxembourg	Luxembourg	2D Polygon
2	Brussels	2D Point	Belgium	Brussels	2D Polygon
3	Amsterdam	2D Point	Netherlands	Netherlands	2D MultiPolygon

The inner join matches Amsterdam, Brussels, and Luxembourg to their countries. Berlin and Paris do not match a Benelux country.

Plot joined and unmatched points

Visualize which points joined to at least one polygon and which did not:

using CairoMakie

matched_cities = unique(joined.city)
matched = subset(points, :city => ByRow(city -> city in matched_cities))
unmatched = subset(points, :city => ByRow(city -> !(city in matched_cities)))

fig = plot(
    zones.geometry;
    color = (:dodgerblue, 0.12),
    strokecolor = :dodgerblue,
    strokewidth = 2,
    axis = (; title = "Cities joined to Benelux countries"),
)
plot!(unmatched.geometry; color = :gray, markersize = 14)
plot!(matched.geometry; color = :tomato, markersize = 14)
fig

Keep unmatched rows

Choose the join shape according to the rows you need to retain:

leftjoin((points, zones), by_pred(:geometry, GeometryOps.within, :geometry))

5×5 DataFrame

Row	city	geometry	country	map_unit	geometry_1
	String	Point…	String?	String?	Abstract…?
1	Luxembourg	2D Point	Luxembourg	Luxembourg	2D Polygon
2	Brussels	2D Point	Belgium	Brussels	2D Polygon
3	Amsterdam	2D Point	Netherlands	Netherlands	2D MultiPolygon
4	Berlin	2D Point	missing	missing	missing
5	Paris	2D Point	missing	missing	missing

leftjoin keeps every point and fills the right-side fields with missing when no zone matches. The same condition also supports:

rightjoin((points, zones), by_pred(:geometry, GeometryOps.within, :geometry))
outerjoin((points, zones), by_pred(:geometry, GeometryOps.within, :geometry))

7×5 DataFrame

Row	city	geometry	country	map_unit	geometry_1
	String?	Point…?	String?	String?	Abstract…?
1	Luxembourg	2D Point	Luxembourg	Luxembourg	2D Polygon
2	Brussels	2D Point	Belgium	Brussels	2D Polygon
3	Amsterdam	2D Point	Netherlands	Netherlands	2D MultiPolygon
4	Berlin	2D Point	missing	missing	missing
5	Paris	2D Point	missing	missing	missing
6	missing	missing	Belgium	Flemish	2D MultiPolygon
7	missing	missing	Belgium	Walloon	2D MultiPolygon

rightjoin keeps every zone; outerjoin keeps rows from both inputs.

Check CRS and candidate filtering

Reproject both tables to the same CRS before joining. GeometryOps and FlexiJoins interpret coordinates as planar x/y values, even when a table records a geographic CRS.

FlexiJoins indexes the second, usually more complex, table with an STR tree. Bounding boxes filter candidate pairs before the exact predicate runs. This reduces candidate checks but does not replace the predicate, and GeoDataFrames' GeometryVector cache is not consumed by current FlexiJoins joins. For join-condition variants, consult the FlexiJoins documentation.