# Spatial join with dataframes

**URL:** <https://discourse.julialang.org/t/spatial-join-with-dataframes/93241>\
**Category:** General Usage\
**Tags:** dataframes, geo\
**Created:** [January 20, 2023, 1:07am UTC](https://discourse.julialang.org/t/spatial-join-with-dataframes/93241 "2023-01-20T01:07:03Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![jade\_mackay](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jade_mackay/32/6659_2.png) [@jade\_mackay](https://discourse.julialang.org/u/jade_mackay)\
**Post date:** [January 20, 2023, 1:07am UTC](https://discourse.julialang.org/t/spatial-join-with-dataframes/93241/1 "2023-01-20T01:07:03Z")

</div>

Hello,

Increasingly I find myself needing to join dataframes on Geographic Information System (GIS) objects and it seems no capability exists to do so in the Julia eco-system. Nor does it seem there is much interest (c.f. [Perform spatial join with GeoDataFrames](https://discourse.julialang.org/t/perform-spatial-join-with-geodataframes/86025)) . Using python I could do this

```plaintext
import geopandas as gpd
import pandas as pd

gdfa = gpd.read_file('./a.csv') # polygons
gdfb = gpd.read_file('./b.csv') # points
gdfc = gdfb.sjoin(gdfa, how="left", predicate="within")

```

Without any specific Julia tooling for this I have tried the following without success

```plaintext
using DataFrames
using ArchGDAL
using CSV

dsa = ArchGDAL.read("./a.csv")
layera = ArchGDAL.getlayer(dsa, 0)
dfa = layera |> DataFrame
rename!(dfa, ""=> :geometry)

dsb = ArchGDAL.read("./b.csv")
layerb = ArchGDAL.getlayer(dsb, 0)
dfb = layerb |> DataFrame
rename!(dfb, ""=> :geometry)

import Base.(==)
import Base.isequal
import Base.isless

function ==(a::ArchGDAL.IGeometry{ArchGDAL.wkbPoint},b::ArchGDAL.IGeometry{ArchGDAL.wkbPolygon})
    ArchGDAL.within(a,b)
end

function ==(a::ArchGDAL.IGeometry{ArchGDAL.wkbPolygon},b::ArchGDAL.IGeometry{ArchGDAL.wkbPoint})
    ArchGDAL.within(b,a)
end

innerjoin(dfa,dfb, on=:geometry, makeunique=true)

```

Thoughts and advice welcome!

The data used in the examples above are here:

> <https://gist.github.com/jademackay/d780816950be775ed7f252bc3010657f>

---

<div class="post-metadata">

**Author:** ![bkamins](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bkamins/32/208538_2.png) [@bkamins](https://discourse.julialang.org/u/bkamins)\
**Post date:** [January 20, 2023, 3:55am UTC](https://discourse.julialang.org/t/spatial-join-with-dataframes/93241/2 "2023-01-20T03:55:11Z")

</div>

I would assume that [[ANN] FlexiJoins.jl: fresh take on joining datasets](https://discourse.julialang.org/t/ann-flexijoins-jl-fresh-take-on-joining-datasets/79655) by @aplavin should be able to do it.

---

<div class="post-metadata">

**Author:** ![evetion](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/evetion/32/22679_2.png) [@evetion](https://discourse.julialang.org/u/evetion)\
**Post date:** [January 20, 2023, 9:46am UTC](https://discourse.julialang.org/t/spatial-join-with-dataframes/93241/3 "2023-01-20T09:46:57Z")

</div>

I fear it’s not yet doable within Julia as I know of no package that allows for completely generic predicates. I couldn’t get FlexiJoins predicates to work with custom methods, and I think both DataFrames as FlexiJoins eventually hash the values to do comparisons, which won’t work on geometries.

For now the workaround is to produce the crossjoin (product of the two tables) and filter them:

```julia
using GeoDataFrames
using DataFrames

dfa = GeoDataFrames.read("a.csv")
dfb = GeoDataFrames.read("b.csv")

X = crossjoin(dfa, dfb, makeunique=true)
subset!(X, [:geometry, :geometry_1] => (a, b) -> intersects.(a, b))
9×6 DataFrame
 Row │ geometry WKT a geometry_1 WKT_1 ⋯
     │ IGeometr… String String IGeometry… Strin ⋯
─────┼─────────────────────────────────────────────────────────────────────────────────────────────
   1 │ Geometry: wkbPolygon POLYGON ((1890386.16948909 55012… 1 Geometry: wkbPoint POINT ⋯
   2 │ Geometry: wkbPolygon POLYGON ((1796386.75606753 55606… 2 Geometry: wkbPoint POINT
   3 │ Geometry: wkbPolygon POLYGON ((1811431.72747276 56325… 3 Geometry: wkbPoint POINT
   4 │ Geometry: wkbPolygon POLYGON ((1776387.29179533 55767… 4 Geometry: wkbPoint POINT
   5 │ Geometry: wkbPolygon POLYGON ((1776318.37508329 55768… 5 Geometry: wkbPoint POINT ⋯
   6 │ Geometry: wkbPolygon POLYGON ((1818262.41701443 55449… 6 Geometry: wkbPoint POINT
   7 │ Geometry: wkbPolygon POLYGON ((1818162.05551015 55440… 7 Geometry: wkbPoint POINT
   8 │ Geometry: wkbPolygon POLYGON ((1818431.18365221 55443… 8 Geometry: wkbPoint POINT
   9 │ Geometry: wkbPolygon POLYGON ((1818162.8216453 554476… 9 Geometry: wkbPoint POINT ⋯

```

---

<div class="post-metadata">

**Author:** ![aplavin](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/aplavin/32/222056_2.png) [@aplavin](https://discourse.julialang.org/u/aplavin)\
**Post date:** [January 20, 2023, 10:24am UTC](https://discourse.julialang.org/t/spatial-join-with-dataframes/93241/4 "2023-01-20T10:24:25Z")

</div>

> [@bkamins](#):
>
> I would assume that [[ANN] FlexiJoins.jl: fresh take on joining datasets ](https://discourse.julialang.org/t/ann-flexijoins-jl-fresh-take-on-joining-datasets/79655) by @aplavin should be able to do it.

FlexiJoins support optimized spatial joins of points - as in “find pairs of points within the specified distance” or “find the closest match in B for each point from A”.  
But no optimization for polygon joins is available there, only the naive `O(n^2)` looping approach. This is potentially in scope for `FlexiJoins`, so feel free to propose an implementation that uses some optimized polygon lookup library. It shouldn’t be difficult, just isn’t in my own plans.

> [@evetion](#):
>
> FlexiJoins eventually hash the values to do comparisons, which won’t work on geometries

FlexiJoins support lots of join predicates that don’t do hashing. Optimized polygon joins just aren’t implemented because I never needed them.

---

<div class="post-metadata">

**Author:** ![rocco\_sprmnt21](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rocco_sprmnt21/32/20127_2.png) [@rocco\_sprmnt21](https://discourse.julialang.org/u/rocco_sprmnt21)\
**Post date:** [January 20, 2023, 10:55am UTC](https://discourse.julialang.org/t/spatial-join-with-dataframes/93241/5 "2023-01-20T10:55:56Z")

</div>

I don’t know if this will help

> [@Efficiently check if points are contained in polygons](https://discourse.julialang.org/t/efficiently-check-if-points-are-contained-in-polygons/74415/15):
>
> To clarify, Rasters.jl polygon ops are also pure julia via PolygonInbounds.jl and other internal code. It’s pretty fast for vectors of points, but not so much for single points. Its very fast for rasterize because the points of a raster are already sorted. (Rasters.jl only uses GDAL for loading some file types, because GDAL has so many and is reliable, and for warping and reprojecting, because projection formats are hard)

---

<div class="post-metadata">

**Author:** ![evetion](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/evetion/32/22679_2.png) [@evetion](https://discourse.julialang.org/u/evetion)\
**Post date:** [January 20, 2023, 11:41am UTC](https://discourse.julialang.org/t/spatial-join-with-dataframes/93241/6 "2023-01-20T11:41:55Z")

</div>

> [@aplavin](#):
>
> But no optimization for polygon joins is available there, only the naive `O(n^2)` looping approach.

How would one execute the naive version with the above example? Can it make use of `intersect`/`within` methods of other packages?

> [@aplavin](#):
>
> FlexiJoins support lots of join predicates that don’t do hashing. Optimized polygon joins just aren’t implemented because I never needed them.

That’s only fair, and I wouldn’t expect you to implement spatial geometry operations. There are many types of geometries and spatial predicates, some of which are implemented in other packages. What could be useful is an interface definition for making optimized lookups for joins, so packages can implement their own algorithms for their specific types.

For example for simple lookups in vectors (no joins), there is [GitHub - andyferris/AcceleratedArrays.jl: Arrays with acceleration indices · GitHub](https://github.com/andyferris/AcceleratedArrays.jl), which speeds up things like `findall` (discussion about it here [DataFrames.jl/issues/2381](https://github.com/JuliaData/DataFrames.jl/issues/2381). Similarly, I’ve made [GitHub - evetion/GeoAcceleratedArrays.jl: AcceleratedArrays with spatial indexing · GitHub](https://github.com/evetion/GeoAcceleratedArrays.jl) in the past, which does the same with a spatial index.

---

<div class="post-metadata">

**Author:** ![evetion](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/evetion/32/22679_2.png) [@evetion](https://discourse.julialang.org/u/evetion)\
**Post date:** [January 20, 2023, 11:47am UTC](https://discourse.julialang.org/t/spatial-join-with-dataframes/93241/7 "2023-01-20T11:47:42Z")

</div>

> [@rocco\_sprmnt21](#):
>
> I don’t know if this will help

PolygonInbounds.jl is a good example of a specific subset of geospatial predicates, written natively in Julia Overall there are more ([DE-9IM - Wikipedia](https://en.wikipedia.org/wiki/DE-9IM#Spatial_predicates)) and these are covered by packages like LibGEOS and ArchGDAL (itself linking to GEOS via GDAL) for the package geometry types. [GitHub - JuliaGeo/GeoInterface.jl: A Julia Protocol for Geospatial Data](https://github.com/JuliaGeo/GeoInterface.jl/) provides an interface for such operations.

Note that GEOS has a thing like _prepared_ geometries, which make repeated spatial operations much faster on them.

---

<div class="post-metadata">

**Author:** ![aplavin](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/aplavin/32/222056_2.png) [@aplavin](https://discourse.julialang.org/u/aplavin)\
**Post date:** [January 20, 2023, 3:15pm UTC](https://discourse.julialang.org/t/spatial-join-with-dataframes/93241/8 "2023-01-20T15:15:47Z")

</div>

> [@evetion](#):
>
> How would one execute the naive version with the above example? Can it make use of `intersect`/`within` methods of other packages?

FlexiJoins can do joins with arbitrary predicates by just looping over all pairs:

```julia
julia> using FlexiJoins

julia> my_predicate(s, r) = occursin(r, s)

julia> ss = ["abc", "def", "xyz"]
julia> rs = [r"a|b", r".{3}"]

julia> innerjoin((;ss, rs), by_pred(identity, my_predicate, identity), mode=FlexiJoins.Mode.NestedLoop(), loop_over_side=1)
4-element StructArray(view(::Vector{String}, [1, 1, 2, 3]), view(::Vector{Regex}, [1, 2, 2, 2])) with eltype NamedTuple{(:ss, :rs), Tuple{String, Regex}}:
 (ss = "abc", rs = r"a|b")
 (ss = "abc", rs = r".{3}")
 (ss = "def", rs = r".{3}")
 (ss = "xyz", rs = r".{3}")

```

It’s just a pretty rare [to me] operation, so maybe requires more boilerplate than strictly needed - eg, `loop_over_side` can be chosen arbitrarily for nested loop joins anyway.

> [@evetion](#):
>
> That’s only fair, and I wouldn’t expect you to implement spatial geometry operations. There are many types of geometries and spatial predicates, some of which are implemented in other packages. What could be useful is an interface definition for making optimized lookups for joins, so packages can implement their own algorithms for their specific types.

That’s exactly how FlexiJoins is designed. As a similar example, I don’t implement spatial point searches myself, but just utilize `NearestNeighbors.jl` for distance queries. Basic support for that is just the first few lines in [src/nearestneighbors.jl · master · Alexander Plavin / FlexiJoins.jl · GitLab](https://gitlab.com/aplavin/FlexiJoins.jl/-/blob/master/src/nearestneighbors.jl#L2-5).

This interface isn’t advertised for now, and maybe not the most optimal one - but it works.

> [@evetion](#):
>
> For example for simple lookups in vectors (no joins), there is [GitHub - andyferris/AcceleratedArrays.jl: Arrays with acceleration indices](https://github.com/andyferris/AcceleratedArrays.jl), which speeds up things like `findall` (discussion about it here [DataFrames.jl/issues/2381](https://github.com/JuliaData/DataFrames.jl/issues/2381). Similarly, I’ve made [GitHub - evetion/GeoAcceleratedArrays.jl: AcceleratedArrays with spatial indexing](https://github.com/evetion/GeoAcceleratedArrays.jl) in the past, which does the same with a spatial index.

Yeah, I’m aware of `AcceleratedArrays`, but it wasn’t really clear on how to utilize them in `FlexiJoins`. Maybe it’s possible though!

---

<div class="post-metadata">

**Author:** ![rocco\_sprmnt21](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rocco_sprmnt21/32/20127_2.png) [@rocco\_sprmnt21](https://discourse.julialang.org/u/rocco_sprmnt21)\
**Post date:** [January 20, 2023, 5:01pm UTC](https://discourse.julialang.org/t/spatial-join-with-dataframes/93241/9 "2023-01-20T17:01:55Z")

</div>

> [@evetion](#):
>
> PolygonInbounds.jl is a good …

I’m not sure I understand the remarks in response to my question.  
Surely also because I misquoted the reference to the discussion “Efficiently check if points are contained in polygons”.  
What I asked to have clarified is if the use of one of the packages mentioned in the post (GMT.jl, rasters.jl, etc.) would not respond to Pavlin’s request

> [@aplavin](#):
>
> so feel free to propose an implementation that uses some optimized polygon lookup library. It shouldn’t be difficult, just isn’t in my own plans.

That is if using FlexiJoins with (for example) GMT doesn’t solve your problem.

---

<div class="post-metadata">

**Author:** ![Dan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dan/32/42581_2.png) [@Dan](https://discourse.julialang.org/u/Dan)\
**Post date:** [January 20, 2023, 5:29pm UTC](https://discourse.julialang.org/t/spatial-join-with-dataframes/93241/10 "2023-01-20T17:29:11Z")

</div>

Here is an attempt to use FlexiJoins by first finding a bounding box for each geometry. This reduces the quadratic complexity of the query (depending on how annoying the geometries are):

```julia
# prep data & environment
using DataFrames, GeoDataFrames, ArchGDAL, FlexiJoins, IntervalSets

dfa = GeoDataFrames.read("a.csv")
dfb = GeoDataFrames.read("b.csv")

getbb(g) = begin
    gg = ArchGDAL.boundingbox(g)
    gg1 = ArchGDAL.getgeom(gg,0)
    x1, y2 = ArchGDAL.getpoint(gg1,0)
    x2, y1 = ArchGDAL.getpoint(gg1,2)
    (x1..x2, y1..y2)
end

getxy(p) = begin
    x,y = ArchGDAL.getpoint(p,0)
    (x, y)
end

```

```julia
# process data
dfa2 = transform(dfa, 
  :geometry => ByRow(getbb) => [:xinterval, :yinterval])
dfb2 = transform(dfb, 
  :geometry => ByRow(getxy) => [:x, :y])

subset!(
  innerjoin(
    (innerjoin((dfb2,dfa2), by_pred(:x, ∈, :xinterval)),dfa2), 
    by_pred(:y, ∈, :yinterval)
  ), [:geometry, :geometry_1] => (a, b) -> intersects.(a, b))

```

The results are similar to the ones from a crossjoin. The bounding boxes and intervals may be of further use possibly.

---

<div class="post-metadata">

**Author:** ![aplavin](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/aplavin/32/222056_2.png) [@aplavin](https://discourse.julialang.org/u/aplavin)\
**Post date:** [January 20, 2023, 11:53pm UTC](https://discourse.julialang.org/t/spatial-join-with-dataframes/93241/11 "2023-01-20T23:53:13Z")

</div>

@Dan nice solution! Should definitely outperform the naive O(n^2) one for large datasets, even with the overhead of two sequential joins.  
Btw, materializing x and y coordinates into separate columns isn’t required: you can pass `r -> first(r.geometry)` instead of `:x` in `by_pred`.

The best interface would be `innerjoin((L, R), by_pred(:geometry, ∈, :geometry))`, and this is potentially possible with FlexiJoins - if someone feels like introducing an optimized implementation for that.

---

<div class="post-metadata">

**Author:** ![jade\_mackay](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jade_mackay/32/6659_2.png) [@jade\_mackay](https://discourse.julialang.org/u/jade_mackay)\
**Post date:** [January 21, 2023, 9:26pm UTC](https://discourse.julialang.org/t/spatial-join-with-dataframes/93241/12 "2023-01-21T21:26:29Z")

</div>

Thanks everyone for engaging in this thread, it is exemplary of a great community. I am especially grateful to @evetion’s workaround and to have be made aware of @aplavin’s FlexiJoin package which I expect to make use of in many other contexts. I also appreciate @Dan’s clever application of FlexiJoins’s interval predicates to solve the problem.

With regard to the particulars that motivated the post, due to the size of the actual problem (3 million polygons and 1 million points), on my computer(32GB mem) the workaround fails out of hand at the crossjoin with `ERROR: OutOfMemoryError()`. And the Flexijoin-with-bounding-boxed-polygons steadily increases memory consumption until killed by the OS memory manager. An optimized FlexiJoin implementation with interface like `innerjoin((L, R), by_pred(:geometry, ∈, :geometry))`, looks a promising way forward.

---

<div class="post-metadata">

**Author:** ![aplavin](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/aplavin/32/222056_2.png) [@aplavin](https://discourse.julialang.org/u/aplavin)\
**Post date:** [January 21, 2023, 10:00pm UTC](https://discourse.julialang.org/t/spatial-join-with-dataframes/93241/13 "2023-01-21T22:00:27Z")

</div>

> [@jade\_mackay](#):
>
> the Flexijoin-with-bounding-boxed-polygons steadily increases memory consumption until killed by the OS memory manager

Joining 1 million numbers with 10 million intervals works just fine on my laptop, and doesn’t noticeably increase memory consumption:

```julia
julia> xs = rand(10^6);
julia> ys = rand(10^7);
julia> ints = Interval.(ys, ys .+ 1e-7);

julia> @time innerjoin((xs, ints), by_pred(identity, ∈, identity)) |> length
 25.102649 seconds (84 allocations: 27.195 MiB)
1001054

```

@Dan’s trick only uses FlexiJoins for this kind of interval joins, so should perform similarly.  
If you can provide a simple reproducer without fancy types (no ArchGDAL, GeoDataFrames, DataFrames - only arrays, numbers, and intervals), I’ll try to look at it to see what goes wrong.

For example, is it possible that just too many point-bbox pairs match, and the resulting array doesn’t fit into memory?

---

<div class="post-metadata">

**Author:** ![jade\_mackay](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jade_mackay/32/6659_2.png) [@jade\_mackay](https://discourse.julialang.org/u/jade_mackay)\
**Post date:** [January 23, 2023, 6:52am UTC](https://discourse.julialang.org/t/spatial-join-with-dataframes/93241/14 "2023-01-23T06:52:40Z")

</div>

> For example, is it possible that just too many point-bbox pairs match, and the resulting array doesn’t fit into memory?

Yes, I think that’s the situation. The result of the geopandas join has size c.a. 60,000 greater than the 3 million starting. That is, about a 1% increase, which I think situation basically corresponds to

```julia
xs = rand(10^6);
ys = rand(10^6);
nts = Interval.(ys, ys .+ 1e-2);
@time innerjoin((xs, ints), by_pred(identity, ∈, identity)) |> length

```

and which, it seems, doesn’t fit in memory.

---

<div class="post-metadata">

**Author:** ![aplavin](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/aplavin/32/222056_2.png) [@aplavin](https://discourse.julialang.org/u/aplavin)\
**Post date:** [January 23, 2023, 11:53am UTC](https://discourse.julialang.org/t/spatial-join-with-dataframes/93241/15 "2023-01-23T11:53:06Z")

</div>

> [@jade\_mackay](#):
>
> ```julia
> xs = rand(10^6);
> ys = rand(10^6);
> nts = Interval.(ys, ys .+ 1e-2);
> 
> ```

This results in about `10^6 * 10^6 * 10^-2 = 10^10` or 10 billion matches. No surprise, they don’t readily fit into memory!

---

<div class="post-metadata">

**Author:** ![Dan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dan/32/42581_2.png) [@Dan](https://discourse.julialang.org/u/Dan)\
**Post date:** [January 23, 2023, 12:07pm UTC](https://discourse.julialang.org/t/spatial-join-with-dataframes/93241/16 "2023-01-23T12:07:17Z")

</div>

When using the BBox to filter the points, the FlexiJoins are nested. So if there are too many records in the first join, we can:

1. switch the order of X and Y join - in case one axis is a better filter than the other and should go first.

2. chunk the points into batches, and perform the double filter on chunks - after two joins there may be fewer points left which would fit memory. When we get to chunks of size 1, this is like a big `for` loop on the points, and if result still doesn’t fit in memory…

3. Instead of materializing any more big DataFrames, just indices to point and polygon frames can be stored. This should drop the memory requirements to about two Ints per candidate intersection.

---

<div class="post-metadata">

**Author:** ![aplavin](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/aplavin/32/222056_2.png) [@aplavin](https://discourse.julialang.org/u/aplavin)\
**Post date:** [January 23, 2023, 12:40pm UTC](https://discourse.julialang.org/t/spatial-join-with-dataframes/93241/17 "2023-01-23T12:40:11Z")

</div>

> [@Dan](#):
>
> 1. Instead of materializing any more big DataFrames, just indices to point and polygon frames can be stored. This should drop the memory requirements to about two Ints per candidate intersection.

That’s how FlexiJoins work already: it computes indices, and the output is a view of the input.  
DataFrames are special though - they don’t follow the collection interface, and require special handling [here](https://gitlab.com/aplavin/FlexiJoins.jl/-/blob/master/src/FlexiJoins.jl#L40-57). This integration is less tested and maybe makes some unnecessary copies?..

---

<div class="post-metadata">

**Author:** ![evetion](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/evetion/32/22679_2.png) [@evetion](https://discourse.julialang.org/u/evetion)\
**Post date:** [January 25, 2023, 8:16pm UTC](https://discourse.julialang.org/t/spatial-join-with-dataframes/93241/18 "2023-01-25T20:16:11Z")

</div>

> [@aplavin](#):
>
> FlexiJoins can do joins with arbitrary predicates by just looping over all pairs [..]  
> It’s just a pretty rare [to me] operation, so maybe requires more boilerplate than strictly needed - eg, `loop_over_side` can be chosen arbitrarily for nested loop joins anyway.

Great stuff! I didn’t know about these keyword arguments, so the MWE would now be:

```julia
using GeoDataFrames
using FlexiJoins

dfa = GeoDataFrames.read("a.csv")
dfb = GeoDataFrames.read("b.csv")

# Either a custom predicate that takes two rows of the DataFrame
my_predicate(s, r) = GeoDataFrames.within(r.geometry, s.geometry) # reversed
innerjoin((dfa, dfb), by_pred(identity, my_predicate, identity), mode=FlexiJoins.Mode.NestedLoop(), loop_over_side=1)

# Or a version that works directly on the fields
innerjoin((dfa, dfb), by_pred(:geometry, GeoDataFrames.contains, :geometry), mode=FlexiJoins.Mode.NestedLoop(), loop_over_side=1)

```

Is the `loop_over_side=1` documented somewhere? Without it you get an error that there’s no method `swap_sides(f)`.

> [@rocco\_sprmnt21](#):
>
> I’m not sure I understand the remarks in response to my question.  
> Surely also because I misquoted the reference to the discussion “Efficiently check if points are contained in polygons”.  
> What I asked to have clarified is if the use of one of the packages mentioned in the post (GMT.jl, rasters.jl, etc.) would not respond to Pavlin’s request

Ah, I indeed didn’t fully understand your question. I took it as a suggestion for spatial intersections/geometry, and I suggested more. I think these can’t do this spatial join themselves, but they can be a provider for the predicate function in a naive, non-optimized (n^2) way. The only fast way is doing spatial indexing, which is essentially an optimized version of the boundingbox approach by [Dan](https://discourse.julialang.org/t/spatial-join-with-dataframes/93241/10) later in this topic.

---

<div class="post-metadata">

**Author:** ![evetion](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/evetion/32/22679_2.png) [@evetion](https://discourse.julialang.org/u/evetion)\
**Post date:** [January 25, 2023, 8:22pm UTC](https://discourse.julialang.org/t/spatial-join-with-dataframes/93241/19 "2023-01-25T20:22:14Z")

</div>

> [@aplavin](#):
>
> Yeah, I’m aware of `AcceleratedArrays`, but it wasn’t really clear on how to utilize them in `FlexiJoins`. Maybe it’s possible though!

That would be great! I guess it could be another Mode? Instead of the predicate getting `element \* element in the Nested Mode, it would require a mode that does element \* vector. The vector being the AcceleratedArray column in a table.

---

<div class="post-metadata">

**Author:** ![jar1](https://avatars.discourse-cdn.com/v4/letter/j/c0e974/32.png) [@jar1](https://discourse.julialang.org/u/jar1)\
**Post date:** [January 25, 2023, 8:31pm UTC](https://discourse.julialang.org/t/spatial-join-with-dataframes/93241/20 "2023-01-25T20:31:21Z")

</div>

There is a literature on [spatial join techniques](http://www.cs.umd.edu/~hjs/pubs/jacoxtrjoin07.pdf). I’m not sure how far along Julia’s packages have gotten.

[Next page](https://discourse.julialang.org/t/spatial-join-with-dataframes/93241.md?page=2)