Updates to pyldas

functions.py

@@ -0,0 +1,35 @@
+import os
+import logging
+
+import numpy as np
+
+logging.basicConfig(format='%(levelname)s: %(message)s', level=logging.INFO)
+
+def walk_up_folder(path, depth=1):
+    """ Walk up a specific number of sub-directories """
+    _cur_depth = 0
+    while _cur_depth < depth:
+        path = os.path.dirname(path)
+        _cur_depth += 1
+    return path
+
+def find_files(path,searchstr):
+    """ Recursive file search with a given search string """
+
+    res = []
+    for root, dirs, files in os.walk(path):
+        for f in files:
+            if f.find(searchstr) != -1:
+                res.append(os.path.join(root, f))
+
+    # multiple files in rc_out possible for
+    multilist = ['catbias','driver','mwRTMparam','RTMparam','catparam','ensprop','obslog','domdecomp','ensupd','obsparam']
+    if len(res) == 0:
+        logging.warning('No files found which contain: "' + searchstr + '".')
+        return None
+    elif len(res) == 1:
+        return res[0]
+    elif len(res) > 1 and multilist.count(searchstr)>=1:
+        return res[0]
+    else:
+        return np.array(res)

grids.py

@@ -33,7 +33,6 @@ class EASE2(object):
     """
 
     def __init__(self, tilecoord=None, tilegrids=None, gtype=None):
-
         self.tilecoord = tilecoord
         self.tilegrids = tilegrids
 
@@ -42,27 +41,36 @@ def __init__(self, tilecoord=None, tilegrids=None, gtype=None):
                 logging.error('No grid type specified and no tilegrids file available.')
                 return
             else:
-                gtype = self.tilegrids.gridtype['global'].decode('utf-8').strip().split('_')[1]
+                grid_str = self.tilegrids.gridtype['global'].decode('utf-8').strip()
+                print(f"DEBUG: Processed gridtype['global']: {grid_str}")
+
+                if 'EASE' in grid_str:
+                    gtype = grid_str.split('_')[1]
+                elif 'LatLon' in grid_str:
+                    gtype = 'LatLon'  # Custom handling for LatLon grids
+                else:
+                    raise ValueError(f"Unsupported gridtype format: {grid_str}")
 
         if (tilecoord is None) or (tilegrids is None):
             logging.warning('tilecoord and/or tilegrids not specified! LatLon - ColRow conversions will not work properly!')
 
         self.setup_grid(gtype)
 
-
     def setup_grid(self, gridtype):
-
         if gridtype == 'M36':
             map_scale = 36032.220840584
         elif gridtype == 'M25':
-            map_scale = 25025.26 # fewer digits are defined by NSIDC to avoid projection issues at +-180 degree
+            map_scale = 25025.26
         elif gridtype == 'M09':
             map_scale = 9008.055210146
         elif gridtype == 'M03':
             map_scale = 3002.6850700487
+        elif gridtype == 'LatLon':
+            self.ease_lons = np.arange(-180, 180, self.tilegrids.loc['global', 'dlon'])
+            self.ease_lats = np.arange(-90, 90, self.tilegrids.loc['global', 'dlat'])
+            return  # Skip EASE2-specific setup
         else:
-            raise NotImplementedError(
-                "Only M03, M09 and M36 grids supported .")
+            raise NotImplementedError("Only M03, M09, M36, and LatLon grids are supported.")
 
         # Upper boundary of north(south)-most M36 pixel. Used for cutting all valid M03/M09 pixels above (below).
         latmax = 85.04457

interface.py

@@ -103,13 +103,14 @@ def __init__(self,
 
             if param == 'ObsFcstAna':
                 self.files = np.sort(list(path.glob(f'**/*{param}.*.bin')))
+            elif param == 'smapL4SMaup':
+                self.files = np.sort(list(path.glob(f'**/*{param}.*.bin')))
             else:
                 if self.mode == 'GEOSldas':
                     self.files = np.sort(list(path.glob(f'**/*{param}.*.nc4')))
                 else:
                     self.files = np.sort(list(path.glob(f'**/*{param}*.bin')))
 
-
             if param == 'hscale':
                 self.pentads = np.array([f.name[-6:-4] for f in self.files]).astype('int')
                 self.orbits = np.array([f.name[-9:-8] for f in self.files])
@@ -918,6 +919,17 @@ def bin2netcdf(self,
         dataset.close()
         self.images = xr.open_dataset(out_file)
 
+    def read_nc4_file(fname):
+        # Open the NetCDF file
+        with Dataset(fname, mode='r') as nc_file:
+            # Create a dictionary to store all variables
+            data = {}
+
+            # Iterate over all variables in the file and read them
+            for var_name in nc_file.variables.keys():
+                data[var_name] = nc_file.variables[var_name][:]
+
+        return data
 
     def mergenc4files(self,
                    overwrite=False,
@@ -984,8 +996,14 @@ def mergenc4files(self,
 
         # Use grid lon lat to avoid rounding issues
         tmp_tilecoord = self.grid.tilecoord.copy()
-        tmp_tilecoord['com_lon'] = self.grid.ease_lons[self.grid.tilecoord.i_indg]
-        tmp_tilecoord['com_lat'] = self.grid.ease_lats[self.grid.tilecoord.j_indg]
+        if 'LatLon' in self.grid.tilegrids.gridtype['global'].decode('utf-8'):
+            tmp_tilecoord['com_lon'] = self.grid.ease_lons[
+                self.grid.tilecoord.i_indg.clip(upper=len(self.grid.ease_lons) - 1)]
+            tmp_tilecoord['com_lat'] = self.grid.ease_lats[
+                self.grid.tilecoord.j_indg.clip(upper=len(self.grid.ease_lats) - 1)]
+        else:
+            tmp_tilecoord['com_lon'] = self.grid.ease_lons[self.grid.tilecoord.i_indg]
+            tmp_tilecoord['com_lat'] = self.grid.ease_lats[self.grid.tilecoord.j_indg]
 
         # Clip region based on specified coordinate boundaries
         ind_img = self.grid.tilecoord[(tmp_tilecoord['com_lon']>=lonmin)&(tmp_tilecoord['com_lon']<=lonmax)&
@@ -1004,18 +1022,28 @@ def mergenc4files(self,
                     if attr[0] != '_':
                         dataset[var].setncattr(attr, ds[var].getncattr(attr))
 
-        for i, (f, dt) in enumerate(zip(files,dates)):
+        for i, (f, dt) in enumerate(zip(files, dates)):
             logging.info('%d / %d' % (i, len(dates)))
 
             with Dataset(f) as data:
-                img = np.full((len(lats),len(lons)), -9999., dtype='float32')
-                ind_lat = self.grid.tilecoord.loc[ind_img, 'j_indg'].values - self.grid.tilegrids.loc['domain','j_offg'] - j_offg_2
-                ind_lon = self.grid.tilecoord.loc[ind_img, 'i_indg'].values - self.grid.tilegrids.loc['domain','i_offg'] - i_offg_2
+                img = np.full((len(lats), len(lons)), -9999., dtype='float32')
+
+                # Extract indices
+                ind_lat = self.grid.tilecoord.loc[ind_img, 'j_indg'].values - self.grid.tilegrids.loc[
+                    'domain', 'j_offg'] - j_offg_2
+                ind_lon = self.grid.tilecoord.loc[ind_img, 'i_indg'].values - self.grid.tilegrids.loc[
+                    'domain', 'i_offg'] - i_offg_2
+
+                # Check if grid type is LatLon, and clip indices if necessary
+                if 'LatLon' in self.grid.tilegrids.gridtype['global'].decode('utf-8'):
+                    ind_lat = np.clip(ind_lat, 0, len(lats) - 1)
+                    ind_lon = np.clip(ind_lon, 0, len(lons) - 1)
 
+                # Assign values
                 for var in variables:
-                    tmp_img = data[var][0,ind_img-1].data
-                    img[ind_lat,ind_lon] = tmp_img
-                    dataset.variables[var][i,:,:] = img
+                    tmp_img = data[var][0, ind_img - 1].data
+                    img[ind_lat, ind_lon] = tmp_img
+                    dataset.variables[var][i, :, :] = img
 
         # Save file to disk and loat it as xarray Dataset into the class variable space
         dataset.close()

paths.py

@@ -7,7 +7,6 @@
 class paths(object):
     """
     Class holding the most important LDAS path definitions
