diff --git a/nes/methods/spatial_join.py b/nes/methods/spatial_join.py
index 20dcd06c4a6cac472014bc9cee85844176f8613b..4de3eff4abe3ea082131d84c8809e467c58166af 100644
--- a/nes/methods/spatial_join.py
+++ b/nes/methods/spatial_join.py
@@ -10,7 +10,7 @@ import pandas as pd
 from shapely.geos import TopologicalError
 
 
-def spatial_join(self, ext_shp, method=None, var_list=None, info=False):
+def spatial_join(self, ext_shp, method=None, var_list=None, info=False, apply_bbox=True):
     """
     Compute overlay intersection of two GeoPandasDataFrames.
 
@@ -24,7 +24,9 @@ def spatial_join(self, ext_shp, method=None, var_list=None, info=False):
     var_list : List or None or str
         Variables that will be included in the resulting shapefile.
     info : bool
-        Indicates if you want to print the process info or no
+        Indicates if you want to print the process info
+    apply_bbox : bool
+        Indicates if you want to reduce the shapefile to a bbox
     """
     
     if self.master and info:
@@ -40,7 +42,8 @@ def spatial_join(self, ext_shp, method=None, var_list=None, info=False):
             sys.stdout.flush()
         self.create_shapefile()
 
-    ext_shp = prepare_external_shapefile(self, ext_shp=ext_shp, var_list=var_list, info=info)
+    ext_shp = prepare_external_shapefile(self, ext_shp=ext_shp, var_list=var_list, info=info, 
+                                         apply_bbox=apply_bbox)
 
     if method == 'nearest':
         # Nearest centroids to the shapefile polygons
@@ -59,7 +62,7 @@ def spatial_join(self, ext_shp, method=None, var_list=None, info=False):
     return None
 
 
-def prepare_external_shapefile(self, ext_shp, var_list, info=False):
+def prepare_external_shapefile(self, ext_shp, var_list, info=False, apply_bbox=True):
     """
     Prepare the external shapefile.
 
@@ -78,20 +81,24 @@ def prepare_external_shapefile(self, ext_shp, var_list, info=False):
         External shapefile variables to be computed
     info : bool
         Indicates if you want to print the information
+    apply_bbox : bool
+        Indicates if you want to reduce the shapefile to a bbox
 
     Returns
     -------
     GeoDataFrame
         External shapefile
     """
-
+    
     if isinstance(ext_shp, str):
         # Reading external shapefile
         if self.master and info:
             print("\tReading external shapefile")
         # ext_shp = gpd.read_file(ext_shp, include_fields=var_list, mask=self.shapefile.geometry)
-        ext_shp = gpd.read_file(ext_shp, include_fields=var_list, bbox=get_bbox(self))
-    
+        if apply_bbox:
+            ext_shp = gpd.read_file(ext_shp, include_fields=var_list, bbox=get_bbox(self))
+        else:
+            ext_shp = gpd.read_file(ext_shp, include_fields=var_list)
     else:
         msg = "WARNING!!! "
         msg += "External shapefile already read. If you pass the path to the shapefile instead of the opened shapefile "
diff --git a/nes/nc_projections/default_nes.py b/nes/nc_projections/default_nes.py
index a7c0fd4544f2359c12602f510cbf90213cf5c843..0c64e6ba2f0d3e7def86e154c74d191b213055d6 100644
--- a/nes/nc_projections/default_nes.py
+++ b/nes/nc_projections/default_nes.py
@@ -3647,7 +3647,7 @@ class Nes(object):
             self, dst_grid, weight_matrix_path=weight_matrix_path, kind=kind, n_neighbours=n_neighbours, info=info, 
             to_providentia=to_providentia, only_create_wm=only_create_wm, wm=wm, flux=flux)
 
-    def spatial_join(self, ext_shp, method=None, var_list=None, info=False):
+    def spatial_join(self, ext_shp, method=None, var_list=None, info=False, apply_bbox=True):
         """
         Compute overlay intersection of two GeoPandasDataFrames.
 
@@ -3660,10 +3660,13 @@ class Nes(object):
         var_list : List or None
             Variables that will be included in the resulting shapefile.
         info : bool
-            Indicates if you want to print the process info or no
+            Indicates if you want to print the process info
+        apply_bbox : bool
+            Indicates if you want to reduce the shapefile to a bbox
         """
 
-        return spatial_join(self, ext_shp=ext_shp, method=method, var_list=var_list, info=info)
+        return spatial_join(self, ext_shp=ext_shp, method=method, var_list=var_list, info=info,
+                            apply_bbox=apply_bbox)
 
     def calculate_grid_area(self, overwrite=True):
         """