diff --git a/Project.toml b/Project.toml
index 56ab6801b..d5b72827d 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "SpeciesDistributionToolkit"
 uuid = "72b53823-5c0b-4575-ad0e-8e97227ad13b"
 authors = ["Timothée Poisot <timothee.poisot@umontreal.ca>"]
-version = "1.2.0"
+version = "1.2.1"
 
 [deps]
 Distances = "b4f34e82-e78d-54a5-968a-f98e89d6e8f7"
diff --git a/docs/src/reference/CHANGELOG.md b/docs/src/reference/CHANGELOG.md
index 36b611a91..7e95b4d76 100644
--- a/docs/src/reference/CHANGELOG.md
+++ b/docs/src/reference/CHANGELOG.md
@@ -5,6 +5,10 @@ All notable changes to this project will be documented in this file.
 The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.1.0/),
 and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
 
+## `v1.2.1`
+
+- **improved** the performance of masking a vector of layers using polygons
+
 ## `v1.2.0`
 
 - **added** `SimpleSDMDatasets` version 1
diff --git a/src/polygons/polygons.jl b/src/polygons/polygons.jl
index 94b6a0397..60050a47d 100644
--- a/src/polygons/polygons.jl
+++ b/src/polygons/polygons.jl
@@ -32,6 +32,7 @@ Return a trimmed version of a layer, according to the feature defined a
 trim(layer::SDMLayer, feature::T) where {T <: GeoJSON.GeoJSONT} =
     trim(mask!(copy(layer), feature))
 
+
 function change_inclusion!(inclusion, layer, polygon, op)
     xytrans = SimpleSDMLayers.Proj.Transformation(
         "+proj=longlat +datum=WGS84 +no_defs",
@@ -68,25 +69,47 @@ function change_inclusion!(inclusion, layer, polygon, op)
     return inclusion
 end
 
+function _get_inclusion_from_polygon!(inclusion, layer, multipolygon::GeoJSON.MultiPolygon)
+    for element in multipolygon
+        for i in eachindex(element)
+            change_inclusion!(inclusion, layer, element[i], true)
+        end
+    end 
+end 
+
+function SimpleSDMLayers.mask!(layers::Vector{SDMLayer}, multipolygon::GeoJSON.MultiPolygon) 
+    inclusion = .!reduce(.|, [l.indices for l in layers])
+    _get_inclusion_from_polygon!(inclusion, first(layers), multipolygon)
+    for layer in layers
+        layer.indices .&= inclusion
+    end
+    return layers
+end
+
 """
     SimpleSDMLayers.mask!(layer::SDMLayer, multipolygon::GeoJSON.MultiPolygon)
 
-Turns of fall the cells outside the polygon (or within holes in the polygon). This modifies the object.
+Turns off all the cells outside the polygon (or within holes in the polygon). This modifies the object.
 """
 function SimpleSDMLayers.mask!(layer::SDMLayer, multipolygon::GeoJSON.MultiPolygon)
     inclusion = zeros(eltype(layer.indices), size(layer))
-    for element in multipolygon
-        change_inclusion!(inclusion, layer, element[1], true)
-        if length(element) > 2
-            for i in 2:length(element)
-                change_inclusion!(inclusion, layer, element[i], false)
-            end
-        end
-    end
+    _get_inclusion_from_polygon!(inclusion, layer, multipolygon)
     layer.indices .&= inclusion
     return layer
 end
 
+SimpleSDMLayers.mask!(layer::SDMLayer, features::GeoJSON.FeatureCollection, feature = 1) =
+    mask!(layer, features[feature])
+
+SimpleSDMLayers.mask!(layers::Vector{SDMLayer}, features::GeoJSON.FeatureCollection, feature = 1) =
+    mask!(layers, features[feature])
+    
+SimpleSDMLayers.mask!(layer::SDMLayer, feature::GeoJSON.Feature) =
+    mask!(layer, feature.geometry)
+
+SimpleSDMLayers.mask!(layers::Vector{SDMLayer}, feature::GeoJSON.Feature) =
+    mask!(layers, feature.geometry)
+
 """
     SimpleSDMLayers.mask(occ::T, multipolygon::GeoJSON.MultiPolygon) where {T <: AsbtractOccurrenceCollection}
 
@@ -120,18 +143,12 @@ function SimpleSDMLayers.mask(
     return elements(occ)[findall(inclusion)]
 end
 
-SimpleSDMLayers.mask!(layer::SDMLayer, features::GeoJSON.FeatureCollection, feature = 1) =
-    mask!(layer, features[feature])
-
 SimpleSDMLayers.mask(
     occ::T,
     features::GeoJSON.FeatureCollection,
     feature = 1,
 ) where {T <: AbstractOccurrenceCollection} = mask(occ, features[feature])
 
-SimpleSDMLayers.mask!(layer::SDMLayer, feature::GeoJSON.Feature) =
-    mask!(layer, feature.geometry)
-
 SimpleSDMLayers.mask(
     occ::T,
     feature::GeoJSON.Feature,