Further optimize extract() for regularly spaced arrays

[alex-s-gardner]

Manually calculating indices or using Interpolations is 3-5x faster than using extract()

Using extract()

using Rasters
using GeoInterface
dd = 0.001
x = X(0:dd:90);
y = Y(0:dd:90);
A = Raster(zeros(UInt8, y, x));

#number of points
n = 1e7;
f = GeoInterface.PointTuple.([(X=rand()*90, Y = rand()*90) for i in 1:n]);

@time pts = extract(A, f, atol=dd/2; index=true, geometry=false);
8.664601 seconds (50 allocations: 228.884 MiB, 0.93% gc time, 0.06% compilation time)

Using calculated index

@time begin
  r = round.(Int64,(pt_x .- x[1])./step(x)).+1;
  c = round.(Int64, (pt_y .- y[1]) ./ step(y)).+1;
  
  idx = [CartesianIndex(a,b) for (a,b) in  zip(r,c)];
  
  A[idx];
end;
1.690515 seconds (29.37 k allocations: 316.762 MiB, 11.33% gc time, 2.62% compilation time)

Using Interpolations

using Interpolations
@time begin    
    interp = Constant();
    itp = Interpolations.extrapolate(scale(Interpolations.interpolate(A.data, BSpline(interp)), A.dims[1].val.data, A.dims[2].val.data), NaN);
    val = itp.(pt_x, pt_y);
end;
2.974507 seconds (21 allocations: 694.414 MiB, 2.11% gc time)

[rafaqz]

On main? And why not Intervals??

Ok I'm getting similar timings. The last PR was just general performance fixes not the Regular case, seems that is still needed.

(also note round is incorrect here as it rounds towards zero, floor is what you want...)

[rafaqz]

I think really this is a DimensionalData.jl selector performance issue. We shouldn't have to think about Contains being slow on Regular here or write any custom code, it should just work and be as fast as just diving by step and calling floor.

[alex-s-gardner]

100% agree... I don't think Rasters should need a specialized extract function... all of that should be handled by DD

[rafaqz]

Ok I made a DD issue to track this and fix at the source.