Using `reproject` in a threaded context makes Julia crash

[mathieu17g]

This issue may be linked to #88

Using reproject in a threaded context makes Julia crash.

julia> versioninfo()
Julia Version 1.11.4
Commit 8561cc3d68 (2025-03-10 11:36 UTC)
Build Info:
  Official https://julialang.org/ release
Platform Info:
  OS: Windows (x86_64-w64-mingw32)
  CPU: 20 × 13th Gen Intel(R) Core(TM) i7-13700H
  WORD_SIZE: 64
  LLVM: libLLVM-16.0.6 (ORCJIT, goldmont)
Threads: 20 default, 0 interactive, 10 GC (on 20 virtual cores)
Environment:
  JULIA_NUM_THREADS = 20
  JULIA_PKG_USE_CLI_GIT = false
  JULIA_SSL_CA_ROOTS_PATH = 
  JULIA_EDITOR = code

(jl_tbk2nN) pkg> st
Status `C:\...\AppData\Local\Temp\jl_tbk2nN\Project.toml`
  [150eb455] CoordinateTransformations v0.6.4
  [e9467ef8] GLMakie v0.11.3
  [68eda718] GeoFormatTypes v0.4.4
  [cf35fbd7] GeoInterface v1.4.1
  [db073c08] GeoMakie v0.7.11
  [3251bfac] GeometryOps v0.1.16
  [c94c279d] Proj v1.8.1

Below is some code to reproduce the crash

1. Create a vector of MuliPolygon

import GeoInterface as GI
import CoordinateTransformations as CT
import GeometryOps as GO

# Lambert-93 (EPSG:2154) center (https://epsg.io/2154)
center_l93 = (489_353.59, 6_587_552.2)

# MultiPolygon pattern creation
exterior = GI.LinearRing(GI.Point.([(0, 0), (40_000, 0), (40_000, 40_000), (0, 40_000), (0, 0)]))
hole = GI.LinearRing(GI.Point.([(10_000, 10_000), (20_000, 10_000), (20_000, 20_000), (10_000, 20_000), (10_000, 10_000)]))
poly1 = GI.Polygon([exterior, hole])
poly2 = GO.transform(CT.Translation(50_000, 0), poly1) # Translate the polygon to the right by 50km
mpoly = GI.MultiPolygon([poly1, poly2])

# Create a 50x50 grid of MultiPolygons
N = 50
mpolys = [GO.transform(CT.Translation(((i - 1) % N) * 100_000, ((i - 1) ÷ N) * 50_000), mpoly) for i = 1:N^2]
# Translate the grid to the center of Lambert-93
mpolys = GO.transform.(Ref(CT.Translation(center_l93...)), mpolys)

2. Define CRS and PROJ String for reprojection

# EPSG:2154 (Lambert-93) to EPSG:3857 (Web Mercator)
import Proj
import GeoFormatTypes as GFT
src_crs = Proj.CRS("EPSG:2154")
src_proj_string = GFT.ProjString(src_crs).val
dst_crs = Proj.CRS("EPSG:3857")
dst_proj_string = GFT.ProjString(dst_crs).val;

3. Broadcasted GeometryOps.reproject works

# Define a plotting function to visualize the MultiPolygons
using GLMakie, GeoMakie
function show_mpolys(mpolys, src_proj_string, rmpolys, dst_proj_string)
    fig = Figure(; size = (1000, 300))
    ga1 = GeoAxis(fig[1, 1], source = src_proj_string; title = "Lambert-93 data source")
    ga2 = GeoAxis(fig[1, 2], source = dst_proj_string; title = "Data source converted to Web Mercator")
    for i = 1:N^2
        plot!(ga1, mpolys[i])
        plot!(ga2, rmpolys[i])
    end
    resize_to_layout!(fig)
    return fig
end

GLMakie.activate!()

rmpolys1 = GO.reproject.(mpolys, Ref(src_crs), Ref(dst_crs))
fig1 = show_mpolys(mpolys, src_proj_string, rmpolys1, dst_proj_string)
display(GLMakie.Screen(title="GeometryOps.jl reproject"), fig1)

