Spatial joins

Project.toml

@@ -8,17 +8,21 @@ CoordinateTransformations = "150eb455-5306-5404-9cee-2592286d6298"
 GeoInterface = "cf35fbd7-0cd7-5166-be24-54bfbe79505f"
 GeometryBasics = "5c1252a2-5f33-56bf-86c9-59e7332b4326"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+SortTileRecursiveTree = "746ee33f-1797-42c2-866d-db2fce69d14d"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 Tables = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
 
 [weakdeps]
+FlexiJoins = "e37f2e79-19fa-4eb7-8510-b63b51fe0a37"
 Proj = "c94c279d-25a6-4763-9509-64d165bea63e"
 
 [extensions]
+GeometryOpsFlexiJoinsExt = "FlexiJoins"
 GeometryOpsProjExt = "Proj"
 
 [compat]
 CoordinateTransformations = "0.5, 0.6"
+FlexiJoins = "0.1.30"
 GeoInterface = "1.2"
 GeometryBasics = "0.4.7"
 LinearAlgebra = "1"
@@ -32,6 +36,7 @@ ArchGDAL = "c9ce4bd3-c3d5-55b8-8973-c0e20141b8c3"
 CoordinateTransformations = "150eb455-5306-5404-9cee-2592286d6298"
 DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
+FlexiJoins = "e37f2e79-19fa-4eb7-8510-b63b51fe0a37"
 GeoFormatTypes = "68eda718-8dee-11e9-39e7-89f7f65f511f"
 GeoJSON = "61d90e0f-e114-555e-ac52-39dfb47a3ef9"
 JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819"
@@ -42,4 +47,4 @@ Shapefile = "8e980c4a-a4fe-5da2-b3a7-4b4b0353a2f4"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["ArchGDAL", "CoordinateTransformations", "DataFrames", "Distributions", "GeoFormatTypes", "GeoJSON", "Proj", "JLD2", "LibGEOS", "Random", "Shapefile", "Test"]
+test = ["ArchGDAL", "CoordinateTransformations", "DataFrames", "Distributions", "FlexiJoins", "GeoFormatTypes", "GeoJSON", "Proj", "JLD2", "LibGEOS", "Random", "Shapefile", "Test"]

docs/Project.toml

@@ -10,6 +10,9 @@ DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 DocumenterVitepress = "4710194d-e776-4893-9690-8d956a29c365"
+DoubleFloats = "497a8b3b-efae-58df-a0af-a86822472b78"
+FlexiJoins = "e37f2e79-19fa-4eb7-8510-b63b51fe0a37"
+GADM = "a8dd9ffe-31dc-4cf5-a379-ea69100a8233"
 GeoDatasets = "ddc7317b-88db-5cb5-a849-8449e5df04f9"
 GeoInterface = "cf35fbd7-0cd7-5166-be24-54bfbe79505f"
 GeoInterfaceMakie = "0edc0954-3250-4c18-859d-ec71c1660c08"
@@ -20,6 +23,8 @@ LibGEOS = "a90b1aa1-3769-5649-ba7e-abc5a9d163eb"
 Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306"
 Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a"
 MakieThemes = "e296ed71-da82-5faf-88ab-0034a9761098"
+Measurements = "eff96d63-e80a-5855-80a2-b1b0885c5ab7"
+MonteCarloMeasurements = "0987c9cc-fe09-11e8-30f0-b96dd679fdca"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 Proj = "c94c279d-25a6-4763-9509-64d165bea63e"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

docs/make.jl

