Merge pull request #21 from CMCC-Foundation/gebco

Changes on gebco handling and plotting

Author: santibravocmcc

config.toml

@@ -13,6 +13,7 @@
     advanced_parameters_path = "WITOIL_iMagine/src/parameters.toml" # this path shuld be provided only if "advanced_parameters" is = true
 [download]
     download_data = true # = true if data are not provided by the user
+    download_bath = true # = true if user wants to crop gebco from ERDDAP
     download_curr = true # = true : OCE data are downloaded from Copernicus Marine Service
     download_wind = true # = true : MET data are downloaded from ECMWF ERA5 product
     copernicus_user = "none"
@@ -25,7 +26,7 @@
     delta = [0.75] # default domain length in degrees (applied to both lon/lat), to download or crop data
     # note: delta is used only if set_domain = false
     [input_files.dtm]
-        bathymetry_path = "WITOIL_iMagine/data/gebco/GEBCO.nc"
+        bathymetry_path = false
         coastline_path = "WITOIL_iMagine/data/gshhs/GSHHS_shp/f/GSHHS_f_L1.shp"
     [input_files.metoce]
         oce_data_path = false # to provide if dowload_curr = false

main.py

@@ -159,6 +159,28 @@ def data_download_medslik(
         else:
             down = "global"
 
+        if config["download"]["download_bath"]:
+            output_path = "data/gebco/"
+            output_name = output_path + "GEBCO_slice_{}_{}_mdk.nc".format(
+                identifier, config["simulation"]["name"]
+            )
+            logger.info("Obtaining GEBCO bathymetry from ERDDAP")
+            gebco_errdap(            
+                lat_min,
+                lat_max,
+                lon_min,
+                lon_max,
+                output_name=output_name
+            )
+            # Copy downloaded GEBCO file(s) to the bnc_files directory
+            source_files = gg(f"{output_path}*{identifier}*{config['simulation']['name']}*.nc")
+            destination = os.path.join(root_directory, "bnc_files")
+            for file in source_files:
+                shutil.copy(file, destination)
+            # Remove the temporary GEBCO files from the output path
+            for file in source_files:
+                os.remove(file)
+
         if config["download"]["download_curr"]:
 
             output_path = "data/COPERNICUS/"

src/download/__init__.py

@@ -1,2 +1,3 @@
 from .download_copernicus_parser import *
 from .download_era5_parser import *
+from .download_gebco_ERDDAP import *

src/download/download_gebco_ERDDAP.py

@@ -0,0 +1,21 @@
+import xarray as xr
+
+# Functions outside this script
+from WITOIL_iMagine.src.utils.utils import *
+
+def gebco_errdap(
+    min_lat,
+    max_lat,
+    min_lon,
+    max_lon,
+    output_name=str,
+    ):
+    
+    ds = xr.open_dataset('https://erddap.cmcc-opa.eu/erddap/griddap/Surf_f204_4c2a_5962')
+
+    ds = ds.sel(latitude = slice(min_lat,max_lat), longitude = slice(min_lon,max_lon))
+
+    # Rename variables only if they exist in the dataset
+    ds = Utils.rename_netcdf_variables_mdk3(ds)
+
+    ds.to_netcdf(output_name)

src/plot/plot_mdk3.py

@@ -97,8 +97,28 @@ def plot_matplotlib(self,lon_min,lon_max,lat_min,lat_max):
         # cropping coastline to area of interest
         rec = land.cx[lon_min:lon_max, lat_min:lat_max]
 
-        # selecting simulation domain
-        curr = curr.sel(lon=slice(lon_min, lon_max), lat=slice(lat_min, lat_max))
+        # Get the full coordinate arrays from the dataset
+        lon_arr = curr.lon.values
+        lat_arr = curr.lat.values
+
+        # Find indices for the outward slice
+        # For the lower bound - the grid value that is <= the lon_min or lat_min.
+        i_min = np.searchsorted(lon_arr, lon_min, side="right") - 1
+        j_min = np.searchsorted(lat_arr, lat_min, side="right") - 1
+
+        # For the upper bound - the grid value that is >= the lon_max or lat_max.
+        i_max = np.searchsorted(lon_arr, lon_max, side="left")
+        j_max = np.searchsorted(lat_arr, lat_max, side="left")
+
+        # Make sure indices are within bounds
+        i_min = max(i_min, 0)
+        j_min = max(j_min, 0)
+        i_max = min(i_max, len(lon_arr) - 1)
+        j_max = min(j_max, len(lat_arr) - 1)
+
+        ### Subset the current data to the plot boundaries ###
+        # Slice the dataset using these indices
+        curr = curr.isel(lon=slice(i_min, i_max + 1), lat=slice(j_min, j_max + 1))
 
         min_concentration = 0.0
         max_concentration = (ds_particles.concentration * 1000).max().values
@@ -782,9 +802,9 @@ def __init__(self, x, y):
         self.width = x[-1] - x[0]
         self.height = y[-1] - y[0]
 
-        if not np.allclose(np.diff(x), self.width / (self.nx - 1), atol=1e-5):
+        if not np.allclose(np.diff(x), self.width / (self.nx - 1), atol=1e-4):
             raise ValueError("'x' values must be equally spaced")
-        if not np.allclose(np.diff(y), self.height / (self.ny - 1), atol=1e-5):
+        if not np.allclose(np.diff(y), self.height / (self.ny - 1), atol=1e-4):
             raise ValueError("'y' values must be equally spaced")
 
     @property

src/preprocessing/preprocessing_mdk3.py

@@ -130,15 +130,20 @@ def _copy_nc_files(self, src_files: str, dst_dir: str) -> None:
     def process_bathymetry(self,gebco):
 
         grid = self.grid
-        gebco = xr.open_dataset(gebco)
-
-        if 'nav_lat' in grid and 'nav_lon' in grid:
-            grid = grid.rename({'nav_lat': 'lat', 'nav_lon': 'lon'})
+        
+        if gebco is None or not gebco:
+            bnc_path = f"{self.exp_folder}/bnc_files/*.nc"
+        else:
+            bnc_path = gebco
 
-        # try:
-        #     grid = grid.rename({'nav_lat':'lat','nav_lon':'lon'})
-        # except:
-        #     pass
+        # Use open_mfdataset with the netcdf4 engine (instead of open_dataset)
+        gebco = xr.open_mfdataset(bnc_path, engine="netcdf4")
+        
+        # --- NEW: Rename variables if they are 'latitude'/'longitude'
+        try:
+            gebco = gebco.rename({"latitude": "lat", "longitude": "lon"})
+        except Exception:
+            pass
 
         # interpolation on medslik grid
         med = gebco.interp(lon=grid.lon.values.tolist(),lat=grid.lat.values.tolist())