diff --git a/src/lib/overloads.jl b/src/lib/overloads.jl index 0cb06b58..5abea4d5 100644 --- a/src/lib/overloads.jl +++ b/src/lib/overloads.jl @@ -21,10 +21,9 @@ Shows a textual representation of the layer. function Base.show(io::IO, layer::T) where {T <: SimpleSDMLayer} itype = eltype(layer) otype = T <: SimpleSDMPredictor ? "predictor" : "response" - print(io, """SDM $(otype) with $(itype) values - $(size(layer,1)) × $(size(layer,2)) - lat.: $(extrema(latitudes(layer))) - lon.: $(extrema(longitudes(layer)))""") + print(io, """SDM $(otype) → $(size(layer,1))×$(size(layer,2)) grid with $(length(layer)) $(itype)-valued cells + Latitudes\t$(extrema(latitudes(layer))) + Longitudes\t$(extrema(longitudes(layer)))""") end """ diff --git a/src/operations/sliding.jl b/src/operations/sliding.jl index e359e2f4..6a1f03ea 100644 --- a/src/operations/sliding.jl +++ b/src/operations/sliding.jl @@ -1,6 +1,4 @@ -function haversine(p1, p2; R=6371.0) - lon1, lat1 = p1 - lon2, lat2 = p2 +function haversine(lon1, lat1, lon2, lat2; R=6371.0) φ1 = lat1 * π/180.0 φ2 = lat2 * π/180.0 Δφ = (lat2-lat1) * π/180.0 @@ -10,6 +8,8 @@ function haversine(p1, p2; R=6371.0) return R*c end +haversine(p1, p2; R=6371.0) = haversine(p1..., p2...; R=R) + """ slidingwindow(L::LT, f::FT, d::IT) where {LT <: SimpleSDMLayer, FT <: Function, IT <: Number} @@ -30,27 +30,24 @@ This function is currently relatively slow. Performance improvements will arrive at some point. """ function slidingwindow(layer::LT, f::FT, d::IT) where {LT <: SimpleSDMLayer, FT <: Function, IT <: Number} + # We infer the return type from a call to the function on the first three elements return_type = typeof(f(collect(layer)[1:min(3, length(layer))])) - newgrid = convert(Matrix{Union{Nothing,return_type}}, fill(nothing, size(layer))) - N = SimpleSDMResponse(newgrid, layer) - pixels = [] - for lat in latitudes(layer) - for lon in longitudes(layer) - if !isnothing(layer[lon,lat]) - push!(pixels, (lon, lat) => layer[lon,lat]) - end - end - end - for p1 in pixels - ok = filter(p2 -> haversine(p2.first, p1.first) < d, pixels) - val = [p2.second for p2 in ok] - N[p1.first...] = f(val) - end + # New layer using typed similar + N = similar(layer, return_type) - #internal_types = unique(typeof.(N.grid)) - #N.grid = convert(Matrix{Union{internal_types...}}, N.grid) - N = typeof(layer) <: SimpleSDMPredictor ? convert(SimpleSDMPredictor, N) : N + # Store latitudes and longitudes + _lat, _lon = latitudes(layer), longitudes(layer) + + # Pre-allocation of a vector of pairs containing the pixel and the + filled_positions = CartesianIndices(layer.grid)[findall(!isnothing, layer.grid)] + + # We then filter in the occupied positions + for pos in filled_positions + neighbors = filter(p -> haversine((_lon[Tuple(pos)[2]], _lat[Tuple(pos)[1]]), (_lon[Tuple(p)[2]], _lat[Tuple(p)[1]])) < d, filled_positions) + N.grid[pos] = f(layer.grid[neighbors]) + end + # And we return the object return N end