I have two arrays of polygons, one with 9,000,000 polygons and the other with 60,000 polygons. I would like to find the intersecting area of each of these. In the MWE example below where I have reduced the example to 201 polygons x 201 polygons, I can only achieve this in 35.522 ms. I have tried using `Threads.@threads`

in the outer loop of `poly_int`

, but this does not work with either sparse arrays (which I need to form a 9,000,000 by 60,000 array) or LibGEOS. Does anyone have any suggestions for improving the performance (e.g., should I focus on finding a way for Threads to work)?

Note: I deleted a previous post about this, as I believe I have asked the question more directly now.

```
using LibGEOS, SparseArrays, BenchmarkTools
function poly_int(grid_1,grid_2)
polygon_intersection = spzeros(length(grid_1["polygon"]),length(grid_2["polygon"]));
for j = 1:length(grid_2["polygon"]) # accessing in column order
for i = 1:length(grid_1["polygon"])
polygon_intersection[i,j] = area(intersection(grid_1["polygon"][i],grid_2["polygon"][j]));
end
end
return polygon_intersection
end
@btime poly_int(grid_1,grid_2);
```

The full code for the MWE example (that defines grid_1 and grid_2) is here.

```
using LibGEOS, SparseArrays, BenchmarkTools
grid_1 = Dict();
# c = center, ll = lower left corner, ul = upper left corner
# ur = upper right corner, lr = lower right corner
grid_1["x_c"] = collect(0:0.5:100);
grid_1["y_c"] = collect(0:0.5:100);
grid_1["x_ll"] = grid_1["x_c"] .- (0.5/2);
grid_1["x_ul"] = grid_1["x_c"] .- (0.5/2);
grid_1["x_ur"] = grid_1["x_c"] .+ (0.5/2);
grid_1["x_lr"] = grid_1["x_c"] .+ (0.5/2);
grid_1["y_ll"] = grid_1["y_c"] .- (0.5/2);
grid_1["y_ul"] = grid_1["y_c"] .+ (0.5/2);
grid_1["y_ur"] = grid_1["y_c"] .+ (0.5/2);
grid_1["y_lr"] = grid_1["y_c"] .- (0.5/2);
grid_1["value"] = rand(length(grid_1["x_c"]));
grid_1["uncertainty"] = rand(length(grid_1["x_c"]));
grid_2 = Dict();
grid_2["x_c"] = collect(1:0.5:101);
grid_2["y_c"] = collect(1:0.5:101);
grid_2["x_ll"] = grid_2["x_c"] .- (0.5/2);
grid_2["x_ul"] = grid_2["x_c"] .- (0.5/2);
grid_2["x_ur"] = grid_2["x_c"] .+ (0.5/2);
grid_2["x_lr"] = grid_2["x_c"] .+ (0.5/2);
grid_2["y_ll"] = grid_2["y_c"] .- (0.5/2);
grid_2["y_ul"] = grid_2["y_c"] .+ (0.5/2);
grid_2["y_ur"] = grid_2["y_c"] .+ (0.5/2);
grid_2["y_lr"] = grid_2["y_c"] .- (0.5/2);
grid_2["value"] = rand(length(grid_2["x_c"]));
grid_2["uncertainty"] = rand(length(grid_2["x_c"]));
for grid in (grid_1,grid_2)
grid["polygon"] = Vector{Polygon}(undef,length(grid["x_c"]));
# makes polygons from the pixcorners for both grids
for idx = 1:length(grid["polygon"])
grid["polygon"][idx] = Polygon([[[grid["x_ll"][idx],grid["y_ll"][idx]],[grid["x_ul"][idx],grid["y_ul"][idx]],[grid["x_ur"][idx],grid["y_ur"][idx]],[grid["x_lr"][idx],grid["y_lr"][idx]],[grid["x_ll"][idx],grid["y_ll"][idx]]]]);
end
grid["polygon_area"] = area.(grid["polygon"]);
end
function poly_int(grid_1,grid_2)
polygon_intersection = spzeros(length(grid_1["polygon"]),length(grid_2["polygon"]));
for j = 1:length(grid_2["polygon"]) # accessing in column order
for i = 1:length(grid_1["polygon"])
polygon_intersection[i,j] = area(intersection(grid_1["polygon"][i],grid_2["polygon"][j]));
end
end
return polygon_intersection
end
@btime poly_int(grid_1,grid_2);
```