diff --git a/batchscripts/create_netcdf_images_GEOSldas_isimip.py b/batchscripts/create_netcdf_images_GEOSldas_isimip.py new file mode 100644 index 0000000..e69de29 diff --git a/functions.py b/functions.py new file mode 100644 index 0000000..749dde3 --- /dev/null +++ b/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) \ No newline at end of file diff --git a/grids.py b/grids.py index bec4373..dfcb12a 100644 --- a/grids.py +++ b/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 diff --git a/interface.py b/interface.py index c05ef00..f1d6615 100644 --- a/interface.py +++ b/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() diff --git a/paths.py b/paths.py index 20452a9..0588e10 100644 --- a/paths.py +++ b/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): diff --git a/templates.py b/templates.py index e10a34b..cbf4648 100644 --- a/templates.py +++ b/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