-
     Parameters
     ----------
     mode : string
@@ -18,7 +17,6 @@ class paths(object):
         experiment name (appended to root path)
     domain : string
         domain name (appended to experiment path)
-
     Attributes
     ----------
     exp_root : string
@@ -35,7 +33,6 @@ class paths(object):
         Path to restart files (for continuing processing or spin-up)
     plots : string
         Path to plots
-
     """
 
     def __init__(self, mode, root=None, exp=None, domain=None):

templates.py

@@ -64,6 +64,9 @@ def get_template(param, **kwargs):
     elif param == 'smosL4SMaup':
         dtype, hdr, length = template_smosL4SMaup()
 
+    elif param == 'smapL4SMaup':
+        dtype, hdr, length = template_smapL4SMaup()
+
     else:
         logging.warning('No template found for "' + param + '".')
         dtype, hdr, length = (None, None, None)
@@ -474,6 +477,48 @@ def template_smosL4SMaup():
     return dtype, hdr, length
 
 
+def template_smapL4SMaup():
+    """"
+    Template for reading smos L4 soil moisture DA output
+
+    """
+
+    hdr = None
+    length = None
+    dtype = np.dtype([('obs_h_time', '>f8'),
+                      ('obs_v_time', '>f8'),
+                      ('res_h_flag', '>i4'),
+                      ('res_v_flag', '>i4'),
+                      ('orbit_h_flag', '>i4'),
+                      ('orbit_v_flag', '>i4'),
+                      ('obs_h', '>f4'),
+                      ('obs_v', '>f4'),
+                      ('obs_h_assim', '>f4'),
+                      ('obs_v_assim', '>f4'),
+                      ('obs_h_errstd', '>f4'),
+                      ('obs_v_errstd', '>f4'),
+                      ('fcst_h', '>f4'),
+                      ('fcst_v', '>f4'),
+                      ('fcst_h_ensstd', '>f4'),
+                      ('fcst_v_ensstd', '>f4'),
+                      ('sm_srfc_fcst', '>f4'),
+                      ('sm_rz_fcst', '>f4'),
+                      ('sm_prof_fcst', '>f4'),
+                      ('surface_temp_fcst', '>f4'),
+                      ('soil_temp_fcst', '>f4'),
+                      ('sm_srfc_ana', '>f4'),
+                      ('sm_rz_ana', '>f4'),
+                      ('sm_prof_ana', '>f4'),
+                      ('surface_temp_ana', '>f4'),
+                      ('soil_temp_ana', '>f4'),
+                      ('sm_srfc_ana_ensstd', '>f4'),
+                      ('sm_rz_ana_ensstd', '>f4'),
+                      ('sm_prof_ana_ensstd', '>f4'),
+                      ('surface_temp_ana_ensstd', '>f4'),
+                      ('soil_temp_ana_ensstd', '>f4'),
+                      ])
+
+    return dtype, hdr, length