@@ -91,6 +91,9 @@ makedocs(;
     pages=[
         "Introduction" => "introduction.md",
         "API Reference" => "api.md",
+        "Tutorials" => [
+            "Spatial Joins" => "tutorials/spatial_joins.md",
+        ],
         "Source code" => literate_pages,
     ],
     warnonly = true,

docs/src/tutorials/spatial_joins.md

@@ -0,0 +1,111 @@
+# Spatial joins
+
+Spatial joins are joins which are based not on equality, but on some predicate ``p(x, y)``, which takes two geometries, and returns a value of either `true` or `false`.  For geometries, the [`DE-9IM`](https://en.wikipedia.org/wiki/DE-9IM) spatial relationship model is used to determine the spatial relationship between two geometries.  
+
+In this tutorial, we will show how to perform a spatial join on first a toy dataset and then two Natural Earth datasets, to show how this can be used in the real world.
+
+In order to perform the spatial join, we use [FlexiJoins.jl](https://github.com/JuliaAPlavin/FlexiJoins.jl) to perform the join, specifically using its `by_pred` joining method.  This allows the user to specify a predicate in the following manner:
+```julia
+[inner/left/right/outer/...]join((table1, table1),
+    by_pred(:table1_column, predicate_function, :table2_column) # & add other conditions here
+)
+```
+
+We have enabled the use of all of GeometryOps' boolean comparisons here.  These are:
+
+```julia
+GO.contains, GO.within, GO.intersects, GO.touches, GO.crosses, GO.disjoint, GO.overlaps, GO.covers, GO.coveredby, GO.equals
+```
+
+!!! tip
+    Always place the dataframe with more complex geometries second, as that is the one which will be sorted into a tree.
+
+## Simple example
+
+This example demonstrates how to perform a spatial join between two datasets: a set of polygons and a set of randomly generated points. 
+
+The polygons are represented as a DataFrame with geometries and colors, while the points are stored in a separate DataFrame. 
+
+The spatial join is performed using the `contains` predicate from GeometryOps, which checks if each point is contained within any of the polygons. The resulting joined DataFrame is then used to plot the points, colored according to the containing polygon.
+
+First, we generate our data.  We create two triangle polygons which, together, span the rectangle (0, 0, 1, 1), and a set of points which are randomly distributed within this rectangle.
+
+```@example spatialjoins
+import GeoInterface as GI, GeometryOps as GO
+using FlexiJoins, DataFrames
+
+using CairoMakie, GeoInterfaceMakie
+
+pl = GI.Polygon([GI.LinearRing([(0, 0), (1, 0), (1, 1), (0, 0)])])
+pu = GI.Polygon([GI.LinearRing([(0, 0), (0, 1), (1, 1), (0, 0)])])
+poly_df = DataFrame(geometry = [pl, pu], color = [:red, :blue])
+f, a, p = Makie.with_theme(Attributes(; Axis = (; aspect = DataAspect()))) do # hide
+f, a, p = poly(poly_df.geometry; color = tuple.(poly_df.color, 0.3))
+end # hide
+```
+
+Here, the upper polygon is blue, and the lower polygon is red.  Keep this in mind!
+
+Now, we generate the points.
+
+```@example spatialjoins
+points = tuple.(rand(1000), rand(1000))
+points_df = DataFrame(geometry = points)
+scatter!(points_df.geometry)
+f
+```
+
+You can see that they are evenly distributed around the box.  But how do we know which points are in which polygons?
+
+We have to join the two dataframes based on which polygon (if any) each point lies within.
+
+Now, we can perform the "spatial join" using FlexiJoins.  We are performing an outer join here
+
+```@example spatialjoins
+@time joined_df = FlexiJoins.innerjoin(
+    (points_df, poly_df), 
+    by_pred(:geometry, GO.within, :geometry)
+)
+```
+
+```@example spatialjoins
+scatter!(a, joined_df.geometry; color = joined_df.color)
+f
+```
+
+Here, you can see that the colors were assigned appropriately to the scattered points!
+
+## Real-world example
+
+Suppose I have a list of polygons representing administrative regions (or mining sites, or what have you), and I have a list of polygons for each country.  I want to find the country each region is in.
+
+```julia real
+import GeoInterface as GI, GeometryOps as GO
+using FlexiJoins, DataFrames, GADM # GADM gives us country and sublevel geometry
+
+using CairoMakie, GeoInterfaceMakie
+
+country_df = GADM.get.(["JPN", "USA", "IND", "DEU", "FRA"]) |> DataFrame
+country_df.geometry = GI.GeometryCollection.(GO.tuples.(country_df.geom))
+
+state_doublets = [
+    ("USA", "New York"),
+    ("USA", "California"),
+    ("IND", "Karnataka"),
+    ("DEU", "Berlin"),
+    ("FRA", "Grand Est"),
+    ("JPN", "Tokyo"),
+]
+
+state_full_df = (x -> GADM.get(x...)).(state_doublets) |> DataFrame
+state_full_df.geom = GO.tuples.(only.(state_full_df.geom))
+state_compact_df = state_full_df[:, [:geom, :NAME_1]]
+```
+
+```julia real
+innerjoin((state_compact_df, country_df), by_pred(:geom, GO.within, :geometry))
+innerjoin((state_compact_df,  view(country_df, 1:1, :)), by_pred(:geom, GO.within, :geometry))
+```
+
+!!! warning
+    This is how you would do this, but it doesn't work yet, since the GeometryOps predicates are quite slow on large polygons.  If you try this, the code will continue to run for a very, very long time (it took 12 hours on my laptop, but with minimal CPU usage).

ext/GeometryOpsFlexiJoinsExt/GeometryOpsFlexiJoinsExt.jl

@@ -0,0 +1,48 @@
+module GeometryOpsFlexiJoinsExt
+
+using GeometryOps
+using FlexiJoins
+
+import GeometryOps as GO, GeoInterface as GI
+using SortTileRecursiveTree, Tables
+
+
+# This module defines the FlexiJoins APIs for GeometryOps' boolean comparison functions, taken from DE-9IM.
+
+# First, we define the joining modes (Tree, NestedLoopFast) that the GO DE-9IM functions support.
+const GO_DE9IM_FUNCS = Union{typeof(GO.contains), typeof(GO.within), typeof(GO.intersects), typeof(GO.disjoint), typeof(GO.touches), typeof(GO.crosses), typeof(GO.overlaps), typeof(GO.covers), typeof(GO.coveredby), typeof(GO.equals)}
+# NestedLoopFast is the naive fallback method
+FlexiJoins.supports_mode(::FlexiJoins.Mode.NestedLoopFast, ::FlexiJoins.ByPred{F}, datas) where F <: GO_DE9IM_FUNCS = true
+# This method allows you to cache a tree, which we do by using an STRtree.
+# TODO: wrap GO predicate functions in a `TreeJoiner` struct or something, to indicate that we want to use trees,
+# since they can be slower in some situations.
+FlexiJoins.supports_mode(::FlexiJoins.Mode.Tree, ::FlexiJoins.ByPred{F}, datas) where F <: GO_DE9IM_FUNCS = true
+
+# Nested loop support is simple, and needs no further support.  
+# However, for trees, we need to define how the tree is prepared and how it is used.
+# This is done by defining the `prepare_for_join` function to return an STRTree,
+# and by defining the `findmatchix` function as querying that tree before checking
+# intersections.
+
+# In theory, one could extract the tree from e.g a GeoPackage or some future GeoDataFrame.
+
+FlexiJoins.prepare_for_join(::FlexiJoins.Mode.Tree, X, cond::FlexiJoins.ByPred{<: GO_DE9IM_FUNCS}) = (X, SortTileRecursiveTree.STRtree(map(cond.Rf, X)))
+function FlexiJoins.findmatchix(::FlexiJoins.Mode.Tree, cond::FlexiJoins.ByPred{F}, ix_a, a, (B, tree)::Tuple, multi::typeof(identity)) where F <: GO_DE9IM_FUNCS
+    idxs = SortTileRecursiveTree.query(tree, cond.Lf(a))
+    intersecting_idxs = filter!(idxs) do idx
+        cond.pred(a, cond.Rf(B[idx]))
+    end
+    return intersecting_idxs
+end
+
+# Finally, for completeness, we define the `swap_sides` function for those predicates which are defined as inversions.
+
+FlexiJoins.swap_sides(::typeof(GO.contains)) = GO.within
+FlexiJoins.swap_sides(::typeof(GO.within)) = GO.contains
+FlexiJoins.swap_sides(::typeof(GO.coveredby)) = GO.covers
+FlexiJoins.swap_sides(::typeof(GO.covers)) = GO.coveredby
+
+# That's a wrap, folks!
+
+end
+

test/extensions/flexijoins.jl

@@ -0,0 +1,22 @@
+import GeometryOps as GO, GeoInterface as GI
+using FlexiJoins, DataFrames
+
+pl = GI.Polygon([GI.LinearRing([(0, 0), (1, 0), (1, 1), (0, 0)])])
+pu = GI.Polygon([GI.LinearRing([(0, 0), (0, 1), (1, 1), (0, 0)])])
+poly_df = DataFrame(geometry = [pl, pu], color = [:red, :blue])
+
+points = tuple.(rand(100), rand(100))
+points_df = DataFrame(geometry = points)
+
+@testset "Basic integration" begin
+
+    @test_nowarn joined_df = FlexiJoins.innerjoin((poly_df, points_df), by_pred(:geometry, GO.contains, :geometry))
+    # Test that the join happened correctly
+    joined_df = FlexiJoins.innerjoin((poly_df, points_df), by_pred(:geometry, GO.contains, :geometry))
+    @test all(GO.contains.((pl,), joined_df.geometry_1[joined_df.color .== :red]))
+    @test all(GO.contains.((pu,), joined_df.geometry_1[joined_df.color .== :blue]))
+    # Test that within also works
+    @test_nowarn joined_df = FlexiJoins.innerjoin((points_df, poly_df), by_pred(:geometry, GO.within, :geometry))
+
+end
+

test/runtests.jl

@@ -41,4 +41,6 @@ const GO = GeometryOps
         include("transformations/correction/closed_ring.jl") 
         include("transformations/correction/intersecting_polygons.jl")
     end
+    # # Extensions
+    @testset "FlexiJoins" begin include("extensions/flexijoins.jl") end
 end