From 0ccf8551dc108a5f02d96eb5d0bd4380db103d19 Mon Sep 17 00:00:00 2001 From: Negin Sobhani Date: Thu, 28 Apr 2022 09:13:40 -0600 Subject: [PATCH 1/9] a first draft of implementing xarray mesh --- .../ctsm/site_and_regional/regional_case.py | 102 ++++++++++++++++++ 1 file changed, 102 insertions(+) diff --git a/python/ctsm/site_and_regional/regional_case.py b/python/ctsm/site_and_regional/regional_case.py index 57bb8474f4..40b96d6d5f 100644 --- a/python/ctsm/site_and_regional/regional_case.py +++ b/python/ctsm/site_and_regional/regional_case.py @@ -8,6 +8,8 @@ # -- 3rd party libraries import numpy as np +import xarray as xr +import tqdm # -- import local classes for this script from ctsm.site_and_regional.base_case import BaseCase, USRDAT_DIR @@ -98,6 +100,7 @@ def __init__( self.reg_name = reg_name self.out_dir = out_dir self.create_tag() + self.create_mesh_at_reg() def create_tag(self): """ @@ -226,3 +229,102 @@ def create_landuse_at_reg(self, indir, file, user_mods_dir): os.path.join(USRDAT_DIR, fluse_out) ) self.write_to_file(line, nl_clm) + + + def create_mesh_at_reg (self): + """ + Create a mesh subsetted for the RegionalCase class. + """ + print ("test") + mesh_in = '/glade/p/cesmdata/cseg/inputdata/share/meshes/fv4x5_050615_polemod_ESMFmesh.nc' + node_coords, subset_element, subset_node, conn_dict = self.subset_mesh_at_reg(mesh_in) + + + def subset_mesh_at_reg (self, mesh_in): + """ + This function subsets the mesh based on lat and lon bounds given by RegionalCase class. + """ + f_in = xr.open_dataset (mesh_in) + elem_count = len (f_in['elementCount']) + elem_conn = f_in['elementConn'] + num_elem_conn = f_in['numElementConn'] + center_coords = f_in['centerCoords'] + node_count = len (f_in['nodeCount']) + node_coords = f_in['nodeCoords'] + + subset_element = [] + cnt = 0 + + for n in tqdm.tqdm(range(elem_count)): + #print (elem_conn.shape) + #print (num_elem_conn.shape) + #print ('-----') + #print (numElementConn[n]) + endx = elem_conn[n,:num_elem_conn[n].values].values + #print ('endx:', endx) + #print ('-----') + endx[:,] -= 1# convert to zero based index + endx = [int(xi) for xi in endx] + #print ('-----') + + #endx = endx.values + #print ('endx:', endx) + #print (node_coords.shape) + #print (node_coords) + nlon = node_coords[endx,0] + nlat = node_coords[endx,1] + #print ('-----') + + l1 = np.logical_or(nlon <= self.lon1,nlon >= self.lon2) + l2 = np.logical_or(nlat <= self.lat1,nlat >= self.lat2) + if np.any(np.logical_or(l1,l2)): + #print(nlon,nlat) + pass + else: + subset_element.append(n) + cnt+=1 + + print (subset_element) + + subset_node = [] + conn_dict = {} + cnt = 1 + + for n in range(node_count): + nlon = node_coords[n,0] + nlat = node_coords[n,1] + + l1 = np.logical_or(nlon <= self.lon1,nlon >= self.lon2) + l2 = np.logical_or(nlat <= self.lat1,nlat >= self.lat2) + if np.logical_or(l1,l2): + conn_dict[n+1] = -9999 + else: + subset_node.append(n) + conn_dict[n+1] = cnt + cnt+=1 + return node_coords, subset_element, subset_node, conn_dict + + + + def write_mesh (self, f_in, node_coords, subset_element, subset_node, conn_dict): + corner_pair_uniq = f_in.variables['nodeCoords'][subset_node,] + + dimensions = f1.dimensions + variables = f1.variables + global_attributes = f1.ncattrs() + var = 'elementConn' + elem_conn_out = np.empty(shape=[elementCount, len(f1.dimensions['maxNodePElement'])]) + elem_conn_index = f1.variables[var][subsetElement,] + + print (elem_conn_out.shape) + for n in range(elementCount): + print ('n:',n) + for m in range(len(f1.dimensions['maxNodePElement'])): + print ('m:',m) + ndx = int (elem_conn_index[n,m]) + print ('ndx:', ndx) + elem_conn_out[n,m] = connDict[ndx] + + + + From c3d981d3f807469d09b7dcf8b3c59e67d5dce220 Mon Sep 17 00:00:00 2001 From: Negin Sobhani Date: Tue, 3 May 2022 14:11:45 -0600 Subject: [PATCH 2/9] subset_mesh capabilities --- .../ctsm/site_and_regional/regional_case.py | 123 ++++++++++++++---- python/ctsm/subset_data.py | 8 +- tools/site_and_regional/default_data.cfg | 2 + 3 files changed, 106 insertions(+), 27 deletions(-) diff --git a/python/ctsm/site_and_regional/regional_case.py b/python/ctsm/site_and_regional/regional_case.py index 40b96d6d5f..b497b49d45 100644 --- a/python/ctsm/site_and_regional/regional_case.py +++ b/python/ctsm/site_and_regional/regional_case.py @@ -9,7 +9,8 @@ # -- 3rd party libraries import numpy as np import xarray as xr -import tqdm +from tqdm import tqdm +from datetime import datetime # -- import local classes for this script from ctsm.site_and_regional.base_case import BaseCase, USRDAT_DIR @@ -79,6 +80,7 @@ def __init__( create_landuse, create_datm, create_user_mods, + create_mesh, out_dir, overwrite, ): @@ -98,9 +100,9 @@ def __init__( self.lon1 = lon1 self.lon2 = lon2 self.reg_name = reg_name + self.create_mesh = create_mesh self.out_dir = out_dir self.create_tag() - self.create_mesh_at_reg() def create_tag(self): """ @@ -157,6 +159,7 @@ def create_surfdata_at_reg(self, indir, file, user_mods_dir): # specify files fsurf_in = os.path.join(indir, file) + print (self.tag) fsurf_out = add_tag_to_filename(fsurf_in, self.tag) logger.info("fsurf_in: %s", fsurf_in) logger.info("fsurf_out: %s", os.path.join(self.out_dir, fsurf_out)) @@ -231,13 +234,18 @@ def create_landuse_at_reg(self, indir, file, user_mods_dir): self.write_to_file(line, nl_clm) - def create_mesh_at_reg (self): + def create_mesh_at_reg(self, mesh_dir, mesh_surf): """ Create a mesh subsetted for the RegionalCase class. """ print ("test") - mesh_in = '/glade/p/cesmdata/cseg/inputdata/share/meshes/fv4x5_050615_polemod_ESMFmesh.nc' + mesh_in = os.path.join(mesh_dir, mesh_surf) + mesh_out = os.path.join(self.out_dir, os.path.splitext(mesh_surf)[0]+'_'+self.tag+'.nc') + print (mesh_out) + print (mesh_in) node_coords, subset_element, subset_node, conn_dict = self.subset_mesh_at_reg(mesh_in) + f_in = xr.open_dataset (mesh_in) + self.write_mesh (f_in, node_coords, subset_element, subset_node, conn_dict) def subset_mesh_at_reg (self, mesh_in): @@ -255,7 +263,7 @@ def subset_mesh_at_reg (self, mesh_in): subset_element = [] cnt = 0 - for n in tqdm.tqdm(range(elem_count)): + for n in tqdm(range(elem_count)): #print (elem_conn.shape) #print (num_elem_conn.shape) #print ('-----') @@ -271,8 +279,8 @@ def subset_mesh_at_reg (self, mesh_in): #print ('endx:', endx) #print (node_coords.shape) #print (node_coords) - nlon = node_coords[endx,0] - nlat = node_coords[endx,1] + nlon = node_coords[endx,0].values + nlat = node_coords[endx,1].values #print ('-----') l1 = np.logical_or(nlon <= self.lon1,nlon >= self.lon2) @@ -291,8 +299,8 @@ def subset_mesh_at_reg (self, mesh_in): cnt = 1 for n in range(node_count): - nlon = node_coords[n,0] - nlat = node_coords[n,1] + nlon = node_coords[n,0].values + nlat = node_coords[n,1].values l1 = np.logical_or(nlon <= self.lon1,nlon >= self.lon2) l2 = np.logical_or(nlat <= self.lat1,nlat >= self.lat2) @@ -307,24 +315,87 @@ def subset_mesh_at_reg (self, mesh_in): def write_mesh (self, f_in, node_coords, subset_element, subset_node, conn_dict): - corner_pair_uniq = f_in.variables['nodeCoords'][subset_node,] - - dimensions = f1.dimensions - variables = f1.variables - global_attributes = f1.ncattrs() - var = 'elementConn' - elem_conn_out = np.empty(shape=[elementCount, len(f1.dimensions['maxNodePElement'])]) - elem_conn_index = f1.variables[var][subsetElement,] - - print (elem_conn_out.shape) - for n in range(elementCount): - print ('n:',n) - for m in range(len(f1.dimensions['maxNodePElement'])): - print ('m:',m) - ndx = int (elem_conn_index[n,m]) - print ('ndx:', ndx) - elem_conn_out[n,m] = connDict[ndx] + corner_pairs = f_in.variables['nodeCoords'][subset_node,] + + dimensions = f_in.dims + variables = f_in.variables + global_attributes = f_in.attrs + + + max_node_dim = len(f_in['maxNodePElement']) + + elem_count = len(subset_element) + elem_conn_out = np.empty(shape=[elem_count, max_node_dim]) + elem_conn_index = f_in.variables['elementConn'][subset_element,] + + for ni in range(elem_count): + for mi in range(max_node_dim): + ndx = int (elem_conn_index[ni,mi]) + elem_conn_out[ni,mi] = conn_dict[ndx] + + + num_elem_conn_out = np.empty(shape=[elem_count,]) + num_elem_conn_out[:] = f_in.variables['numElementConn'][subset_element,] + + center_coords_out = np.empty(shape=[elem_count,2]) + center_coords_out[:,:]=f_in.variables['centerCoords'][subset_element,:] + + if 'elementMask' in variables: + elem_mask_out = np.empty(shape=[elem_count,]) + elem_mask_out[:]=f_in.variables['elementMask'][subset_element,] + + if 'elementArea' in variables: + elem_area_out = np.empty(shape=[elem_count,]) + elem_area_out[:]=f_in.variables['elementArea'][subset_element,] + + # -- create output dataset + f_out = xr.Dataset() + + f_out['nodeCoords'] = xr.DataArray(corner_pairs, + dims=('nodeCount', 'coordDim'), + attrs={'units': 'degrees'}) + + f_out['elementConn'] = xr.DataArray(elem_conn_out, + dims=('elementCount', 'maxNodePElement'), + attrs={'long_name': 'Node indices that define the element connectivity'}) + f_out.elementConn.encoding = {'dtype': np.int32} + + f_out['numElementConn'] = xr.DataArray(num_elem_conn_out, + dims=('elementCount'), + attrs={'long_name': 'Number of nodes per element'}) + + f_out['centerCoords'] = xr.DataArray(center_coords_out, + dims=('elementCount', 'coordDim'), + attrs={'units': 'degrees'}) + + + #-- add mask + f_out['elementMask'] = xr.DataArray(elem_mask_out, + dims=('elementCount'), + attrs={'units': 'unitless'}) + f_out.elementMask.encoding = {'dtype': np.int32} + + #-- setting fill values + for var in variables: + if '_FillValue' in f_in[var].encoding: + f_out[var].encoding['_FillValue'] = f_in[var].encoding['_FillValue'] + else: + f_out[var].encoding['_FillValue'] = None + + #-- add global attributes + + for attr in global_attributes: + if attr != 'timeGenerated': + f_out.attrs[attr] = global_attributes[attr] + f_out.attrs['timeGenerated'] = datetime.datetime.now() + #out.attrs = {'title': 'ESMF unstructured grid file for rectangular grid', + # 'created_by': os.path.basename(__file__), + # 'date_created': '{}'.format(datetime.now()), + # 'conventions': 'ESMFMESH', + # } + mesh_out = '/glade/scratch/negins/ctsm-mesh/mesh_out_test.nc' + out.to_netcdf(mesh_out) diff --git a/python/ctsm/subset_data.py b/python/ctsm/subset_data.py index a1b8d21c15..88147c9301 100644 --- a/python/ctsm/subset_data.py +++ b/python/ctsm/subset_data.py @@ -398,10 +398,13 @@ def setup_files(args, defaults, cesmroot): 'fdomain_in': defaults.get("domain", "file"), 'fsurf_dir': os.path.join(defaults.get("main", "clmforcingindir"), os.path.join(defaults.get("surfdat", "dir"))), + 'mesh_dir': os.path.join(defaults.get("main", "clmforcingindir"), + defaults.get("surfdat", "mesh_dir")), 'fluse_dir': os.path.join(defaults.get("main", "clmforcingindir"), os.path.join(defaults.get("landuse", "dir"))), 'fsurf_in': fsurf_in, 'fluse_in': fluse_in, + 'mesh_surf' : defaults.get("surfdat","mesh_surf"), 'datm_tuple': DatmFiles(dir_input_datm, dir_output_datm, defaults.get(datm_type, "domain"), @@ -415,7 +418,6 @@ def setup_files(args, defaults, cesmroot): defaults.get(datm_type, 'precname'), defaults.get(datm_type, 'tpqwname')) } - return file_dict @@ -502,6 +504,7 @@ def subset_region(args, file_dict: dict): create_landuse = args.create_landuse, create_datm = args.create_datm, create_user_mods = args.create_user_mods, + create_mesh = args.create_mesh, out_dir = args.out_dir, overwrite = args.overwrite, ) @@ -517,6 +520,9 @@ def subset_region(args, file_dict: dict): region.create_surfdata_at_reg(file_dict["fsurf_dir"], file_dict["fsurf_in"], args.user_mods_dir) + if region.create_mesh: + region.create_mesh_at_reg (file_dict["mesh_dir"], file_dict["mesh_surf"]) + # -- Create CTSM transient landuse data file if region.create_landuse: region.create_landuse_at_reg(file_dict["fluse_dir"], file_dict["fluse_in"], diff --git a/tools/site_and_regional/default_data.cfg b/tools/site_and_regional/default_data.cfg index f689c99044..630013d502 100644 --- a/tools/site_and_regional/default_data.cfg +++ b/tools/site_and_regional/default_data.cfg @@ -18,6 +18,8 @@ tpqwname = CLMGSWP3v1.TPQW dir = lnd/clm2/surfdata_map/release-clm5.0.18 surfdat_16pft = surfdata_0.9x1.25_hist_16pfts_Irrig_CMIP6_simyr2000_c190214.nc surfdat_78pft = surfdata_0.9x1.25_hist_78pfts_CMIP6_simyr2000_c190214.nc +mesh_dir = share/meshes/ +mesh_surf = fv0.9x1.25_141008_ESMFmesh.nc [landuse] dir = lnd/clm2/surfdata_map/release-clm5.0.18 From 66f46c393c4c945b6427991ac44c82b347b6c0fc Mon Sep 17 00:00:00 2001 From: Negin Sobhani Date: Fri, 3 Jun 2022 10:53:27 -0600 Subject: [PATCH 3/9] a few changes --- .../ctsm/site_and_regional/regional_case.py | 39 +++++++++++-------- 1 file changed, 23 insertions(+), 16 deletions(-) diff --git a/python/ctsm/site_and_regional/regional_case.py b/python/ctsm/site_and_regional/regional_case.py index b497b49d45..40fad6efee 100644 --- a/python/ctsm/site_and_regional/regional_case.py +++ b/python/ctsm/site_and_regional/regional_case.py @@ -160,8 +160,10 @@ def create_surfdata_at_reg(self, indir, file, user_mods_dir): # specify files fsurf_in = os.path.join(indir, file) print (self.tag) + print (fsurf_in) fsurf_out = add_tag_to_filename(fsurf_in, self.tag) logger.info("fsurf_in: %s", fsurf_in) + print (fsurf_out) logger.info("fsurf_out: %s", os.path.join(self.out_dir, fsurf_out)) # create 1d coordinate variables to enable sel() method @@ -238,14 +240,14 @@ def create_mesh_at_reg(self, mesh_dir, mesh_surf): """ Create a mesh subsetted for the RegionalCase class. """ - print ("test") + print ("Creating mesh for this: .....") mesh_in = os.path.join(mesh_dir, mesh_surf) mesh_out = os.path.join(self.out_dir, os.path.splitext(mesh_surf)[0]+'_'+self.tag+'.nc') print (mesh_out) print (mesh_in) node_coords, subset_element, subset_node, conn_dict = self.subset_mesh_at_reg(mesh_in) f_in = xr.open_dataset (mesh_in) - self.write_mesh (f_in, node_coords, subset_element, subset_node, conn_dict) + self.write_mesh (f_in, node_coords, subset_element, subset_node, conn_dict, mesh_out) def subset_mesh_at_reg (self, mesh_in): @@ -314,7 +316,7 @@ def subset_mesh_at_reg (self, mesh_in): - def write_mesh (self, f_in, node_coords, subset_element, subset_node, conn_dict): + def write_mesh (self, f_in, node_coords, subset_element, subset_node, conn_dict, mesh_out): corner_pairs = f_in.variables['nodeCoords'][subset_node,] dimensions = f_in.dims @@ -370,32 +372,37 @@ def write_mesh (self, f_in, node_coords, subset_element, subset_node, conn_dict) #-- add mask - f_out['elementMask'] = xr.DataArray(elem_mask_out, - dims=('elementCount'), - attrs={'units': 'unitless'}) - f_out.elementMask.encoding = {'dtype': np.int32} + if 'elementMask' in variables: + f_out['elementMask'] = xr.DataArray(elem_mask_out, + dims=('elementCount'), + attrs={'units': 'unitless'}) + print (elem_mask_out) + f_out.elementMask.encoding = {'dtype': np.int32} + + if 'elementArea' in variables: + f_out['elementArea'] = xr.DataArray(elem_area_out, + dims=('elementCount'), + attrs={'units': 'unitless'}) #-- setting fill values for var in variables: + print (var) if '_FillValue' in f_in[var].encoding: f_out[var].encoding['_FillValue'] = f_in[var].encoding['_FillValue'] + print ('FillValue') else: f_out[var].encoding['_FillValue'] = None #-- add global attributes - for attr in global_attributes: if attr != 'timeGenerated': f_out.attrs[attr] = global_attributes[attr] - f_out.attrs['timeGenerated'] = datetime.datetime.now() - #out.attrs = {'title': 'ESMF unstructured grid file for rectangular grid', - # 'created_by': os.path.basename(__file__), - # 'date_created': '{}'.format(datetime.now()), - # 'conventions': 'ESMFMESH', - # } + f_out.attrs = {'title': 'ESMF unstructured grid file for a region', + 'created_by': 'subset_data', + 'date_created': '{}'.format(datetime.now()), + } - mesh_out = '/glade/scratch/negins/ctsm-mesh/mesh_out_test.nc' - out.to_netcdf(mesh_out) + f_out.to_netcdf(mesh_out) From 29ed6fab54e3dba0a3cd3d67872fb37e874655f2 Mon Sep 17 00:00:00 2001 From: Negin Sobhani Date: Wed, 8 Jun 2022 00:19:01 -0600 Subject: [PATCH 4/9] clean ups --- .../ctsm/site_and_regional/regional_case.py | 46 +++++++++---------- 1 file changed, 22 insertions(+), 24 deletions(-) diff --git a/python/ctsm/site_and_regional/regional_case.py b/python/ctsm/site_and_regional/regional_case.py index 40fad6efee..97eb17ee12 100644 --- a/python/ctsm/site_and_regional/regional_case.py +++ b/python/ctsm/site_and_regional/regional_case.py @@ -159,11 +159,8 @@ def create_surfdata_at_reg(self, indir, file, user_mods_dir): # specify files fsurf_in = os.path.join(indir, file) - print (self.tag) - print (fsurf_in) fsurf_out = add_tag_to_filename(fsurf_in, self.tag) logger.info("fsurf_in: %s", fsurf_in) - print (fsurf_out) logger.info("fsurf_out: %s", os.path.join(self.out_dir, fsurf_out)) # create 1d coordinate variables to enable sel() method @@ -240,12 +237,26 @@ def create_mesh_at_reg(self, mesh_dir, mesh_surf): """ Create a mesh subsetted for the RegionalCase class. """ - print ("Creating mesh for this: .....") + logger.info( + "----------------------------------------------------------------------" + ) + logger.info( + "Subsetting mesh file for region: %s", + self.tag + ) + + today = datetime.today() + today_string = today.strftime("%y%m%d") + + mesh_in = os.path.join(mesh_dir, mesh_surf) - mesh_out = os.path.join(self.out_dir, os.path.splitext(mesh_surf)[0]+'_'+self.tag+'.nc') - print (mesh_out) - print (mesh_in) + mesh_out = os.path.join(self.out_dir, os.path.splitext(mesh_surf)[0]+'_'+self.tag+'_c'+today_string+'.nc') + + logger.info("mesh_in : %s", mesh_in) + logger.info("mesh_out : %s", mesh_out) + node_coords, subset_element, subset_node, conn_dict = self.subset_mesh_at_reg(mesh_in) + f_in = xr.open_dataset (mesh_in) self.write_mesh (f_in, node_coords, subset_element, subset_node, conn_dict, mesh_out) @@ -266,35 +277,21 @@ def subset_mesh_at_reg (self, mesh_in): cnt = 0 for n in tqdm(range(elem_count)): - #print (elem_conn.shape) - #print (num_elem_conn.shape) - #print ('-----') - #print (numElementConn[n]) endx = elem_conn[n,:num_elem_conn[n].values].values - #print ('endx:', endx) - #print ('-----') endx[:,] -= 1# convert to zero based index endx = [int(xi) for xi in endx] - #print ('-----') - #endx = endx.values - #print ('endx:', endx) - #print (node_coords.shape) - #print (node_coords) nlon = node_coords[endx,0].values nlat = node_coords[endx,1].values - #print ('-----') l1 = np.logical_or(nlon <= self.lon1,nlon >= self.lon2) l2 = np.logical_or(nlat <= self.lat1,nlat >= self.lat2) if np.any(np.logical_or(l1,l2)): - #print(nlon,nlat) pass else: subset_element.append(n) cnt+=1 - print (subset_element) subset_node = [] conn_dict = {} @@ -317,6 +314,9 @@ def subset_mesh_at_reg (self, mesh_in): def write_mesh (self, f_in, node_coords, subset_element, subset_node, conn_dict, mesh_out): + """ + This function writes out the subsetted mesh file. + """ corner_pairs = f_in.variables['nodeCoords'][subset_node,] dimensions = f_in.dims @@ -376,7 +376,6 @@ def write_mesh (self, f_in, node_coords, subset_element, subset_node, conn_dict, f_out['elementMask'] = xr.DataArray(elem_mask_out, dims=('elementCount'), attrs={'units': 'unitless'}) - print (elem_mask_out) f_out.elementMask.encoding = {'dtype': np.int32} if 'elementArea' in variables: @@ -386,10 +385,8 @@ def write_mesh (self, f_in, node_coords, subset_element, subset_node, conn_dict, #-- setting fill values for var in variables: - print (var) if '_FillValue' in f_in[var].encoding: f_out[var].encoding['_FillValue'] = f_in[var].encoding['_FillValue'] - print ('FillValue') else: f_out[var].encoding['_FillValue'] = None @@ -404,5 +401,6 @@ def write_mesh (self, f_in, node_coords, subset_element, subset_node, conn_dict, } f_out.to_netcdf(mesh_out) + logger.info("Successfully created file (mesh_out) %s", mesh_out) From 9155ec8552991bff85efec80701070a6a8d2a7aa Mon Sep 17 00:00:00 2001 From: Negin Sobhani Date: Wed, 8 Jun 2022 00:48:11 -0600 Subject: [PATCH 5/9] adding write shell commands --- python/ctsm/site_and_regional/regional_case.py | 18 +++++++++++++++++- python/ctsm/subset_data.py | 15 +++++++++++++++ 2 files changed, 32 insertions(+), 1 deletion(-) diff --git a/python/ctsm/site_and_regional/regional_case.py b/python/ctsm/site_and_regional/regional_case.py index 97eb17ee12..f1d1d40552 100644 --- a/python/ctsm/site_and_regional/regional_case.py +++ b/python/ctsm/site_and_regional/regional_case.py @@ -255,6 +255,9 @@ def create_mesh_at_reg(self, mesh_dir, mesh_surf): logger.info("mesh_in : %s", mesh_in) logger.info("mesh_out : %s", mesh_out) + self.mesh = mesh_out + print (self) + node_coords, subset_element, subset_node, conn_dict = self.subset_mesh_at_reg(mesh_in) f_in = xr.open_dataset (mesh_in) @@ -403,4 +406,17 @@ def write_mesh (self, f_in, node_coords, subset_element, subset_node, conn_dict, f_out.to_netcdf(mesh_out) logger.info("Successfully created file (mesh_out) %s", mesh_out) - + def write_shell_commands(self, namelist): + """ + writes out xml commands commands to a file (i.e. shell_commands) for single-point runs + """ + # write_to_file surrounds text with newlines + with open(namelist, "w") as nl_file: + self.write_to_file( + "# Change below line if you move the subset data directory", nl_file + ) + self.write_to_file( + "./xmlchange {}={}".format(USRDAT_DIR, self.out_dir), nl_file + ) + self.write_to_file("./xmlchange ATM_DOMAIN_MESH={}".format(str(self.mesh)), nl_file) + self.write_to_file("./xmlchange LAND_DOMAIN_MESH={}".format(str(self.mesh)), nl_file) diff --git a/python/ctsm/subset_data.py b/python/ctsm/subset_data.py index 88147c9301..9e1a704788 100644 --- a/python/ctsm/subset_data.py +++ b/python/ctsm/subset_data.py @@ -528,6 +528,21 @@ def subset_region(args, file_dict: dict): region.create_landuse_at_reg(file_dict["fluse_dir"], file_dict["fluse_in"], args.user_mods_dir) + + # -- Write shell commands + if region.create_user_mods: + if not region.create_mesh: + err_msg = """ + \n + ERROR: For regional cases, you can not create user_mods + without creating the mesh file. + + Please rerun the script adding --create-mesh to subset the mesh file. + """ + raise argparse.ArgumentTypeError(err_msg) + + region.write_shell_commands(os.path.join(args.user_mods_dir, "shell_commands")) + logger.info("Successfully ran script for a regional case.") From 18c9ca7ced8ba942b2517a457fee60e9cff37017 Mon Sep 17 00:00:00 2001 From: Negin Sobhani Date: Wed, 8 Jun 2022 12:15:56 -0600 Subject: [PATCH 6/9] clean ups --- python/ctsm/site_and_regional/regional_case.py | 1 - 1 file changed, 1 deletion(-) diff --git a/python/ctsm/site_and_regional/regional_case.py b/python/ctsm/site_and_regional/regional_case.py index f1d1d40552..ee7c0cd7b5 100644 --- a/python/ctsm/site_and_regional/regional_case.py +++ b/python/ctsm/site_and_regional/regional_case.py @@ -256,7 +256,6 @@ def create_mesh_at_reg(self, mesh_dir, mesh_surf): logger.info("mesh_out : %s", mesh_out) self.mesh = mesh_out - print (self) node_coords, subset_element, subset_node, conn_dict = self.subset_mesh_at_reg(mesh_in) From b8c91151c8ca934101ed4f7753060f170c1d6849 Mon Sep 17 00:00:00 2001 From: Negin Sobhani Date: Tue, 28 Jun 2022 19:17:07 -0600 Subject: [PATCH 7/9] adding bound checks --- .../ctsm/site_and_regional/regional_case.py | 55 +++++++++++++++++++ 1 file changed, 55 insertions(+) diff --git a/python/ctsm/site_and_regional/regional_case.py b/python/ctsm/site_and_regional/regional_case.py index ee7c0cd7b5..89751fec55 100644 --- a/python/ctsm/site_and_regional/regional_case.py +++ b/python/ctsm/site_and_regional/regional_case.py @@ -5,6 +5,7 @@ # -- Import Python Standard Libraries import logging import os +import argparse # -- 3rd party libraries import numpy as np @@ -52,6 +53,15 @@ class RegionalCase(BaseCase): Methods ------- + check_region_bounds + Check for the regional bounds + + check_region_lons + Check for the regional lons + + check_region_lats + Check for the regional lats + create_tag Create a tag for this region which is either region's name or a combination of bounds of this @@ -102,6 +112,7 @@ def __init__( self.reg_name = reg_name self.create_mesh = create_mesh self.out_dir = out_dir + self.check_region_bounds() self.create_tag() def create_tag(self): @@ -117,6 +128,50 @@ def create_tag(self): str(self.lon1), str(self.lon2), str(self.lat1), str(self.lat2) ) + def check_region_bounds (self): + """ + Check for the regional bounds + """ + self.check_region_lons() + self.check_region_lats() + + def check_region_lons (self): + """ + Check for the regional lon bounds + """ + if self.lon1 >= self.lon2: + err_msg = """ + \n + ERROR: lon1 is bigger than lon2. + lon1 points to the westernmost longitude of the region. {} + lon2 points to the easternmost longitude of the region. {} + Please make sure lon1 is smaller than lon2. + + Please note that if longitude in -180-0, the code automatically + convert it to 0-360. + """.format( + self.lon1, self.lon2 + ) + raise argparse.ArgumentTypeError(err_msg) + + def check_region_lats (self): + """ + Check for the regional lat bound + """ + if self.lat1 >= self.lat2: + err_msg = """ + \n + ERROR: lat1 is bigger than lat2. + lat1 points to the westernmost longitude of the region. {} + lat2 points to the easternmost longitude of the region. {} + Please make sure lat1 is smaller than lat2. + + """.format( + self.lat1, self.lat2 + ) + raise argparse.ArgumentTypeError(err_msg) + + def create_domain_at_reg(self, indir, file): """ Create domain file for this RegionalCase class. From 76821859f888763d7a37761cc472fd8e3f5a2d58 Mon Sep 17 00:00:00 2001 From: Negin Sobhani Date: Thu, 30 Jun 2022 11:29:05 -0600 Subject: [PATCH 8/9] adding w/e s/n instead of start/end to help --- python/ctsm/subset_data.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/python/ctsm/subset_data.py b/python/ctsm/subset_data.py index 9e1a704788..7230a72a35 100644 --- a/python/ctsm/subset_data.py +++ b/python/ctsm/subset_data.py @@ -170,7 +170,7 @@ def get_parser(): # -- region-specific parser options rg_parser.add_argument( "--lat1", - help="Region start latitude. [default: %(default)s]", + help="Region southernmost latitude. [default: %(default)s]", action="store", dest="lat1", required=False, @@ -179,7 +179,7 @@ def get_parser(): ) rg_parser.add_argument( "--lat2", - help="Region end latitude. [default: %(default)s]", + help="Region northernmost latitude. [default: %(default)s]", action="store", dest="lat2", required=False, @@ -188,7 +188,7 @@ def get_parser(): ) rg_parser.add_argument( "--lon1", - help="Region start longitude. [default: %(default)s]", + help="Region westernmost longitude. [default: %(default)s]", action="store", dest="lon1", required=False, @@ -197,7 +197,7 @@ def get_parser(): ) rg_parser.add_argument( "--lon2", - help="Region end longitude. [default: %(default)s]", + help="Region easternmost longitude. [default: %(default)s]", action="store", dest="lon2", required=False, From 60b232f03c99189d4fb30c4f5cc69e2056a24410 Mon Sep 17 00:00:00 2001 From: Negin Sobhani Date: Thu, 30 Jun 2022 11:32:22 -0600 Subject: [PATCH 9/9] lnd for land_domain_mesh --- python/ctsm/site_and_regional/regional_case.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/python/ctsm/site_and_regional/regional_case.py b/python/ctsm/site_and_regional/regional_case.py index 89751fec55..92892b9340 100644 --- a/python/ctsm/site_and_regional/regional_case.py +++ b/python/ctsm/site_and_regional/regional_case.py @@ -473,4 +473,4 @@ def write_shell_commands(self, namelist): "./xmlchange {}={}".format(USRDAT_DIR, self.out_dir), nl_file ) self.write_to_file("./xmlchange ATM_DOMAIN_MESH={}".format(str(self.mesh)), nl_file) - self.write_to_file("./xmlchange LAND_DOMAIN_MESH={}".format(str(self.mesh)), nl_file) + self.write_to_file("./xmlchange LND_DOMAIN_MESH={}".format(str(self.mesh)), nl_file)