Modified plot_sectormaps and K16 definition

src/tbxsectors.py

@@ -15,7 +15,7 @@ def define_sectors(ref="RH"):
         - ref="Z83" or "ZwallyEtAl1983": Zwally et al. (1983, see also doi:10.1029/2000JC000733), based on geography (same definition as Cavalieri & Parkinson [2008], e.g.)
     - for Northern Hemisphere:
         - ref="A21" or "ArthunEtAl2021": Arthun et al. (2021, doi:10.1029/2020GL090825), only covers high Arctic
-        - ref="K16" or "KoenigkEtAl2016": Koenigk et al. (2016, doi:10.1007/s00382-015-2586-1), covers sub-Arctic seas
+        - ref="K16" or "KoenigkEtAl2016": based on Koenigk et al. (2016, doi:10.1007/s00382-015-2586-1), covers sub-Arctic seas but modified from original publication to remove gaps and overlaps.
         - ref="NSIDC" (recommended): NSIDC new (as of 2023) definition (NSIDC Special report 25, figure 3)
     Returns a dictionary with names of the sectors as keys and a 2-items list of the bounding longitudes (E and W) for each sector.
 
@@ -59,7 +59,7 @@ def define_sectors(ref="RH"):
         sectors['GreenlandSea'] = [[50,80],[-40, 15]]
         sectors['BarentsKara']  = [[65,80],[15, 100]]
         sectors['SiberianSeas'] = [[65,80],[100, 180]]
-        sectors['Chukchi']      = [[65,80],[-180,-160]]
+        sectors['Chukchi']      = [[60,80],[-180,-160]]
         sectors['Beaufort']     = [[65,80],[-160,-100]]
    # NSIDC's 2023 definition: covers the whole Arctic, many sectors (18!)
     elif ref.lower() in ["nsidc"]:
@@ -73,7 +73,7 @@ def define_sectors(ref="RH"):
     return sectors
 
 def interp_mask2grid(target_grid=None,
-                    maskFile=Path.cwd().joinpath('../data/NSIDCRegions_N3.125km_v1.1_wLatLon_df1.nc')):
+                    maskFile=Path.home().joinpath('toolbox/data/NSIDCRegions_N3.125km_v1.1_wLatLon_df1.nc')):
     """
     This function creates a region mask based on the NSIDC definition for the Arctic Ocean.
     It loads the NSIDC mask and interpolates it onto the provided grid.
@@ -96,9 +96,9 @@ def interp_mask2grid(target_grid=None,
     # Before calculating weights, check if this is one of the pre-calculated situations
     if isinstance(target_grid,str):
         if target_grid == 'ORCA1_nh':
-            ds_mask_interp = xr.open_dataset(Path.cwd().joinpath('../data/NSIDCRegionsMask_ORCA1_nh.nc'))
+            ds_mask_interp = xr.open_dataset(Path.home().joinpath('toolbox/data/NSIDCRegionsMask_ORCA1_nh.nc'))
         elif target_grid == 'ORCA025_nh':
-            ds_mask_interp = xr.open_dataset(Path.cwd().joinpath('../data/NSIDCRegionsMask_ORCA025_nh.nc'))
+            ds_mask_interp = xr.open_dataset(Path.home().joinpath('toolbox/data/NSIDCRegionsMask_ORCA025_nh.nc'))
         else:
             ds_mask_interp = xr.open_dataset(target_grid)
     else: # If the target grid is a dataset or dataarray, interpolate the mask onto it.
@@ -117,7 +117,7 @@ def interp_mask2grid(target_grid=None,
 
 def groupby_sectors(ds, ref=None,
                     target_gd=None,
-                    maskFile=Path.cwd().joinpath('../data/NSIDCRegions_N3.125km_v1.1_wLatLon_df1.nc')):
+                    maskFile=Path.home().joinpath('toolbox/data/NSIDCRegions_N3.125km_v1.1_wLatLon_df1.nc')):
     """
     Define sectors and Groupby a dataset or dataArray into those sectors. 
     This can then be used to perform calculations such as ice extent or ice area for each sector.
@@ -127,7 +127,7 @@ def groupby_sectors(ds, ref=None,
                    'Z83' or 'ZwallyEtAl1983' (for Zwally et al. (1983), AA),
                    'RH' or 'RaphaelHobbs' (for Raphael and Hobbs (2014), AA),
                    'A21' or 'Arthun2021' (for Arthun et al. (2021), Arctic),
-                   'K16' or 'Koenigk2016' (for Koenigk et al. (2016), Arctic),
+                   'K16' or 'Koenigk2016' (modified Koenigk et al. (2016), Arctic),
                    'NSIDC' or 'nsidc' (for NSIDC mask, new as of 2023, Arctic).
                    If None (not recommended), then attempts to figure out if `ds` covers Arctic or AA,
                    then use either the NSIDC mask for Arctic or the RH ref for the AA.
@@ -226,11 +226,17 @@ def groupby_sectors(ds, ref=None,
     return ds_out.groupby('sector')
 
 
-def plot_sectormaps():
+def plot_sectormaps(listRefs = ['K16', 'A21', 'NSIDC', 'Z83', 'RH']):
     """ 
     Plot a 5-panel figure which illustrates the different sector definitions available. 
-    This is called by default when executing the module.
+    Input: -listRefs: list of definitions to plot. By default, include all availables, 
+                      but the user can ask for a specific one or a subset. 
+    Outputs: return the figure instance.
+
     """
+    # Quick test to see if listRefs is iterable
+    if not isinstance(listRefs,list):
+        listRefs=[listRefs]
     ### Some preparatory stuff: define sectors, prepare the circular frame, etc.
     # Prep circular boundary
     r_extent = 4651194.319 * 1.005
@@ -249,9 +255,8 @@ def plot_sectormaps():
 
     # Second, apply the sector definitions to get a groupby, then ungroup everything to get one grid
     # This is a bit circonvoluted, but hey it works!
-    listRefs = ['K16', 'A21', 'NSIDC', 'Z83', 'RH']
-    listFullNames = ['Koenigk et al. (2016)', 'Arthun et al. (2021)', 'NSIDC (2023)',
-                     'Zwally et al. (1983)', 'Raphael and Hobbs (2014)']
+    dictFullNames = {'K16': 'Koenigk et al. (2016)', 'A21': 'Arthun et al. (2021)', 'NSIDC': 'NSIDC (2023)',
+                     'Z83': 'Zwally et al. (1983)', 'RH': 'Raphael and Hobbs (2014)'}
     for n, ref in enumerate(listRefs):
         # Get a groupby dataArray, and ungroup it to retrieve a mask
         try:  # If xesmf exists, create or load the regridder and regrid
@@ -304,12 +309,12 @@ def plot_sectormaps():
         ax.set_boundary(circle_path) # Create a round boundary instead of a square one. Just aesthetic.
         ax.set_frame_on(False)  # hide the boundary frame
         # Now can put the title on top of the panel
-        ax.set_title(listFullNames[n])
-        # Restrain the layou to make it fit tighter.
+        ax.set_title(dictFullNames[ref])
+        # Restrain the layout to make it fit tighter.
     fig.tight_layout()
     return fig
 
 
 if __name__ == "__main__":
     print("Available sector definitions:")
-    plot_sectormaps()
+    plot_sectormaps()
\ No newline at end of file