|
| 1 | +# Spatial joins |
| 2 | + |
| 3 | +Spatial joins are [table joins](https://www.geeksforgeeks.org/sql-join-set-1-inner-left-right-and-full-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. |
| 4 | + |
| 5 | +Spatial joins can be done between any geometry types (from geometrycollections to points), just as geometrical predicates can be evaluated on any geometries. |
| 6 | + |
| 7 | +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. |
| 8 | + |
| 9 | +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: |
| 10 | +```julia |
| 11 | +[inner/left/right/outer/...]join((table1, table1), |
| 12 | + by_pred(:table1_column, predicate_function, :table2_column) # & add other conditions here |
| 13 | +) |
| 14 | +``` |
| 15 | + |
| 16 | +We have enabled the use of all of GeometryOps' boolean comparisons here. These are: |
| 17 | + |
| 18 | +```julia |
| 19 | +GO.contains, GO.within, GO.intersects, GO.touches, GO.crosses, GO.disjoint, GO.overlaps, GO.covers, GO.coveredby, GO.equals |
| 20 | +``` |
| 21 | + |
| 22 | +!!! tip |
| 23 | + Always place the dataframe with more complex geometries second, as that is the one which will be sorted into a tree. |
| 24 | + |
| 25 | +## Simple example |
| 26 | + |
| 27 | +This example demonstrates how to perform a spatial join between two datasets: a set of polygons and a set of randomly generated points. |
| 28 | + |
| 29 | +The polygons are represented as a DataFrame with geometries and colors, while the points are stored in a separate DataFrame. |
| 30 | + |
| 31 | +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. |
| 32 | + |
| 33 | +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. |
| 34 | + |
| 35 | +```@example spatialjoins |
| 36 | +import GeoInterface as GI, GeometryOps as GO |
| 37 | +using FlexiJoins, DataFrames |
| 38 | +
|
| 39 | +using CairoMakie, GeoInterfaceMakie |
| 40 | +
|
| 41 | +pl = GI.Polygon([GI.LinearRing([(0, 0), (1, 0), (1, 1), (0, 0)])]) |
| 42 | +pu = GI.Polygon([GI.LinearRing([(0, 0), (0, 1), (1, 1), (0, 0)])]) |
| 43 | +poly_df = DataFrame(geometry = [pl, pu], color = [:red, :blue]) |
| 44 | +f, a, p = Makie.with_theme(Attributes(; Axis = (; aspect = DataAspect()))) do # hide |
| 45 | +f, a, p = poly(poly_df.geometry; color = tuple.(poly_df.color, 0.3)) |
| 46 | +end # hide |
| 47 | +``` |
| 48 | + |
| 49 | +Here, the upper polygon is blue, and the lower polygon is red. Keep this in mind! |
| 50 | + |
| 51 | +Now, we generate the points. |
| 52 | + |
| 53 | +```@example spatialjoins |
| 54 | +points = tuple.(rand(1000), rand(1000)) |
| 55 | +points_df = DataFrame(geometry = points) |
| 56 | +scatter!(points_df.geometry) |
| 57 | +f |
| 58 | +``` |
| 59 | + |
| 60 | +You can see that they are evenly distributed around the box. But how do we know which points are in which polygons? |
| 61 | + |
| 62 | +We have to join the two dataframes based on which polygon (if any) each point lies within. |
| 63 | + |
| 64 | +Now, we can perform the "spatial join" using FlexiJoins. We are performing an outer join here |
| 65 | + |
| 66 | +```@example spatialjoins |
| 67 | +@time joined_df = FlexiJoins.innerjoin( |
| 68 | + (points_df, poly_df), |
| 69 | + by_pred(:geometry, GO.within, :geometry) |
| 70 | +) |
| 71 | +``` |
| 72 | + |
| 73 | +```@example spatialjoins |
| 74 | +scatter!(a, joined_df.geometry; color = joined_df.color) |
| 75 | +f |
| 76 | +``` |
| 77 | + |
| 78 | +Here, you can see that the colors were assigned appropriately to the scattered points! |
| 79 | + |
| 80 | +## Real-world example |
| 81 | + |
| 82 | +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. |
| 83 | + |
| 84 | +```julia real |
| 85 | +import GeoInterface as GI, GeometryOps as GO |
| 86 | +using FlexiJoins, DataFrames, GADM # GADM gives us country and sublevel geometry |
| 87 | + |
| 88 | +using CairoMakie, GeoInterfaceMakie |
| 89 | + |
| 90 | +country_df = GADM.get.(["JPN", "USA", "IND", "DEU", "FRA"]) |> DataFrame |
| 91 | +country_df.geometry = GI.GeometryCollection.(GO.tuples.(country_df.geom)) |
| 92 | + |
| 93 | +state_doublets = [ |
| 94 | + ("USA", "New York"), |
| 95 | + ("USA", "California"), |
| 96 | + ("IND", "Karnataka"), |
| 97 | + ("DEU", "Berlin"), |
| 98 | + ("FRA", "Grand Est"), |
| 99 | + ("JPN", "Tokyo"), |
| 100 | +] |
| 101 | + |
| 102 | +state_full_df = (x -> GADM.get(x...)).(state_doublets) |> DataFrame |
| 103 | +state_full_df.geom = GO.tuples.(only.(state_full_df.geom)) |
| 104 | +state_compact_df = state_full_df[:, [:geom, :NAME_1]] |
| 105 | +``` |
| 106 | + |
| 107 | +```julia real |
| 108 | +innerjoin((state_compact_df, country_df), by_pred(:geom, GO.within, :geometry)) |
| 109 | +innerjoin((state_compact_df, view(country_df, 1:1, :)), by_pred(:geom, GO.within, :geometry)) |
| 110 | +``` |
| 111 | + |
| 112 | +!!! warning |
| 113 | + 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). |
| 114 | + |
| 115 | +## Enabling custom predicates |
| 116 | + |
| 117 | +In case you want to use a custom predicate, you only need to define a method to tell FlexiJoins how to use it. |
| 118 | + |
| 119 | +For example, let's suppose you wanted to perform a spatial join on geometries which are some distance away from each other: |
| 120 | + |
| 121 | +```julia |
| 122 | +my_predicate_function = <(5) ∘ abs ∘ GO.distance |
| 123 | +``` |
| 124 | + |
| 125 | +You would need to define `FlexiJoins.supports_mode` on your predicate: |
| 126 | + |
| 127 | +```julia{3} |
| 128 | +FlexiJoins.supports_mode( |
| 129 | + ::FlexiJoins.Mode.NestedLoopFast, |
| 130 | + ::FlexiJoins.ByPred{typeof(my_predicate_function)}, |
| 131 | + datas |
| 132 | +) = true |
| 133 | +``` |
| 134 | + |
| 135 | +This will enable FlexiJoins to support your custom function, when it's passed to `by_pred(:geometry, my_predicate_function, :geometry)`. |
0 commit comments