diff --git a/nes/nc_projections/default_nes.py b/nes/nc_projections/default_nes.py
index d716b35918de93cc494265f957b476f7ac7d591a..964f8dac35ea20c19074fdccb6da245775dbd731 100644
--- a/nes/nc_projections/default_nes.py
+++ b/nes/nc_projections/default_nes.py
@@ -2112,7 +2112,9 @@ class Nes(object):
         fids = fids.reshape((len(self._lat['data']), len(self._lon['data'])))
         fids = fids[self.read_axis_limits['y_min']:self.read_axis_limits['y_max'],
                     self.read_axis_limits['x_min']:self.read_axis_limits['x_max']]
-        gdf = gpd.GeoDataFrame(index=pd.Index(name='FID', data=fids.ravel()), geometry=geometry, crs="EPSG:4326")
+        gdf = gpd.GeoDataFrame(index=pd.Index(name='FID', data=fids.ravel()), 
+                               geometry=geometry, 
+                               crs="EPSG:4326")
         self.shapefile = gdf
 
         return gdf
@@ -2141,6 +2143,127 @@ class Nes(object):
         
         return None
 
+    def spatial_join(self, mask, method=None):
+        """
+        Compute overlay intersection of two GeoPandasDataFrames
+
+        Parameters
+        ----------
+        mask : GeoPandasDataFrame
+            File from where the data will be obtained on the intersection.
+        method : str
+            Overlay method. Accepted values: ['nearest', 'intersection', None].
+        """
+
+        # Nearest centroids to the mask polygons
+        if method == 'nearest':
+            
+            # Get centroids of shapefile to mask
+            shapefile_aux = deepcopy(self.shapefile)
+            shapefile_aux.geometry = self.shapefile.centroid
+
+            # Calculate spatial joint by distance
+            shapefile_aux = gpd.sjoin_nearest(shapefile_aux, mask.to_crs(self.shapefile.crs), distance_col='distance')
+            
+            # Get data from closest shapes to centroids
+            del shapefile_aux['geometry'], shapefile_aux['index_right']
+            self.shapefile.loc[shapefile_aux.index, shapefile_aux.columns] = shapefile_aux
+        
+        # Intersect the areas of the mask polygons, outside of the mask there will be NaN
+        elif method == 'intersection':
+        
+            # Get intersected areas
+            inp, res = mask.sindex.query_bulk(self.shapefile.geometry, predicate='intersects')
+            print('Rank {0:03d}: {1} intersected areas found'.format(self.rank, len(inp)))
+            
+            # Calculate intersected areas and fractions
+            intersection = pd.concat([self.shapefile.geometry[inp].reset_index(), mask.geometry[res].reset_index()], 
+                                     axis=1, ignore_index=True)
+            intersection.columns = (list(self.shapefile.geometry[inp].reset_index().columns) + 
+                                    list(mask.geometry[res].reset_index().rename(columns={'geometry': 'geometry_mask', 
+                                                                                          'index': 'index_mask'}).columns))
+            intersection['area'] = intersection.apply(lambda x: x['geometry'].intersection(x['geometry_mask']).buffer(0).area, 
+                                                      axis=1)
+            intersection['fraction'] = intersection.apply(lambda x: x['area'] / x['geometry'].area, axis=1)
+            
+            # Choose biggest area from intersected areas with multiple options
+            intersection.sort_values('fraction', ascending=False, inplace=True)
+            intersection = intersection.drop_duplicates(subset='FID', keep="first")
+            intersection = intersection.sort_values('FID').set_index('FID')
+
+            # Get data from mask
+            del mask['geometry']
+            self.shapefile.loc[intersection.index, mask.columns] = np.array(mask.loc[intersection.index_mask, :])
+
+        # Centroids that fall on the mask polygons, outside of the mask there will be NaN
+        elif method is None:
+
+            # Get centroids of shapefile to mask
+            shapefile_aux = deepcopy(self.shapefile)
+            shapefile_aux.geometry = self.shapefile.centroid
+
+            # Calculate spatial joint
+            shapefile_aux = gpd.sjoin(shapefile_aux, mask.to_crs(self.shapefile.crs))
+            
+            # Get data from shapes where there are centroids, rest will be NaN
+            del shapefile_aux['geometry'], shapefile_aux['index_right']
+            self.shapefile.loc[shapefile_aux.index, shapefile_aux.columns] = shapefile_aux
+            
+        return None
+
+    @staticmethod
+    def spatial_overlays(df1, df2, how='intersection'):
+        """
+        Compute overlay intersection of two GeoPandasDataFrames df1 and df2
+
+        https://github.com/geopandas/geopandas/issues/400
+
+        :param df1: GeoDataFrame
+        :param df2: GeoDataFrame
+        :param how: Operation to do
+        :return: GeoDataFrame
+        """
+        from functools import reduce
+
+        df1 = df1.copy()
+        df2 = df2.copy()
+        df1['geometry'] = df1.geometry.buffer(0)
+        df2['geometry'] = df2.geometry.buffer(0)
+        if how == 'intersection':
+            # Spatial Index to create intersections
+            spatial_index = df2.sindex
+            df1['bbox'] = df1.geometry.apply(lambda x: x.bounds)
+            df1['histreg'] = df1.bbox.apply(lambda x: list(spatial_index.intersection(x)))
+            pairs = df1['histreg'].to_dict()
+            nei = []
+            for i, j in pairs.items():
+                for k in j:
+                    nei.append([i, k])
+            pairs = pd.DataFrame(nei, columns=['idx1', 'idx2'])
+            pairs = pairs.merge(df1, left_on='idx1', right_index=True)
+            pairs = pairs.merge(df2, left_on='idx2', right_index=True, suffixes=['_1', '_2'])
+            pairs['geometry'] = pairs.apply(lambda x: (x['geometry_1'].intersection(x['geometry_2'])).buffer(0), axis=1)
+
+            pairs.drop(columns=['geometry_1', 'geometry_2', 'histreg', 'bbox'], inplace=True)
+            pairs = gpd.GeoDataFrame(pairs, columns=pairs.columns, crs=df1.crs)
+            pairs = pairs.loc[~pairs.geometry.is_empty]
+
+            return_value = pairs
+        elif how == 'difference':
+            spatial_index = df2.sindex
+            df1['bbox'] = df1.geometry.apply(lambda x: x.bounds)
+            df1['histreg'] = df1.bbox.apply(lambda x: list(spatial_index.intersection(x)))
+            df1['new_g'] = df1.apply(lambda x: reduce(lambda x, y: x.difference(y).buffer(0),
+                                                      [x.geometry] + list(df2.iloc[x.histreg].geometry)), axis=1)
+            df1.geometry = df1.new_g
+            df1 = df1.loc[~df1.geometry.is_empty].copy()
+            df1.drop(['bbox', 'histreg', 'new_g'], axis=1, inplace=True)
+            return_value = df1
+        else:
+            raise NotImplementedError(how)
+
+        return return_value
+
     def __gather_data_py_object(self):
         """
         Gather all the variable data into the MPI rank 0 to perform a serial write.
diff --git a/nes/nc_projections/lcc_nes.py b/nes/nc_projections/lcc_nes.py
index e8048044774c87488130d87ade90335717c5663d..f9aa7e282507ca4b3e9436929df02ed40adc77fe 100644
--- a/nes/nc_projections/lcc_nes.py
+++ b/nes/nc_projections/lcc_nes.py
@@ -455,7 +455,9 @@ class LCCNes(Nes):
         fids = fids.reshape((self._lat['data'].shape[0],self._lat['data'].shape[1]))
         fids = fids[self.read_axis_limits['y_min']:self.read_axis_limits['y_max'],
                     self.read_axis_limits['x_min']:self.read_axis_limits['x_max']]
-        gdf = gpd.GeoDataFrame(index=pd.Index(name='FID', data=fids.ravel()), geometry=geometry, crs="EPSG:4326")
+        gdf = gpd.GeoDataFrame(index=pd.Index(name='FID', data=fids.ravel()), 
+                               geometry=geometry, 
+                               crs="EPSG:4326")
         self.shapefile = gdf
 
         return gdf
diff --git a/nes/nc_projections/mercator_nes.py b/nes/nc_projections/mercator_nes.py
index 4ca7d4e3cb91e18d0acc67a7ab784700377c06bc..7c42bd6d52424002aac54b1abbf5bd44160df811 100644
--- a/nes/nc_projections/mercator_nes.py
+++ b/nes/nc_projections/mercator_nes.py
@@ -432,7 +432,9 @@ class MercatorNes(Nes):
         fids = fids.reshape((self._lat['data'].shape[0],self._lat['data'].shape[1]))
         fids = fids[self.read_axis_limits['y_min']:self.read_axis_limits['y_max'],
                     self.read_axis_limits['x_min']:self.read_axis_limits['x_max']]
-        gdf = gpd.GeoDataFrame(index=pd.Index(name='FID', data=fids.ravel()), geometry=geometry, crs="EPSG:4326")
+        gdf = gpd.GeoDataFrame(index=pd.Index(name='FID', data=fids.ravel()), 
+                               geometry=geometry, 
+                               crs="EPSG:4326")
         self.shapefile = gdf
 
         return gdf
diff --git a/nes/nc_projections/points_nes.py b/nes/nc_projections/points_nes.py
index f42934c5862cdb2836521dbdfb5a6ac9495bd524..018bd79a351780dfc41c54be17503d972be9c798 100644
--- a/nes/nc_projections/points_nes.py
+++ b/nes/nc_projections/points_nes.py
@@ -293,7 +293,7 @@ class PointsNes(Nes):
         elif len(var_dims) == 2:
             if 'strlen' in nc_var.dimensions:
                 data = nc_var[self.read_axis_limits['x_min']:self.read_axis_limits['x_max'], :]
-                data = np.array([''.join(i) for i in np.char.decode(data)])
+                data = np.array([''.join(i) for i in data])
             else:
                 data = nc_var[self.read_axis_limits['t_min']:self.read_axis_limits['t_max'],
                               self.read_axis_limits['x_min']:self.read_axis_limits['x_max']]
@@ -605,7 +605,6 @@ class PointsNes(Nes):
         variables = {}
         interpolated_variables = deepcopy(self.variables)
         for var_name, var_info in interpolated_variables.items():
-            print(var_name)
             variables[var_name] = {}
             if var_info['dimensions'] == ('time', 'lat', 'lon') and len(var_info['data'].shape) == 2:
                 variables[var_name]['data'] = var_info['data'].T
@@ -646,7 +645,9 @@ class PointsNes(Nes):
         geometry = gpd.points_from_xy(self.lon['data'], self.lat['data'])
         fids = np.arange(len(self._lon['data']))
         fids = fids[self.read_axis_limits['x_min']:self.read_axis_limits['x_max']]
-        gdf = gpd.GeoDataFrame(index=pd.Index(name='FID', data=fids.ravel()), geometry=geometry, crs="EPSG:4326")
+        gdf = gpd.GeoDataFrame(index=pd.Index(name='FID', data=fids), 
+                               geometry=geometry, 
+                               crs="EPSG:4326")
         self.shapefile = gdf
 
         return gdf
diff --git a/nes/nc_projections/rotated_nes.py b/nes/nc_projections/rotated_nes.py
index 54f9e44d2a4ef4068e0303878a4bf8a27a07d0f9..f933c184b41351536c371c5af5aed292636d607c 100644
--- a/nes/nc_projections/rotated_nes.py
+++ b/nes/nc_projections/rotated_nes.py
@@ -481,7 +481,9 @@ class RotatedNes(Nes):
         fids = fids.reshape((self._lat['data'].shape[0],self._lat['data'].shape[1]))
         fids = fids[self.read_axis_limits['y_min']:self.read_axis_limits['y_max'],
                     self.read_axis_limits['x_min']:self.read_axis_limits['x_max']]
-        gdf = gpd.GeoDataFrame(index=pd.Index(name='FID', data=fids.ravel()), geometry=geometry, crs="EPSG:4326")
+        gdf = gpd.GeoDataFrame(index=pd.Index(name='FID', data=fids.ravel()), 
+                               geometry=geometry, 
+                               crs="EPSG:4326")
         self.shapefile = gdf
 
         return gdf
diff --git a/tests/1-nes_tests_by_size.py b/tests/1-test_read_write_size.py
similarity index 100%
rename from tests/1-nes_tests_by_size.py
rename to tests/1-test_read_write_size.py
diff --git a/tests/2-nes_tests_by_projection.py b/tests/2-test_read_write_projection.py
similarity index 98%
rename from tests/2-nes_tests_by_projection.py
rename to tests/2-test_read_write_projection.py
index 12ef38e30ef6c7a2e903f945554e3efff247b8dd..5f4590c1482f5e2690643caeaaf765ab9db2cbb2 100644
--- a/tests/2-nes_tests_by_projection.py
+++ b/tests/2-test_read_write_projection.py
@@ -25,7 +25,7 @@ paths = {'regular_file': {'path': '/gpfs/scratch/bsc32/bsc32538/mr_multiplyby/or
                       'projection': 'lcc',
                       'variables': [], # all
                       'parallel_methods': ['X', 'Y', 'T']},
-         'mercator_file': {'path': '/esarchive/scratch/avilanova/software/NES/Jupyter_notebooks/input/mercator_grid_example.nc',
+         'mercator_file': {'path': '/esarchive/scratch/avilanova/software/NES/tutorials/data/mercator_grid_example.nc',
                            'projection': 'mercator',
                            'variables': [], # all
                            'parallel_methods': ['X', 'Y', 'T']}
diff --git a/tests/scalability_test_nord3v2.bash b/tests/2-test_read_write_projection_nord3v2.bash
similarity index 96%
rename from tests/scalability_test_nord3v2.bash
rename to tests/2-test_read_write_projection_nord3v2.bash
index c7256757d51f3597aa9da92d265745e507d8a1d3..87fb5c63555df7c92d09891d4d5b43dbf201ceb7 100644
--- a/tests/scalability_test_nord3v2.bash
+++ b/tests/2-test_read_write_projection_nord3v2.bash
@@ -3,7 +3,7 @@
 #EXPORTPATH="/esarchive/scratch/avilanova/software/NES"
 EXPORTPATH="/gpfs/projects/bsc32/models/NES"
 SRCPATH="/esarchive/scratch/avilanova/software/NES/tests"
-EXE="2-nes_tests_by_projection.py"
+EXE="2-test_read_write_projection.py"
 
 module purge
 module load Python/3.7.4-GCCcore-8.3.0
diff --git a/tests/3-test_spatial_join.py b/tests/3-test_spatial_join.py
new file mode 100644
index 0000000000000000000000000000000000000000..e456c64ee2cd1ea24e63e864f851678d37cc5936
--- /dev/null
+++ b/tests/3-test_spatial_join.py
@@ -0,0 +1,124 @@
+#!/usr/bin/env python
+import geopandas as gpd
+import pandas as pd
+import timeit
+import sys
+from mpi4py import MPI
+from nes import *
+
+# Hide warning
+pd.options.mode.chained_assignment = None
+
+comm = MPI.COMM_WORLD
+rank = comm.Get_rank()
+size = comm.Get_size()
+
+results = []
+for method in ['spatial_overlay', 'spatial_join']:
+    for projection in ['regular', 'rotated']:
+        for projection_type in ['created', 'read']:
+
+            # Regular projection
+            if projection == 'regular':
+                # Create dataset and get shapefile
+                if projection_type == 'created':
+                    lat_orig = 41.1
+                    lon_orig = 1.8
+                    inc_lat = 0.2
+                    inc_lon = 0.2
+                    n_lat = 100
+                    n_lon = 100
+                    coordinates = create_nes(comm=None, info=True, 
+                                             projection='regular',
+                                             lat_orig=lat_orig, 
+                                             lon_orig=lon_orig, 
+                                             inc_lat=inc_lat, 
+                                             inc_lon=inc_lon, 
+                                             n_lat=n_lat, 
+                                             n_lon=n_lon)
+                    coordinates.create_shapefile()
+
+                # Open dataset and get shapefile
+                elif projection_type == 'read':
+                    coordinates_path = '/gpfs/scratch/bsc32/bsc32538/mr_multiplyby/original_file/MONARCH_d01_2008123100.nc'
+                    coordinates = open_netcdf(path=coordinates_path, info=True)
+                    coordinates.create_shapefile()
+                    coordinates.keep_vars(['O3'])
+                    coordinates.load()
+                    coordinates.shapefile['O3'] = coordinates.variables['O3']['data'][-1, -1, :].ravel()
+
+            # Rotated projection
+            elif projection == 'rotated':
+                # Create dataset and get shapefile
+                if projection_type == 'created':
+                    centre_lat = 51
+                    centre_lon = 10
+                    west_boundary = -35
+                    south_boundary = -27
+                    inc_rlat = 0.2
+                    inc_rlon = 0.2
+                    coordinates = create_nes(comm=None, info=True, 
+                                             projection='rotated',
+                                             centre_lat=centre_lat, 
+                                             centre_lon=centre_lon,
+                                             west_boundary=west_boundary, 
+                                             south_boundary=south_boundary,
+                                             inc_rlat=inc_rlat, 
+                                             inc_rlon=inc_rlon)
+                    coordinates.create_shapefile()
+
+                # Open dataset and get shapefile
+                elif projection_type == 'read':
+                    coordinates_path = '/gpfs/scratch/bsc32/bsc32538/original_files/CAMS_MONARCH_d01_2022070412.nc'
+                    coordinates = open_netcdf(path=coordinates_path, info=True)
+                    coordinates.create_shapefile()
+                    coordinates.keep_vars(['O3'])
+                    coordinates.load()
+                    coordinates.shapefile['O3'] = coordinates.variables['O3']['data'][-1, -1, :].ravel()
+
+            coordinates.write_shapefile('coordinates_{0}_{1}'.format(projection, 
+                                                                     projection_type))
+
+            mask_path = '/esarchive/scratch/avilanova/software/NES/tutorials/data/timezones_2021c/timezones_2021c.shp'
+            mask = gpd.read_file(mask_path)
+
+            # Spatial overlay (old method)
+            if method == 'spatial_overlay':
+                start_time = timeit.default_timer()
+                intersection = coordinates.shapefile.copy()
+                #intersection['area'] = intersection.geometry.area
+                intersection = coordinates.spatial_overlays(intersection, mask)
+                #intersection.rename(columns={'idx1': 'FID', 'idx2': 'shp_id'}, inplace=True)
+                #intersection['fraction'] = intersection.geometry.area / intersection['area']
+                #intersection.sort_values('fraction', ascending=False, inplace=True)
+                #intersection = intersection.drop_duplicates(subset='FID', keep="first")
+                #intersection.set_index('FID', inplace=True)
+                #coordinates.loc[intersection.index, coordinates.shp_colname] = intersection[coordinates.shp_colname]
+                time = timeit.default_timer() - start_time
+
+            # Spatial join (new method)
+            elif method == 'spatial_join':
+                start_time = timeit.default_timer()
+                coordinates.spatial_join(mask, method='intersection')
+                time = timeit.default_timer() - start_time
+
+            coordinates.write_shapefile('masked_coordinates_{0}_{1}_{2}'.format(projection, 
+                                                                                projection_type, 
+                                                                                method))
+
+        results = []
+        results.append({'Projection': projection,
+                        'Projection type': projection_type,
+                        'Method': '{min:02d}:{sec:02.3f}'.format(
+                            min=int(time // 60), sec=time - (int(time // 60) * 60))
+                        })
+
+        comm.Barrier()
+
+comm.Barrier()
+
+if rank == 0:
+    table = pd.DataFrame(results)
+    print('RESULTS TABLE')
+    print(table)
+    sys.stdout.flush()
diff --git a/tests/3-test_spatial_join_nord3v2.bash b/tests/3-test_spatial_join_nord3v2.bash
new file mode 100644
index 0000000000000000000000000000000000000000..079aa7928696599a280f2c555704dbd543b65385
--- /dev/null
+++ b/tests/3-test_spatial_join_nord3v2.bash
@@ -0,0 +1,23 @@
+#!/bin/bash
+
+EXPORTPATH="/esarchive/scratch/avilanova/software/NES"
+SRCPATH="/esarchive/scratch/avilanova/software/NES/tests"
+EXE="3-test_spatial_join.py"
+
+module purge
+module load Python/3.7.4-GCCcore-8.3.0
+module load netcdf4-python/1.5.3-foss-2019b-Python-3.7.4
+module load cfunits/1.8-foss-2019b-Python-3.7.4
+module load xarray/0.17.0-foss-2019b-Python-3.7.4
+module load pandas/1.2.4-foss-2019b-Python-3.7.4
+module load mpi4py/3.0.3-foss-2019b-Python-3.7.4
+module load filelock/3.7.1-foss-2019b-Python-3.7.4
+module load pyproj/2.5.0-foss-2019b-Python-3.7.4
+module load eccodes-python/0.9.5-foss-2019b-Python-3.7.4
+module load geopandas/0.8.1-foss-2019b-Python-3.7.4
+module load Shapely/1.7.1-foss-2019b-Python-3.7.4
+
+for nprocs in 1
+do
+    JOB_ID=`sbatch --ntasks=${nprocs} --exclusive --job-name=nes_${nprocs} --output=./log_nord3v2_NES_${nprocs}_%J.out --error=./log_nord3v2_NES_${nprocs}_%J.err -D . --time=02:00:00 --wrap="export PYTHONPATH=${EXPORTPATH}:${PYTHONPATH}; cd ${SRCPATH}; mpirun --mca mpi_warn_on_fork 0 -np ${nprocs} python ${SRCPATH}/${EXE}"`
+done
\ No newline at end of file
diff --git a/tests/4-test_bounds.py b/tests/4-test_bounds.py
new file mode 100644
index 0000000000000000000000000000000000000000..e3d3063fddd238f0dc601a6a1c2837e1bcd98e94
--- /dev/null
+++ b/tests/4-test_bounds.py
@@ -0,0 +1,41 @@
+#!/usr/bin/env python
+import sys
+import timeit
+import pandas as pd
+from mpi4py import MPI
+from nes import *
+
+comm = MPI.COMM_WORLD
+rank = comm.Get_rank()
+
+for projection_type in ['read', 'created']:
+
+    # Open dataset
+    if projection_type == 'read':
+        test_path = "/gpfs/scratch/bsc32/bsc32538/mr_multiplyby/OUT/stats_bnds/monarch/a45g/regional/daily_max/O3_all/O3_all-000_2021080300.nc"
+        nessy = open_netcdf(path=test_path, info=True)
+
+    # Create dataset
+    elif projection_type == 'created':
+        lat_orig = 41.1
+        lon_orig = 1.8
+        inc_lat = 0.2
+        inc_lon = 0.2
+        n_lat = 100
+        n_lon = 100
+        nessy = create_nes(comm=None, info=True, projection='regular',
+                           lat_orig=lat_orig, lon_orig=lon_orig, 
+                           inc_lat=inc_lat, inc_lon=inc_lon, 
+                           n_lat=n_lat, n_lon=n_lon)
+    
+    # Add bounds
+    nessy.create_spatial_bounds()
+
+    print('NES', projection_type, '-', 'Rank', rank, '-', nessy)
+    print('NES', projection_type, '-', 'Rank', rank, '-', 'Lat bounds', 
+          nessy.lat_bnds)
+    print('NES', projection_type, '-', 'Rank', rank, '-', 'Lon bounds', 
+          nessy.lon_bnds)
+
+    comm.Barrier()
+    sys.stdout.flush()
\ No newline at end of file
diff --git a/tests/4-test_bounds_nord3v2.bash b/tests/4-test_bounds_nord3v2.bash
new file mode 100644
index 0000000000000000000000000000000000000000..fc6e06b0824d8682c80372867da0d429df1a80b6
--- /dev/null
+++ b/tests/4-test_bounds_nord3v2.bash
@@ -0,0 +1,23 @@
+#!/bin/bash
+
+EXPORTPATH="/esarchive/scratch/avilanova/software/NES"
+SRCPATH="/esarchive/scratch/avilanova/software/NES/tests"
+EXE="4-test_bounds.py"
+
+module purge
+module load Python/3.7.4-GCCcore-8.3.0
+module load netcdf4-python/1.5.3-foss-2019b-Python-3.7.4
+module load cfunits/1.8-foss-2019b-Python-3.7.4
+module load xarray/0.17.0-foss-2019b-Python-3.7.4
+module load pandas/1.2.4-foss-2019b-Python-3.7.4
+module load mpi4py/3.0.3-foss-2019b-Python-3.7.4
+module load filelock/3.7.1-foss-2019b-Python-3.7.4
+module load pyproj/2.5.0-foss-2019b-Python-3.7.4
+module load eccodes-python/0.9.5-foss-2019b-Python-3.7.4
+module load geopandas/0.8.1-foss-2019b-Python-3.7.4
+module load Shapely/1.7.1-foss-2019b-Python-3.7.4
+
+for nprocs in 1 2
+do
+    JOB_ID=`sbatch --ntasks=${nprocs} --exclusive --job-name=nes_${nprocs} --output=./log_nord3v2_NES_${nprocs}_%J.out --error=./log_nord3v2_NES_${nprocs}_%J.err -D . --time=02:00:00 --wrap="export PYTHONPATH=${EXPORTPATH}:${PYTHONPATH}; cd ${SRCPATH}; mpirun --mca mpi_warn_on_fork 0 -np ${nprocs} python ${SRCPATH}/${EXE}"`
+done
\ No newline at end of file
diff --git a/tutorials/2.Creation/2.5.Create_Points_CSIC.ipynb b/tutorials/2.Creation/2.5.Create_Points_CSIC.ipynb
index 48cd49d6883628b607be18fed904d3f40ceb9282..65d623e8c48565ecd1279c7d57168bbcac3e44c2 100644
--- a/tutorials/2.Creation/2.5.Create_Points_CSIC.ipynb
+++ b/tutorials/2.Creation/2.5.Create_Points_CSIC.ipynb
@@ -1869,7 +1869,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.7.4"
+   "version": "3.6.8"
   },
   "vscode": {
    "interpreter": {
diff --git a/tutorials/5.Others/5.2.Add_Coordinates_Bounds.ipynb b/tutorials/5.Others/5.2.Add_Coordinates_Bounds.ipynb
index 29ca013c2b0ef269dec76430d975e01c6b5c8a97..537349001ae305b0d7ac343a269a5bb13b18ad22 100644
--- a/tutorials/5.Others/5.2.Add_Coordinates_Bounds.ipynb
+++ b/tutorials/5.Others/5.2.Add_Coordinates_Bounds.ipynb
@@ -209,7 +209,7 @@
  ],
  "metadata": {
   "kernelspec": {
-   "display_name": "Python 3",
+   "display_name": "Python 3.6.8 64-bit",
    "language": "python",
    "name": "python3"
   },
@@ -223,7 +223,12 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.7.4"
+   "version": "3.6.8"
+  },
+  "vscode": {
+   "interpreter": {
+    "hash": "31f2aee4e71d21fbe5cf8b01ff0e069b9275f58929596ceb00d14d90e3e16cd6"
+   }
   }
  },
  "nbformat": 4,
diff --git a/tutorials/5.Others/5.4.Spatial_Join.ipynb b/tutorials/5.Others/5.4.Spatial_Join.ipynb
new file mode 100644
index 0000000000000000000000000000000000000000..b00e62288dc5f117c373579760892b9c5d86452a
--- /dev/null
+++ b/tutorials/5.Others/5.4.Spatial_Join.ipynb
@@ -0,0 +1,892 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# How to make spatial joins"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from nes import *\n",
+    "import xarray as xr\n",
+    "import geopandas as gpd\n",
+    "import pandas as pd\n",
+    "\n",
+    "pd.options.mode.chained_assignment = None"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Method 1: Centroids"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/html": [
+       "<div>\n",
+       "<style scoped>\n",
+       "    .dataframe tbody tr th:only-of-type {\n",
+       "        vertical-align: middle;\n",
+       "    }\n",
+       "\n",
+       "    .dataframe tbody tr th {\n",
+       "        vertical-align: top;\n",
+       "    }\n",
+       "\n",
+       "    .dataframe thead th {\n",
+       "        text-align: right;\n",
+       "    }\n",
+       "</style>\n",
+       "<table border=\"1\" class=\"dataframe\">\n",
+       "  <thead>\n",
+       "    <tr style=\"text-align: right;\">\n",
+       "      <th></th>\n",
+       "      <th>geometry</th>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>FID</th>\n",
+       "      <th></th>\n",
+       "    </tr>\n",
+       "  </thead>\n",
+       "  <tbody>\n",
+       "    <tr>\n",
+       "      <th>0</th>\n",
+       "      <td>POLYGON ((1.80000 41.10000, 2.00000 41.10000, ...</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>1</th>\n",
+       "      <td>POLYGON ((2.00000 41.10000, 2.20000 41.10000, ...</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>2</th>\n",
+       "      <td>POLYGON ((2.20000 41.10000, 2.40000 41.10000, ...</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>3</th>\n",
+       "      <td>POLYGON ((2.40000 41.10000, 2.60000 41.10000, ...</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>4</th>\n",
+       "      <td>POLYGON ((2.60000 41.10000, 2.80000 41.10000, ...</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>...</th>\n",
+       "      <td>...</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>9995</th>\n",
+       "      <td>POLYGON ((20.80000 60.90000, 21.00000 60.90000...</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>9996</th>\n",
+       "      <td>POLYGON ((21.00000 60.90000, 21.20000 60.90000...</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>9997</th>\n",
+       "      <td>POLYGON ((21.20000 60.90000, 21.40000 60.90000...</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>9998</th>\n",
+       "      <td>POLYGON ((21.40000 60.90000, 21.60000 60.90000...</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>9999</th>\n",
+       "      <td>POLYGON ((21.60000 60.90000, 21.80000 60.90000...</td>\n",
+       "    </tr>\n",
+       "  </tbody>\n",
+       "</table>\n",
+       "<p>10000 rows × 1 columns</p>\n",
+       "</div>"
+      ],
+      "text/plain": [
+       "                                               geometry\n",
+       "FID                                                    \n",
+       "0     POLYGON ((1.80000 41.10000, 2.00000 41.10000, ...\n",
+       "1     POLYGON ((2.00000 41.10000, 2.20000 41.10000, ...\n",
+       "2     POLYGON ((2.20000 41.10000, 2.40000 41.10000, ...\n",
+       "3     POLYGON ((2.40000 41.10000, 2.60000 41.10000, ...\n",
+       "4     POLYGON ((2.60000 41.10000, 2.80000 41.10000, ...\n",
+       "...                                                 ...\n",
+       "9995  POLYGON ((20.80000 60.90000, 21.00000 60.90000...\n",
+       "9996  POLYGON ((21.00000 60.90000, 21.20000 60.90000...\n",
+       "9997  POLYGON ((21.20000 60.90000, 21.40000 60.90000...\n",
+       "9998  POLYGON ((21.40000 60.90000, 21.60000 60.90000...\n",
+       "9999  POLYGON ((21.60000 60.90000, 21.80000 60.90000...\n",
+       "\n",
+       "[10000 rows x 1 columns]"
+      ]
+     },
+     "execution_count": 2,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "lat_orig = 41.1\n",
+    "lon_orig = 1.8\n",
+    "inc_lat = 0.2\n",
+    "inc_lon = 0.2\n",
+    "n_lat = 100\n",
+    "n_lon = 100\n",
+    "coordinates = create_nes(comm=None, info=False, projection='regular',\n",
+    "                         lat_orig=lat_orig, lon_orig=lon_orig, inc_lat=inc_lat, inc_lon=inc_lon, \n",
+    "                         n_lat=n_lat, n_lon=n_lon)\n",
+    "coordinates.create_shapefile()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 3,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "mask_path = '/esarchive/scratch/avilanova/software/NES/tutorials/data/timezones_2021c/timezones_2021c.shp'\n",
+    "mask = gpd.read_file(mask_path)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 4,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stderr",
+     "output_type": "stream",
+     "text": [
+      "/esarchive/scratch/avilanova/software/NES/nes/nc_projections/default_nes.py:2190: UserWarning: Geometry is in a geographic CRS. Results from 'centroid' are likely incorrect. Use 'GeoSeries.to_crs()' to re-project geometries to a projected CRS before this operation.\n",
+      "\n",
+      "  shapefile_aux.geometry = self.shapefile.centroid\n"
+     ]
+    },
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "CPU times: user 15.2 s, sys: 680 ms, total: 15.9 s\n",
+      "Wall time: 15.9 s\n"
+     ]
+    }
+   ],
+   "source": [
+    "%time coordinates.spatial_join(mask)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 5,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/html": [
+       "<div>\n",
+       "<style scoped>\n",
+       "    .dataframe tbody tr th:only-of-type {\n",
+       "        vertical-align: middle;\n",
+       "    }\n",
+       "\n",
+       "    .dataframe tbody tr th {\n",
+       "        vertical-align: top;\n",
+       "    }\n",
+       "\n",
+       "    .dataframe thead th {\n",
+       "        text-align: right;\n",
+       "    }\n",
+       "</style>\n",
+       "<table border=\"1\" class=\"dataframe\">\n",
+       "  <thead>\n",
+       "    <tr style=\"text-align: right;\">\n",
+       "      <th></th>\n",
+       "      <th>geometry</th>\n",
+       "      <th>tzid</th>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>FID</th>\n",
+       "      <th></th>\n",
+       "      <th></th>\n",
+       "    </tr>\n",
+       "  </thead>\n",
+       "  <tbody>\n",
+       "    <tr>\n",
+       "      <th>0</th>\n",
+       "      <td>POLYGON ((1.80000 41.10000, 2.00000 41.10000, ...</td>\n",
+       "      <td>Europe/Madrid</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>1</th>\n",
+       "      <td>POLYGON ((2.00000 41.10000, 2.20000 41.10000, ...</td>\n",
+       "      <td>Europe/Madrid</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>2</th>\n",
+       "      <td>POLYGON ((2.20000 41.10000, 2.40000 41.10000, ...</td>\n",
+       "      <td>Europe/Madrid</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>3</th>\n",
+       "      <td>POLYGON ((2.40000 41.10000, 2.60000 41.10000, ...</td>\n",
+       "      <td>NaN</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>4</th>\n",
+       "      <td>POLYGON ((2.60000 41.10000, 2.80000 41.10000, ...</td>\n",
+       "      <td>NaN</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>...</th>\n",
+       "      <td>...</td>\n",
+       "      <td>...</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>9995</th>\n",
+       "      <td>POLYGON ((20.80000 60.90000, 21.00000 60.90000...</td>\n",
+       "      <td>Europe/Helsinki</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>9996</th>\n",
+       "      <td>POLYGON ((21.00000 60.90000, 21.20000 60.90000...</td>\n",
+       "      <td>Europe/Helsinki</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>9997</th>\n",
+       "      <td>POLYGON ((21.20000 60.90000, 21.40000 60.90000...</td>\n",
+       "      <td>Europe/Helsinki</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>9998</th>\n",
+       "      <td>POLYGON ((21.40000 60.90000, 21.60000 60.90000...</td>\n",
+       "      <td>Europe/Helsinki</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>9999</th>\n",
+       "      <td>POLYGON ((21.60000 60.90000, 21.80000 60.90000...</td>\n",
+       "      <td>Europe/Helsinki</td>\n",
+       "    </tr>\n",
+       "  </tbody>\n",
+       "</table>\n",
+       "<p>10000 rows × 2 columns</p>\n",
+       "</div>"
+      ],
+      "text/plain": [
+       "                                               geometry             tzid\n",
+       "FID                                                                     \n",
+       "0     POLYGON ((1.80000 41.10000, 2.00000 41.10000, ...    Europe/Madrid\n",
+       "1     POLYGON ((2.00000 41.10000, 2.20000 41.10000, ...    Europe/Madrid\n",
+       "2     POLYGON ((2.20000 41.10000, 2.40000 41.10000, ...    Europe/Madrid\n",
+       "3     POLYGON ((2.40000 41.10000, 2.60000 41.10000, ...              NaN\n",
+       "4     POLYGON ((2.60000 41.10000, 2.80000 41.10000, ...              NaN\n",
+       "...                                                 ...              ...\n",
+       "9995  POLYGON ((20.80000 60.90000, 21.00000 60.90000...  Europe/Helsinki\n",
+       "9996  POLYGON ((21.00000 60.90000, 21.20000 60.90000...  Europe/Helsinki\n",
+       "9997  POLYGON ((21.20000 60.90000, 21.40000 60.90000...  Europe/Helsinki\n",
+       "9998  POLYGON ((21.40000 60.90000, 21.60000 60.90000...  Europe/Helsinki\n",
+       "9999  POLYGON ((21.60000 60.90000, 21.80000 60.90000...  Europe/Helsinki\n",
+       "\n",
+       "[10000 rows x 2 columns]"
+      ]
+     },
+     "execution_count": 5,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "coordinates.shapefile"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "coordinates.shapefile.to_file('spatial_join_method_1')"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Method 2: Nearest centroids"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 6,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/html": [
+       "<div>\n",
+       "<style scoped>\n",
+       "    .dataframe tbody tr th:only-of-type {\n",
+       "        vertical-align: middle;\n",
+       "    }\n",
+       "\n",
+       "    .dataframe tbody tr th {\n",
+       "        vertical-align: top;\n",
+       "    }\n",
+       "\n",
+       "    .dataframe thead th {\n",
+       "        text-align: right;\n",
+       "    }\n",
+       "</style>\n",
+       "<table border=\"1\" class=\"dataframe\">\n",
+       "  <thead>\n",
+       "    <tr style=\"text-align: right;\">\n",
+       "      <th></th>\n",
+       "      <th>geometry</th>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>FID</th>\n",
+       "      <th></th>\n",
+       "    </tr>\n",
+       "  </thead>\n",
+       "  <tbody>\n",
+       "    <tr>\n",
+       "      <th>0</th>\n",
+       "      <td>POLYGON ((1.80000 41.10000, 2.00000 41.10000, ...</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>1</th>\n",
+       "      <td>POLYGON ((2.00000 41.10000, 2.20000 41.10000, ...</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>2</th>\n",
+       "      <td>POLYGON ((2.20000 41.10000, 2.40000 41.10000, ...</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>3</th>\n",
+       "      <td>POLYGON ((2.40000 41.10000, 2.60000 41.10000, ...</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>4</th>\n",
+       "      <td>POLYGON ((2.60000 41.10000, 2.80000 41.10000, ...</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>...</th>\n",
+       "      <td>...</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>95</th>\n",
+       "      <td>POLYGON ((2.80000 42.90000, 3.00000 42.90000, ...</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>96</th>\n",
+       "      <td>POLYGON ((3.00000 42.90000, 3.20000 42.90000, ...</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>97</th>\n",
+       "      <td>POLYGON ((3.20000 42.90000, 3.40000 42.90000, ...</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>98</th>\n",
+       "      <td>POLYGON ((3.40000 42.90000, 3.60000 42.90000, ...</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>99</th>\n",
+       "      <td>POLYGON ((3.60000 42.90000, 3.80000 42.90000, ...</td>\n",
+       "    </tr>\n",
+       "  </tbody>\n",
+       "</table>\n",
+       "<p>100 rows × 1 columns</p>\n",
+       "</div>"
+      ],
+      "text/plain": [
+       "                                              geometry\n",
+       "FID                                                   \n",
+       "0    POLYGON ((1.80000 41.10000, 2.00000 41.10000, ...\n",
+       "1    POLYGON ((2.00000 41.10000, 2.20000 41.10000, ...\n",
+       "2    POLYGON ((2.20000 41.10000, 2.40000 41.10000, ...\n",
+       "3    POLYGON ((2.40000 41.10000, 2.60000 41.10000, ...\n",
+       "4    POLYGON ((2.60000 41.10000, 2.80000 41.10000, ...\n",
+       "..                                                 ...\n",
+       "95   POLYGON ((2.80000 42.90000, 3.00000 42.90000, ...\n",
+       "96   POLYGON ((3.00000 42.90000, 3.20000 42.90000, ...\n",
+       "97   POLYGON ((3.20000 42.90000, 3.40000 42.90000, ...\n",
+       "98   POLYGON ((3.40000 42.90000, 3.60000 42.90000, ...\n",
+       "99   POLYGON ((3.60000 42.90000, 3.80000 42.90000, ...\n",
+       "\n",
+       "[100 rows x 1 columns]"
+      ]
+     },
+     "execution_count": 6,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "lat_orig = 41.1\n",
+    "lon_orig = 1.8\n",
+    "inc_lat = 0.2\n",
+    "inc_lon = 0.2\n",
+    "n_lat = 10\n",
+    "n_lon = 10\n",
+    "coordinates = create_nes(comm=None, info=False, projection='regular',\n",
+    "                         lat_orig=lat_orig, lon_orig=lon_orig, inc_lat=inc_lat, inc_lon=inc_lon, \n",
+    "                         n_lat=n_lat, n_lon=n_lon)\n",
+    "coordinates.create_shapefile()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 22,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "mask_path = '/esarchive/scratch/avilanova/software/NES/tutorials/data/timezones_2021c/timezones_2021c.shp'\n",
+    "mask = gpd.read_file(mask_path)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 8,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "#%time coordinates.spatial_join(mask, method='nearest')"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 9,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/html": [
+       "<div>\n",
+       "<style scoped>\n",
+       "    .dataframe tbody tr th:only-of-type {\n",
+       "        vertical-align: middle;\n",
+       "    }\n",
+       "\n",
+       "    .dataframe tbody tr th {\n",
+       "        vertical-align: top;\n",
+       "    }\n",
+       "\n",
+       "    .dataframe thead th {\n",
+       "        text-align: right;\n",
+       "    }\n",
+       "</style>\n",
+       "<table border=\"1\" class=\"dataframe\">\n",
+       "  <thead>\n",
+       "    <tr style=\"text-align: right;\">\n",
+       "      <th></th>\n",
+       "      <th>geometry</th>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>FID</th>\n",
+       "      <th></th>\n",
+       "    </tr>\n",
+       "  </thead>\n",
+       "  <tbody>\n",
+       "    <tr>\n",
+       "      <th>0</th>\n",
+       "      <td>POLYGON ((1.80000 41.10000, 2.00000 41.10000, ...</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>1</th>\n",
+       "      <td>POLYGON ((2.00000 41.10000, 2.20000 41.10000, ...</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>2</th>\n",
+       "      <td>POLYGON ((2.20000 41.10000, 2.40000 41.10000, ...</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>3</th>\n",
+       "      <td>POLYGON ((2.40000 41.10000, 2.60000 41.10000, ...</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>4</th>\n",
+       "      <td>POLYGON ((2.60000 41.10000, 2.80000 41.10000, ...</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>...</th>\n",
+       "      <td>...</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>95</th>\n",
+       "      <td>POLYGON ((2.80000 42.90000, 3.00000 42.90000, ...</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>96</th>\n",
+       "      <td>POLYGON ((3.00000 42.90000, 3.20000 42.90000, ...</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>97</th>\n",
+       "      <td>POLYGON ((3.20000 42.90000, 3.40000 42.90000, ...</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>98</th>\n",
+       "      <td>POLYGON ((3.40000 42.90000, 3.60000 42.90000, ...</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>99</th>\n",
+       "      <td>POLYGON ((3.60000 42.90000, 3.80000 42.90000, ...</td>\n",
+       "    </tr>\n",
+       "  </tbody>\n",
+       "</table>\n",
+       "<p>100 rows × 1 columns</p>\n",
+       "</div>"
+      ],
+      "text/plain": [
+       "                                              geometry\n",
+       "FID                                                   \n",
+       "0    POLYGON ((1.80000 41.10000, 2.00000 41.10000, ...\n",
+       "1    POLYGON ((2.00000 41.10000, 2.20000 41.10000, ...\n",
+       "2    POLYGON ((2.20000 41.10000, 2.40000 41.10000, ...\n",
+       "3    POLYGON ((2.40000 41.10000, 2.60000 41.10000, ...\n",
+       "4    POLYGON ((2.60000 41.10000, 2.80000 41.10000, ...\n",
+       "..                                                 ...\n",
+       "95   POLYGON ((2.80000 42.90000, 3.00000 42.90000, ...\n",
+       "96   POLYGON ((3.00000 42.90000, 3.20000 42.90000, ...\n",
+       "97   POLYGON ((3.20000 42.90000, 3.40000 42.90000, ...\n",
+       "98   POLYGON ((3.40000 42.90000, 3.60000 42.90000, ...\n",
+       "99   POLYGON ((3.60000 42.90000, 3.80000 42.90000, ...\n",
+       "\n",
+       "[100 rows x 1 columns]"
+      ]
+     },
+     "execution_count": 9,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "coordinates.shapefile"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "coordinates.shapefile.to_file('spatial_join_method_2')"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Method 3: Areas intersection"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 12,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/html": [
+       "<div>\n",
+       "<style scoped>\n",
+       "    .dataframe tbody tr th:only-of-type {\n",
+       "        vertical-align: middle;\n",
+       "    }\n",
+       "\n",
+       "    .dataframe tbody tr th {\n",
+       "        vertical-align: top;\n",
+       "    }\n",
+       "\n",
+       "    .dataframe thead th {\n",
+       "        text-align: right;\n",
+       "    }\n",
+       "</style>\n",
+       "<table border=\"1\" class=\"dataframe\">\n",
+       "  <thead>\n",
+       "    <tr style=\"text-align: right;\">\n",
+       "      <th></th>\n",
+       "      <th>geometry</th>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>FID</th>\n",
+       "      <th></th>\n",
+       "    </tr>\n",
+       "  </thead>\n",
+       "  <tbody>\n",
+       "    <tr>\n",
+       "      <th>0</th>\n",
+       "      <td>POLYGON ((1.80000 41.10000, 2.00000 41.10000, ...</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>1</th>\n",
+       "      <td>POLYGON ((2.00000 41.10000, 2.20000 41.10000, ...</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>2</th>\n",
+       "      <td>POLYGON ((2.20000 41.10000, 2.40000 41.10000, ...</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>3</th>\n",
+       "      <td>POLYGON ((2.40000 41.10000, 2.60000 41.10000, ...</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>4</th>\n",
+       "      <td>POLYGON ((2.60000 41.10000, 2.80000 41.10000, ...</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>...</th>\n",
+       "      <td>...</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>9995</th>\n",
+       "      <td>POLYGON ((20.80000 60.90000, 21.00000 60.90000...</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>9996</th>\n",
+       "      <td>POLYGON ((21.00000 60.90000, 21.20000 60.90000...</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>9997</th>\n",
+       "      <td>POLYGON ((21.20000 60.90000, 21.40000 60.90000...</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>9998</th>\n",
+       "      <td>POLYGON ((21.40000 60.90000, 21.60000 60.90000...</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>9999</th>\n",
+       "      <td>POLYGON ((21.60000 60.90000, 21.80000 60.90000...</td>\n",
+       "    </tr>\n",
+       "  </tbody>\n",
+       "</table>\n",
+       "<p>10000 rows × 1 columns</p>\n",
+       "</div>"
+      ],
+      "text/plain": [
+       "                                               geometry\n",
+       "FID                                                    \n",
+       "0     POLYGON ((1.80000 41.10000, 2.00000 41.10000, ...\n",
+       "1     POLYGON ((2.00000 41.10000, 2.20000 41.10000, ...\n",
+       "2     POLYGON ((2.20000 41.10000, 2.40000 41.10000, ...\n",
+       "3     POLYGON ((2.40000 41.10000, 2.60000 41.10000, ...\n",
+       "4     POLYGON ((2.60000 41.10000, 2.80000 41.10000, ...\n",
+       "...                                                 ...\n",
+       "9995  POLYGON ((20.80000 60.90000, 21.00000 60.90000...\n",
+       "9996  POLYGON ((21.00000 60.90000, 21.20000 60.90000...\n",
+       "9997  POLYGON ((21.20000 60.90000, 21.40000 60.90000...\n",
+       "9998  POLYGON ((21.40000 60.90000, 21.60000 60.90000...\n",
+       "9999  POLYGON ((21.60000 60.90000, 21.80000 60.90000...\n",
+       "\n",
+       "[10000 rows x 1 columns]"
+      ]
+     },
+     "execution_count": 12,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "lat_orig = 41.1\n",
+    "lon_orig = 1.8\n",
+    "inc_lat = 0.2\n",
+    "inc_lon = 0.2\n",
+    "n_lat = 100\n",
+    "n_lon = 100\n",
+    "coordinates = create_nes(comm=None, info=False, projection='regular',\n",
+    "                         lat_orig=lat_orig, lon_orig=lon_orig, inc_lat=inc_lat, inc_lon=inc_lon, \n",
+    "                         n_lat=n_lat, n_lon=n_lon)\n",
+    "coordinates.create_shapefile()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 13,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "mask_path = '/esarchive/scratch/avilanova/software/NES/tutorials/data/timezones_2021c/timezones_2021c.shp'\n",
+    "mask = gpd.read_file(mask_path)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 14,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "Rank 000: 9270 intersected areas found\n",
+      "CPU times: user 8min 19s, sys: 9.29 s, total: 8min 29s\n",
+      "Wall time: 8min 30s\n"
+     ]
+    }
+   ],
+   "source": [
+    "%time coordinates.spatial_join(mask, method='intersection')"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 17,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/html": [
+       "<div>\n",
+       "<style scoped>\n",
+       "    .dataframe tbody tr th:only-of-type {\n",
+       "        vertical-align: middle;\n",
+       "    }\n",
+       "\n",
+       "    .dataframe tbody tr th {\n",
+       "        vertical-align: top;\n",
+       "    }\n",
+       "\n",
+       "    .dataframe thead th {\n",
+       "        text-align: right;\n",
+       "    }\n",
+       "</style>\n",
+       "<table border=\"1\" class=\"dataframe\">\n",
+       "  <thead>\n",
+       "    <tr style=\"text-align: right;\">\n",
+       "      <th></th>\n",
+       "      <th>geometry</th>\n",
+       "      <th>tzid</th>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>FID</th>\n",
+       "      <th></th>\n",
+       "      <th></th>\n",
+       "    </tr>\n",
+       "  </thead>\n",
+       "  <tbody>\n",
+       "    <tr>\n",
+       "      <th>0</th>\n",
+       "      <td>POLYGON ((1.80000 41.10000, 2.00000 41.10000, ...</td>\n",
+       "      <td>Europe/Madrid</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>1</th>\n",
+       "      <td>POLYGON ((2.00000 41.10000, 2.20000 41.10000, ...</td>\n",
+       "      <td>Europe/Madrid</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>2</th>\n",
+       "      <td>POLYGON ((2.20000 41.10000, 2.40000 41.10000, ...</td>\n",
+       "      <td>Europe/Madrid</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>3</th>\n",
+       "      <td>POLYGON ((2.40000 41.10000, 2.60000 41.10000, ...</td>\n",
+       "      <td>Europe/Madrid</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>4</th>\n",
+       "      <td>POLYGON ((2.60000 41.10000, 2.80000 41.10000, ...</td>\n",
+       "      <td>NaN</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>...</th>\n",
+       "      <td>...</td>\n",
+       "      <td>...</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>9995</th>\n",
+       "      <td>POLYGON ((20.80000 60.90000, 21.00000 60.90000...</td>\n",
+       "      <td>Europe/Helsinki</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>9996</th>\n",
+       "      <td>POLYGON ((21.00000 60.90000, 21.20000 60.90000...</td>\n",
+       "      <td>Europe/Helsinki</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>9997</th>\n",
+       "      <td>POLYGON ((21.20000 60.90000, 21.40000 60.90000...</td>\n",
+       "      <td>Europe/Helsinki</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>9998</th>\n",
+       "      <td>POLYGON ((21.40000 60.90000, 21.60000 60.90000...</td>\n",
+       "      <td>Europe/Helsinki</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>9999</th>\n",
+       "      <td>POLYGON ((21.60000 60.90000, 21.80000 60.90000...</td>\n",
+       "      <td>Europe/Helsinki</td>\n",
+       "    </tr>\n",
+       "  </tbody>\n",
+       "</table>\n",
+       "<p>10000 rows × 2 columns</p>\n",
+       "</div>"
+      ],
+      "text/plain": [
+       "                                               geometry             tzid\n",
+       "FID                                                                     \n",
+       "0     POLYGON ((1.80000 41.10000, 2.00000 41.10000, ...    Europe/Madrid\n",
+       "1     POLYGON ((2.00000 41.10000, 2.20000 41.10000, ...    Europe/Madrid\n",
+       "2     POLYGON ((2.20000 41.10000, 2.40000 41.10000, ...    Europe/Madrid\n",
+       "3     POLYGON ((2.40000 41.10000, 2.60000 41.10000, ...    Europe/Madrid\n",
+       "4     POLYGON ((2.60000 41.10000, 2.80000 41.10000, ...              NaN\n",
+       "...                                                 ...              ...\n",
+       "9995  POLYGON ((20.80000 60.90000, 21.00000 60.90000...  Europe/Helsinki\n",
+       "9996  POLYGON ((21.00000 60.90000, 21.20000 60.90000...  Europe/Helsinki\n",
+       "9997  POLYGON ((21.20000 60.90000, 21.40000 60.90000...  Europe/Helsinki\n",
+       "9998  POLYGON ((21.40000 60.90000, 21.60000 60.90000...  Europe/Helsinki\n",
+       "9999  POLYGON ((21.60000 60.90000, 21.80000 60.90000...  Europe/Helsinki\n",
+       "\n",
+       "[10000 rows x 2 columns]"
+      ]
+     },
+     "execution_count": 17,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "coordinates.shapefile"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 18,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "coordinates.shapefile.to_file('spatial_join_method_3')"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.7.4"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 4
+}
diff --git a/tutorials/Jupyter_bash_nord3v2.cmd b/tutorials/Jupyter_bash_nord3v2.cmd
index 8138ab9c270ee96ba8c47d17fc94151a77bdaf36..eefd0f7219bcc337065fb5eef61df954e423919e 100644
--- a/tutorials/Jupyter_bash_nord3v2.cmd
+++ b/tutorials/Jupyter_bash_nord3v2.cmd
@@ -1,10 +1,11 @@
 #!/bin/bash
 #SBATCH --ntasks 1
-#SBATCH --time 03:00:00
+#SBATCH --time 02:00:00
 #SBATCH --job-name NES
 #SBATCH --output log_jupyter-notebook-%J.out
 #SBATCH --error log_jupyter-notebook-%J.err
 #SBATCH --exclusive
+#SBATCH --qos debug
 
 # get tunneling info
 XDG_RUNTIME_DIR=""
@@ -32,7 +33,7 @@ module load cfunits/1.8-foss-2019b-Python-3.7.4
 module load filelock/3.7.1-foss-2019b-Python-3.7.4
 module load pyproj/2.5.0-foss-2019b-Python-3.7.4
 module load eccodes-python/0.9.5-foss-2019b-Python-3.7.4
-module load geopandas/0.8.1-foss-2019b-Python-3.7.4
+module load geopandas/0.10.2-foss-2019b-Python-3.7.4
 module load Shapely/1.7.1-foss-2019b-Python-3.7.4
 #module load NES/0.9.0-foss-2019b-Python-3.7.4
 
diff --git a/tutorials/data/NH3_barcelona_2019_csic.csv b/tutorials/data/NH3_barcelona_2019_csic.csv
new file mode 100644
index 0000000000000000000000000000000000000000..0b1fdd941b1d3dc57aea5ca3431abd63be90191c
--- /dev/null
+++ b/tutorials/data/NH3_barcelona_2019_csic.csv
@@ -0,0 +1,13 @@
+Date-hour in,traffic_site,urban_site
+01/01/2019 00:00:00,4.98898814531849,2.55323513371909
+02/01/2019 00:00:00,3.42253492604592,1.55622607730934
+03/01/2019 00:00:00,2.67506499076435,1.68635511757375
+04/01/2019 00:00:00,3.42552223044561,1.97548631953393
+05/01/2019 00:00:00,5.31480909108289,1.11924536797087
+06/01/2019 00:00:00,3.13949500437943,1.62656669999897
+07/01/2019 00:00:00,0,2.22685564636695
+08/01/2019 00:00:00,0,2.4696380137646
+09/01/2019 00:00:00,0,3.72735535650666
+10/01/2019 00:00:00,0,1.53505610127377
+11/01/2019 00:00:00,0,2.51115189477121
+12/01/2019 00:00:00,0,0
diff --git a/tutorials/data/NH3_stations_CSIC.csv b/tutorials/data/NH3_stations_CSIC.csv
new file mode 100644
index 0000000000000000000000000000000000000000..cf9a54d4dd99109a2aee009c61f72ea663f03574
--- /dev/null
+++ b/tutorials/data/NH3_stations_CSIC.csv
@@ -0,0 +1,3 @@
+station,Lon,Lat
+urban_site,2.1151,41.3875
+traffic_site,2.1534,41.3987
diff --git a/tutorials/data/timezones_2021c/timezones_2021c.dbf b/tutorials/data/timezones_2021c/timezones_2021c.dbf
new file mode 100644
index 0000000000000000000000000000000000000000..aed52eac47a3f420d8ecd782863d6145e1a7f300
Binary files /dev/null and b/tutorials/data/timezones_2021c/timezones_2021c.dbf differ
diff --git a/tutorials/data/timezones_2021c/timezones_2021c.prj b/tutorials/data/timezones_2021c/timezones_2021c.prj
new file mode 100644
index 0000000000000000000000000000000000000000..f45cbadf0074d8b7b2669559a93bc50bb95f82d4
--- /dev/null
+++ b/tutorials/data/timezones_2021c/timezones_2021c.prj
@@ -0,0 +1 @@
+GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]]
\ No newline at end of file
diff --git a/tutorials/data/timezones_2021c/timezones_2021c.shp b/tutorials/data/timezones_2021c/timezones_2021c.shp
new file mode 100644
index 0000000000000000000000000000000000000000..5dbd64b964295a7689fb70646caf4db8f0745ddb
Binary files /dev/null and b/tutorials/data/timezones_2021c/timezones_2021c.shp differ
diff --git a/tutorials/data/timezones_2021c/timezones_2021c.shx b/tutorials/data/timezones_2021c/timezones_2021c.shx
new file mode 100644
index 0000000000000000000000000000000000000000..dda8f5adb453c5460d2c39d7741a2c62d71c2153
Binary files /dev/null and b/tutorials/data/timezones_2021c/timezones_2021c.shx differ