Some cleaning of the code

src/tbxsectors.py

@@ -5,7 +5,6 @@ import matplotlib.path as mpath
 import cartopy.crs as ccrs
 import cartopy.feature as cfeat
 from pathlib import Path
-from pyproj import CRS, Transformer
 
 def define_sectors(ref="RH"):
     """
@@ -86,8 +85,9 @@ def interp_mask2grid(target_grid=None,
     Inputs: - target_grid: the target grid on which to interpolate the NSIDC mask. 
                            It can be a xr.Dataset containing lat and lon variables defining the grid,
                            or a string mentioning the grid: 
-                               if weights='ORCA1_nh', use the weights for Northern Hemisphere ORCA1 config; 
-                               if 'ORCA025_nh', use the weights for NorthH for ORCA025
+                            - if 'ORCA1_nh', use the mask for Northern Hemisphere ORCA1 config; 
+                            - if 'ORCA025_nh', use the mask for Northern Hemisphere for ORCA025;
+                            - any other string: a user-calculated mask already interpolated on the target grid.
                            #TODO: need to add full ORCA1 and ORCA025, as well as ORCA12 and ORCA12_nh
             -maskFile: path and name of file containing the NSIDC mask.
 
@@ -101,41 +101,22 @@ def interp_mask2grid(target_grid=None,
             ds_mask_interp = xr.open_dataset(Path.cwd().joinpath('../data/NSIDCRegionsMask_ORCA025_nh.nc'))
         else:
             ds_mask_interp = xr.open_dataset(target_grid)
-    else:
+    else: # If the target grid is a dataset or dataarray, interpolate the mask onto it.
         # First of all, load the mask
         ds_mask = xr.open_dataset(maskFile)
-        # ------------- If no lat-lon in the ds_mask, calculate and add them -------------
-        if 'lat' not in list(ds_mask.coords):
-            # Convert from projected coordinates to lat-lon coordinates 
-            # A first issue is that the mask is provided in North Polar Stereographic projection
-            # Extract the x and y (projection) coordinates and grid them (from 1D to 2D)
-            x_mask, y_mask = np.meshgrid(ds_mask.x, ds_mask.y)
-            # Define the projection and geographic systems according to CRS conventions.
-            orig_crs = CRS.from_epsg(3413) # North Polar Stereographic from NSIDC
-            target_crs = CRS.from_epsg(4326) # Geographic coordinate system (i.e. lat-lon)
-            # Use pyproj.Transformer to convert from NSIDC projection to lat-lon coordinates
-            transformer = Transformer.from_crs(orig_crs, target_crs)
-            lat_src, lon_src = transformer.transform(x_mask, y_mask)
-            # Assign the new coordinates to the dataset
-            ds_mask = ds_mask.assign_coords({'lon': (['y', 'x'], lon_src), 'lat': (['y', 'x'], lat_src)})
-
         # ------------- Bulk of the function: interpolate the mask to the ds grid -------------
-        # If 'weights' are given, load them directly. If not, calculate them
-        # But since xesmf does not exist yet on `cyclone`, need to deal with it.
-        # try:  # If xesmf exists, create or load the regridder and regrid
-        import xesmf as xe # type: ignore
+        # Import xesmf. Will throw an error if not available
+        import xesmf as xe
+        # Create a regridder.
         Regridder = xe.Regridder(ds_mask, target_grid, method='bilinear')
         # Now apply this regridder to the mask
         ds_mask_interp = Regridder(ds_mask)
-        # except ImportError: # If the import failed, send an error
-            # print("ERROR: xesmf could not be loaded, go to `coriolis`.")
-            # return
-    # return the interpolated mask.
+    # Return the interpolated mask.
     return ds_mask_interp
 
 
 def groupby_sectors(ds, ref=None,
-                    target_gd=None, 
+                    target_gd=None,
                     maskFile=Path.cwd().joinpath('../data/NSIDCRegions_N3.125km_v1.1_wLatLon_df1.nc')):
     """
     Define sectors and Groupby a dataset or dataArray into those sectors.