diff --git a/README.md b/README.md
index 937266f..1345f98 100644
--- a/README.md
+++ b/README.md
@@ -36,13 +36,12 @@ This tool is developed as part of a [Python Toolbox](https://pro.arcgis.com/en/p
 ## Getting Started
 
 ### Prerequisites
-The following is a list of the programs and libraries used in the tool, with their respective versions:
 
-* ArcGIS Pro - ```arcpy``` (version 3.1)
-* ```pandas``` (version 1.4)
-* ```numpy``` (version 1.20)
-* ```shapely``` (version 2.0)
-* ```statsmodels``` (version 0.13)
+Check that you have installed all the required libraries used in the toolbox. All packages with their tested versions are listed in [requirements.txt](https://github.com/AlbertGallegoJimenez/shoreline-evolution-tool/tree/main/requirements.txt). Note that when working with a cloned version of the ArcGIS ```anaconda``` environment there are already pre-installed libraries, those that are not and need to be installed manually are the following:
+
+* ```shapely```
+* ```statsmodels``` (from ArcGIS version 3.2 this package is included in the base ArcGIS ```anaconda``` environment)
+* ```cartopy```
 
 In terms of data, this tool relies on the use of the following two files:
 * <ins>**Baseline**</ins> (Vector - Polyline). This is the reference line used to assess the evolution of the coastal stretch. It can be digitized manually by the user with the help of a background orthophoto, it is recommended to place the baseline **inland**. The baseline must capture the general orientation of the coast.
@@ -51,8 +50,8 @@ In terms of data, this tool relies on the use of the following two files:
 
 ### Installation
 
-0. Make sure you have cloned the base ArcGIS' ```anaconda``` environment so you can install more packages. More info [here](https://pro.arcgis.com/en/pro-app/latest/arcpy/get-started/clone-an-environment.htm).
-1. Install both ```shapely``` and ```statsmodels``` packages.
+0. Make sure you have cloned the base ArcGIS ```anaconda``` environment so you can install more packages. More info [here](https://pro.arcgis.com/en/pro-app/latest/arcpy/get-started/clone-an-environment.htm).
+1. Install ```shapely```, ```statsmodels``` and ```cartopy``` packages if you do not have it installed yet.
 2. Download the content in the [src](https://github.com/AlbertGallegoJimenez/shoreline-evolution-tool/tree/main/src) folder.
 3. Open the Catalog Pane in ArcGIS Pro and open the downloaded Toolbox (.pyt) to see the tools.
 <div align="center">
diff --git a/requirements.txt b/requirements.txt
new file mode 100644
index 0000000..4e34f5d
--- /dev/null
+++ b/requirements.txt
@@ -0,0 +1,11 @@
+# Name                    Version
+arcgispro (software)      3.1
+arcpy                     3.1
+cartopy                   0.21.1
+matplotlib                3.6.0
+numpy                     1.20.1
+pandas                    1.4.4
+regex                     2022.7.9
+seaborn                   0.12.1
+shapely                   2.0.1
+statsmodels               0.13.5
\ No newline at end of file
diff --git a/src/tools/correctTransects.py b/src/tools/correctTransects.py
index de24aa9..c718fd9 100644
--- a/src/tools/correctTransects.py
+++ b/src/tools/correctTransects.py
@@ -74,7 +74,6 @@ def execute(self, parameters, messages):
 
         # Add Bearing attribute
         arcpy.management.CalculateGeometryAttributes(transectsFeature, "Bearing LINE_BEARING")
-
         cursor = arcpy.da.SearchCursor(transectsFeature, ["transect_id", "Bearing"])
         df = pd.DataFrame(data=[row for row in cursor], columns=["transect_id", "Bearing"])
 
@@ -87,7 +86,7 @@ def execute(self, parameters, messages):
         # Make the changes in the feature class
         RotateFeatures(df, transectsFeature)
 
-        # Recalulate the bearing with the transects inverted
+        # Recalculate the bearing with the transects inverted
         arcpy.management.CalculateGeometryAttributes(transectsFeature, "Bearing LINE_BEARING")
         cursor = arcpy.da.SearchCursor(transectsFeature, ["transect_id", "Bearing"])
         df = pd.DataFrame(data=[row for row in cursor], columns=["transect_id", "Bearing"])
@@ -111,6 +110,7 @@ def execute(self, parameters, messages):
     def postExecute(self, parameters):
         """This method takes place after outputs are processed and
         added to the display."""
+        # Update the labels of the transect feature
         aprx = arcpy.mp.ArcGISProject('CURRENT')
         aprxMap = aprx.activeMap
         transectsFeature = aprxMap.listLayers(os.path.basename(parameters[0].valueAsText))[0]
diff --git a/src/tools/generateTransects.py b/src/tools/generateTransects.py
index ca09340..ed2f7aa 100644
--- a/src/tools/generateTransects.py
+++ b/src/tools/generateTransects.py
@@ -108,6 +108,7 @@ def execute(self, parameters, messages):
     def postExecute(self, parameters):
         """This method takes place after outputs are processed and
         added to the display."""
+        # Update the labels of the transect feature
         aprx = arcpy.mp.ArcGISProject('CURRENT')
         aprxMap = aprx.activeMap
         transectsFeature = aprxMap.listLayers(os.path.basename(parameters[0].valueAsText))[0]
diff --git a/src/tools/performAnalysis.py b/src/tools/performAnalysis.py
index fdfb923..432ad2e 100644
--- a/src/tools/performAnalysis.py
+++ b/src/tools/performAnalysis.py
@@ -81,11 +81,9 @@ def execute(self, parameters, messages):
         for field in metrics_fields:
             arcpy.management.AddField(transectsFeature, field, "DOUBLE")
 
-        count = 0
         with arcpy.da.UpdateCursor(transectsFeature, metrics_fields) as cursor:
-            for row in cursor:
-                cursor.updateRow(shore_metrics.loc[count, metrics_fields].tolist())
-                count += 1
+            for i, _ in enumerate(cursor):
+                cursor.updateRow(shore_metrics.loc[i, metrics_fields].tolist())
 
         return
 
@@ -102,8 +100,6 @@ def postExecute(self, parameters):
 
         # Get the current symbology settings of the layer
         sym = transectsLayerObj.symbology
-        # Get the layer's CIM definition
-        cim = transectsLayerObj.getDefinition('V3')
 
         # Set the renderer to Graduated Colors
         sym.updateRenderer('GraduatedColorsRenderer')
@@ -113,7 +109,6 @@ def postExecute(self, parameters):
 
         # Specify the field used for classification
         sym.renderer.classificationField = 'LRR'
-        cim.renderer.heading = 'LRR (m/year)'
 
         # Set the number of class breaks
         sym.renderer.breakCount = 6
@@ -125,7 +120,7 @@ def postExecute(self, parameters):
         labels = ["-50.0 - -4.0", "-4.0 - -2.0", "-2.0 - 0.0", "0.0 - 2.0", "2.0 - 4.0", "4.0 - 50.0"]
 
         # Define sizes for each class
-        sizes = [4, 2, 1, 1, 2, 4]
+        sizes = [6, 3, 1.5, 1.5, 3, 6]
 
         # Update values for each class
         for i, brk in enumerate(sym.renderer.classBreaks):
@@ -135,6 +130,21 @@ def postExecute(self, parameters):
 
         # Apply the updated symbology settings to the layer
         transectsLayerObj.symbology = sym
-        transectsLayerObj.setDefinition(cim)
 
+        # Get the layer's CIM definition
+        cim = transectsLayerObj.getDefinition('V3')
+
+        # Set the label for the symbolised field
+        cim.renderer.heading = 'LRR (m/year)'
+
+        # Exclude non-significant transects and apply unique symbology
+        cim.renderer.useExclusionSymbol = True
+        cim.renderer.exclusionClause = 'Pvalue > 0.05'
+        cim.renderer.exclusionLabel = 'Non-significant transect'
+        cim.renderer.exclusionSymbol.symbol.symbolLayers[0].color.values = [130, 130, 130, 100]
+        cim.renderer.exclusionSymbol.symbol.symbolLayers[0].width = 1.5
+        
+        # Update the CIM
+        transectsLayerObj.setDefinition(cim)
+        
         return
\ No newline at end of file
diff --git a/src/tools/plotResults.py b/src/tools/plotResults.py
new file mode 100644
index 0000000..b695d67
--- /dev/null
+++ b/src/tools/plotResults.py
@@ -0,0 +1,122 @@
+import arcpy
+from tools.utils.plot_results import PlottingUtils
+
+class PlotResults(object):
+    def __init__(self):
+        """Define the tool (tool name is the name of the class)."""
+        self.label = "5. Plot The Analysis Results"
+        self.description = ""
+        self.canRunInBackground = False
+
+    def getParameterInfo(self):
+        """Define parameter definitions"""
+
+        shoreline_param = arcpy.Parameter(
+            displayName="Input Shorelines Intersection Points Feature",
+            name="shore_features",
+            datatype="GPFeatureLayer",
+            parameterType="Required",
+            direction="Input")
+        shoreline_param.filter.list = ["Point"]
+
+        transects_param = arcpy.Parameter(
+            displayName="Input Transects Feature",
+            name="transects_features",
+            datatype="GPFeatureLayer",
+            parameterType="Required",
+            direction="Input")
+        transects_param.filter.list = ["Polyline"]
+
+        transects_ID_2plot_param = arcpy.Parameter(
+            displayName="Input Transects ID To Plot",
+            name="transects_ID_2plot",
+            datatype="GPValueTable",
+            parameterType="Required",
+            direction="Input")
+        transects_ID_2plot_param.columns = [['GPLong', 'Transects ID']]
+
+        parameters = [shoreline_param, transects_param, transects_ID_2plot_param]
+
+        return parameters
+
+    def isLicensed(self):
+        """Set whether tool is licensed to execute."""
+        return True
+
+    def updateParameters(self, parameters):
+        """Modify the values and properties of parameters before internal
+        validation is performed.  This method is called whenever a parameter
+        has been changed."""
+        return
+
+    def updateMessages(self, parameters):
+        """Modify the messages created by internal validation for each tool
+        parameter.  This method is called after internal validation."""
+        parameters[2].setWarningMessage(
+            "Check significance before selecting transects (Pvalue < 0.05).")
+        return
+
+    def execute(self, parameters, messages):
+        """The source code of the tool."""
+        shoreFeatures = parameters[0].valueAsText
+        transectsFeature = parameters[1].valueAsText
+        transectsID_2plot = parameters[2].value
+        transectsID_2plot = [id[0] for id in transectsID_2plot] # As the parameter value is passed as a list of lists, instead of a list of integers
+
+        # Check for errors in the transect IDs selected
+        self._check_transects_id(transectsID_2plot, transectsFeature)
+
+        # Initialize the class
+        plotter = PlottingUtils(transects=transectsFeature,
+                                shore_intersections=shoreFeatures)
+        
+        # Plot the spatial evolution
+        plotter.plot_spatial_evolution()
+
+        # Plot the time series for the transects selected
+        plotter.plot_time_series(transects2plot=transectsID_2plot)
+
+        # Plot the seasonality for the transects selected. Set a minimum of 2 years of data to plot the seasonality.
+        if plotter.shore_intersections_df['date'].dt.year.max() - plotter.shore_intersections_df['date'].dt.year.min() >= 2:
+            plotter.plot_seasonality(transects2plot=transectsID_2plot)
+
+        # Plot the LRR map
+        plotter.plot_map('LRR')
+
+        # Plot the SCE map
+        plotter.plot_map('SCE')
+
+        # Plot the NSM map
+        plotter.plot_map('NSM')
+
+        return
+
+    def postExecute(self, parameters):
+        """This method takes place after outputs are processed and
+        added to the display."""
+        return
+    
+    def _check_transects_id(self, transectsID_2plot, transectsFeature):
+        """This private method checks that all transects 
+        passed by the user are valid."""
+
+        # Check that all IDs are numerical. This code may be unnecessary since the parameter is defined as 'GPLong' and this means that ArcGIS will only accept numerical values.
+        if not all(isinstance(id, (int, float)) for id in transectsID_2plot):
+            arcpy.AddError('Invalid transect ID, check that all IDs are numeric.')
+            raise Exception('Invalid transect ID, check that all IDs are numeric.')
+        
+        # Check for duplicates
+        if len(transectsID_2plot) != len(set(transectsID_2plot)):
+            arcpy.AddError('Some transects are repeated.')
+            raise Exception('Some transects are repeated.')
+        
+        # Check if the user only has selected one transect
+        if len(transectsID_2plot) == 1:
+            arcpy.AddError('Please, select more than one transect to plot.')
+            raise Exception('Please, select more than one transect to plot.')
+        
+        # Check for ID out of range
+        list_transect_id = [row[0] for row in arcpy.da.SearchCursor(transectsFeature, ['transect_id'])]
+        if set(transectsID_2plot) - set(list_transect_id):
+            arcpy.AddError('Invalid transect ID, check that all IDs are within the valid range.')
+            raise Exception('Invalid transect ID, check that all IDs are within the valid range.')
diff --git a/src/tools/utils/__pycache__/intersect_lines.cpython-39.pyc b/src/tools/utils/__pycache__/intersect_lines.cpython-39.pyc
index cf46271..11a22c9 100644
Binary files a/src/tools/utils/__pycache__/intersect_lines.cpython-39.pyc and b/src/tools/utils/__pycache__/intersect_lines.cpython-39.pyc differ
diff --git a/src/tools/utils/__pycache__/plot_results.cpython-39.pyc b/src/tools/utils/__pycache__/plot_results.cpython-39.pyc
new file mode 100644
index 0000000..f0f12c5
Binary files /dev/null and b/src/tools/utils/__pycache__/plot_results.cpython-39.pyc differ
diff --git a/src/tools/utils/__pycache__/shoreline_evolution.cpython-39.pyc b/src/tools/utils/__pycache__/shoreline_evolution.cpython-39.pyc
index 7b91ffa..e230549 100644
Binary files a/src/tools/utils/__pycache__/shoreline_evolution.cpython-39.pyc and b/src/tools/utils/__pycache__/shoreline_evolution.cpython-39.pyc differ
diff --git a/src/tools/utils/__pycache__/transect_processor.cpython-39.pyc b/src/tools/utils/__pycache__/transect_processor.cpython-39.pyc
index b4739e0..3dbc1ea 100644
Binary files a/src/tools/utils/__pycache__/transect_processor.cpython-39.pyc and b/src/tools/utils/__pycache__/transect_processor.cpython-39.pyc differ
diff --git a/src/tools/utils/plot_results.py b/src/tools/utils/plot_results.py
new file mode 100644
index 0000000..4b60cd1
--- /dev/null
+++ b/src/tools/utils/plot_results.py
@@ -0,0 +1,316 @@
+import arcpy
+import seaborn as sns
+import matplotlib.pyplot as plt
+import numpy as np
+import pandas as pd
+import os
+import re
+import cartopy.crs as ccrs
+from tools.utils.intersect_lines import *
+from matplotlib.colors import Normalize, TwoSlopeNorm
+from matplotlib.cm import ScalarMappable
+import matplotlib.lines as mlines
+from matplotlib.gridspec import GridSpec
+
+# Define the settings for the plots
+plt.style.use('classic')
+plt.rcParams['axes.grid'] = True
+plt.rcParams['grid.linestyle'] = '--'
+plt.rcParams['grid.linewidth'] = 0.5
+plt.rcParams['axes.prop_cycle'] = plt.cycler(color=['#5B9CFD', '#FECB31', '#bfbbd9', '#fa8174', '#81b1d2',
+                                                    '#fdb462', '#b3de69', '#bc82bd', '#ccebc4', '#ffed6f'])
+
+class PlottingUtils():
+    def __init__(self, transects, shore_intersections):
+        """
+        Constructor method for initializing PlottingUtils class.
+
+        Parameters:
+            transects (arcpy.FeatureClass): Feature Class object for transects (Line geometries).
+            shore_intersections (arcpy.FeatureClass): Feature Class object for shoreline intersections (Point geometries).
+        """
+        self.transects = transects
+
+        # Set the directory where plots will be stored
+        aprx = arcpy.mp.ArcGISProject('CURRENT')
+        self.out_dir = os.path.join(aprx.homeFolder, 'Plots results')
+        if not os.path.exists(self.out_dir):
+            os.mkdir(self.out_dir)
+
+        # Convert Feature Classes to Pandas DataFrames
+        self.transects_df = self._feature_class_to_dataframe(transects)
+        self.shore_intersections_df = self._feature_class_to_dataframe(shore_intersections)
+        self.shore_intersections_df['date'] = pd.to_datetime(self.shore_intersections_df['date'])
+
+        # Create Shapely Object for transects geometries
+        self.transects_shapely = line_arcgis2shapely(feature=transects, id='transect_id')
+    
+
+    def _feature_class_to_dataframe(self, feature_class):
+        """
+        Private method to convert an ArcGIS Feature Class to a Pandas DataFrame.
+
+        Parameters:
+            feature_class (arcpy.FeatureClass): Feature Class object.
+
+        Returns:
+            pd.DataFrame: Pandas DataFrame containing feature class data.
+        """
+        fields = [field.name for field in arcpy.ListFields(feature_class)]
+
+        return pd.DataFrame([row for row in arcpy.da.SearchCursor(feature_class, fields)], columns=fields)
+    
+
+    def _get_UTM_projection(self):
+        """
+        Private method to get UTM projection from transects Feature Class.
+
+        Returns:
+            UTM_number (int): UTM number.
+            southern_hemisphere (bool): Whether or not the Feature Class is in the southern hemisphere.
+        """
+        # CRS string (e.g. 'WGS_1984_UTM_Zone_19N')
+        crs = arcpy.Describe(self.transects).SpatialReference.name
+
+        # Regex codes to extract the UTM number and the hemisphere
+        num_code = r'\b(\d+)[A-Za-z]?\b'
+        hemisphere_code = r'\b\d+([A-Za-z])\b'
+
+        UTM_number = int(re.findall(num_code, crs.split('_')[-1])[0])
+        hemisphere_UTM = re.findall(hemisphere_code, crs.split('_')[-1])[0]
+
+        if hemisphere_UTM == 'N':
+            southern_hemisphere = False
+        else:
+            southern_hemisphere = True
+        
+        return UTM_number, southern_hemisphere
+    
+
+    def _set_map_configuration(self, metric):
+        """
+        Private method to set the map configuration for plotting LRR, SCE and NSM.
+        That is, the colormap, the type of normalization scale and the type of colorbar according to the metric.
+
+        Parameters:
+            metric (string): Name of the feature to plot.
+
+        Returns:
+            cmap (matplotlib.colors.LinearSegmentedColormap): The colormap.
+            norm (matplotlib.colors.TwoSlopeNorm or matplotlib.colors.Normalize): The type of the normalization.
+            extend_cbar (string): The type of the colorbar according to cmap and norm.
+        """
+
+        metric_min, metric_max = self.transects_df[metric].describe()[['min', 'max']]
+
+        if metric_min < 0 and metric_max > 0:
+            cmap = plt.get_cmap('RdBu')
+            norm = TwoSlopeNorm(vmin=metric_min, vcenter=0, vmax=metric_max)
+            extend_cbar = 'both'
+        elif metric_min < 0 and metric_max < 0:
+            cmap = plt.get_cmap('Reds')
+            norm = Normalize(vmin=metric_min, vmax=0)
+            extend_cbar = 'min'
+        elif metric_min > 0 and metric_max > 0:
+            cmap = plt.get_cmap('Blues')
+            norm = Normalize(vmin=0, vmax=metric_max)
+            extend_cbar = 'max'
+        
+        return cmap, norm, extend_cbar
+    
+     #  ===== FROM HERE, ALL THE FUNCTIONS TO CREATE THE PLOTS ARE DEFINED =====
+  
+    def plot_spatial_evolution(self):
+        # Plot the distances between all shorelines and the baseline on each transect
+        
+        fig, ax = plt.subplots(figsize=(12, 5))
+
+        # Plot all values as scatter
+        ax.scatter(self.shore_intersections_df['transect_id'],
+                   self.shore_intersections_df['distance_from_base'],
+                   alpha=.1, lw=0, zorder=1)
+
+        # Plot average width value for each transect
+        ax.plot(self.shore_intersections_df['transect_id'].unique(),
+                self.shore_intersections_df.groupby('transect_id')['distance_from_base'].mean(),
+                label='avg', color='#fa8174', lw=2, zorder=2)
+        
+        # Plot mean +-2*std width value for each transect
+        ax.fill_between(self.shore_intersections_df['transect_id'].unique(),
+                        self.shore_intersections_df.groupby('transect_id')['distance_from_base'].mean()+\
+                        2 * self.shore_intersections_df.groupby('transect_id')['distance_from_base'].std(),
+                        self.shore_intersections_df.groupby('transect_id')['distance_from_base'].mean()-\
+                        2 * self.shore_intersections_df.groupby('transect_id')['distance_from_base'].std(),
+                        color='#bfbbd9', alpha=.5, lw=0, label='avg\u00B12std', zorder=0)
+        
+        # Plot settings
+        ax.set_xticks((np.arange(0, max(self.shore_intersections_df['transect_id']) + 2, 2)).tolist())
+        ax.set_xlabel('transect_id')
+        ax.set_ylabel('distance from baseline (m)')
+        ax.set_xlim([-1, max(self.shore_intersections_df['transect_id']) + 2])
+        plt.text(-0.05, 0.9, 'seaward', transform=plt.gca().transAxes, horizontalalignment='right', fontstyle='italic')
+        plt.text(-0.05, .1, 'landward', transform=plt.gca().transAxes, horizontalalignment='right', fontstyle='italic')
+        plt.grid(linestyle='--', alpha=0.3)
+        ax.legend()
+        ax.set_title('Spatial shoreline evolution (%.f - %.f)' % (self.shore_intersections_df['date'].min().year,
+                                                                  self.shore_intersections_df['date'].max().year))
+        fig.savefig(os.path.join(self.out_dir, 'Spatial shoreline evolution.png'), dpi=300, bbox_inches='tight')
+
+
+    def plot_time_series(self, transects2plot):
+        # Plot time series for selected transects
+
+        fig = plt.figure(figsize=(12, 2 * len(transects2plot)))
+        gs = GridSpec(len(transects2plot), 7, figure=fig)
+
+        for i, t in enumerate(transects2plot):
+            ax = fig.add_subplot(gs[i, :-1])
+
+            # Prepare the data to plot
+            data_transect = self.shore_intersections_df.loc[self.shore_intersections_df['transect_id'] == t, :]
+            data_transect.index = pd.to_datetime(data_transect['date'])
+            data_transect = data_transect.sort_index()
+            data_transect['Time'] = np.arange(len(data_transect.index))
+            X = data_transect.index.map(pd.Timestamp.toordinal)  # features
+            y = data_transect.loc[:, 'distance_from_base']  # target
+
+            # Plot time series
+            #y.plot(ax=ax, color='#5B9CFD', linestyle='-', marker='o', markersize=2, label='shoreline positions')
+            ax.plot(X, y.values, linestyle='-', marker='o', markersize=2, label='shoreline positions')
+            sns.regplot(x=X, y=y, ci=95, scatter=False, label='linear regression fit', ax=ax)
+
+            # Plot settings
+            ax.set_ylabel('distance from\nbaseline (m)')
+            ax.locator_params(axis='y', nbins=4)
+            ax.set_title('transect ' + str(t))
+            ax.grid(alpha=0.3)
+
+            # Plot LRR error plot
+            ax2 = fig.add_subplot(gs[i, -1:])
+            lrr = self.transects_df.loc[self.transects_df['transect_id']==t, 'LRR']
+            lci = self.transects_df.loc[self.transects_df['transect_id']==t, ['LCI_low', 'LCI_upp']].to_numpy()[0]
+            ax2.errorbar(0, lrr, yerr=([abs(lci[0] - lrr.values[0])], [lci[1] - lrr.values[0]]),
+                        fmt='or', markersize=8, capsize=5)
+            # Plot settings
+            ax2.set_xticklabels([])
+            ax2.locator_params(axis='y', nbins=4)
+            ax2.locator_params(axis='x', nbins=1)
+            ax2.grid(axis='y', alpha=0.3)
+            ax2.grid(axis='x', alpha=0)
+            
+            if i == len(transects2plot) - 1:
+                # Time series plot
+                ax.set_xlabel('')
+                xticks = ax.get_xticks()
+                labels = [pd.Timestamp.fromordinal(int(label)).date().strftime("%m-%Y") for label in xticks]
+                ax.set_xticks(xticks)
+                ax.set_xticklabels(labels)
+                ax.legend(fontsize=8)
+                
+                # LRR error plot
+                ax2.set_xlabel('LRR\u00B195%\n(m/year)')
+                
+            else:
+                # Time series plot
+                ax.set_xlabel('')
+                ax.set_xticklabels([])
+            
+            plt.subplots_adjust(hspace=0.4, wspace=1)
+        
+        fig.savefig(os.path.join(self.out_dir, 'Time shoreline evolution.png'), bbox_inches='tight', dpi=300)
+
+
+    def plot_seasonality(self, transects2plot):
+        # Boxplot by months for selected transects
+
+        shore_month = self.shore_intersections_df.copy()
+        shore_month['month'] = shore_month['date'].dt.month
+
+        fig, ax = plt.subplots(len(transects2plot), 1, figsize=(10, 15), sharex=True)
+
+        for i, t in enumerate(transects2plot):
+            # Grab the data
+            data = shore_month.loc[shore_month['transect_id']==t, ['month', 'distance_from_base']]
+            
+            # Median values Line plot
+            ax[i].plot(data.groupby('month')['distance_from_base'].median(),
+                    color='#FECB31', marker='o', markeredgecolor=None, alpha=.7)
+            
+            # Box plot
+            medianprops = dict(linestyle='-', linewidth=2, color='#FECB31')
+            whiskerprops = dict(linestyle='-', color='#5B9CFD')
+            data.boxplot(column='distance_from_base', by='month', ax=ax[i], medianprops=medianprops, whiskerprops=whiskerprops)
+            
+            # Plot settings
+            ax[i].set_xlabel('')
+            ax[i].set_ylabel('distance from\nbaseline (m)')
+            ax[i].locator_params(axis='y', nbins=5)
+            ax[i].set_title('transect ' + str(t))
+            ax[i].grid(linestyle='--', alpha=0.3)
+            ax[i].grid(axis='x', alpha=0)
+            
+            if i == len(transects2plot) - 1:
+                ax[i].set_xticklabels(['Jan', 'Feb', 'Mar', 'Apr', 'May', 'Jun', 'Jul', 'Aug', 'Sep', 'Oct', 'Nov', 'Dec'] * len(transects2plot))
+                
+        fig.suptitle(None)
+        plt.subplots_adjust(hspace=0.2)
+        fig.savefig(os.path.join(self.out_dir, 'Shoreline evolution seasonality.png'), dpi=300, bbox_inches='tight')
+
+
+    def plot_map(self, metric):
+        # Plot a map of the selected metric (LRR, SCE or NSM)
+
+        # Set the cmap, norm and the type of cbar of the plot
+        cmap, norm, extend_cbar = self._set_map_configuration(metric)
+        
+        # Set projection parameter
+        UTM_number, southern_hemisphere = self._get_UTM_projection()
+        proj = ccrs.UTM(UTM_number, southern_hemisphere=southern_hemisphere)
+
+        # Create a subplot with UTM projection
+        fig, ax = plt.subplots(layout='compressed', subplot_kw={'projection':proj})
+
+        # Iterate through transects and plot lines
+        for i, t in self.transects_shapely.items():
+            
+            # Check if Pvalue is less than 0.05 for significance
+            if self.transects_df.loc[self.transects_df['transect_id'] == i, 'Pvalue'].values <= 0.05:
+                color = cmap(norm(self.transects_df.loc[self.transects_df['transect_id'] == i, metric]))
+                ls = '-'
+            else:
+                color = 'gray'
+                ls = '--'
+            
+            # Plot the transect line
+            ax.plot(*t.xy,
+                    color=color,
+                    transform=proj,
+                    lw=3,
+                    ls=ls)
+
+        # Create a legend entry for non-significant transects
+        legend_entry = mlines.Line2D([], [], color='gray', lw=2, ls='--', label='Non-significant transect')
+
+        # Customize gridlines, ticks, and appearance
+        ax.gridlines(crs=proj, linewidth=1, color='black', alpha=0, linestyle="--")
+        ax.set_xticks(ax.get_xticks(), crs=proj)
+        ax.set_yticks(ax.get_yticks(), crs=proj)
+        ax.grid(alpha=.3)
+
+        # Add colorbar and set title
+        if self.transects_df['Pvalue'].min() <= 0.05:
+            fig.colorbar(ScalarMappable(norm=norm, cmap=cmap), extend=extend_cbar, ax=ax)
+        if metric == 'LRR':
+            ax.set_title('Linear Regression Rate, LRR (m/year)', y=1.05)
+        elif metric == 'SCE':
+            ax.set_title('Shoreline Change Envelope, SCE (m)', y=1.05)
+        elif metric == 'NSM':
+            ax.set_title('Net Shoreline Movement, NSM (m)', y=1.05)
+
+        # Set limits, labels, legend, and save the figure
+        lons = [x for t in self.transects_shapely.values() for x in t.xy[0]] # Extract the longitudes for plotting purposes (Sometimes the plot is centered to the west and a blank space is left to the east of the study area)
+        ax.set_xlim([ax.get_xlim()[0], max(lons)])
+        ax.set_xlabel('Eastings (m)')
+        ax.set_ylabel('Northings (m)')
+        ax.legend(handles=[legend_entry], fontsize='small')
+        fig.savefig(os.path.join(self.out_dir, '{0}_transects.png'.format(metric)), dpi=300, bbox_inches='tight')
diff --git a/src/tools/utils/transect_processor.py b/src/tools/utils/transect_processor.py
index d8de6da..ae487d1 100644
--- a/src/tools/utils/transect_processor.py
+++ b/src/tools/utils/transect_processor.py
@@ -17,6 +17,7 @@ def invert_angles(self):
         # Detect the first transects that are inverted (if there are more than one change). To do so, find a difference angle around ~180º (it has been selected a range of 180 +-50).
         start_change_transects = self.df[(self.df['diffBear'].abs() >= 130) & (self.df['diffBear'].abs() < 230)]['transect_id'].to_list()
 
+        # Empty list of the transects that need to be inverted
         transects2correct = []
         for i, _ in enumerate(start_change_transects):
             if (i + 1) % 2 != 0: # Odd number
@@ -30,7 +31,7 @@ def invert_angles(self):
                 pass
         
         self.df['Angle'] = 0
-        self.df.loc[self.df['transect_id'].isin(transects2correct), 'Angle'] = 180 # Rotate 180 degrees the bearing anle
+        self.df.loc[self.df['transect_id'].isin(transects2correct), 'Angle'] = 180 # Rotate 180 degrees the bearing angle
 
     def classify_transects(self):
         # Classify transects with large differences using the corrFactor value. In addition, try not to take into account the differences in the 360-0 sector.
@@ -49,11 +50,9 @@ def __init__(self, df, fclass):
         # Initialize the class by adding an angle field with the values calculated above
         arcpy.management.AddField(fclass, 'Angle', 'DOUBLE')
         
-        count = 0
         with arcpy.da.UpdateCursor(fclass, 'Angle') as cursor:
-            for row in cursor:
-                cursor.updateRow([df.loc[count, 'Angle']])
-                count += 1
+            for i, row in enumerate(cursor):
+                cursor.updateRow([df.loc[i, 'Angle']])
 
         # Rebuild each polyine with rotated vertices
         with arcpy.da.UpdateCursor(fclass, ['SHAPE@', 'Angle']) as cursor: