diff --git a/Project.toml b/Project.toml
index e9e6b0b3..8624f823 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "SimpleSDMLayers"
 uuid = "2c645270-77db-11e9-22c3-0f302a89c64c"
 authors = ["Timothée Poisot <timothee.poisot@umontreal.ca>", "Gabriel Dansereau <gabriel.dansereau@umontreal.ca>"]
-version = "0.4.6"
+version = "0.4.7"
 
 [deps]
 ArchGDAL = "c9ce4bd3-c3d5-55b8-8973-c0e20141b8c3"
diff --git a/src/integrations/GBIF.jl b/src/integrations/GBIF.jl
index 4066e86b..ec296fa7 100644
--- a/src/integrations/GBIF.jl
+++ b/src/integrations/GBIF.jl
@@ -86,6 +86,15 @@ function Base.getindex(layer::T, records::GBIF.GBIFRecords) where {T <: SimpleSD
    return convert(Vector{eltype(layer)}, filter(!isnothing, [layer[records[i]] for i in 1:length(records)]))
 end
 
+"""
+    Base.getindex(layer::T, records::Vector{GBIF.GBIFRecord}) where {T <: SimpleSDMLayer}
+
+Returns the values of a layer at all occurrences in a `GBIFRecord` array.
+"""
+function Base.getindex(layer::T, records::Vector{GBIF.GBIFRecord}) where {T <: SimpleSDMLayer}
+   return [layer[record] for record in records]
+end
+
 """
     mask!(layer::SimpleSDMResponse{T}, records::GBIF.GBIFRecords) where {T <: AbstractBool}
 
@@ -94,7 +103,9 @@ if an occurrence is found in the cell, `false` if not.
 """
 function mask!(layer::SimpleSDMResponse{T}, records::GBIF.GBIFRecords) where {T <: Bool}
     for record in records
-        layer[record] = true
+        if !isnothing(layer[record])
+            layer[record] = true
+        end
     end
     return layer
 end
@@ -107,7 +118,9 @@ the number of occurrences in the cell.
 """
 function mask!(layer::SimpleSDMResponse{T}, records::GBIF.GBIFRecords) where {T <: Number}
     for record in records
-        layer[record] = layer[record] + one(T)
+        if !isnothing(layer[record])
+            layer[record] = layer[record] + one(T)
+        end
     end
     return layer
 end
diff --git a/src/operations/sliding.jl b/src/operations/sliding.jl
index 6a1f03ea..441d3acb 100644
--- a/src/operations/sliding.jl
+++ b/src/operations/sliding.jl
@@ -29,7 +29,7 @@ have data.
 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}
+function slidingwindow(layer::LT, f::FT, d::IT; threaded::Bool=Threads.nthreads()>1) 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))]))
 
@@ -39,13 +39,20 @@ function slidingwindow(layer::LT, f::FT, d::IT) where {LT <: SimpleSDMLayer, FT
     # Store latitudes and longitudes
     _lat, _lon = latitudes(layer), longitudes(layer)
 
-    # Pre-allocation of a vector of pairs containing the pixel and the 
+    # Vector of all positions with a value
     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])
+    if threaded
+        Threads.@threads 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
+    else
+        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
     end
 
     # And we return the object
diff --git a/test/gbif.jl b/test/gbif.jl
index c74f1187..33fbffae 100644
--- a/test/gbif.jl
+++ b/test/gbif.jl
@@ -9,6 +9,9 @@ kingfisher = GBIF.taxon("Megaceryle alcyon", strict=true)
 
 o = GBIF.occurrences(kingfisher, "hasCoordinate" => "true")
 
+# Array of occurrences is the same thing as occurrences collection
+@test temperature[o] == temperature[o.occurrences[1:length(o)]]
+
 # Extract from a single record
 for oc in o
     @test typeof(temperature[oc]) <: Number