4. Threaded GeometryOps.reproject on MultiPolygons leads to a Julia crash

rmpolys2 = fetch.([Threads.@spawn(GO.reproject(mpoly, src_proj_string, dst_proj_string)) for mpoly in mpolys])

---

Side issue, maybe worth a dedicated issue

While looking for a workaround to speed up my pipeline, I found that the performance of the broadcasted version of GeometryOps.reproject could be enhanced when compared to a broadcasted MultiPolygon conversion with a Proj transformation

using BenchmarkTools
reproject_mp(T::Proj.Transformation, mp::GI.MultiPolygon) =
    GI.MultiPolygon([[[[T(pt)...] for pt in pts] for pts in subgeom] for subgeom in GI.coordinates(mp)])
trans = Proj.Transformation(src_proj_string, dst_proj_string)

julia> @btime GO.reproject.($mpolys, Ref($src_crs), Ref($dst_crs));
  489.643 ms (7662 allocations: 334.42 KiB)

julia> @btime reproject_mp.(Ref($trans), $mpolys);
  2.640 ms (45460 allocations: 1.48 MiB)

[asinghvi17]

The thread crashing thing is likely because we need to assign Proj a different context for every Julia task that we spawn. Thanks for the detailed bug report!

That performance issue was unexpected - will look into it as well.

[asinghvi17]

So the issue with this benchmark is that you're broadcasting over the vector of multipolys, meaning that a new Proj transformation is created for each multipolygon, whereas in the manual method you only create it once outside of the function.

Apply will automatically broadcast over any vector or table, transform all geometries, then reconstruct - so you shouldn't have to broadcast over apply at all unless you have a vector of vectors of geometry, or something like that!

Changing this helps the runtime on my machine:
Manual method:

julia> @be reproject_mp.($((trans,)), $mpolys)
Benchmark: 6 samples with 1 evaluation
 min    17.645 ms (505010 allocs: 16.461 MiB)
 median 18.770 ms (505010 allocs: 16.461 MiB, 3.50% gc time)
 mean   18.606 ms (505010 allocs: 16.461 MiB, 3.88% gc time)
 max    19.448 ms (505010 allocs: 16.461 MiB, 9.23% gc time)

Apply with pre-constructed transform

julia> @be GO.reproject($mpolys, $trans)
Benchmark: 12 samples with 1 evaluation
 min    7.500 ms (75009 allocs: 3.109 MiB)
 median 7.573 ms (75009 allocs: 3.109 MiB)
 mean   8.385 ms (75009.08 allocs: 3.109 MiB, 3.13% gc time)
 max    10.729 ms (75010 allocs: 3.109 MiB, 20.59% gc time)

Apply with CRSs passed in

@be GO.reproject($mpolys, $src_crs, $dst_crs)
Benchmark: 11 samples with 1 evaluation
 min    8.724 ms (75011 allocs: 3.357 MiB)
 median 8.882 ms (75011 allocs: 3.357 MiB)
 mean   9.761 ms (75011.82 allocs: 3.357 MiB, 2.89% gc time)
 max    12.690 ms (75017 allocs: 3.357 MiB, 16.02% gc time)

Notice that if you only pass in the source and destination CRSs, it takes longer - this is because
Proj has to construct the transformation. Let's test this:

julia> @be Proj.Transformation($src_crs, $dst_crs; always_xy = true)
Benchmark: 79 samples with 1 evaluation
 min    1.154 ms (1 allocs: 32 bytes)
 median 1.194 ms (1 allocs: 32 bytes)
 mean   1.265 ms (1 allocs: 32 bytes)
 max    1.863 ms (1 allocs: 32 bytes)

Note the always_xy = true here: that forces the coordinate order to be (x, y) instead of (y, x). If you are creating Proj transforms yourself you will generally always need to have this, otherwise coordinate order will be flipped and GeometryOps won't know about it.

[asinghvi17]

You can also use multithreading that is built in to GeometryOps and all apply and applyreduce based functions, using the threaded = true or threaded = GO.True() (for maximum stability and speed) keyword:

julia> @be GO.reproject($mpolys, $trans; threaded = $true)
Benchmark: 57 samples with 1 evaluation
 min    1.507 ms (75275 allocs: 3.306 MiB)
 median 1.560 ms (75275 allocs: 3.306 MiB)
 mean   1.765 ms (75275 allocs: 3.306 MiB, 8.67% gc time)
 max    3.548 ms (75275 allocs: 3.306 MiB, 53.24% gc time)

that's quite the speedup from 7 ms! I would recommend, if you do this, that you pass the CRSs directly and not the transform - in future we will make this more threadsafe if you pass the CRS.

[asinghvi17]

In summary - Proj can be thread unsafe unless you construct one transformation per thread with unique Proj contexts. multithreading is available through apply and you should use that! And apply will do the broadcasting for you in addition to moving a lot of work out of the broadcast, so you don't have to broadcast when using it (usually :D).

[rafaqz]

I suspected this might just be really fast with threaded=true GG
Good to see it is!

But yes we should make sure it's always threadsafe.

[asinghvi17]

This is another case for an overridable apply pipeline - it would be great if reproject had access to the low-level maptasks call somehow...

[rafaqz]

Ah we need a dispatch object to pass through apply. Pretty doable

[rafaqz]

Just a thought... can we just make a channel at the top level, fill it with however many Transformations we need (nthreads) and then access them in the closure that runs in maptasks ?

Then it can all happen in reproject

Untested attempt:

function reproject(geom, transforms::Channel{Proj.Transformation}; time=Inf, target_crs=nothing, kw...)
     if _is3d(geom)
        return apply(GI.PointTrait(), geom; crs=target_crs, kw...) do p
            transform = take!(transforms)
            result = transform(GI.x(p), GI.y(p), GI.z(p))
            transform = put!(transforms, transform)
            result
        end
    else
        return apply(GI.PointTrait(), geom; crs=target_crs, kw...) do p
            transform = take!(transforms)
            result = transform(GI.x(p), GI.y(p))
            transform = put!(transforms, transform)
            result
        end
    end
end

And somewhere above add:

    if threaded
        channel = Channel{Proj.Transformation}()
        for _ in 1:Threads.nthreads() 
            put!(channel, Proj.Transformation(s, t; always_xy))
        end
        return reproject(geom, channel; time, target_crs, kw...)
    end

Possibly the work done by each apply is too small to make up for using the channel

[asinghvi17]

work done by each apply is too small to make up for using the channel

Exactly - channel put and take would be 99.999% of the runtime.

But that approach could work if we can somehow use task local values...again the cost is the repeated lookup per point, though. That's why I want to drill down to maptasks...

[rafaqz]

Yeah I didn't look closely enough at the target of apply...

So we need to insert this at the task partitioning level, we can just make a vector of functors in a wrapper

struct ThreadFunctors{F} 
    functors::F
    tasks_per_thread::Int
end

@inline function _maptasks(tf::ThreadFunctors, taskrange, threaded::True)::Vector
    ntasks = length(taskrange)
    chunk_size = max(1, ntasks ÷ (tf.tasks_per_thread * nthreads()))
    # partition the range into chunks
    task_chunks = Iterators.partition(taskrange, chunk_size)
    # Map over the chunks
    tasks = map(task_chunks, view(tf.functors, eachindex(task_chunks))) do chunk, f
        # Spawn a task to process this chunk
        StableTasks.@spawn begin
            # Where we map `f` over the chunk indices
            map(f, chunk)
        end
    end

    # Finally we join the results into a new vector
    return mapreduce(fetch, vcat, tasks)
end

Then we define them in reproject:

tasks_per_thread = 2
tf = ThreadFunctors([Proj.Transformation(s, t; always_xy) for _ in 1:Threads.nthreads() * tasks_per_thread], tasks_per_thread)

[asinghvi17]

That could work, let me prototype

[rafaqz]

On it

[asinghvi17]

Assuming you ran into similar issues I ran into