Merge pull request #1375 from DestinE-Climate-DT/icon_healpix

Icon healpix grids update

Author: mnurisso

CHANGELOG.md

@@ -9,15 +9,15 @@ Unreleased in the current development version:
 
 AQUA core complete list:
 
+- Mixed updates to support data for NextGEMS cycle4 hackathon (#1375)
 - Preprocess functionality added to the `Reader` class (#1298)
 - The AQUAthon material has been moved under the `notebooks` folder (#1342)
 - `aqua` source code has been moved to the `src` folder (#1332)
 - A diagnostic module, called `aqua.diagnostics`, has been created under the `src` folder (#1332, #1341)
-- LRA generator tool support for multiple relizations (#1357)
+- LRA generator tool support for multiple relizations (#1357, #1375)
 - LRA generator requires `catalog` as a mandatory argument (#1357)
 - AQUA console revisiting, adding `avail` method and `update` method (#1346)
 - AQUA install now requires mandatory machine name (#1346)
-
 - Fix to make keyword step optional in request (#1360)
 
 AQUA diagnostic complete list:

cli/grid-from-data/config-hpx-icon.yaml

@@ -1,17 +1,21 @@
 # These are needed for the AQUA retrieve
+catalog: nextgems3
 model: ICON
-exp: historical-1990
-source: daily-hpz10-oce3d
-var: avg_so
+exp: ngc3028
+source: P1D
+var: tke
 
 # These are needed for grid conversion with CDO
-zoom: 10
+zoom: 5
 nested: True # otherwise ring
 
 # These are extra info:
 # model name for saving the grid
 model_name: icon
+resolution: R02B09
+extra: null
+
 # tmp folder for the grid generation
-tmp: /pfs/lustrep3/scratch/project_465000454/padavini
+tmp: /home/b/b382076/scratch/healpix-icon-ng3/tmp
 # target folder for the grid generated
-tgt: /pfs/lustrep3/projappl/project_465000454/data/AQUA/grids/HealPix
+tgt: /home/b/b382076/scratch/healpix-icon-ng3

cli/grid-from-data/hpx-from-sources.py

@@ -11,17 +11,21 @@
 from aqua.util import load_yaml, get_arg, create_folder
 from aqua.logger import log_configure
 
+reasonable_vert_coords = ['depth_full', 'depth_half', 'level']
+
 
 def parse_arguments(args):
     """
     Parse command line arguments
     """
 
-    parser = argparse.ArgumentParser(description='AQUA gribber')
+    parser = argparse.ArgumentParser(description='AQUA grids from data')
     parser.add_argument('-c', '--config', type=str, required=True,
                         help='yaml file with exp information and directories')
     parser.add_argument('-l', '--loglevel', type=str,
                         help='log level [default: WARNING]')
+    parser.add_argument('-z', '--zoom', type=str,
+                        help='override zoom: convenient for loop on many zoom levels')
 
     return parser.parse_args(args)
 
@@ -42,28 +46,44 @@ def parse_arguments(args):
     source = config['source']
     var = config['var']
     resolution = config.get('resolution')
+    extra = config.get('extra')
     logger.info('Retrieving %s from %s %s %s', var, model, exp, source)
 
     # Configuration needed to save the file
     tmp = config['tmp']
     model_name = config['model_name']
-    zoom = config['zoom']
+    zoom = int(get_arg(args, 'zoom', config.get('zoom'))) #HACK for zoom
+    catalog = config.get('catalog')
     nested = config['nested']
 
+    # adding zoom if found
+    mykwargs = {}
+    if zoom:
+        mykwargs = {**mykwargs, **{'zoom': zoom}}
+    if catalog:
+        mykwargs = {**mykwargs, **{'catalog': catalog}}
+        
+
     # Create Reader object
     reader = Reader(model=model, exp=exp, source=source,
-                    areas=False, fix=False, loglevel=loglevel)
+                    areas=False, fix=False, loglevel=loglevel, **mykwargs)
     data = reader.retrieve(var=var)
 
+    # detect vertical coordinate from a set of plausible values
+    vert_coord = list(set(reasonable_vert_coords) & set(data.coords))
+
     #automatic detection of 3d
-    if data['level'] is not None:
+    if vert_coord:
         model3d = True
+        if len(vert_coord)>1:
+            raise KeyError("Too many vertical coordinates identified, check the data manually")
+        vert_coord = vert_coord[0]
     else: 
         model3d = False
 
     if model3d:
         logger.debug("Modifying level axis attributes as Z")
-        data['level'].attrs['axis'] = 'Z'
+        data[vert_coord].attrs['axis'] = 'Z'
 
     # Save data in a netcdf file on the temporary directory
     create_folder(tmp)
@@ -78,28 +98,38 @@ def parse_arguments(args):
 
     # cdo command setup:
     nside = 2**zoom
+   
+    # Setting grid name
     if nested:
-        grid_name = 'hp' + str(nside) + '_nested'
+        grid_name = f'hp{nside}_nested'
     else:  # ring
-        grid_name = 'hp' + str(nside) + '_ring'
+        grid_name = f'hp{nside}_ring'
 
-    # setting output filename
+    # Setting output filename
     tgt = config['tgt']
     create_folder(tgt)
-    filename_tgt = tgt + '/' + model_name
+    filename_tgt = os.path.join(tgt, model_name)
+
+    # add info on original data resolution if available
     if resolution:
-        filename_tgt = filename_tgt + "-" + resolution
+        filename_tgt = f"{filename_tgt}-{resolution}"
 
-    filename_tgt = filename_tgt + '_hpz' + str(zoom)
+    # add info on hpz zoom
+    filename_tgt = f"{filename_tgt}_hpz{zoom}"
+
+    # options for nested and ringed
     if nested:
-        filename_tgt = filename_tgt + '_nested_oce'
+        filename_tgt = f"{filename_tgt}_nested_oce"
     else:
-        filename_tgt = filename_tgt + '_ring_oce'
+        filename_tgt = f"{filename_tgt}_ring_oce"
 
     if model3d:
-        filename_tgt = filename_tgt + '_level.nc'
-    else:
-        filename_tgt = filename_tgt + '.nc'
+        filename_tgt = f"{filename_tgt}_{vert_coord}"
+
+    if extra:
+        filename_tgt = f"{filename_tgt}_{extra}"
+
+    filename_tgt = f'{filename_tgt}.nc'
 
     logger.info('Setting grid %s', grid_name)
     logger.info('Saving data in %s', filename_tgt)

cli/lra/auto_LRA_overnight.yaml

@@ -2,14 +2,18 @@
 target:
   resolution: r100
   frequency: monthly
+  catalog: climatedt-phase1
+
+path:
   outdir: /pfs/lustrep3/projappl/project_465000454/data/AQUA/LRA
   tmpdir: /pfs/lustrep3/projappl/project_465000454/padavini/tmp
 
-loglevel: INFO
 
+options:
 # Set to True if lra-r100-monthly-zarr should be created and checked
-zarr: False
-verify_zarr: False
+  zarr: False
+  verify_zarr: False
+  loglevel: INFO
 
 slurm:
   partition: small

cli/lra/cli_lra_parallel_slurm.py

@@ -66,10 +66,12 @@ def submit_sbatch(model, exp, source, varname, slurm_dict, yaml_file,
     ]
 
     if dependency is not None:
+        print(dependency)
         sbatch_cmd.append('--dependency=afterany:'+ str(dependency))
 
     # Add script command
-    sbatch_cmd.append('aqua lra')
+    sbatch_cmd.append('aqua')
+    sbatch_cmd.append('lra')
     sbatch_cmd.append('--config')
     sbatch_cmd.append(yaml_file)
     sbatch_cmd.append('--model')
@@ -92,11 +94,17 @@ def submit_sbatch(model, exp, source, varname, slurm_dict, yaml_file,
         if overwrite:
             sbatch_cmd.append('-o')
         sbatch_cmd.append('-d')
-        result = subprocess.run(sbatch_cmd, capture_output = True, check=True).stdout.decode('utf-8')
-        jobid = re.findall(r'\b\d+\b', result)[-1]
-        return jobid
-
-    #print(sbatch_cmd)
+        try:
+            result = subprocess.run(sbatch_cmd, capture_output = True, check=True).stdout.decode('utf-8')
+            jobid = re.findall(r'\b\d+\b', result)[-1]
+            return jobid
+        except subprocess.CalledProcessError as e:
+            # Print the error message and stderr if the command fails
+            print(f"Command failed with return code {e.returncode}")
+            print(f"stdout: {e.stdout}")
+            print(f"stderr: {e.stderr}")
+
+    print(sbatch_cmd)
     return 0
 
 
@@ -147,7 +155,7 @@ def parse_arguments(arguments):
                 varnames = config['data'][model][exp][source]['vars']
                 for varname in varnames:
                     if (COUNT % int(parallel)) == 0 and COUNT != 0:
-                        print('Updating parent job to' + jobid)
+                        print('Updating parent job to' + str(jobid))
                         PARENT_JOB = str(jobid)
                     COUNT = COUNT + 1
                     print(' '.join(['Submitting', model, exp, source, varname]))

config/fixes/ICON-default.yaml

@@ -109,7 +109,9 @@ fixer_name:
             avg_von:
                 source: v
                 grib: true
-
+            avg_so:
+                source: so
+                grib: true
             #fixes for NextGems data    
             cfh_lnd: 
                 source: cfh_lnd

config/fixes/OBS-default.yaml

@@ -80,3 +80,11 @@ fixer_name:
           source: sftof
           src_units: frac
           units: frac
+  
+  era5-arco:
+    data_model: false
+    coords:
+      lat:
+        source: latitude
+      lon:
+        source: longitude

config/grids/ICON.yaml

@@ -1,3 +1,4 @@
+# More info on ICON here: https://easy.gems.dkrz.de/Processing/healpix/index.html
 grids:
   # ICON oceanic healpix
   icon-hpz10-nested:
@@ -8,11 +9,11 @@ grids:
       level: '{{ grids }}/HealPix/icon_hpz10_nested_oce_level.nc'
 
   #ICON
-  R02B08:
+  R02B08: # ICON grid R02B08, corresponding to 10km resolution
     path: '{{ grids }}/ICON/icon_grid_0033_R02B08_G.nc'
     space_coord: ["cell"]
     extra: -selname,cell_index
-  R02B09:
+  R02B09: # ICON grid R02B09, corresponding to 5km resolution
     path: '{{ grids }}/ICON/icon_grid_0015_R02B09_G.nc'
     space_coord: ["cell"]
     extra: -selname,cell_index
@@ -25,12 +26,22 @@ grids:
       path: '{{ grids }}/ICON/icon_grid_0015_R02B09_G.nc'
       extra: -selname,cell_index
       space_coord: ["ncells"]
-  icon-healpix:
+  icon-R02B08-hp-nested: # NG4, from ICON grid R02B08, corresponding to 10km resolution, 72 depth_full levels
     path:
-      2d: '{{ grids }}/HealPix_old/icon_hpx{zoom}_atm_2d.nc'   # this is the default 2d grid
-      2dm: '{{ grids }}/HealPix_old/icon_hpx{zoom}_oce_2d.nc'  # this is an additional and optional 2d grid used if data are masked
-      depth_full: '{{ grids }}/HealPix_old/icon_hpx{zoom}_oce_depth_full.nc'
-      depth_half: '{{ grids }}/HealPix_old/icon_hpx{zoom}_oce_depth_half.nc'
+      2d: '{{ grids }}/HealPix/hpz{zoom}_nested_atm.nc'   # this is the default 2d grid
+      2dm: '{{ grids }}/HealPix/icon-R02B08_hpz{zoom}_nested_oce.nc'  # this is an additional and optional 2d grid used if data are masked
+      depth_full: '{{ grids }}/HealPix/icon-R02B08_hpz{zoom}_nested_oce_depth_full.nc'
+      depth_half: '{{ grids }}/HealPix/icon-R02B08_hpz{zoom}_nested_oce_depth_half.nc'
+    masked:   # This is the attribute used to distinguish variables which should go into the masked category
+      component: ocean
+    space_coord: ["cell"]
+    vert_coord: ["depth_half", "depth_full"]
+  icon-R02B09-hp-nested: # NG3, from ICON grid R02B09, corresponding to 5km resolution, 128 depth_full levels
+    path:
+      2d: '{{ grids }}/HealPix/hpz{zoom}_nested_atm.nc'   # this is the default 2d grid
+      2dm: '{{ grids }}/HealPix/icon-R02B09_hpz{zoom}_nested_oce.nc'  # this is an additional and optional 2d grid used if data are masked
+      depth_full: '{{ grids }}/HealPix/icon-R02B09_hpz{zoom}_nested_oce_depth_full.nc'
+      depth_half: '{{ grids }}/HealPix/icon-R02B09_hpz{zoom}_nested_oce_depth_half.nc'
     masked:   # This is the attribute used to distinguish variables which should go into the masked category
       component: ocean
     space_coord: ["cell"]

config/grids/default.yaml

@@ -22,13 +22,25 @@ grids:
   lon-lat-depth:
     vert_coord: ["depth"]
     space_coord: [lon, lat]
+  longitude-latitude:
+    vert_coord: ["2d"]
+    space_coord: [longitude, latitude]
 
   # target regrid
+  # rxxxs are staggered grids
+  r005s: r7200x3601
   r005: r7200x3600
+  r010s: r3600x1801
   r010: r3600x1800
+  r020s: r1800x901
   r020: r1800x900
+  r025s: r1440x721
   r025: r1440x720
+  r050s: r720x361
   r050: r720x360
+  r100s: r360x181
   r100: r360x180
+  r200s: r180x91
   r200: r180x90
+  r250s: r144x73
   r250: r144x72

pyproject.toml

@@ -40,6 +40,7 @@ dependencies = [
     "dask_jobqueue",
     "eccodes==2.38.0",
     "docker",
+    "gcsfs", # Google Cloud Storage File System
     "gribscan",
     "pandas",
     "intake==0.7.0", # Version 2 is not yet compatible with intake-xarray

src/aqua/cli/lra.py

@@ -55,6 +55,8 @@ def lra_parser(parser = None):
                         help='source to be processed. Use with coherence with --exp and --var')
     parser.add_argument('-v', '--var', type=str,
                         help='var to be processed. Use with coherence with --source')
+    #parser.add_argument('-r', '--realization', type=str,
+    #                    help="realization to be processed. Use with coherence with --var")
 
     #return parser.parse_args(arguments)
     return parser
@@ -149,7 +151,7 @@ def lra_cli(args, config, catalog, resolution, frequency, fix, outdir, tmpdir, l
             sources = to_list(get_arg(args, 'source', config['data'][model][exp].keys()))
             for source in sources:
                 # get info on potential realizations
-                realizations = get_arg(args, 'var', config['data'][model][exp][source].get('realizations'))
+                realizations = config['data'][model][exp][source].get('realizations')
                 loop_realizations = to_list(realizations) if realizations is not None else [1]
 
                 # get info on varlist and workers

src/aqua/exceptions.py

@@ -31,3 +31,10 @@ class NoObservationError(Exception):
     def __init__(self, message="No observation available"):
         self.message = message
         super().__init__(self.message)
+
+class NoRegridError(Exception):
+    """Exception raised when no regrid is available"""
+
+    def __init__(self, message="No regrid available"):
+        self.message = message
+        super().__init__(self.message)