V2 environment update

ProFSea-environment.yml

@@ -3,13 +3,13 @@ channels:
   - conda-forge
 dependencies:
   - numpy
-  - pandas
-  - python=3.7.7
-  - scipy
+  - pandas<=2.3.1
+  - python=3.13.5
+  - scipy<=1.16.1
   - pyyaml
   - netcdf4
   - matplotlib
-  - libwebp>=1.3.2
+  - libwebp
   - cartopy
   - iris
   - tqdm

README.md

@@ -11,6 +11,9 @@ Users are advised to read this guide before using any products as it describes t
 
 The ProFSea tool user guide should be cited as: Perks, R. J. (2023). Met Office Projecting Future Sea Level (ProFSea) tool - User Guide (1.1). Met Office. https://doi.org/10.5281/zenodo.10245409
 
+---
+Once the `ProFSea-env` conda environment has been set up, users must also run `pip install -e .` from the root (top-most) directory in order to run the tool.
+
 ## Contributors
 Several people have contributed to the development of the ProFSea tool and User Guide documentation, namely: Rachel Perks, Jacob Cheung, Benjamin Harrison, Katie Hodge, Mathew Palmer, Michael Sanderson, Hamish Steptoe, Jennifer Weeks and Gregory Munday.
 
@@ -20,7 +23,7 @@ This work was supported by the UK Research & Innovation (UKRI) Strategic Priorit
 ## Licence
 ProFSea is licensed under the [Open Government Licence 3.0](https://www.nationalarchives.gov.uk/doc/open-government-licence/version/3/).
 
-If you have any queries about this tool please contact: enquiries@metoffice.gov.uk
+If you have any queries or feedback about this tool please contact Rachel Perks at rachel.perks@metoffice.gov.uk or the Met Office Service Desk at enquiries@metoffice.gov.uk.
 
 <h5 align="center">
 <img src="https://www.metoffice.gov.uk/binaries/content/gallery/metofficegovuk/images/about-us/website/mo_master_black_mono_for_light_backg_rbg.png" width="200" alt="Met Office"> <br>

profsea/slr_pkg/__init__.py

@@ -42,7 +42,11 @@ def extract_dyn_steric_regression(models, df, scenarios):
     :param scenarios: list of RCP scenarios
     """
     # Base directory for CMIP "zos" and "zostoga" data
-    datadir = settings["cmipinfo"]["sealevelbasedir"]
+    if settings["datalocation"] != "":
+        datadir = os.path.join(settings["datalocation"],"cmip5/")
+    else:
+        datadir = settings["cmipinfo"]["sealevelbasedir"]
+
     # Dictionary of CMIP models and experiments
     zos_dict = cmip.zos_dictionary()
 
@@ -88,21 +92,21 @@ def extract_dyn_steric_regression(models, df, scenarios):
                 try:
                     # dynamic sea level (zos)
                     zos_date = zos_dict[model][scenario]['driftcorr']
-                    zos_file = f'{datadir}normalized_zos_Omon_{model}_' \
-                               f'{scenario}_{zos_date}_driftcorr.nc'
+                    zos_file = os.path.join(datadir, f'normalized_zos_Omon_{model}_' \
+                               f'{scenario}_{zos_date}_driftcorr.nc')
                     zos = cubedata.read_zos_cube(zos_file)[0][:, j, i]
                     # --------------------------------------------------------
                     # global mean thermosteric (zostoga)
                     # Extract, and drift-correct CMIP "zostoga" data
                     # Normal (concatenated)
                     zostoga_date = zos_dict[model][scenario]['zostoga']
-                    zostoga_file = f'{datadir}zostoga_Omon_{model}_' \
-                                   f'{scenario}_{zostoga_date}.nc'
+                    zostoga_file = os.path.join(datadir, f'zostoga_Omon_{model}_' \
+                                   f'{scenario}_{zostoga_date}.nc')
                     zostoga_raw = cubedata.read_zos_cube(zostoga_file)[0]
                     # piControl (concatenated)
                     piControl_date = zos_dict[model][scenario]['piControl']
-                    zostoga_pic_file = f'{datadir}zostoga_Omon_' \
-                                       f'{model}_piControl_{piControl_date}.nc'
+                    zostoga_pic_file = os.path.join(datadir, f'zostoga_Omon_' \
+                                       f'{model}_piControl_{piControl_date}.nc')
                     zostoga_pic = cubedata.read_zos_cube(zostoga_pic_file)[0]
 
                     regr = process.Regress('linear')
@@ -374,9 +378,18 @@ def choose_montecarlo_dir():
     """
     end_yr = settings["projection_end_year"]
     if (end_yr >= 2050) & (end_yr <= 2100):
-        mcdir = settings["short_montecarlodir"]
+        if settings["datalocation"] != "":
+            mcdir = os.path.join(settings["datalocation"],
+                                 "monte_carlo_timeseries")
+        else:
+            mcdir = settings["short_montecarlodir"]
+
     elif (end_yr > 2100) & (end_yr <= 2300):
-        mcdir = settings["long_montecarlodir"]
+        if settings["datalocation"] != "":
+            mcdir = os.path.join(settings["datalocation"], "slr")
+        else:
+            mcdir = settings["long_montecarlodir"]       
+
     else:
         raise ValueError('Projection end year must be between 2050 and 2300')
 

profsea/step1_extract_cmip.py

@@ -72,7 +72,11 @@ def extract_ssh_data(cmip_sea):
     :param cmip_sea: variable to distinguish between which CMIP models to use
     :return: CMIP model names, and associated SSH data cubes
     """
-    cmip_dir = settings["cmipinfo"]["sealevelbasedir"]
+    if settings["datalocation"] != "":
+        cmip_dir = os.path.join(settings["datalocation"],"cmip5/")
+    else:
+        cmip_dir = settings["cmipinfo"]["sealevelbasedir"]
+
     # Select CMIP models to use depending on whether location is within a
     # marginal sea
     if cmip_sea == 'all':
@@ -89,8 +93,14 @@ def extract_ssh_data(cmip_sea):
     for model in model_names:
         print(f'Getting data for {model} model')
         cmip_date = cmip_dict[model]['historical']
-        cmip_file = f'{cmip_dir}zos_Omon_{model}_historical_{cmip_date}.nc'
-        cube = cubeutils.loadcube(cmip_file, ncvar='zos')[0]
+        cmip_file = os.path.join(cmip_dir, f'zos_Omon_{model}_historical_{cmip_date}.nc')
+        try:
+            cube = cubeutils.loadcube(cmip_file, ncvar='zos')[0]
+        except IOError:
+            raise FileNotFoundError(os.path.join(cmip_file),
+                                    '- CMIP5 sea level data not found, please' \
+                                    ' check file path')            
+
         cubes.append(cube.slices(['latitude', 'longitude']).next())
 
     return model_names, cubes
@@ -299,9 +309,9 @@ def main():
         print(f'    No lat lon specified - use tide gauge metadata if '
               f'available')
     print(f'User specified science method is: {settings["sciencemethod"]}')
-    if {settings["cmipinfo"]["cmip_sea"]} == {'all'}:
+    if {settings["cmip_sea"]} == {'all'}:
         print('User specified all CMIP models')
-    elif {settings["cmipinfo"]["cmip_sea"]} == {'marginal'}:
+    elif {settings["cmip_sea"]} == {'marginal'}:
         print('User specified CMIP models for marginal seas only')
 
     # Extract site data from station list (e.g. tide gauge location) or
@@ -314,7 +324,7 @@ def main():
 
     # Find the nearest, appropriate ocean point in CMIP models to specified
     # site location
-    cmip_models, ssh_cubes = extract_ssh_data(settings["cmipinfo"]["cmip_sea"])
+    cmip_models, ssh_cubes = extract_ssh_data(settings["cmip_sea"])
     ocean_point_wrapper(df_site_data, cmip_models, ssh_cubes)
 
 

profsea/step2_extract_steric_dyn_regression.py

@@ -20,14 +20,14 @@ def extract_cmip5_steric_dyn_regression(df):
     scenarios = ['rcp26', 'rcp45', 'rcp85']
 
     # Select CMIP5 models to use
-    if settings["cmipinfo"]["cmip_sea"] == 'all':
+    if settings["cmip_sea"] == 'all':
         model_names = models.cmip5_names()
-    elif settings["cmipinfo"]["cmip_sea"] == 'marginal':
+    elif settings["cmip_sea"] == 'marginal':
         model_names = models.cmip5_names_marginal()
     else:
         raise UnboundLocalError(
             'The selected CMIP5 models to use - cmip_sea = ' +
-            f'{settings["cmipinfo"]["cmip_sea"]} - ' +
+            f'{settings["cmip_sea"]} - ' +
             'is not recognised')
 
     # Calculate the regression parameters and plot the results
@@ -51,9 +51,9 @@ def main():
         print(f'    No lat lon specified - use tide gauge metadata if '
               f'available')
     print(f'User specified science method is: {settings["sciencemethod"]}')
-    if {settings["cmipinfo"]["cmip_sea"]} == {'all'}:
+    if {settings["cmip_sea"]} == {'all'}:
         print('User specified all CMIP models')
-    elif {settings["cmipinfo"]["cmip_sea"]} == {'marginal'}:
+    elif {settings["cmip_sea"]} == {'marginal'}:
         print('User specified CMIP models for marginal seas only')
 
     # Extract site data from station list (e.g. tide gauge location) or

profsea/step3_process_regional_sealevel_projections.py

@@ -203,7 +203,14 @@ def calculate_sl_components(mcdir, components, scenario, site_loc, loc_coords,
 
         offset = G_offset * offset_slopes[comp]
 
-        cube = iris.load_cube(os.path.join(mcdir, f'{scenario}_{comp}.nc'))
+        try:
+            cube = iris.load_cube(os.path.join(mcdir, f'{scenario}_{comp}.nc'))
+        except IOError:
+            raise FileNotFoundError(os.path.join(mcdir, 
+                                                f'{scenario}_{comp}.nc'),
+                                    '- monte carlo file not found, please ' \
+                                    'check file path')            
+
         montecarlo_G[cc, :, :] = cube.data[:nyrs, resamples] + offset
 
         if comp == 'exp':
@@ -322,7 +329,13 @@ def create_FP_interpolator(datadir, dfile, method='linear'):
     :param method: interpolation type --> 'linear' or 'nearest'
     :return: 2D Interpolator object
     """
-    cube = iris.load_cube(os.path.join(datadir, dfile))
+    try:
+        cube = iris.load_cube(os.path.join(datadir, dfile))
+    except IOError:
+        raise FileNotFoundError(os.path.join(datadir, dfile),
+                                '- grd fingerprint file not found, ' \
+                                'please check file path')
+
     lon = cube.coord('longitude').points
     lat = cube.coord('latitude').points
 
@@ -395,15 +408,21 @@ def load_CMIP5_slope_coeffs_UK(scenario):
     :return: 1D array of slope coefficients and 1D array of weights
     """
     print('running function load_CMIP5_slope_coeffs_UK')
-    in_zosdir_uk = settings["cmipinfo"]["slopecoeffsuk"]
+    if settings["datalocation"] != "":
+        in_zosdir_uk = os.path.join(settings["datalocation"],
+                                    "uk_cmip_slope_coefficients")
+    else:
+        in_zosdir_uk = settings["cmipinfo"]["slopecoeffsuk"]
+
     filename_uk = f'{scenario}_CMIP5_regress_coeffs_uk_mask_1.pickle'
 
     try:
         with open(os.path.join(in_zosdir_uk, filename_uk), 'rb') as f:
             data = pickle.load(f, encoding='latin1')['uk_mask_1']
     except FileNotFoundError:
-        raise FileNotFoundError(filename_uk,
-                                '- scenario selected does not exist')
+        raise FileNotFoundError(os.path.join(in_zosdir_uk, filename_uk),
+                                '- CMIP5 UK Slope Coefficients not found, ' \
+                                'please check file path')
 
     # Keys are: 'coeffs', 'models', 'weights'
     coeffs = data['coeffs']
@@ -425,21 +444,39 @@ def read_gia_estimates(sci_method, coords):
     print('running function read_gia_estimates')
     # Directories containing GIA data (independent of scenario)
     if sci_method == 'global':
-        gia_file = settings["giaestimates"]["global"]
+        if settings["datalocation"] != "":
+            gia_file = os.path.join(settings["datalocation"], "gia_estimates",
+                                    "global_GIA_interpolators.pickle")
+        else:
+            gia_file = settings["giaestimates"]["global"]
+
     elif sci_method == 'UK':
-        gia_file = settings["giaestimates"]["uk"]
+        if settings["datalocation"] != "":
+            gia_file = os.path.join(settings["datalocation"], "gia_estimates",
+                                    "Bradley_GIA_interpolator.pickle")
+        else:
+            gia_file = settings["giaestimates"]["uk"]
     else:
         raise UnboundLocalError('The selected GIA estimate - ' +
                                 f'{sci_method} - is not available')
 
-    with open(gia_file, "rb") as ifp:
-        GIA_dict = pickle.load(ifp, encoding='latin1')
+    try:
+        with open(gia_file, "rb") as ifp:
+            GIA_dict = pickle.load(ifp, encoding='latin1')
+    except FileNotFoundError:
+        raise FileNotFoundError(os.path.join("gia_estimates", gia_file),
+                                '- gia estimates file not found, please ' \
+                                'check file path')
 
     GIA_vals = []
     lat, lon = coords
+
     # The GIA_dict contains interpolator objects
     for key in list(GIA_dict.keys()):
-        val = GIA_dict[key]([lat, lon])[0]
+        old_grid = GIA_dict[key].grid
+        old_vals = GIA_dict[key].values
+        new_rgi = RegularGridInterpolator(old_grid, old_vals)
+        val = np.squeeze(new_rgi([lat, lon]))
         GIA_vals.append(val)
 
     nGIA = len(GIA_vals)
@@ -461,9 +498,18 @@ def setup_FP_interpolators(components, sci_method):
     print('running function setup_FP_interpolators')
 
     # Directories for the Slangen, Spada and Klemann fingerprints
-    slangendir = settings["fingerprints"]["slangendir"]
-    spadadir = settings["fingerprints"]["spadadir"]
-    klemanndir = settings["fingerprints"]["klemanndir"]
+    if settings["datalocation"] != "":
+        slangendir = os.path.join(settings["datalocation"],
+                                  "grd_fingerprints")
+        spadadir = os.path.join(settings["datalocation"],
+                                "grd_fingerprints")
+        klemanndir = os.path.join(settings["datalocation"],
+                                  "grd_fingerprints")
+    else:
+        slangendir = settings["fingerprints"]["slangendir"]
+        spadadir = settings["fingerprints"]["spadadir"]
+        klemanndir = settings["fingerprints"]["klemanndir"]
+
 
     # Create empty dictionaries for the Slangen, Spada and Klemann fingerprints
     # interpolator objects.
@@ -516,9 +562,9 @@ def main():
         print(f'    No lat lon specified - use tide gauge metadata if '
               f'available')
     print(f'User specified science method is: {settings["sciencemethod"]}')
-    if {settings["cmipinfo"]["cmip_sea"]} == {'all'}:
+    if {settings["cmip_sea"]} == {'all'}:
         print('User specified all CMIP models')
-    elif {settings["cmipinfo"]["cmip_sea"]} == {'marginal'}:
+    elif {settings["cmip_sea"]} == {'marginal'}:
         print('User specified CMIP models for marginal seas only')
 
     print(f'\nProjecting out to: {settings["projection_end_year"]}\n')

profsea/step4_plot_regional_sealevel.py

@@ -912,9 +912,9 @@ def main():
         print(f'    No lat lon specified - use tide gauge metadata if '
               f'available')
     print(f'User specified science method is: {settings["sciencemethod"]}')
-    if {settings["cmipinfo"]["cmip_sea"]} == {'all'}:
+    if {settings["cmip_sea"]} == {'all'}:
         print('User specified all CMIP models')
-    elif {settings["cmipinfo"]["cmip_sea"]} == {'marginal'}:
+    elif {settings["cmip_sea"]} == {'marginal'}:
         print('User specified CMIP models for marginal seas only')
 
     # Extract site data from station list (e.g. tide gauge location) or

profsea/tide_gauge_locations.py

@@ -95,7 +95,7 @@ def extract_site_info(data_source, data_type, region, site_name, latlon):
                             keep='first')].transpose()
 
             df_temp = pd.DataFrame(site_data).transpose()
-            dfObj = dfObj.append(df_temp)
+            dfObj = dfObj._append(df_temp)
 
         elif latlon != [[]]:
             print(f'{site_to_check} - Site metadata taken from user input')
@@ -110,7 +110,7 @@ def extract_site_info(data_source, data_type, region, site_name, latlon):
                 'Location', 'Dataset type', 'Station ID', 'Latitude',
                 'Longitude'])
             df_temp = df_temp.set_index('Location')
-            dfObj = dfObj.append(df_temp)
+            dfObj = dfObj._append(df_temp)
 
         else:
             raise IndexError(f'No site metadata for this site: {site_to_check}, have you spelled it correctly?')