From 0ccf8551dc108a5f02d96eb5d0bd4380db103d19 Mon Sep 17 00:00:00 2001 From: Negin Sobhani Date: Thu, 28 Apr 2022 09:13:40 -0600 Subject: [PATCH 001/595] 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 002/595] 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 003/595] 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 004/595] 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 005/595] 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 006/595] 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 007/595] 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 008/595] 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 009/595] 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) From 6aef680d6084cd18beb314b3aaa24a2525bcb445 Mon Sep 17 00:00:00 2001 From: Negin Sobhani Date: Thu, 27 Oct 2022 13:12:15 -0600 Subject: [PATCH 010/595] saving progress for old debugging --- .../ctsm/site_and_regional/regional_case.py | 47 +++++++++++++------ 1 file changed, 33 insertions(+), 14 deletions(-) diff --git a/python/ctsm/site_and_regional/regional_case.py b/python/ctsm/site_and_regional/regional_case.py index 92892b9340..16669ffb44 100644 --- a/python/ctsm/site_and_regional/regional_case.py +++ b/python/ctsm/site_and_regional/regional_case.py @@ -337,35 +337,50 @@ def subset_mesh_at_reg (self, mesh_in): endx = elem_conn[n,:num_elem_conn[n].values].values endx[:,] -= 1# convert to zero based index endx = [int(xi) for xi in endx] - + #print (endx) nlon = node_coords[endx,0].values nlat = node_coords[endx,1].values - 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)): - pass - else: + l1 = np.logical_and(nlon > self.lon1,nlon < self.lon2) + l2 = np.logical_and(nlat > self.lat1,nlat < self.lat2) + if np.any(l1) and np.any(l2): + #if np.any(np.logical_or(l1,l2)): + # pass + #else: subset_element.append(n) cnt+=1 - + print (cnt) subset_node = [] conn_dict = {} cnt = 1 - + print (node_count) + print (elem_count) + print (node_coords) for n in range(node_count): + #n = 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) - if np.logical_or(l1,l2): - conn_dict[n+1] = -9999 - else: + #l1 = np.logical_or(nlon <= self.lon1,nlon >= self.lon2) + #l2 = np.logical_or(nlat <= self.lat1,nlat >= self.lat2) + l1 = np.logical_and(nlon >= self.lon1,nlon <= self.lon2) + l2 = np.logical_and(nlat >= self.lat1,nlat <= self.lat2) + if np.any(l1) and np.any(l2): subset_node.append(n) conn_dict[n+1] = cnt cnt+=1 + else: + conn_dict[n+1] = -9999 + + #if np.logical_or(l1,l2): + # conn_dict[n+1] = -9999 + #if np.any(l1) and np.any(l2): + #else: + # subset_node.append(n) + # conn_dict[n+1] = cnt + # cnt+=1 + print ('cnt:',cnt) return node_coords, subset_element, subset_node, conn_dict @@ -374,8 +389,10 @@ def write_mesh (self, f_in, node_coords, subset_element, subset_node, conn_dict, """ This function writes out the subsetted mesh file. """ + print (len(subset_node)) corner_pairs = f_in.variables['nodeCoords'][subset_node,] - + print (corner_pairs) + #print (corner_pairs.shape()) dimensions = f_in.dims variables = f_in.variables global_attributes = f_in.attrs @@ -390,6 +407,7 @@ def write_mesh (self, f_in, node_coords, subset_element, subset_node, conn_dict, for ni in range(elem_count): for mi in range(max_node_dim): ndx = int (elem_conn_index[ni,mi]) + #print ('mi',mi, 'ndx',ndx) elem_conn_out[ni,mi] = conn_dict[ndx] @@ -422,6 +440,7 @@ def write_mesh (self, f_in, node_coords, subset_element, subset_node, conn_dict, f_out['numElementConn'] = xr.DataArray(num_elem_conn_out, dims=('elementCount'), attrs={'long_name': 'Number of nodes per element'}) + f_out.numElementConn.encoding = {'dtype': np.int32} f_out['centerCoords'] = xr.DataArray(center_coords_out, dims=('elementCount', 'coordDim'), From 88ed680536a59bc13d18ff6046d5b50e6d328ff2 Mon Sep 17 00:00:00 2001 From: Negin Sobhani Date: Fri, 28 Oct 2022 18:35:46 -0600 Subject: [PATCH 011/595] working mesh type class --- python/ctsm/site_and_regional/mesh_type.py | 370 +++++++++++++++++++++ 1 file changed, 370 insertions(+) create mode 100644 python/ctsm/site_and_regional/mesh_type.py diff --git a/python/ctsm/site_and_regional/mesh_type.py b/python/ctsm/site_and_regional/mesh_type.py new file mode 100644 index 0000000000..640f5548b5 --- /dev/null +++ b/python/ctsm/site_and_regional/mesh_type.py @@ -0,0 +1,370 @@ +""" +This module includes the definition and functions for defining a Mesh type. +""" + +import logging +import datetime +import numpy as np +import xarray as xr +import pandas as pd +import dask.array as da +import dask.dataframe as dd +from dask.diagnostics import ProgressBar + +logger = logging.getLogger(__name__) + +class MeshType: + """ + An object for describing mesh or grid files. + """ + + def __init__(self, center_lats, center_lons, mesh_name=None, mask=None): + """ + Construct a mesh object + + center_lats : + latitudes (either 1D or 2D) + center_lons : + longitudes (either 1D or 2D) + mesh_name : str or None, optional + Name of the mesh + mask : np array or None + numpy array that include landmask + """ + + self.mesh_name = mesh_name + self.center_lats = center_lats + self.center_lons = center_lons + + # -- dims of lat and lon (1d or 2d) + self.lat_dims = len(self.center_lats.dims) + self.lon_dims = len(self.center_lons.dims) + self.check_lat_lon_dims() + + if mask is None: + self.create_artificial_mask() + else: + self.mask = mask + + def __str__(self): + """ + Converts ingredients of this class to string for printing. + """ + return "{}\n{}".format( + str(self.__class__), + "\n".join( + ( + "{} = {}".format(str(key), str(self.__dict__[key])) + for key in sorted(self.__dict__) + ) + ), + ) + + def check_lat_lon_dims(self): + """ + Check latitude and longitude dimensions to make sure they are either 1 or 2. + """ + if self.lat_dims not in [1, 2]: + print( + "Unrecognized grid! The dimension of latitude should be either 1 or 2 but it is {}.".format( + self.lat_dims + ) + ) + if self.lon_dims not in [1, 2]: + print( + "Unrecognized grid! The dimension of longitude should be either 1 or 2 but it is {}.".format( + self.lon_dims + ) + ) + + def create_artificial_mask(self): + logger.info("Creating an artificial mask for this region...") + + if self.lat_dims == 1: + # -- 1D mask (lat x lon) + lats_size = self.center_lats.size + lons_size = self.center_lons.size + mask = np.ones([lons_size, lats_size], dtype=np.int8) + elif self.lat_dims == 2: + # -- 2D mask + mask = np.ones( + center_lats.shape, dtype=np.int8 + ) # np.ones(tuple(self.center_lats.sizes.values()), dtype=np.int8) # + mask_da = da.from_array(mask) + self.mask = mask_da + + def create_2d_coords(self): + """ + Create 2d center points for our mesh + and convert them to Dask Array. + """ + self.center_lats = self.center_lats.astype(np.float64, copy=False) + self.center_lons = self.center_lons.astype(np.float64, copy=False) + + if self.lat_dims == 1: + # -- 1D lats and lons + lats_size = self.center_lats.size + lons_size = self.center_lons.size + + # -- convert center points from 1d to 2d + self.center_lat2d = da.broadcast_to( + self.center_lats.values[None, :], (lons_size, lats_size) + ) + self.center_lon2d = da.broadcast_to( + self.center_lons.values[:, None], (lons_size, lats_size) + ) + elif self.lat_dims == 2: + # -- 2D lats and lons + dims = self.center_lons.shape + + # -- convert 2D lats and lons to number x and y + lons_size = dims[0] + lats_size = dims[1] + + # -- convert to dask array + self.center_lat2d = da.from_array(self.center_lats) + self.center_lon2d = da.from_array(self.center_lons) + + def calculate_corners(self, unit="degrees"): + """ + calculate corner coordinates by averaging adjacent cells + + Parameters + ---------- + unit : {'degrees', 'radians'}, optional + The unit of corner coordinates. + """ + + self.create_2d_coords() + # -- pad center_lats for calculating edge gridpoints + # -- otherwise we cannot calculate the corner coords + # -- for the edge rows/columns. + + padded_lat2d = da.from_array( + np.pad( + self.center_lat2d.compute(), (1, 1), mode="reflect", reflect_type="odd" + ) + ) + + # -- pad center_lons for calculating edge grids + padded_lon2d = da.from_array( + np.pad( + self.center_lon2d.compute(), (1, 1), mode="reflect", reflect_type="odd" + ) + ) + + # -- calculate corner lats for each grid + north_east = ( + padded_lat2d[1:-1, 1:-1] + + padded_lat2d[0:-2, 1:-1] + + padded_lat2d[1:-1, 2:] + + padded_lat2d[0:-2, 2:] + ) / 4.0 + north_west = ( + padded_lat2d[1:-1, 1:-1] + + padded_lat2d[0:-2, 1:-1] + + padded_lat2d[1:-1, 0:-2] + + padded_lat2d[0:-2, 0:-2] + ) / 4.0 + south_west = ( + padded_lat2d[1:-1, 1:-1] + + padded_lat2d[1:-1, 0:-2] + + padded_lat2d[2:, 1:-1] + + padded_lat2d[2:, 0:-2] + ) / 4.0 + south_east = ( + padded_lat2d[1:-1, 1:-1] + + padded_lat2d[1:-1, 2:] + + padded_lat2d[2:, 1:-1] + + padded_lat2d[2:, 2:] + ) / 4.0 + + # -- order counter-clockwise + self.corner_lats = da.stack( + [ + north_west.T.reshape((-1,)).T, + south_west.T.reshape((-1,)).T, + south_east.T.reshape((-1,)).T, + north_east.T.reshape((-1,)).T, + ], + axis=1, + ) + + # -- calculate corner lons for each grid + north_east = ( + padded_lon2d[1:-1, 1:-1] + + padded_lon2d[0:-2, 1:-1] + + padded_lon2d[1:-1, 2:] + + padded_lon2d[0:-2, 2:] + ) / 4.0 + north_west = ( + padded_lon2d[1:-1, 1:-1] + + padded_lon2d[0:-2, 1:-1] + + padded_lon2d[1:-1, 0:-2] + + padded_lon2d[0:-2, 0:-2] + ) / 4.0 + south_west = ( + padded_lon2d[1:-1, 1:-1] + + padded_lon2d[1:-1, 0:-2] + + padded_lon2d[2:, 1:-1] + + padded_lon2d[2:, 0:-2] + ) / 4.0 + south_east = ( + padded_lon2d[1:-1, 1:-1] + + padded_lon2d[1:-1, 2:] + + padded_lon2d[2:, 1:-1] + + padded_lon2d[2:, 2:] + ) / 4.0 + + # -- order counter-clockwise + self.corner_lons = da.stack( + [ + north_west.T.reshape((-1,)).T, + south_west.T.reshape((-1,)).T, + south_east.T.reshape((-1,)).T, + north_east.T.reshape((-1,)).T, + ], + axis=1, + ) + self.unit = unit + + def calculate_node_coords(self): + """ + Calculates coordinates of each node (for 'nodeCoords' in ESMF mesh). + In ESMF mesh, 'nodeCoords' is a two-dimensional array with dimension ('nodeCount','coordDim') + """ + # -- create an array of corner pairs + corner_pairs = da.stack( + [self.corner_lons.T.reshape((-1,)).T, self.corner_lats.T.reshape((-1,)).T], + axis=1, + ) + + # -- remove coordinates that are shared between the elements + node_coords = dd.from_dask_array(corner_pairs).drop_duplicates().values + node_coords.compute_chunk_sizes() + # -- check size of unique coordinate pairs + dims = self.mask.shape + nlon = dims[0] + nlat = dims[1] + elem_conn_size = nlon * nlat + nlon + nlat + 1 + self.node_coords = node_coords + + if self.node_coords.shape[0] != elem_conn_size: + logger.warning( + "The size of unique coordinate pairs is {} but expected size is {}!".format( + self.node_coords.shape[0], elem_conn_size + ) + ) + sys.exit(2) + + def calculate_elem_conn(self): + """ + Calculate element connectivity (for 'elementConn' in ESMF mesh). + In ESMF mesh, 'elementConn' describes how the nodes are connected together. + """ + corners = dd.concat( + [ + dd.from_dask_array(corner) + for corner in [ + self.corner_lons.T.reshape((-1,)).T, + self.corner_lats.T.reshape((-1,)).T, + ] + ], + axis=1, + ) + corners.columns = ["lon", "lat"] + + elem_conn = corners.compute().groupby(["lon", "lat"], sort=False).ngroup() + 1 + elem_conn = da.from_array(elem_conn.to_numpy()) + # -- reshape to write to ESMF + self.elem_conn = elem_conn.T.reshape((4, -1)).T + + def create_esmf(self, mesh_fname, area=None): + """ + Create an ESMF mesh file for the mesh + + Parameters + ---------- + mesh_fname : str + The path to write the ESMF meshfile + + area : numpy.ndarray or None + Array containing element areas for the ESMF mesh file + If None, ESMF calculates element areas internally. + """ + # -- calculate node coordinates + self.calculate_node_coords() + + # -- calculate element connections + self.calculate_elem_conn() + + center_coords = da.stack( + [ + self.center_lon2d.T.reshape((-1,)).T, + self.center_lat2d.T.reshape((-1,)).T, + ], + axis=1, + ) + # create output Xarray dataset + ds_out = xr.Dataset() + + ds_out["origGridDims"] = xr.DataArray( + np.array(self.center_lon2d.shape, dtype=np.int32), dims=("origGridRank") + ) + ds_out["nodeCoords"] = xr.DataArray( + self.node_coords, dims=("nodeCount", "coordDim"), attrs={"units": self.unit} + ) + ds_out["elementConn"] = xr.DataArray( + self.elem_conn, + dims=("elementCount", "maxNodePElement"), + attrs={ + "long_name": "Node indices that define the element connectivity", + "_FillValue": -1, + }, + ) + ds_out.elementConn.encoding = {"dtype": np.int32} + + ds_out["numElementConn"] = xr.DataArray( + 4 * np.ones(self.center_lon2d.size, dtype=np.int32), + dims=("elementCount"), + attrs={"long_name": "Number of nodes per element"}, + ) + ds_out["centerCoords"] = xr.DataArray( + center_coords, dims=("elementCount", "coordDim"), attrs={"units": self.unit} + ) + + # -- add mask + ds_out["elementMask"] = xr.DataArray( + self.mask.T.reshape((-1,)).T, + dims=("elementCount"), + attrs={"units": "unitless", "_FillValue": -9999.0}, + ) + ds_out.elementMask.encoding = {"dtype": np.int32} + + # -- add area if provided + if area: + ds_out["elementArea"] = xr.DataArray( + area.T.reshape((-1,)).T, + dims=("elementCount"), + attrs={"units": "radians^2", "long_name": "area weights"}, + ) + + # -- force no '_FillValue' if not specified (default Nan) + for var in ds_out.variables: + if "_FillValue" not in ds_out[var].encoding: + ds_out[var].encoding["_FillValue"] = None + + # -- add global attributes + ds_out.attrs["title"] = "ESMF unstructured grid file" + ds_out.attrs["gridType"] = "unstructured mesh" + ds_out.attrs["version"] = "0.9" + ds_out.attrs["conventions"] = "ESMFMESH" + ds_out.attrs["date_created"] = datetime.datetime.now().strftime( + "%Y-%m-%d %H:%M:%S" + ) + + # -- write Xarray dataset to file + if mesh_fname is not None: + logger.info("Writing ESMF Mesh file to : %s", mesh_fname) + ds_out.to_netcdf(mesh_fname) + logger.info("Successfully created ESMF Mesh file : %s", mesh_fname) From 8e88f66fc065bb9605cb177bb51e550f8677fd7c Mon Sep 17 00:00:00 2001 From: Negin Sobhani Date: Fri, 28 Oct 2022 18:39:55 -0600 Subject: [PATCH 012/595] mesh_maker --- python/ctsm/site_and_regional/mesh_maker.py | 52 +++++++++++++++++++++ 1 file changed, 52 insertions(+) create mode 100644 python/ctsm/site_and_regional/mesh_maker.py diff --git a/python/ctsm/site_and_regional/mesh_maker.py b/python/ctsm/site_and_regional/mesh_maker.py new file mode 100644 index 0000000000..4922f17b42 --- /dev/null +++ b/python/ctsm/site_and_regional/mesh_maker.py @@ -0,0 +1,52 @@ +from mesh_type import MeshType +import xarray as xr +import os +oned = True +if oned: + ifile = "/glade/scratch/negins/this_region_4x5_new_2/surfdata_4x5_hist_78pfts_CMIP6_simyr1850_275.0-330.0_-40-15_c220705.nc" +else: + ifile = "/glade/scratch/negins/wrf-ctsm_production/WRF/test/em_real_sim1_noahmp/wrfinput_d01" + +import os, sys, getopt + +if os.path.isfile(ifile): + ds = xr.open_dataset(ifile, mask_and_scale=False, decode_times=False).transpose() +else: + print("Input file could not find!") + sys.exit(2) + + import os, sys, getopt + + if os.path.isfile(ifile): + ds = xr.open_dataset( + ifile, mask_and_scale=False, decode_times=False + ).transpose() + else: + print("Input file could not find!") + sys.exit(2) + +if oned: + lat_name = "lsmlat" + lon_name = "lsmlon" +else: + lat_name = "XLAT" + lon_name = "XLONG" + +lats = ds[lat_name] +lons = ds[lon_name] + +hasTime = "Time" in ds[lat_name].dims +print(hasTime) + +if hasTime: + lats = lats[:, :, 0] + lons = lons[:, :, 0] + + +this_mesh = MeshType(lats, lons) + +this_mesh.calculate_corners() +filename = "/glade/scratch/negins/this_region_4x5_new_2/khar2.nc" + +this_mesh.create_esmf(filename, area=None) +print("Done") From dd7f5b817e3b2e17d8625bd9fa8a8f7678257690 Mon Sep 17 00:00:00 2001 From: Negin Sobhani Date: Tue, 1 Nov 2022 17:12:22 -0600 Subject: [PATCH 013/595] mesh_maker --- python/ctsm/site_and_regional/mesh_maker.py | 187 +++++++++++++++----- 1 file changed, 147 insertions(+), 40 deletions(-) mode change 100644 => 100755 python/ctsm/site_and_regional/mesh_maker.py diff --git a/python/ctsm/site_and_regional/mesh_maker.py b/python/ctsm/site_and_regional/mesh_maker.py old mode 100644 new mode 100755 index 4922f17b42..a9449ddb16 --- a/python/ctsm/site_and_regional/mesh_maker.py +++ b/python/ctsm/site_and_regional/mesh_maker.py @@ -1,52 +1,159 @@ +#!/usr/bin/env python3 + from mesh_type import MeshType import xarray as xr import os -oned = True -if oned: - ifile = "/glade/scratch/negins/this_region_4x5_new_2/surfdata_4x5_hist_78pfts_CMIP6_simyr1850_275.0-330.0_-40-15_c220705.nc" -else: - ifile = "/glade/scratch/negins/wrf-ctsm_production/WRF/test/em_real_sim1_noahmp/wrfinput_d01" - -import os, sys, getopt - -if os.path.isfile(ifile): - ds = xr.open_dataset(ifile, mask_and_scale=False, decode_times=False).transpose() -else: - print("Input file could not find!") - sys.exit(2) - - import os, sys, getopt - - if os.path.isfile(ifile): - ds = xr.open_dataset( - ifile, mask_and_scale=False, decode_times=False - ).transpose() +import argparse +import sys +import textwrap + +def get_parser(): + """ + Get the parser object for mesh_maker.py script. + + Returns: + parser (ArgumentParser): + ArgumentParser which includes all the parser information. + """ + parser = argparse.ArgumentParser( + description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter + ) + + parser.print_usage = parser.print_help + + parser.add_argument( + "--input", + help="Netcdf input file for creating ESMF mesh.", + action="store", + dest="input", + required=True, + ) + parser.add_argument( + "--output", + help="Name of the ESMF mesh created.", + action="store", + dest="output", + required=False, + ) + parser.add_argument( + "--lat", + help="Name of latitude varibale on netcdf input file. If none given, looks to find variables that include 'lat'.", + action="store", + dest="lat_name", + type = str, + required=True, + ) + parser.add_argument( + "--lon", + help="Name of latitude varibale on netcdf input file. If none given, looks to find variables that include 'lon'.", + action="store", + dest="lon_name", + type = str, + required=True, + ) + parser.add_argument( + "--mask", + help="Name of mask varibale on netcdf input file. If none given, create a fake mask with values of 1.", + action="store", + dest="mask_name", + type = str, + required=False, + ) + parser.add_argument( + "--area", + help="Name of area variable on netcdf input file. If none given, ESMF calculates element areas automatically. ", + action="store", + dest="area_name", + type = str, + required=False, + ) + parser.add_argument( + "--overwrite", + help="If meshfile exists, overwrite the meshfile.", + action="store_true", + dest="overwrite", + required=False, + ) + return parser + + +def main (): + oned = True + if oned: + ifile = "/glade/scratch/negins/this_region_4x5_new_2/surfdata_4x5_hist_78pfts_CMIP6_simyr1850_275.0-330.0_-40-15_c220705.nc" + else: + ifile = "/glade/scratch/negins/wrf-ctsm_production/WRF/test/em_real_sim1_noahmp/wrfinput_d01" + + if oned: + lat_name = "lsmlat" + lon_name = "lsmlon" + else: + lat_name = "XLAT" + lon_name = "XLONG" + + parser = get_parser() + args = parser.parse_args() + + nc_file = args.input + lat_name = args.lat_name + lon_name = args.lon_name + mesh_out = args.output + overwrite = args.overwrite + mask_name = args.mask_name + area_name = args.area_name + + if os.path.isfile(nc_file): + ds = xr.open_dataset(nc_file, mask_and_scale=False, decode_times=False).transpose() + else: + err_msg = textwrap.dedent( + """\ + \n ------------------------------------ + \n Input file not found. Please make sure to provide the full path of Netcdf input file for making mesh. + \n ------------------------------------ + """ + ) + raise parser.error(err_msg) + + if lat_name not in ds.coords and lat_name not in ds.variables : + print('Input file does not have variable named {}.'.format(lat_name)) + else: + print (lat_name, "exist in the provided netcdf file with dimension of ", len(ds[lat_name].dims)) + + if lon_name not in ds.coords and lat_name not in ds.variables : + print('Input file does not have variable named {}.'.format(lon_name)) else: - print("Input file could not find!") - sys.exit(2) + print (lat_name, "exist in the provided netcdf file with dimension of ", len(ds[lat_name].dims)) -if oned: - lat_name = "lsmlat" - lon_name = "lsmlon" -else: - lat_name = "XLAT" - lon_name = "XLONG" -lats = ds[lat_name] -lons = ds[lon_name] + lats = ds[lat_name] + lons = ds[lon_name] -hasTime = "Time" in ds[lat_name].dims -print(hasTime) + if (len(lats.dims)>3) or (len(lons.dims)>3): + time_dims = [dim for dim in lats.dims if 'time' in dim.lower()] + if time_dims: + print ('time dimension found on lat', time_dims) + lats = lats [:,:,0] + lats = lons [:,:,0] + else: + print ('latitude or longitude has more than 2 dimensions and the third dimension cannot be detected as time.') -if hasTime: - lats = lats[:, :, 0] - lons = lons[:, :, 0] + if mesh_out: + if os.path.exists(mesh_out): + if overwrite: + os.remove(mesh_out) + else: + print ('output meshfile exists, please choose --overwrite to overwrite the mesh file.') + if mask_name is not None: + mask = ds[mask_name].values() -this_mesh = MeshType(lats, lons) + if area_name is not None: + area = ds[mask_name].values() -this_mesh.calculate_corners() -filename = "/glade/scratch/negins/this_region_4x5_new_2/khar2.nc" + this_mesh = MeshType(lats, lons, mask=None) + this_mesh.calculate_corners() + this_mesh.create_esmf(mesh_out, area=None) + print ('DONE!') -this_mesh.create_esmf(filename, area=None) -print("Done") +if __name__ == "__main__": + main() From 1a3115321671048b231d245ee10768a46f8ebdda Mon Sep 17 00:00:00 2001 From: Negin Sobhani Date: Tue, 1 Nov 2022 19:08:31 -0600 Subject: [PATCH 014/595] update subset_data to use the new mesh type --- python/ctsm/site_and_regional/mesh_type.py | 2 +- .../ctsm/site_and_regional/regional_case.py | 32 +++++++++++++++++++ python/ctsm/subset_data.py | 8 +++-- 3 files changed, 39 insertions(+), 3 deletions(-) diff --git a/python/ctsm/site_and_regional/mesh_type.py b/python/ctsm/site_and_regional/mesh_type.py index 640f5548b5..fddb0f226f 100644 --- a/python/ctsm/site_and_regional/mesh_type.py +++ b/python/ctsm/site_and_regional/mesh_type.py @@ -18,7 +18,7 @@ class MeshType: An object for describing mesh or grid files. """ - def __init__(self, center_lats, center_lons, mesh_name=None, mask=None): + def __init__(self, center_lats, center_lons, mask=None, mesh_name=None): """ Construct a mesh object diff --git a/python/ctsm/site_and_regional/regional_case.py b/python/ctsm/site_and_regional/regional_case.py index 16669ffb44..d7d0fa7c91 100644 --- a/python/ctsm/site_and_regional/regional_case.py +++ b/python/ctsm/site_and_regional/regional_case.py @@ -15,6 +15,7 @@ # -- import local classes for this script from ctsm.site_and_regional.base_case import BaseCase, USRDAT_DIR +from ctsm.site_and_regional.mesh_type import MeshType from ctsm.utils import add_tag_to_filename logger = logging.getLogger(__name__) @@ -218,6 +219,13 @@ def create_surfdata_at_reg(self, indir, file, user_mods_dir): logger.info("fsurf_in: %s", fsurf_in) logger.info("fsurf_out: %s", os.path.join(self.out_dir, fsurf_out)) + self.fsurf_out = os.path.join(self.out_dir, fsurf_out) + self.fsurf_out = fsurf_out + + if self.create_mesh: + mesh_out = os.path.join(self.out_dir, os.path.splitext(fsurf_out)[0]+'_ESMF_UNSTRUCTURED_MESH.nc') + self.mesh = mesh_out + # create 1d coordinate variables to enable sel() method f_in = self.create_1d_coord(fsurf_in, "LONGXY", "LATIXY", "lsmlon", "lsmlat") lat = f_in["lat"] @@ -244,6 +252,29 @@ def create_surfdata_at_reg(self, indir, file, user_mods_dir): with open(os.path.join(user_mods_dir, "user_nl_clm"), "a") as nl_clm: line = "fsurdat = '${}'".format(os.path.join(USRDAT_DIR, fsurf_out)) self.write_to_file(line, nl_clm) + if self.create_mesh: + print ('creating mesh_file for this surface dataset.') + self.extract_mesh_at_reg (f_out) + + def extract_mesh_at_reg (self, ds): + """ + Create Mesh from Surface dataset netcdf file. + """ + logger.info("Creating meshfile for at region: %s", self.tag) + + lat_name = "lsmlat" + lon_name = "lsmlon" + + lats = ds[lat_name] + lons = ds[lon_name] + + this_mesh = MeshType(lats, lons) + this_mesh.calculate_corners() + this_mesh.create_esmf(self.mesh) + + print ("DONE") + + def create_landuse_at_reg(self, indir, file, user_mods_dir): """ @@ -493,3 +524,4 @@ def write_shell_commands(self, namelist): ) self.write_to_file("./xmlchange ATM_DOMAIN_MESH={}".format(str(self.mesh)), nl_file) self.write_to_file("./xmlchange LND_DOMAIN_MESH={}".format(str(self.mesh)), nl_file) + self.write_to_file("./xmlchange MASK_MESH={}".format(str('/glade/p/cesmdata/cseg/inputdata/share/meshes/gx1v6_090205_ESMFmesh.nc')), nl_file) diff --git a/python/ctsm/subset_data.py b/python/ctsm/subset_data.py index 7230a72a35..d149267304 100644 --- a/python/ctsm/subset_data.py +++ b/python/ctsm/subset_data.py @@ -520,8 +520,8 @@ 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"]) + #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: @@ -542,6 +542,10 @@ def subset_region(args, file_dict: dict): raise argparse.ArgumentTypeError(err_msg) region.write_shell_commands(os.path.join(args.user_mods_dir, "shell_commands")) + print ('-----') + print (args.user_mods_dir) + print ('-----') + logger.info("Successfully ran script for a regional case.") From d315d2b4667fd0bd507122839df64e82a4c711db Mon Sep 17 00:00:00 2001 From: Negin Sobhani Date: Wed, 2 Nov 2022 20:05:09 -0600 Subject: [PATCH 015/595] clean ups... --- python/ctsm/site_and_regional/mesh_type.py | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/python/ctsm/site_and_regional/mesh_type.py b/python/ctsm/site_and_regional/mesh_type.py index fddb0f226f..d4e9451ae1 100644 --- a/python/ctsm/site_and_regional/mesh_type.py +++ b/python/ctsm/site_and_regional/mesh_type.py @@ -1,7 +1,7 @@ """ -This module includes the definition and functions for defining a Mesh type. +This module includes the definition and functions for defining a Mesh type. """ - +import sys import logging import datetime import numpy as np @@ -238,6 +238,7 @@ def calculate_node_coords(self): [self.corner_lons.T.reshape((-1,)).T, self.corner_lats.T.reshape((-1,)).T], axis=1, ) + corner_pairs = corner_pairs.astype(np.float32, copy=False) # -- remove coordinates that are shared between the elements node_coords = dd.from_dask_array(corner_pairs).drop_duplicates().values @@ -247,15 +248,17 @@ def calculate_node_coords(self): nlon = dims[0] nlat = dims[1] elem_conn_size = nlon * nlat + nlon + nlat + 1 + self.node_coords = node_coords + # -- error check to avoid issues later if self.node_coords.shape[0] != elem_conn_size: logger.warning( "The size of unique coordinate pairs is {} but expected size is {}!".format( self.node_coords.shape[0], elem_conn_size ) ) - sys.exit(2) + logger.warning ("This may result your simulation to crash later.") def calculate_elem_conn(self): """ From 527cdc036bedc9e1f281500f4b6d35dc8ef3f97a Mon Sep 17 00:00:00 2001 From: Negin Sobhani Date: Thu, 3 Nov 2022 00:42:03 -0600 Subject: [PATCH 016/595] update namelist --- tools/site_and_regional/default_data.cfg | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/tools/site_and_regional/default_data.cfg b/tools/site_and_regional/default_data.cfg index 630013d502..c38a00adcc 100644 --- a/tools/site_and_regional/default_data.cfg +++ b/tools/site_and_regional/default_data.cfg @@ -16,10 +16,10 @@ tpqwname = CLMGSWP3v1.TPQW [surfdat] 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 +surfdat_16pft = surfdata_4x5_hist_16pfts_Irrig_CMIP6_simyr2000_c190214.nc +surfdat_78pft = surfdata_4x5_hist_78pfts_CMIP6_simyr2000_c190214.nc mesh_dir = share/meshes/ -mesh_surf = fv0.9x1.25_141008_ESMFmesh.nc +mesh_surf = fv4x5_050615_polemod_ESMFmesh.nc [landuse] dir = lnd/clm2/surfdata_map/release-clm5.0.18 From 947b7681618040ea9c34f12d0ed1a039e423af80 Mon Sep 17 00:00:00 2001 From: Negin Sobhani Date: Thu, 3 Nov 2022 05:46:56 -0600 Subject: [PATCH 017/595] fixing a few issues --- python/ctsm/site_and_regional/mesh_type.py | 87 +++++++++++++++------- 1 file changed, 61 insertions(+), 26 deletions(-) diff --git a/python/ctsm/site_and_regional/mesh_type.py b/python/ctsm/site_and_regional/mesh_type.py index d4e9451ae1..0b04829ec1 100644 --- a/python/ctsm/site_and_regional/mesh_type.py +++ b/python/ctsm/site_and_regional/mesh_type.py @@ -1,9 +1,12 @@ """ -This module includes the definition and functions for defining a Mesh type. +This module includes the definition and functions for defining a Grid or Mesh. +This module enables creating ESMF mesh file (unstructured grid file)for valid 1D or 2D lats and lons. """ import sys import logging +import argparse import datetime + import numpy as np import xarray as xr import pandas as pd @@ -22,14 +25,39 @@ def __init__(self, center_lats, center_lons, mask=None, mesh_name=None): """ Construct a mesh object - center_lats : - latitudes (either 1D or 2D) - center_lons : - longitudes (either 1D or 2D) + Attributes + ---------- + center_lats : xarray DataArray + array of latitudes (either 1D or 2D) + center_lons : xarray DataArray + array longitudes (either 1D or 2D) mesh_name : str or None, optional Name of the mesh mask : np array or None numpy array that include landmask + + Methods + ------- + check_lat_lon_dims: + check if dimensions of latitude and longitude is valid. + + create_artificial_mask: + create an artificial mask of ones if mask does not exits. + + create_2d_coords: + if 1d coords is provided, this method creates 2d coords from 1d coords. + + calculate_corners: + calculate corner coordinates for each polygon in the grid + + calculate_node_coords: + extract node coordiantes for nodeCoords variable on mesh file + + calculate_elem_conn + calculate element connectivity for elementConn variable on mesh file + + create_esmf + write mesh file to netcdf file. """ self.mesh_name = mesh_name @@ -46,38 +74,41 @@ def __init__(self, center_lats, center_lons, mask=None, mesh_name=None): else: self.mask = mask - def __str__(self): - """ - Converts ingredients of this class to string for printing. - """ - return "{}\n{}".format( - str(self.__class__), - "\n".join( - ( - "{} = {}".format(str(key), str(self.__dict__[key])) - for key in sorted(self.__dict__) - ) - ), - ) - def check_lat_lon_dims(self): """ - Check latitude and longitude dimensions to make sure they are either 1 or 2. + Check latitude and longitude dimensions to make sure they are valid. + + ------------- + Raises: + Error (ArgumentTypeError): + If the provided latitude has dimension >2. + Error (ArgumentTypeError): + If the provided longitude has dimension >2. """ + if self.lat_dims not in [1, 2]: - print( - "Unrecognized grid! The dimension of latitude should be either 1 or 2 but it is {}.".format( + err_msg = """ + Unrecognized grid! \n + The dimension of latitude should be either 1 or 2 but it is {}. + """.format( self.lat_dims ) - ) + raise argparse.ArgumentTypeError(err_msg) + if self.lon_dims not in [1, 2]: - print( - "Unrecognized grid! The dimension of longitude should be either 1 or 2 but it is {}.".format( + err_msg = """ + Unrecognized grid! \n + The dimension of longitude should be either 1 or 2 but it is {}. + """.format( self.lon_dims ) - ) + raise argparse.ArgumentTypeError(err_msg) def create_artificial_mask(self): + """ + create an artificial mask of 1 if no land mask is provided. + """ + logger.info("Creating an artificial mask for this region...") if self.lat_dims == 1: @@ -238,11 +269,14 @@ def calculate_node_coords(self): [self.corner_lons.T.reshape((-1,)).T, self.corner_lats.T.reshape((-1,)).T], axis=1, ) + + # -- convert to float32 to find duplicates corner_pairs = corner_pairs.astype(np.float32, copy=False) # -- remove coordinates that are shared between the elements node_coords = dd.from_dask_array(corner_pairs).drop_duplicates().values node_coords.compute_chunk_sizes() + # -- check size of unique coordinate pairs dims = self.mask.shape nlon = dims[0] @@ -279,6 +313,7 @@ def calculate_elem_conn(self): elem_conn = corners.compute().groupby(["lon", "lat"], sort=False).ngroup() + 1 elem_conn = da.from_array(elem_conn.to_numpy()) + # -- reshape to write to ESMF self.elem_conn = elem_conn.T.reshape((4, -1)).T From c0c6aa9e216f686eff088b3c589b956b7564d809 Mon Sep 17 00:00:00 2001 From: Negin Sobhani Date: Thu, 3 Nov 2022 06:07:07 -0600 Subject: [PATCH 018/595] updates --- .../ctsm/site_and_regional/regional_case.py | 82 +++++++++---------- 1 file changed, 37 insertions(+), 45 deletions(-) diff --git a/python/ctsm/site_and_regional/regional_case.py b/python/ctsm/site_and_regional/regional_case.py index d7d0fa7c91..d77a14a7d0 100644 --- a/python/ctsm/site_and_regional/regional_case.py +++ b/python/ctsm/site_and_regional/regional_case.py @@ -10,7 +10,6 @@ # -- 3rd party libraries import numpy as np import xarray as xr -from tqdm import tqdm from datetime import datetime # -- import local classes for this script @@ -54,6 +53,11 @@ class RegionalCase(BaseCase): Methods ------- + create_tag + Create a tag for this region which is either + region's name or a combination of bounds of this + region lat1-lat2_lon1-lon2 + check_region_bounds Check for the regional bounds @@ -63,20 +67,20 @@ class RegionalCase(BaseCase): 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 - region lat1-lat2_lon1-lon2 - create_domain_at_reg Create domain file at this region create_surfdata_at_reg Create surface dataset at this region + extract_mesh_at_reg + Extract mesh from surface dataset created by create_surfdata_at_reg + create_landuse_at_reg Create landuse file at this region + write_shell_commands(namelist) + write out xml commands to a file for usermods (i.e. shell_commands) for regional settings. """ def __init__( @@ -253,7 +257,7 @@ def create_surfdata_at_reg(self, indir, file, user_mods_dir): line = "fsurdat = '${}'".format(os.path.join(USRDAT_DIR, fsurf_out)) self.write_to_file(line, nl_clm) if self.create_mesh: - print ('creating mesh_file for this surface dataset.') + logger.info('creating mesh file from surface_dataset={}'.format(str(str(self.mesh))) self.extract_mesh_at_reg (f_out) def extract_mesh_at_reg (self, ds): @@ -272,9 +276,6 @@ def extract_mesh_at_reg (self, ds): this_mesh.calculate_corners() this_mesh.create_esmf(self.mesh) - print ("DONE") - - def create_landuse_at_reg(self, indir, file, user_mods_dir): """ @@ -364,66 +365,58 @@ def subset_mesh_at_reg (self, mesh_in): subset_element = [] cnt = 0 - for n in tqdm(range(elem_count)): + for n in range(elem_count): endx = elem_conn[n,:num_elem_conn[n].values].values endx[:,] -= 1# convert to zero based index endx = [int(xi) for xi in endx] - #print (endx) + nlon = node_coords[endx,0].values nlat = node_coords[endx,1].values - - l1 = np.logical_and(nlon > self.lon1,nlon < self.lon2) - l2 = np.logical_and(nlat > self.lat1,nlat < self.lat2) - if np.any(l1) and np.any(l2): - #if np.any(np.logical_or(l1,l2)): - # pass - #else: + + 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)): + pass + else: subset_element.append(n) cnt+=1 - print (cnt) subset_node = [] conn_dict = {} - cnt = 1 - print (node_count) - print (elem_count) - print (node_coords) + cnt = 1 for n in range(node_count): - #n = 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) - l1 = np.logical_and(nlon >= self.lon1,nlon <= self.lon2) - l2 = np.logical_and(nlat >= self.lat1,nlat <= self.lat2) - if np.any(l1) and np.any(l2): + + 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 + conn_dict[n+1] = cnt cnt+=1 - else: - conn_dict[n+1] = -9999 - - #if np.logical_or(l1,l2): - # conn_dict[n+1] = -9999 + + # -- reverse logic + #l1 = np.logical_and(nlon >= self.lon1,nlon <= self.lon2) + #l2 = np.logical_and(nlat >= self.lat1,nlat <= self.lat2) #if np.any(l1) and np.any(l2): - #else: # subset_node.append(n) # conn_dict[n+1] = cnt # cnt+=1 - print ('cnt:',cnt) - return node_coords, subset_element, subset_node, conn_dict + #else: + # conn_dict[n+1] = -9999 + return node_coords, subset_element, subset_node, conn_dict def write_mesh (self, f_in, node_coords, subset_element, subset_node, conn_dict, mesh_out): """ This function writes out the subsetted mesh file. """ - print (len(subset_node)) corner_pairs = f_in.variables['nodeCoords'][subset_node,] - print (corner_pairs) - #print (corner_pairs.shape()) dimensions = f_in.dims variables = f_in.variables global_attributes = f_in.attrs @@ -438,7 +431,6 @@ def write_mesh (self, f_in, node_coords, subset_element, subset_node, conn_dict, for ni in range(elem_count): for mi in range(max_node_dim): ndx = int (elem_conn_index[ni,mi]) - #print ('mi',mi, 'ndx',ndx) elem_conn_out[ni,mi] = conn_dict[ndx] @@ -524,4 +516,4 @@ def write_shell_commands(self, namelist): ) self.write_to_file("./xmlchange ATM_DOMAIN_MESH={}".format(str(self.mesh)), nl_file) self.write_to_file("./xmlchange LND_DOMAIN_MESH={}".format(str(self.mesh)), nl_file) - self.write_to_file("./xmlchange MASK_MESH={}".format(str('/glade/p/cesmdata/cseg/inputdata/share/meshes/gx1v6_090205_ESMFmesh.nc')), nl_file) + self.write_to_file("./xmlchange MASK_MESH={}".format(str(str(self.mesh))), nl_file) From 3f47819b4e78c09b6156535074da095b27668cc5 Mon Sep 17 00:00:00 2001 From: Negin Sobhani Date: Thu, 3 Nov 2022 06:08:46 -0600 Subject: [PATCH 019/595] black-ify the code --- .../ctsm/site_and_regional/regional_case.py | 276 ++++++++++-------- 1 file changed, 159 insertions(+), 117 deletions(-) diff --git a/python/ctsm/site_and_regional/regional_case.py b/python/ctsm/site_and_regional/regional_case.py index d77a14a7d0..61d0fdf78e 100644 --- a/python/ctsm/site_and_regional/regional_case.py +++ b/python/ctsm/site_and_regional/regional_case.py @@ -133,14 +133,14 @@ def create_tag(self): str(self.lon1), str(self.lon2), str(self.lat1), str(self.lat2) ) - def check_region_bounds (self): + def check_region_bounds(self): """ Check for the regional bounds """ self.check_region_lons() self.check_region_lats() - def check_region_lons (self): + def check_region_lons(self): """ Check for the regional lon bounds """ @@ -159,7 +159,7 @@ def check_region_lons (self): ) raise argparse.ArgumentTypeError(err_msg) - def check_region_lats (self): + def check_region_lats(self): """ Check for the regional lat bound """ @@ -176,7 +176,6 @@ def check_region_lats (self): ) raise argparse.ArgumentTypeError(err_msg) - def create_domain_at_reg(self, indir, file): """ Create domain file for this RegionalCase class. @@ -227,7 +226,10 @@ def create_surfdata_at_reg(self, indir, file, user_mods_dir): self.fsurf_out = fsurf_out if self.create_mesh: - mesh_out = os.path.join(self.out_dir, os.path.splitext(fsurf_out)[0]+'_ESMF_UNSTRUCTURED_MESH.nc') + mesh_out = os.path.join( + self.out_dir, + os.path.splitext(fsurf_out)[0] + "_ESMF_UNSTRUCTURED_MESH.nc", + ) self.mesh = mesh_out # create 1d coordinate variables to enable sel() method @@ -257,10 +259,10 @@ def create_surfdata_at_reg(self, indir, file, user_mods_dir): line = "fsurdat = '${}'".format(os.path.join(USRDAT_DIR, fsurf_out)) self.write_to_file(line, nl_clm) if self.create_mesh: - logger.info('creating mesh file from surface_dataset={}'.format(str(str(self.mesh))) - self.extract_mesh_at_reg (f_out) + logger.info("creating mesh file from surface_dataset: %s", wfile) + self.extract_mesh_at_reg(f_out) - def extract_mesh_at_reg (self, ds): + def extract_mesh_at_reg(self, ds): """ Create Mesh from Surface dataset netcdf file. """ @@ -276,7 +278,6 @@ def extract_mesh_at_reg (self, ds): this_mesh.calculate_corners() this_mesh.create_esmf(self.mesh) - def create_landuse_at_reg(self, indir, file, user_mods_dir): """ Create land use data file for this RegionalCase class. @@ -319,7 +320,6 @@ def create_landuse_at_reg(self, indir, file, user_mods_dir): ) self.write_to_file(line, nl_clm) - def create_mesh_at_reg(self, mesh_dir, mesh_surf): """ Create a mesh subsetted for the RegionalCase class. @@ -327,177 +327,213 @@ def create_mesh_at_reg(self, mesh_dir, mesh_surf): logger.info( "----------------------------------------------------------------------" ) - logger.info( - "Subsetting mesh file for region: %s", - self.tag - ) + 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+'_c'+today_string+'.nc') + 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) self.mesh = 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) + 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 + ) - def subset_mesh_at_reg (self, 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'] + 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 range(elem_count): - endx = elem_conn[n,:num_elem_conn[n].values].values - endx[:,] -= 1# convert to zero based index + endx = elem_conn[n, : num_elem_conn[n].values].values + endx[ + :, + ] -= 1 # convert to zero based index endx = [int(xi) for xi in endx] - - nlon = node_coords[endx,0].values - nlat = node_coords[endx,1].values - - 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)): + + nlon = node_coords[endx, 0].values + nlat = node_coords[endx, 1].values + + 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)): pass else: subset_element.append(n) - cnt+=1 + cnt += 1 subset_node = [] conn_dict = {} - cnt = 1 + cnt = 1 for n in range(node_count): - 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) - - if np.logical_or(l1,l2): - conn_dict[n+1] = -9999 + 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) + + if np.logical_or(l1, l2): + conn_dict[n + 1] = -9999 else: subset_node.append(n) - conn_dict[n+1] = cnt - cnt+=1 - + conn_dict[n + 1] = cnt + cnt += 1 + # -- reverse logic - #l1 = np.logical_and(nlon >= self.lon1,nlon <= self.lon2) - #l2 = np.logical_and(nlat >= self.lat1,nlat <= self.lat2) - #if np.any(l1) and np.any(l2): + # l1 = np.logical_and(nlon >= self.lon1,nlon <= self.lon2) + # l2 = np.logical_and(nlat >= self.lat1,nlat <= self.lat2) + # if np.any(l1) and np.any(l2): # subset_node.append(n) # conn_dict[n+1] = cnt # cnt+=1 - #else: + # else: # conn_dict[n+1] = -9999 return node_coords, subset_element, subset_node, conn_dict - - def write_mesh (self, f_in, node_coords, subset_element, subset_node, conn_dict, mesh_out): + 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,] + corner_pairs = f_in.variables["nodeCoords"][ + subset_node, + ] dimensions = f_in.dims - variables = f_in.variables - global_attributes = f_in.attrs - + variables = f_in.variables + global_attributes = f_in.attrs - max_node_dim = len(f_in['maxNodePElement']) + 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,] + 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] + 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,] + 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["nodeCoords"] = xr.DataArray( + corner_pairs, dims=("nodeCount", "coordDim"), attrs={"units": "degrees"} + ) - f_out['numElementConn'] = xr.DataArray(num_elem_conn_out, - dims=('elementCount'), - attrs={'long_name': 'Number of nodes per element'}) - f_out.numElementConn.encoding = {'dtype': np.int32} + 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['centerCoords'] = xr.DataArray(center_coords_out, - dims=('elementCount', 'coordDim'), - attrs={'units': 'degrees'}) + f_out["numElementConn"] = xr.DataArray( + num_elem_conn_out, + dims=("elementCount"), + attrs={"long_name": "Number of nodes per element"}, + ) + f_out.numElementConn.encoding = {"dtype": np.int32} + f_out["centerCoords"] = xr.DataArray( + center_coords_out, + dims=("elementCount", "coordDim"), + attrs={"units": "degrees"}, + ) - #-- add mask - if 'elementMask' in variables: - f_out['elementMask'] = xr.DataArray(elem_mask_out, - dims=('elementCount'), - attrs={'units': 'unitless'}) - f_out.elementMask.encoding = {'dtype': np.int32} + # -- add mask + if "elementMask" in variables: + f_out["elementMask"] = xr.DataArray( + elem_mask_out, dims=("elementCount"), attrs={"units": "unitless"} + ) + f_out.elementMask.encoding = {"dtype": np.int32} - if 'elementArea' in variables: - f_out['elementArea'] = xr.DataArray(elem_area_out, - dims=('elementCount'), - attrs={'units': 'unitless'}) + if "elementArea" in variables: + f_out["elementArea"] = xr.DataArray( + elem_area_out, dims=("elementCount"), attrs={"units": "unitless"} + ) - #-- setting fill values + # -- setting fill values for var in variables: - if '_FillValue' in f_in[var].encoding: - f_out[var].encoding['_FillValue'] = f_in[var].encoding['_FillValue'] + if "_FillValue" in f_in[var].encoding: + f_out[var].encoding["_FillValue"] = f_in[var].encoding["_FillValue"] else: - f_out[var].encoding['_FillValue'] = None + f_out[var].encoding["_FillValue"] = None - #-- add global attributes + # -- add global attributes for attr in global_attributes: - if attr != 'timeGenerated': + if attr != "timeGenerated": f_out.attrs[attr] = global_attributes[attr] - f_out.attrs = {'title': 'ESMF unstructured grid file for a region', - 'created_by': 'subset_data', - 'date_created': '{}'.format(datetime.now()), - } + f_out.attrs = { + "title": "ESMF unstructured grid file for a region", + "created_by": "subset_data", + "date_created": "{}".format(datetime.now()), + } f_out.to_netcdf(mesh_out) logger.info("Successfully created file (mesh_out) %s", mesh_out) @@ -514,6 +550,12 @@ def write_shell_commands(self, namelist): 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 LND_DOMAIN_MESH={}".format(str(self.mesh)), nl_file) - self.write_to_file("./xmlchange MASK_MESH={}".format(str(str(self.mesh))), nl_file) + self.write_to_file( + "./xmlchange ATM_DOMAIN_MESH={}".format(str(self.mesh)), nl_file + ) + self.write_to_file( + "./xmlchange LND_DOMAIN_MESH={}".format(str(self.mesh)), nl_file + ) + self.write_to_file( + "./xmlchange MASK_MESH={}".format(str(str(self.mesh))), nl_file + ) From ca2d2da580a25b0a3c9b0e93ff35068187691820 Mon Sep 17 00:00:00 2001 From: Negin Sobhani Date: Thu, 3 Nov 2022 08:01:44 -0600 Subject: [PATCH 020/595] final edits --- python/ctsm/site_and_regional/mesh_maker.py | 113 ++++++++++++++------ 1 file changed, 81 insertions(+), 32 deletions(-) diff --git a/python/ctsm/site_and_regional/mesh_maker.py b/python/ctsm/site_and_regional/mesh_maker.py index a9449ddb16..8a5b61d5a7 100755 --- a/python/ctsm/site_and_regional/mesh_maker.py +++ b/python/ctsm/site_and_regional/mesh_maker.py @@ -1,11 +1,13 @@ #!/usr/bin/env python3 -from mesh_type import MeshType -import xarray as xr import os -import argparse import sys +import logging +import argparse import textwrap +from datetime import datetime +import xarray as xr +from mesh_type import MeshType def get_parser(): """ @@ -35,6 +37,16 @@ def get_parser(): dest="output", required=False, ) + + parser.add_argument( + "--outdir", + help="Output directory (only if name of output mesh is not defined)", + action="store", + dest="out_dir", + type=str, + ) + + parser.add_argument( "--lat", help="Name of latitude varibale on netcdf input file. If none given, looks to find variables that include 'lat'.", @@ -74,30 +86,32 @@ def get_parser(): dest="overwrite", required=False, ) - return parser + parser.add_argument( + "-v", "--verbose", + help="Increase output verbosity", + action="store_true", + dest = "verbose", + required = False, + ) + return parser def main (): - oned = True - if oned: - ifile = "/glade/scratch/negins/this_region_4x5_new_2/surfdata_4x5_hist_78pfts_CMIP6_simyr1850_275.0-330.0_-40-15_c220705.nc" - else: - ifile = "/glade/scratch/negins/wrf-ctsm_production/WRF/test/em_real_sim1_noahmp/wrfinput_d01" - - if oned: - lat_name = "lsmlat" - lon_name = "lsmlon" - else: - lat_name = "XLAT" - lon_name = "XLONG" parser = get_parser() args = parser.parse_args() + if args.verbose: + logging.basicConfig(level= logging.DEBUG, format=' %(levelname)-8s :: %(message)s') + else: + logging.basicConfig(level= logging.INFO, format=' %(levelname)-8s :: %(message)s') + + nc_file = args.input lat_name = args.lat_name lon_name = args.lon_name mesh_out = args.output + out_dir = args.out_dir overwrite = args.overwrite mask_name = args.mask_name area_name = args.area_name @@ -115,34 +129,70 @@ def main (): raise parser.error(err_msg) if lat_name not in ds.coords and lat_name not in ds.variables : - print('Input file does not have variable named {}.'.format(lat_name)) + logging.error('Input file does not have variable named %s',lat_name) + sys.exit() + else: - print (lat_name, "exist in the provided netcdf file with dimension of ", len(ds[lat_name].dims)) + logging.debug ("- %s exist in the provided netcdf file with dimension of %s.", lat_name, len(ds[lat_name].dims).__str__()) - if lon_name not in ds.coords and lat_name not in ds.variables : - print('Input file does not have variable named {}.'.format(lon_name)) + if lon_name not in ds.coords and lon_name not in ds.variables : + logging.error('Input file does not have variable named %s',lon_name) + sys.exit() else: - print (lat_name, "exist in the provided netcdf file with dimension of ", len(ds[lat_name].dims)) + logging.debug ("- %s exist in the provided netcdf file with dimension of %s.", lon_name, len(ds[lon_name].dims).__str__()) lats = ds[lat_name] lons = ds[lon_name] - if (len(lats.dims)>3) or (len(lons.dims)>3): + if (len(lats.dims)>2) or (len(lons.dims)>2): time_dims = [dim for dim in lats.dims if 'time' in dim.lower()] if time_dims: - print ('time dimension found on lat', time_dims) + logging.debug ('- time dimension found on lat {}'.format(time_dims)) + logging.debug ('- removing time dimensions from lats and lons. ') lats = lats [:,:,0] - lats = lons [:,:,0] + lons = lons [:,:,0] else: - print ('latitude or longitude has more than 2 dimensions and the third dimension cannot be detected as time.') + logging.error ('latitude or longitude has more than 2 dimensions and the third dimension cannot be detected as time.') + sys.exit() + + today = datetime.today() + today_string = today.strftime("%y%m%d") - if mesh_out: - if os.path.exists(mesh_out): - if overwrite: - os.remove(mesh_out) - else: - print ('output meshfile exists, please choose --overwrite to overwrite the mesh file.') + if mesh_out and out_dir: + logging.error(" Both --outdir and --output cannot be provided at the same time.") + err_msg = textwrap.dedent(''' + \n ------------------------------------ + \n You have provided both --outdir and --output. + \n Please provide only one of these options to proceed: + \n --outdir : directory to save mesh file. mesh file name automatically created. + \n --output : Absolute or relative path of the ESMF mesh file created.\n + ''' + ) + raise parser.error(err_msg) + + #-- no file name and output path: + if not mesh_out and not out_dir: + out_dir = os.path.join(os.getcwd(),'meshes') + + if not mesh_out: + #-- make output path if does not exist. + if not os.path.isdir(out_dir): + os.mkdir(out_dir) + + mesh_out = os.path.join(out_dir, + os.path.splitext(nc_file)[0] + "_ESMF_UNSTRUCTURED_MESH"+ "_c"+ today_string+".nc") + + logging.info ('Creating mesh file from : %s', nc_file) + logging.info ('Writing mesh file to : %s', mesh_out) + + # -- exit if mesh_out exists and --overwrite is not specified. + if os.path.exists(mesh_out): + if overwrite: + os.remove(mesh_out) + else: + logging.error ('output meshfile exists, please choose --overwrite to overwrite the mesh file.') + sys.exit() if mask_name is not None: mask = ds[mask_name].values() @@ -153,7 +203,6 @@ def main (): this_mesh = MeshType(lats, lons, mask=None) this_mesh.calculate_corners() this_mesh.create_esmf(mesh_out, area=None) - print ('DONE!') if __name__ == "__main__": main() From 9ec631b30121092f1d2834f4dd433e316e21995b Mon Sep 17 00:00:00 2001 From: Negin Sobhani Date: Thu, 3 Nov 2022 08:04:23 -0600 Subject: [PATCH 021/595] updates --- .../site_and_regional/mesh_maker.py | 119 +++++++++++------- 1 file changed, 77 insertions(+), 42 deletions(-) rename {python/ctsm => tools}/site_and_regional/mesh_maker.py (63%) diff --git a/python/ctsm/site_and_regional/mesh_maker.py b/tools/site_and_regional/mesh_maker.py similarity index 63% rename from python/ctsm/site_and_regional/mesh_maker.py rename to tools/site_and_regional/mesh_maker.py index 8a5b61d5a7..a130042167 100755 --- a/python/ctsm/site_and_regional/mesh_maker.py +++ b/tools/site_and_regional/mesh_maker.py @@ -7,7 +7,15 @@ import textwrap from datetime import datetime import xarray as xr -from mesh_type import MeshType + +# -- add python/ctsm to path +_CTSM_PYTHON = os.path.join( + os.path.dirname(os.path.realpath(__file__)), os.pardir, os.pardir, "python" +) +sys.path.insert(1, _CTSM_PYTHON) + +from ctsm.site_and_regional.mesh_type import MeshType + def get_parser(): """ @@ -44,15 +52,14 @@ def get_parser(): action="store", dest="out_dir", type=str, - ) - + ) parser.add_argument( "--lat", help="Name of latitude varibale on netcdf input file. If none given, looks to find variables that include 'lat'.", action="store", dest="lat_name", - type = str, + type=str, required=True, ) parser.add_argument( @@ -60,7 +67,7 @@ def get_parser(): help="Name of latitude varibale on netcdf input file. If none given, looks to find variables that include 'lon'.", action="store", dest="lon_name", - type = str, + type=str, required=True, ) parser.add_argument( @@ -68,7 +75,7 @@ def get_parser(): help="Name of mask varibale on netcdf input file. If none given, create a fake mask with values of 1.", action="store", dest="mask_name", - type = str, + type=str, required=False, ) parser.add_argument( @@ -76,7 +83,7 @@ def get_parser(): help="Name of area variable on netcdf input file. If none given, ESMF calculates element areas automatically. ", action="store", dest="area_name", - type = str, + type=str, required=False, ) parser.add_argument( @@ -87,25 +94,30 @@ def get_parser(): required=False, ) parser.add_argument( - "-v", "--verbose", + "-v", + "--verbose", help="Increase output verbosity", - action="store_true", - dest = "verbose", - required = False, + action="store_true", + dest="verbose", + required=False, ) return parser -def main (): + +def main(): parser = get_parser() args = parser.parse_args() if args.verbose: - logging.basicConfig(level= logging.DEBUG, format=' %(levelname)-8s :: %(message)s') + logging.basicConfig( + level=logging.DEBUG, format=" %(levelname)-8s :: %(message)s" + ) else: - logging.basicConfig(level= logging.INFO, format=' %(levelname)-8s :: %(message)s') - + logging.basicConfig( + level=logging.INFO, format=" %(levelname)-8s :: %(message)s" + ) nc_file = args.input lat_name = args.lat_name @@ -117,7 +129,9 @@ def main (): area_name = args.area_name if os.path.isfile(nc_file): - ds = xr.open_dataset(nc_file, mask_and_scale=False, decode_times=False).transpose() + ds = xr.open_dataset( + nc_file, mask_and_scale=False, decode_times=False + ).transpose() else: err_msg = textwrap.dedent( """\ @@ -128,70 +142,90 @@ def main (): ) raise parser.error(err_msg) - if lat_name not in ds.coords and lat_name not in ds.variables : - logging.error('Input file does not have variable named %s',lat_name) + if lat_name not in ds.coords and lat_name not in ds.variables: + logging.error("Input file does not have variable named %s", lat_name) sys.exit() else: - logging.debug ("- %s exist in the provided netcdf file with dimension of %s.", lat_name, len(ds[lat_name].dims).__str__()) + logging.debug( + "- %s exist in the provided netcdf file with dimension of %s.", + lat_name, + len(ds[lat_name].dims).__str__(), + ) - if lon_name not in ds.coords and lon_name not in ds.variables : - logging.error('Input file does not have variable named %s',lon_name) + if lon_name not in ds.coords and lon_name not in ds.variables: + logging.error("Input file does not have variable named %s", lon_name) sys.exit() else: - logging.debug ("- %s exist in the provided netcdf file with dimension of %s.", lon_name, len(ds[lon_name].dims).__str__()) - + logging.debug( + "- %s exist in the provided netcdf file with dimension of %s.", + lon_name, + len(ds[lon_name].dims).__str__(), + ) lats = ds[lat_name] lons = ds[lon_name] - if (len(lats.dims)>2) or (len(lons.dims)>2): - time_dims = [dim for dim in lats.dims if 'time' in dim.lower()] + if (len(lats.dims) > 2) or (len(lons.dims) > 2): + time_dims = [dim for dim in lats.dims if "time" in dim.lower()] if time_dims: - logging.debug ('- time dimension found on lat {}'.format(time_dims)) - logging.debug ('- removing time dimensions from lats and lons. ') - lats = lats [:,:,0] - lons = lons [:,:,0] + logging.debug("- time dimension found on lat {}".format(time_dims)) + logging.debug("- removing time dimensions from lats and lons. ") + lats = lats[:, :, 0] + lons = lons[:, :, 0] else: - logging.error ('latitude or longitude has more than 2 dimensions and the third dimension cannot be detected as time.') + logging.error( + "latitude or longitude has more than 2 dimensions and the third dimension cannot be detected as time." + ) sys.exit() today = datetime.today() today_string = today.strftime("%y%m%d") if mesh_out and out_dir: - logging.error(" Both --outdir and --output cannot be provided at the same time.") - err_msg = textwrap.dedent(''' + logging.error( + " Both --outdir and --output cannot be provided at the same time." + ) + err_msg = textwrap.dedent( + """ \n ------------------------------------ \n You have provided both --outdir and --output. \n Please provide only one of these options to proceed: \n --outdir : directory to save mesh file. mesh file name automatically created. \n --output : Absolute or relative path of the ESMF mesh file created.\n - ''' - ) + """ + ) raise parser.error(err_msg) - #-- no file name and output path: + # -- no file name and output path: if not mesh_out and not out_dir: - out_dir = os.path.join(os.getcwd(),'meshes') + out_dir = os.path.join(os.getcwd(), "meshes") if not mesh_out: - #-- make output path if does not exist. + # -- make output path if does not exist. if not os.path.isdir(out_dir): os.mkdir(out_dir) - mesh_out = os.path.join(out_dir, - os.path.splitext(nc_file)[0] + "_ESMF_UNSTRUCTURED_MESH"+ "_c"+ today_string+".nc") + mesh_out = os.path.join( + out_dir, + os.path.splitext(nc_file)[0] + + "_ESMF_UNSTRUCTURED_MESH" + + "_c" + + today_string + + ".nc", + ) - logging.info ('Creating mesh file from : %s', nc_file) - logging.info ('Writing mesh file to : %s', mesh_out) + logging.info("Creating mesh file from : %s", nc_file) + logging.info("Writing mesh file to : %s", mesh_out) # -- exit if mesh_out exists and --overwrite is not specified. if os.path.exists(mesh_out): if overwrite: os.remove(mesh_out) else: - logging.error ('output meshfile exists, please choose --overwrite to overwrite the mesh file.') + logging.error( + "output meshfile exists, please choose --overwrite to overwrite the mesh file." + ) sys.exit() if mask_name is not None: @@ -204,5 +238,6 @@ def main (): this_mesh.calculate_corners() this_mesh.create_esmf(mesh_out, area=None) + if __name__ == "__main__": main() From d9acef07c4c1b00987d421c8c0aedb0bc10aa3b3 Mon Sep 17 00:00:00 2001 From: Negin Sobhani Date: Thu, 3 Nov 2022 08:05:10 -0600 Subject: [PATCH 022/595] updates --- python/ctsm/site_and_regional/mesh_type.py | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/python/ctsm/site_and_regional/mesh_type.py b/python/ctsm/site_and_regional/mesh_type.py index 0b04829ec1..12bbac3e3d 100644 --- a/python/ctsm/site_and_regional/mesh_type.py +++ b/python/ctsm/site_and_regional/mesh_type.py @@ -119,8 +119,8 @@ def create_artificial_mask(self): elif self.lat_dims == 2: # -- 2D mask mask = np.ones( - center_lats.shape, dtype=np.int8 - ) # np.ones(tuple(self.center_lats.sizes.values()), dtype=np.int8) # + self.center_lats.shape, dtype=np.int8 + ) mask_da = da.from_array(mask) self.mask = mask_da @@ -393,7 +393,10 @@ def create_esmf(self, mesh_fname, area=None): ds_out[var].encoding["_FillValue"] = None # -- add global attributes - ds_out.attrs["title"] = "ESMF unstructured grid file" + if self.mesh_name: + ds_out.attrs["title"] = "ESMF unstructured grid file " + self.mesh_name + else: + ds_out.attrs["title"] = "ESMF unstructured grid file " ds_out.attrs["gridType"] = "unstructured mesh" ds_out.attrs["version"] = "0.9" ds_out.attrs["conventions"] = "ESMFMESH" From 07270a934fd5500caee37aaef6a8ab309f21c823 Mon Sep 17 00:00:00 2001 From: Negin Sobhani Date: Thu, 3 Nov 2022 08:12:25 -0600 Subject: [PATCH 023/595] updates to subset_data --- python/ctsm/subset_data.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/python/ctsm/subset_data.py b/python/ctsm/subset_data.py index d149267304..48293b5248 100644 --- a/python/ctsm/subset_data.py +++ b/python/ctsm/subset_data.py @@ -542,10 +542,10 @@ def subset_region(args, file_dict: dict): raise argparse.ArgumentTypeError(err_msg) region.write_shell_commands(os.path.join(args.user_mods_dir, "shell_commands")) - print ('-----') - print (args.user_mods_dir) - print ('-----') + print ('\nFor running this regional case with the created user_mods : ') + print ('./create_newcase --case case --res CLM_USRDAT --compset I2000Clm51BgcCrop', + '--run-unsupported --user-mods-dirs ', args.user_mods_dir, '\n\n') logger.info("Successfully ran script for a regional case.") From 734854a43bb6f3003dd0af225e51069ca5b5e9f3 Mon Sep 17 00:00:00 2001 From: Negin Sobhani Date: Thu, 3 Nov 2022 08:18:15 -0600 Subject: [PATCH 024/595] updating readme. --- tools/site_and_regional/mesh_maker.py | 14 +++++++++++++- 1 file changed, 13 insertions(+), 1 deletion(-) diff --git a/tools/site_and_regional/mesh_maker.py b/tools/site_and_regional/mesh_maker.py index a130042167..3229bea2e9 100755 --- a/tools/site_and_regional/mesh_maker.py +++ b/tools/site_and_regional/mesh_maker.py @@ -1,5 +1,17 @@ #!/usr/bin/env python3 - +""" +|------------------------------------------------------------------| +|--------------------- Instructions -----------------------------| +|------------------------------------------------------------------| +This script creates ESMF unstructured GRID (mesh file) from a netcdf +file with valid lats and lons. Provided lats and lons can be 1D or 2D. + +For example for running WRF-CTSM cases, the user can create a mesh +file for their domain : + ./mesh_maker.py --input wrfinput_d01 --output my_region + --lat XLAT --lon XLONG --verbose + +""" import os import sys import logging From b35d76bbeb21bbc085c2b8dd7b984e2b0e9065d0 Mon Sep 17 00:00:00 2001 From: Negin Sobhani Date: Thu, 3 Nov 2022 08:20:50 -0600 Subject: [PATCH 025/595] update README --- tools/site_and_regional/README | 3 +++ 1 file changed, 3 insertions(+) diff --git a/tools/site_and_regional/README b/tools/site_and_regional/README index 930a07a7ac..ecb7fb8381 100644 --- a/tools/site_and_regional/README +++ b/tools/site_and_regional/README @@ -19,6 +19,9 @@ subset_data For extracting domain files, surface dataset, and DATM files at a region, use: ./subset_data region +mesh_maker.py + This script creates a mesh file from a netcdf file with valid lats and lons. + modify_singlept_site_neon.py After running subset_data.py overwrite some fields with site-specific data for neon sites. From 56cdbae4032aad43a999f3de2d4e6c8eff5d46fd Mon Sep 17 00:00:00 2001 From: Negin Sobhani Date: Thu, 3 Nov 2022 08:41:33 -0600 Subject: [PATCH 026/595] updates to mesh file type to plot --- python/ctsm/site_and_regional/mesh_type.py | 59 +++++++++++++++++++++- 1 file changed, 58 insertions(+), 1 deletion(-) diff --git a/python/ctsm/site_and_regional/mesh_type.py b/python/ctsm/site_and_regional/mesh_type.py index 12bbac3e3d..bcd1199f77 100644 --- a/python/ctsm/site_and_regional/mesh_type.py +++ b/python/ctsm/site_and_regional/mesh_type.py @@ -1,6 +1,6 @@ """ This module includes the definition and functions for defining a Grid or Mesh. -This module enables creating ESMF mesh file (unstructured grid file)for valid 1D or 2D lats and lons. +This enables creating ESMF mesh file (unstructured grid file)for valid 1D or 2D lats and lons. """ import sys import logging @@ -14,6 +14,12 @@ import dask.dataframe as dd from dask.diagnostics import ProgressBar +# -- libraries for plotting mesh (make_mesh_plot) +import matplotlib.pyplot as plt +import cartopy.crs as ccrs +import matplotlib.pyplot as plt +import cartopy.feature as cfeature + logger = logging.getLogger(__name__) class MeshType: @@ -409,3 +415,54 @@ def create_esmf(self, mesh_fname, area=None): logger.info("Writing ESMF Mesh file to : %s", mesh_fname) ds_out.to_netcdf(mesh_fname) logger.info("Successfully created ESMF Mesh file : %s", mesh_fname) + + plot_fname = mesh_fname +'.png' + self.make_plot = False + + if self.make_plot: + self.make_mesh_plot(plot_fname) + + def make_mesh_plot(self, plot_fname): + """ + Create an ESMF mesh file for the mesh + + Parameters + ---------- + plot_fname : str + The path to write the ESMF meshfile plot + """ + plt.figure(num=None, figsize=(25, 11), facecolor='w', edgecolor='k') + ax = plt.axes(projection=ccrs.PlateCarree()) + + + ax.coastlines() + ax.gridlines(color="black", linestyle="dotted") + ax.add_feature(cfeature.OCEAN) + ax.add_feature(cfeature.BORDERS) + ax.add_feature(cfeature.LAND, edgecolor='black') + ax.add_feature(cfeature.LAKES, edgecolor='black') + ax.add_feature(cfeature.RIVERS) + + x, y = self.node_coords.T + print (x.shape) + print (y.shape) + df = pd.DataFrame({'x':x, 'y':y}, index = [0]) + grouped_x = df.groupby('x') + for name, group in grouped_x: + plt.plot (group['x']-180, group['y'],color='black', marker='o') + grouped_y = df.groupby('y') + for name, group in grouped_y: + plt.plot (group['x']-180, group['y'],color='black', marker='o') + + x, y = this_mesh.center_coords.T + + plt.scatter (x-180,y, color='red', marker='x') + plt.savefig (plot_fname, bbox_inches='tight') + + logger.info("Successfully created plots for ESMF Mesh file : %s", plot_fname) + + + + + + From 075e6c0d4d963ac64b5b1164ce147575ce709986 Mon Sep 17 00:00:00 2001 From: Negin Sobhani Date: Thu, 3 Nov 2022 14:23:53 -0600 Subject: [PATCH 027/595] updating mesh_type -- working... --- python/ctsm/site_and_regional/mesh_type.py | 32 ++++++++++++++++++- .../ctsm/site_and_regional/regional_case.py | 4 +-- 2 files changed, 33 insertions(+), 3 deletions(-) diff --git a/python/ctsm/site_and_regional/mesh_type.py b/python/ctsm/site_and_regional/mesh_type.py index bcd1199f77..bd66bff875 100644 --- a/python/ctsm/site_and_regional/mesh_type.py +++ b/python/ctsm/site_and_regional/mesh_type.py @@ -277,7 +277,7 @@ def calculate_node_coords(self): ) # -- convert to float32 to find duplicates - corner_pairs = corner_pairs.astype(np.float32, copy=False) + ##corner_pairs = corner_pairs.astype(np.float32, copy=False) # -- remove coordinates that are shared between the elements node_coords = dd.from_dask_array(corner_pairs).drop_duplicates().values @@ -300,6 +300,36 @@ def calculate_node_coords(self): ) logger.warning ("This may result your simulation to crash later.") + # Identify unique nodes + #dtr = np.pi/180 + #index_radius = 1e-6 + #all_unique_nodes = [] + #num_all_nodes = corner_pairs.shape[0] + + #all_unique_nodes = [] + #for i in range(corner_pairs.shape[0]): + # if np.mod(i,100) == 0: + # print('i in num_all_nodes ',i,corner_pairs) + + # dlon = corner_pairs[i+1:,0] - corner_pairs[i,0] + # dlat = corner_pairs[i+1:,1] - corner_pairs[i,1] + # dlon = dtr*dlon + # dlat = dtr*dlat + # dist = np.power(np.sin(dlat/2),2) + (np.cos(dtr*corner_pairs[i+1:,1]) * np.cos(dtr*corner_pairs[i,1]) * np.power(np.sin(dlon/2),2)) + # dist = np.arctan2( np.sqrt(dist), np.sqrt(1-dist)) + # npoints = np.sum(np.where(dist < index_radius,1,0)) + # if npoints < 1: + # all_unique_nodes.append(corner_pairs[i,:]) + + #node_coords = np.vstack(all_unique_nodes) + #self.node_coords = node_coords + + #if node_coords.shape[0] != elem_conn_size: + # print ("khar!") + # sys.exit() + + + def calculate_elem_conn(self): """ Calculate element connectivity (for 'elementConn' in ESMF mesh). diff --git a/python/ctsm/site_and_regional/regional_case.py b/python/ctsm/site_and_regional/regional_case.py index 61d0fdf78e..ef8a53a575 100644 --- a/python/ctsm/site_and_regional/regional_case.py +++ b/python/ctsm/site_and_regional/regional_case.py @@ -271,8 +271,8 @@ def extract_mesh_at_reg(self, ds): lat_name = "lsmlat" lon_name = "lsmlon" - lats = ds[lat_name] - lons = ds[lon_name] + lats = ds[lat_name].astype(np.float32) + lons = ds[lon_name].astype(np.float32) this_mesh = MeshType(lats, lons) this_mesh.calculate_corners() From a3199f1cc3f7472c08a210d00f219b8ade93383b Mon Sep 17 00:00:00 2001 From: Negin Sobhani Date: Mon, 7 Nov 2022 13:27:38 -0700 Subject: [PATCH 028/595] plotting capabilities --- python/ctsm/site_and_regional/mesh_type.py | 116 ++++++++++++++++----- 1 file changed, 90 insertions(+), 26 deletions(-) diff --git a/python/ctsm/site_and_regional/mesh_type.py b/python/ctsm/site_and_regional/mesh_type.py index bd66bff875..08639afc04 100644 --- a/python/ctsm/site_and_regional/mesh_type.py +++ b/python/ctsm/site_and_regional/mesh_type.py @@ -2,6 +2,7 @@ This module includes the definition and functions for defining a Grid or Mesh. This enables creating ESMF mesh file (unstructured grid file)for valid 1D or 2D lats and lons. """ +import os import sys import logging import argparse @@ -19,6 +20,8 @@ import cartopy.crs as ccrs import matplotlib.pyplot as plt import cartopy.feature as cfeature +import matplotlib.patches as mpatches + logger = logging.getLogger(__name__) @@ -372,7 +375,7 @@ def create_esmf(self, mesh_fname, area=None): # -- calculate element connections self.calculate_elem_conn() - center_coords = da.stack( + self.center_coords = da.stack( [ self.center_lon2d.T.reshape((-1,)).T, self.center_lat2d.T.reshape((-1,)).T, @@ -404,7 +407,7 @@ def create_esmf(self, mesh_fname, area=None): attrs={"long_name": "Number of nodes per element"}, ) ds_out["centerCoords"] = xr.DataArray( - center_coords, dims=("elementCount", "coordDim"), attrs={"units": self.unit} + self.center_coords, dims=("elementCount", "coordDim"), attrs={"units": self.unit} ) # -- add mask @@ -446,53 +449,114 @@ def create_esmf(self, mesh_fname, area=None): ds_out.to_netcdf(mesh_fname) logger.info("Successfully created ESMF Mesh file : %s", mesh_fname) - plot_fname = mesh_fname +'.png' - self.make_plot = False + self.make_plot = True if self.make_plot: - self.make_mesh_plot(plot_fname) + plot_regional = ( + os.path.splitext(mesh_fname)[0] + + '_regional' + + ".png") + + plot_global = ( + os.path.splitext(mesh_fname)[0] + + '_global' + + ".png") + + self.make_mesh_plot(plot_regional, plot_global) - def make_mesh_plot(self, plot_fname): - """ - Create an ESMF mesh file for the mesh + + def make_mesh_plot(self, plot_regional, plot_global): + """ + Create a plot for the ESMF mesh file Parameters ---------- - plot_fname : str - The path to write the ESMF meshfile plot + plot_regional : str + The path to write the ESMF meshfile regional plot + plot_global : str + The path to write the ESMF meshfile global plot """ - plt.figure(num=None, figsize=(25, 11), facecolor='w', edgecolor='k') - ax = plt.axes(projection=ccrs.PlateCarree()) + # -- regional plot + plt.figure(num=None, figsize=(15, 13), facecolor='w', edgecolor='k') + ax = plt.axes(projection=ccrs.PlateCarree()) - ax.coastlines() - ax.gridlines(color="black", linestyle="dotted") + ax.add_feature(cfeature.COASTLINE, edgecolor="black") + ax.add_feature(cfeature.BORDERS, edgecolor="black") ax.add_feature(cfeature.OCEAN) ax.add_feature(cfeature.BORDERS) ax.add_feature(cfeature.LAND, edgecolor='black') ax.add_feature(cfeature.LAKES, edgecolor='black') ax.add_feature(cfeature.RIVERS) + ax.gridlines(color="black", linestyle="dotted", draw_labels=True,) + + #-- plot corner coordinates + clon, clat = self.node_coords.T.compute() + df = pd.DataFrame({'lon':list(clon), 'lat':list(clat)}) - x, y = self.node_coords.T - print (x.shape) - print (y.shape) - df = pd.DataFrame({'x':x, 'y':y}, index = [0]) - grouped_x = df.groupby('x') + grouped_x = df.groupby('lon') for name, group in grouped_x: - plt.plot (group['x']-180, group['y'],color='black', marker='o') - grouped_y = df.groupby('y') + plt.plot (group['lon'], group['lat'],color='black', linewidth = + 0.5, transform = ccrs.PlateCarree()) + grouped_y = df.groupby('lat') for name, group in grouped_y: - plt.plot (group['x']-180, group['y'],color='black', marker='o') + plt.plot (group['lon'], group['lat'],color='black', linewidth = + 0.5, transform = ccrs.PlateCarree()) + + #-- plot center coordinates + clon, clat = self.center_coords.T.compute() + + plt.scatter (clon,clat, color='tomato', marker='x') + lc_colors = { + 'Corner Coordinates': "black", # value=0 + 'Center Coordinates': "tomato", # value=1 + + } + labels, handles = zip(*[(k, mpatches.Rectangle((0, 0), 1, 1, facecolor=v)) for k,v in lc_colors.items()]) + + ax.legend(handles, labels) - x, y = this_mesh.center_coords.T + plt.savefig (plot_regional, bbox_inches='tight') - plt.scatter (x-180,y, color='red', marker='x') - plt.savefig (plot_fname, bbox_inches='tight') + logger.info("Successfully created regional plots for ESMF Mesh file : %s", plot_regional) - logger.info("Successfully created plots for ESMF Mesh file : %s", plot_fname) + # -- global plot + fig = plt.figure(num=None, figsize=(15, 10), facecolor='w', edgecolor='k') + + ax = fig.add_subplot(1,1,1, projection=ccrs.Robinson()) + + ax.add_feature(cfeature.COASTLINE, edgecolor="black") + ax.add_feature(cfeature.BORDERS, edgecolor="black") + ax.add_feature(cfeature.OCEAN) + ax.add_feature(cfeature.BORDERS) + ax.add_feature(cfeature.LAND, edgecolor='black') + ax.add_feature(cfeature.LAKES, edgecolor='black') + ax.add_feature(cfeature.RIVERS) + ax.gridlines(color="black", linestyle="dotted", draw_labels=True,) + ax.set_global() + + #-- plot corner coordinates + clon, clat = self.node_coords.T.compute() + df = pd.DataFrame({'lon':list(clon), 'lat':list(clat)}) + + grouped_x = df.groupby('lon') + for name, group in grouped_x: + plt.plot (group['lon'], group['lat'],color='black', transform=ccrs.PlateCarree(), linewidth = + 0.5) + grouped_y = df.groupby('lat') + for name, group in grouped_y: + plt.plot (group['lon'], group['lat'],color='black', transform=ccrs.PlateCarree(), linewidth = + 0.5) + + #-- plot center coordinates + clon, clat = self.center_coords.T.compute() + plt.scatter (clon,clat, color='tomato', marker='o',s = 1, transform=ccrs.PlateCarree()) + ax.legend(handles, labels) + plt.savefig (plot_global, bbox_inches='tight') + logger.info("Successfully created regional plots for ESMF Mesh file : %s", plot_global) From 8ef1a086e46c5d9a42a3d4960a5fba1923ad518a Mon Sep 17 00:00:00 2001 From: Negin Sobhani Date: Mon, 7 Nov 2022 14:38:41 -0700 Subject: [PATCH 029/595] adding high res as default --- tools/site_and_regional/default_data.cfg | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tools/site_and_regional/default_data.cfg b/tools/site_and_regional/default_data.cfg index c38a00adcc..bc57ace149 100644 --- a/tools/site_and_regional/default_data.cfg +++ b/tools/site_and_regional/default_data.cfg @@ -16,8 +16,8 @@ tpqwname = CLMGSWP3v1.TPQW [surfdat] dir = lnd/clm2/surfdata_map/release-clm5.0.18 -surfdat_16pft = surfdata_4x5_hist_16pfts_Irrig_CMIP6_simyr2000_c190214.nc -surfdat_78pft = surfdata_4x5_hist_78pfts_CMIP6_simyr2000_c190214.nc +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 = fv4x5_050615_polemod_ESMFmesh.nc From e7826c228d9a634936bd0f326f93807a756ab1e5 Mon Sep 17 00:00:00 2001 From: Negin Sobhani Date: Mon, 7 Nov 2022 15:01:25 -0700 Subject: [PATCH 030/595] high res as default --- tools/site_and_regional/default_data.cfg | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tools/site_and_regional/default_data.cfg b/tools/site_and_regional/default_data.cfg index bc57ace149..630013d502 100644 --- a/tools/site_and_regional/default_data.cfg +++ b/tools/site_and_regional/default_data.cfg @@ -19,7 +19,7 @@ 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 = fv4x5_050615_polemod_ESMFmesh.nc +mesh_surf = fv0.9x1.25_141008_ESMFmesh.nc [landuse] dir = lnd/clm2/surfdata_map/release-clm5.0.18 From 129492f828709333fd9d76f964680ee360313257 Mon Sep 17 00:00:00 2001 From: Negin Sobhani Date: Mon, 7 Nov 2022 15:02:54 -0700 Subject: [PATCH 031/595] blackify --- python/ctsm/site_and_regional/mesh_type.py | 182 ++++++++++++--------- 1 file changed, 105 insertions(+), 77 deletions(-) diff --git a/python/ctsm/site_and_regional/mesh_type.py b/python/ctsm/site_and_regional/mesh_type.py index 08639afc04..a330413f02 100644 --- a/python/ctsm/site_and_regional/mesh_type.py +++ b/python/ctsm/site_and_regional/mesh_type.py @@ -25,6 +25,7 @@ logger = logging.getLogger(__name__) + class MeshType: """ An object for describing mesh or grid files. @@ -100,8 +101,8 @@ def check_lat_lon_dims(self): Unrecognized grid! \n The dimension of latitude should be either 1 or 2 but it is {}. """.format( - self.lat_dims - ) + self.lat_dims + ) raise argparse.ArgumentTypeError(err_msg) if self.lon_dims not in [1, 2]: @@ -109,8 +110,8 @@ def check_lat_lon_dims(self): Unrecognized grid! \n The dimension of longitude should be either 1 or 2 but it is {}. """.format( - self.lon_dims - ) + self.lon_dims + ) raise argparse.ArgumentTypeError(err_msg) def create_artificial_mask(self): @@ -127,9 +128,7 @@ def create_artificial_mask(self): mask = np.ones([lons_size, lats_size], dtype=np.int8) elif self.lat_dims == 2: # -- 2D mask - mask = np.ones( - self.center_lats.shape, dtype=np.int8 - ) + mask = np.ones(self.center_lats.shape, dtype=np.int8) mask_da = da.from_array(mask) self.mask = mask_da @@ -294,26 +293,26 @@ def calculate_node_coords(self): self.node_coords = node_coords - # -- error check to avoid issues later + # -- error check to avoid issues later if self.node_coords.shape[0] != elem_conn_size: logger.warning( "The size of unique coordinate pairs is {} but expected size is {}!".format( self.node_coords.shape[0], elem_conn_size ) ) - logger.warning ("This may result your simulation to crash later.") - - # Identify unique nodes - #dtr = np.pi/180 - #index_radius = 1e-6 - #all_unique_nodes = [] - #num_all_nodes = corner_pairs.shape[0] - - #all_unique_nodes = [] - #for i in range(corner_pairs.shape[0]): + logger.warning("This may result your simulation to crash later.") + + # Identify unique nodes + # dtr = np.pi/180 + # index_radius = 1e-6 + # all_unique_nodes = [] + # num_all_nodes = corner_pairs.shape[0] + + # all_unique_nodes = [] + # for i in range(corner_pairs.shape[0]): # if np.mod(i,100) == 0: # print('i in num_all_nodes ',i,corner_pairs) - + # dlon = corner_pairs[i+1:,0] - corner_pairs[i,0] # dlat = corner_pairs[i+1:,1] - corner_pairs[i,1] # dlon = dtr*dlon @@ -323,15 +322,13 @@ def calculate_node_coords(self): # npoints = np.sum(np.where(dist < index_radius,1,0)) # if npoints < 1: # all_unique_nodes.append(corner_pairs[i,:]) - - #node_coords = np.vstack(all_unique_nodes) - #self.node_coords = node_coords - - #if node_coords.shape[0] != elem_conn_size: - # print ("khar!") - # sys.exit() + # node_coords = np.vstack(all_unique_nodes) + # self.node_coords = node_coords + # if node_coords.shape[0] != elem_conn_size: + # print ("khar!") + # sys.exit() def calculate_elem_conn(self): """ @@ -407,7 +404,9 @@ def create_esmf(self, mesh_fname, area=None): attrs={"long_name": "Number of nodes per element"}, ) ds_out["centerCoords"] = xr.DataArray( - self.center_coords, dims=("elementCount", "coordDim"), attrs={"units": self.unit} + self.center_coords, + dims=("elementCount", "coordDim"), + attrs={"units": self.unit}, ) # -- add mask @@ -435,7 +434,7 @@ def create_esmf(self, mesh_fname, area=None): if self.mesh_name: ds_out.attrs["title"] = "ESMF unstructured grid file " + self.mesh_name else: - ds_out.attrs["title"] = "ESMF unstructured grid file " + ds_out.attrs["title"] = "ESMF unstructured grid file " ds_out.attrs["gridType"] = "unstructured mesh" ds_out.attrs["version"] = "0.9" ds_out.attrs["conventions"] = "ESMFMESH" @@ -452,20 +451,13 @@ def create_esmf(self, mesh_fname, area=None): self.make_plot = True if self.make_plot: - plot_regional = ( - os.path.splitext(mesh_fname)[0] - + '_regional' - + ".png") + plot_regional = os.path.splitext(mesh_fname)[0] + "_regional" + ".png" - plot_global = ( - os.path.splitext(mesh_fname)[0] - + '_global' - + ".png") + plot_global = os.path.splitext(mesh_fname)[0] + "_global" + ".png" self.make_mesh_plot(plot_regional, plot_global) - - def make_mesh_plot(self, plot_regional, plot_global): + def make_mesh_plot(self, plot_regional, plot_global): """ Create a plot for the ESMF mesh file @@ -478,85 +470,121 @@ def make_mesh_plot(self, plot_regional, plot_global): """ # -- regional plot - plt.figure(num=None, figsize=(15, 13), facecolor='w', edgecolor='k') + plt.figure(num=None, figsize=(15, 13), facecolor="w", edgecolor="k") ax = plt.axes(projection=ccrs.PlateCarree()) ax.add_feature(cfeature.COASTLINE, edgecolor="black") ax.add_feature(cfeature.BORDERS, edgecolor="black") ax.add_feature(cfeature.OCEAN) ax.add_feature(cfeature.BORDERS) - ax.add_feature(cfeature.LAND, edgecolor='black') - ax.add_feature(cfeature.LAKES, edgecolor='black') + ax.add_feature(cfeature.LAND, edgecolor="black") + ax.add_feature(cfeature.LAKES, edgecolor="black") ax.add_feature(cfeature.RIVERS) - ax.gridlines(color="black", linestyle="dotted", draw_labels=True,) + ax.gridlines( + color="black", + linestyle="dotted", + draw_labels=True, + ) - #-- plot corner coordinates + # -- plot corner coordinates clon, clat = self.node_coords.T.compute() - df = pd.DataFrame({'lon':list(clon), 'lat':list(clat)}) + df = pd.DataFrame({"lon": list(clon), "lat": list(clat)}) - grouped_x = df.groupby('lon') + grouped_x = df.groupby("lon") for name, group in grouped_x: - plt.plot (group['lon'], group['lat'],color='black', linewidth = - 0.5, transform = ccrs.PlateCarree()) - grouped_y = df.groupby('lat') + plt.plot( + group["lon"], + group["lat"], + color="black", + linewidth=0.5, + transform=ccrs.PlateCarree(), + ) + grouped_y = df.groupby("lat") for name, group in grouped_y: - plt.plot (group['lon'], group['lat'],color='black', linewidth = - 0.5, transform = ccrs.PlateCarree()) + plt.plot( + group["lon"], + group["lat"], + color="black", + linewidth=0.5, + transform=ccrs.PlateCarree(), + ) - #-- plot center coordinates + # -- plot center coordinates clon, clat = self.center_coords.T.compute() - plt.scatter (clon,clat, color='tomato', marker='x') + plt.scatter(clon, clat, color="tomato", marker="x") lc_colors = { - 'Corner Coordinates': "black", # value=0 - 'Center Coordinates': "tomato", # value=1 - + "Corner Coordinates": "black", # value=0 + "Center Coordinates": "tomato", # value=1 } - labels, handles = zip(*[(k, mpatches.Rectangle((0, 0), 1, 1, facecolor=v)) for k,v in lc_colors.items()]) + labels, handles = zip( + *[ + (k, mpatches.Rectangle((0, 0), 1, 1, facecolor=v)) + for k, v in lc_colors.items() + ] + ) ax.legend(handles, labels) - plt.savefig (plot_regional, bbox_inches='tight') - - logger.info("Successfully created regional plots for ESMF Mesh file : %s", plot_regional) + plt.savefig(plot_regional, bbox_inches="tight") + logger.info( + "Successfully created regional plots for ESMF Mesh file : %s", plot_regional + ) # -- global plot - fig = plt.figure(num=None, figsize=(15, 10), facecolor='w', edgecolor='k') + fig = plt.figure(num=None, figsize=(15, 10), facecolor="w", edgecolor="k") - ax = fig.add_subplot(1,1,1, projection=ccrs.Robinson()) + ax = fig.add_subplot(1, 1, 1, projection=ccrs.Robinson()) ax.add_feature(cfeature.COASTLINE, edgecolor="black") ax.add_feature(cfeature.BORDERS, edgecolor="black") ax.add_feature(cfeature.OCEAN) ax.add_feature(cfeature.BORDERS) - ax.add_feature(cfeature.LAND, edgecolor='black') - ax.add_feature(cfeature.LAKES, edgecolor='black') + ax.add_feature(cfeature.LAND, edgecolor="black") + ax.add_feature(cfeature.LAKES, edgecolor="black") ax.add_feature(cfeature.RIVERS) - ax.gridlines(color="black", linestyle="dotted", draw_labels=True,) + ax.gridlines( + color="black", + linestyle="dotted", + draw_labels=True, + ) ax.set_global() - #-- plot corner coordinates + # -- plot corner coordinates clon, clat = self.node_coords.T.compute() - df = pd.DataFrame({'lon':list(clon), 'lat':list(clat)}) + df = pd.DataFrame({"lon": list(clon), "lat": list(clat)}) - grouped_x = df.groupby('lon') + grouped_x = df.groupby("lon") for name, group in grouped_x: - plt.plot (group['lon'], group['lat'],color='black', transform=ccrs.PlateCarree(), linewidth = - 0.5) - grouped_y = df.groupby('lat') + plt.plot( + group["lon"], + group["lat"], + color="black", + transform=ccrs.PlateCarree(), + linewidth=0.5, + ) + grouped_y = df.groupby("lat") for name, group in grouped_y: - plt.plot (group['lon'], group['lat'],color='black', transform=ccrs.PlateCarree(), linewidth = - 0.5) + plt.plot( + group["lon"], + group["lat"], + color="black", + transform=ccrs.PlateCarree(), + linewidth=0.5, + ) - #-- plot center coordinates + # -- plot center coordinates clon, clat = self.center_coords.T.compute() - plt.scatter (clon,clat, color='tomato', marker='o',s = 1, transform=ccrs.PlateCarree()) + plt.scatter( + clon, clat, color="tomato", marker="o", s=1, transform=ccrs.PlateCarree() + ) ax.legend(handles, labels) - plt.savefig (plot_global, bbox_inches='tight') - - logger.info("Successfully created regional plots for ESMF Mesh file : %s", plot_global) + plt.savefig(plot_global, bbox_inches="tight") + logger.info( + "Successfully created regional plots for ESMF Mesh file : %s", plot_global + ) From 8c8279b9e75b68d70af16e9583b63fb052d1cfe1 Mon Sep 17 00:00:00 2001 From: Negin Sobhani Date: Mon, 7 Nov 2022 15:08:10 -0700 Subject: [PATCH 032/595] extra comments --- python/ctsm/site_and_regional/mesh_type.py | 28 ---------------------- 1 file changed, 28 deletions(-) diff --git a/python/ctsm/site_and_regional/mesh_type.py b/python/ctsm/site_and_regional/mesh_type.py index a330413f02..84ccecdf83 100644 --- a/python/ctsm/site_and_regional/mesh_type.py +++ b/python/ctsm/site_and_regional/mesh_type.py @@ -302,34 +302,6 @@ def calculate_node_coords(self): ) logger.warning("This may result your simulation to crash later.") - # Identify unique nodes - # dtr = np.pi/180 - # index_radius = 1e-6 - # all_unique_nodes = [] - # num_all_nodes = corner_pairs.shape[0] - - # all_unique_nodes = [] - # for i in range(corner_pairs.shape[0]): - # if np.mod(i,100) == 0: - # print('i in num_all_nodes ',i,corner_pairs) - - # dlon = corner_pairs[i+1:,0] - corner_pairs[i,0] - # dlat = corner_pairs[i+1:,1] - corner_pairs[i,1] - # dlon = dtr*dlon - # dlat = dtr*dlat - # dist = np.power(np.sin(dlat/2),2) + (np.cos(dtr*corner_pairs[i+1:,1]) * np.cos(dtr*corner_pairs[i,1]) * np.power(np.sin(dlon/2),2)) - # dist = np.arctan2( np.sqrt(dist), np.sqrt(1-dist)) - # npoints = np.sum(np.where(dist < index_radius,1,0)) - # if npoints < 1: - # all_unique_nodes.append(corner_pairs[i,:]) - - # node_coords = np.vstack(all_unique_nodes) - # self.node_coords = node_coords - - # if node_coords.shape[0] != elem_conn_size: - # print ("khar!") - # sys.exit() - def calculate_elem_conn(self): """ Calculate element connectivity (for 'elementConn' in ESMF mesh). From 2816f95b936b14f00b138431b6ea4982015deac8 Mon Sep 17 00:00:00 2001 From: Negin Sobhani Date: Mon, 7 Nov 2022 15:52:04 -0700 Subject: [PATCH 033/595] quick updates --- tools/site_and_regional/mesh_maker.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tools/site_and_regional/mesh_maker.py b/tools/site_and_regional/mesh_maker.py index 3229bea2e9..4ce8642f96 100755 --- a/tools/site_and_regional/mesh_maker.py +++ b/tools/site_and_regional/mesh_maker.py @@ -68,7 +68,7 @@ def get_parser(): parser.add_argument( "--lat", - help="Name of latitude varibale on netcdf input file. If none given, looks to find variables that include 'lat'.", + help="Name of latitude varibale on netcdf input file.", action="store", dest="lat_name", type=str, @@ -76,7 +76,7 @@ def get_parser(): ) parser.add_argument( "--lon", - help="Name of latitude varibale on netcdf input file. If none given, looks to find variables that include 'lon'.", + help="Name of latitude varibale on netcdf input file.", action="store", dest="lon_name", type=str, From f614b60bff8f2a01dccd31ea019aa0982d9d2ee8 Mon Sep 17 00:00:00 2001 From: Negin Sobhani Date: Mon, 7 Nov 2022 17:01:09 -0700 Subject: [PATCH 034/595] efficient plotting for regions --- python/ctsm/site_and_regional/mesh_type.py | 78 +++++++++------------- 1 file changed, 33 insertions(+), 45 deletions(-) diff --git a/python/ctsm/site_and_regional/mesh_type.py b/python/ctsm/site_and_regional/mesh_type.py index 84ccecdf83..aa157f30c2 100644 --- a/python/ctsm/site_and_regional/mesh_type.py +++ b/python/ctsm/site_and_regional/mesh_type.py @@ -459,32 +459,26 @@ def make_mesh_plot(self, plot_regional, plot_global): ) # -- plot corner coordinates - clon, clat = self.node_coords.T.compute() - df = pd.DataFrame({"lon": list(clon), "lat": list(clat)}) - - grouped_x = df.groupby("lon") - for name, group in grouped_x: - plt.plot( - group["lon"], - group["lat"], - color="black", - linewidth=0.5, - transform=ccrs.PlateCarree(), - ) - grouped_y = df.groupby("lat") - for name, group in grouped_y: - plt.plot( - group["lon"], - group["lat"], - color="black", - linewidth=0.5, - transform=ccrs.PlateCarree(), - ) + clats, clons = self.node_coords.T.compute() + elem_conn_vals = self.elem_conn.compute() + element_counts = elem_conn_vals.shape[0] + + for index in range(element_counts): + conns = [ int(x)-1 for x in elem_conn_vals[index]] + + lat_corners = clats[conns] + lon_corners = clons[conns] + poly_corners = np.zeros((len(lat_corners), 2)) + poly_corners[:,1] = lon_corners + poly_corners[:,0] = lat_corners + poly = mpatches.Polygon(poly_corners, closed=True, ec='black', lw=1, transform=ccrs.PlateCarree(), zorder = + 10, facecolor='none') + ax.add_patch(poly) # -- plot center coordinates clon, clat = self.center_coords.T.compute() - plt.scatter(clon, clat, color="tomato", marker="x") + ax.scatter(clon, clat, color="tomato", marker="x", transform=ccrs.PlateCarree(), zorder = 11) lc_colors = { "Corner Coordinates": "black", # value=0 "Center Coordinates": "tomato", # value=1 @@ -524,33 +518,27 @@ def make_mesh_plot(self, plot_regional, plot_global): ax.set_global() # -- plot corner coordinates - clon, clat = self.node_coords.T.compute() - df = pd.DataFrame({"lon": list(clon), "lat": list(clat)}) - - grouped_x = df.groupby("lon") - for name, group in grouped_x: - plt.plot( - group["lon"], - group["lat"], - color="black", - transform=ccrs.PlateCarree(), - linewidth=0.5, - ) - grouped_y = df.groupby("lat") - for name, group in grouped_y: - plt.plot( - group["lon"], - group["lat"], - color="black", - transform=ccrs.PlateCarree(), - linewidth=0.5, - ) + clats, clons = self.node_coords.T.compute() + elem_conn_vals = self.elem_conn.compute() + element_counts = elem_conn_vals.shape[0] + + for index in range(element_counts): + conns = [ int(x)-1 for x in elem_conn_vals[index]] + + lat_corners = clats[conns] + lon_corners = clons[conns] + poly_corners = np.zeros((len(lat_corners), 2)) + poly_corners[:,1] = lon_corners + poly_corners[:,0] = lat_corners + poly = mpatches.Polygon(poly_corners, closed=True, ec='black', transform=ccrs.PlateCarree(), zorder = + 10, facecolor='none', linewidth=0.5) + ax.add_patch(poly) # -- plot center coordinates clon, clat = self.center_coords.T.compute() - plt.scatter( - clon, clat, color="tomato", marker="o", s=1, transform=ccrs.PlateCarree() + ax.scatter( + clon, clat, color="tomato", marker="o", s=1, transform=ccrs.PlateCarree(), zorder = 11 ) ax.legend(handles, labels) From 796bdeabba4354b21a37d46274b63d5177f1b95a Mon Sep 17 00:00:00 2001 From: Negin Sobhani Date: Mon, 7 Nov 2022 18:35:42 -0700 Subject: [PATCH 035/595] quick clean ups --- python/ctsm/site_and_regional/mesh_type.py | 53 ++++++++++++++++------ 1 file changed, 39 insertions(+), 14 deletions(-) diff --git a/python/ctsm/site_and_regional/mesh_type.py b/python/ctsm/site_and_regional/mesh_type.py index aa157f30c2..f0bc10dfef 100644 --- a/python/ctsm/site_and_regional/mesh_type.py +++ b/python/ctsm/site_and_regional/mesh_type.py @@ -16,13 +16,11 @@ from dask.diagnostics import ProgressBar # -- libraries for plotting mesh (make_mesh_plot) -import matplotlib.pyplot as plt import cartopy.crs as ccrs import matplotlib.pyplot as plt import cartopy.feature as cfeature import matplotlib.patches as mpatches - logger = logging.getLogger(__name__) @@ -464,21 +462,35 @@ def make_mesh_plot(self, plot_regional, plot_global): element_counts = elem_conn_vals.shape[0] for index in range(element_counts): - conns = [ int(x)-1 for x in elem_conn_vals[index]] + conns = [int(x) - 1 for x in elem_conn_vals[index]] lat_corners = clats[conns] lon_corners = clons[conns] poly_corners = np.zeros((len(lat_corners), 2)) - poly_corners[:,1] = lon_corners - poly_corners[:,0] = lat_corners - poly = mpatches.Polygon(poly_corners, closed=True, ec='black', lw=1, transform=ccrs.PlateCarree(), zorder = - 10, facecolor='none') + poly_corners[:, 1] = lon_corners + poly_corners[:, 0] = lat_corners + poly = mpatches.Polygon( + poly_corners, + closed=True, + ec="black", + lw=1, + transform=ccrs.PlateCarree(), + zorder=10, + facecolor="none", + ) ax.add_patch(poly) # -- plot center coordinates clon, clat = self.center_coords.T.compute() - ax.scatter(clon, clat, color="tomato", marker="x", transform=ccrs.PlateCarree(), zorder = 11) + ax.scatter( + clon, + clat, + color="tomato", + marker="x", + transform=ccrs.PlateCarree(), + zorder=11, + ) lc_colors = { "Corner Coordinates": "black", # value=0 "Center Coordinates": "tomato", # value=1 @@ -523,22 +535,35 @@ def make_mesh_plot(self, plot_regional, plot_global): element_counts = elem_conn_vals.shape[0] for index in range(element_counts): - conns = [ int(x)-1 for x in elem_conn_vals[index]] + conns = [int(x) - 1 for x in elem_conn_vals[index]] lat_corners = clats[conns] lon_corners = clons[conns] poly_corners = np.zeros((len(lat_corners), 2)) - poly_corners[:,1] = lon_corners - poly_corners[:,0] = lat_corners - poly = mpatches.Polygon(poly_corners, closed=True, ec='black', transform=ccrs.PlateCarree(), zorder = - 10, facecolor='none', linewidth=0.5) + poly_corners[:, 1] = lon_corners + poly_corners[:, 0] = lat_corners + poly = mpatches.Polygon( + poly_corners, + closed=True, + ec="black", + transform=ccrs.PlateCarree(), + zorder=10, + facecolor="none", + linewidth=0.5, + ) ax.add_patch(poly) # -- plot center coordinates clon, clat = self.center_coords.T.compute() ax.scatter( - clon, clat, color="tomato", marker="o", s=1, transform=ccrs.PlateCarree(), zorder = 11 + clon, + clat, + color="tomato", + marker="o", + s=1, + transform=ccrs.PlateCarree(), + zorder=11, ) ax.legend(handles, labels) From fcb3405b3092e4aaa6a35771472f43f5174b478f Mon Sep 17 00:00:00 2001 From: Negin Sobhani Date: Mon, 7 Nov 2022 19:03:35 -0700 Subject: [PATCH 036/595] pylint --- python/ctsm/site_and_regional/regional_case.py | 11 +++++++---- 1 file changed, 7 insertions(+), 4 deletions(-) diff --git a/python/ctsm/site_and_regional/regional_case.py b/python/ctsm/site_and_regional/regional_case.py index ef8a53a575..780ecbe866 100644 --- a/python/ctsm/site_and_regional/regional_case.py +++ b/python/ctsm/site_and_regional/regional_case.py @@ -6,11 +6,11 @@ import logging import os import argparse +from datetime import datetime # -- 3rd party libraries import numpy as np import xarray as xr -from datetime import datetime # -- import local classes for this script from ctsm.site_and_regional.base_case import BaseCase, USRDAT_DIR @@ -83,6 +83,9 @@ class RegionalCase(BaseCase): write out xml commands to a file for usermods (i.e. shell_commands) for regional settings. """ + # pylint: disable=too-many-instance-attributes + # the ones we have are useful + def __init__( self, lat1, @@ -262,7 +265,7 @@ def create_surfdata_at_reg(self, indir, file, user_mods_dir): logger.info("creating mesh file from surface_dataset: %s", wfile) self.extract_mesh_at_reg(f_out) - def extract_mesh_at_reg(self, ds): + def extract_mesh_at_reg(self, ds_in): """ Create Mesh from Surface dataset netcdf file. """ @@ -271,8 +274,8 @@ def extract_mesh_at_reg(self, ds): lat_name = "lsmlat" lon_name = "lsmlon" - lats = ds[lat_name].astype(np.float32) - lons = ds[lon_name].astype(np.float32) + lats = ds_in[lat_name].astype(np.float32) + lons = ds_in[lon_name].astype(np.float32) this_mesh = MeshType(lats, lons) this_mesh.calculate_corners() From 0b6dc1d179c3617761e97a65b597a33588456ce4 Mon Sep 17 00:00:00 2001 From: Negin Sobhani Date: Mon, 7 Nov 2022 19:22:13 -0700 Subject: [PATCH 037/595] pylint issues --- .../ctsm/site_and_regional/regional_case.py | 19 +++++++++---------- 1 file changed, 9 insertions(+), 10 deletions(-) diff --git a/python/ctsm/site_and_regional/regional_case.py b/python/ctsm/site_and_regional/regional_case.py index 780ecbe866..16d834c847 100644 --- a/python/ctsm/site_and_regional/regional_case.py +++ b/python/ctsm/site_and_regional/regional_case.py @@ -119,6 +119,7 @@ def __init__( self.lon2 = lon2 self.reg_name = reg_name self.create_mesh = create_mesh + self.mesh = None self.out_dir = out_dir self.check_region_bounds() self.create_tag() @@ -225,8 +226,7 @@ def create_surfdata_at_reg(self, indir, file, user_mods_dir): logger.info("fsurf_in: %s", fsurf_in) logger.info("fsurf_out: %s", os.path.join(self.out_dir, fsurf_out)) - self.fsurf_out = os.path.join(self.out_dir, fsurf_out) - self.fsurf_out = fsurf_out + fsurf_out = os.path.join(self.out_dir, fsurf_out) if self.create_mesh: mesh_out = os.path.join( @@ -357,7 +357,7 @@ def create_mesh_at_reg(self, mesh_dir, mesh_surf): f_in = xr.open_dataset(mesh_in) self.write_mesh( - f_in, node_coords, subset_element, subset_node, conn_dict, mesh_out + f_in, subset_element, subset_node, conn_dict, mesh_out ) def subset_mesh_at_reg(self, mesh_in): @@ -368,7 +368,6 @@ def subset_mesh_at_reg(self, 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"] @@ -423,8 +422,9 @@ def subset_mesh_at_reg(self, mesh_in): return node_coords, subset_element, subset_node, conn_dict + @staticmethod def write_mesh( - self, f_in, node_coords, subset_element, subset_node, conn_dict, mesh_out + f_in, subset_element, subset_node, conn_dict, mesh_out ): """ This function writes out the subsetted mesh file. @@ -432,7 +432,6 @@ def write_mesh( corner_pairs = f_in.variables["nodeCoords"][ subset_node, ] - dimensions = f_in.dims variables = f_in.variables global_attributes = f_in.attrs @@ -444,10 +443,10 @@ def write_mesh( 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] + for n in range(elem_count): + for m in range(max_node_dim): + ndx = int(elem_conn_index[n, m]) + elem_conn_out[n, m] = conn_dict[ndx] num_elem_conn_out = np.empty( shape=[ From ae82c8404c8d44a6d6cff8cbe0370f617b1fd726 Mon Sep 17 00:00:00 2001 From: Negin Sobhani Date: Mon, 7 Nov 2022 19:30:47 -0700 Subject: [PATCH 038/595] pylint complaints --- python/ctsm/site_and_regional/regional_case.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/python/ctsm/site_and_regional/regional_case.py b/python/ctsm/site_and_regional/regional_case.py index 16d834c847..93b31f8102 100644 --- a/python/ctsm/site_and_regional/regional_case.py +++ b/python/ctsm/site_and_regional/regional_case.py @@ -357,7 +357,7 @@ def create_mesh_at_reg(self, mesh_dir, mesh_surf): f_in = xr.open_dataset(mesh_in) self.write_mesh( - f_in, subset_element, subset_node, conn_dict, mesh_out + f_in, node_coords, subset_element, subset_node, conn_dict, mesh_out ) def subset_mesh_at_reg(self, mesh_in): @@ -424,8 +424,9 @@ def subset_mesh_at_reg(self, mesh_in): @staticmethod def write_mesh( - f_in, subset_element, subset_node, conn_dict, mesh_out + f_in, node_coords, subset_element, subset_node, conn_dict, mesh_out ): + # pylint: disable=unused-argument """ This function writes out the subsetted mesh file. """ From 05e8facf86781cc6f2f888f4a2cb75c4a79e9854 Mon Sep 17 00:00:00 2001 From: Negin Sobhani Date: Mon, 7 Nov 2022 19:31:05 -0700 Subject: [PATCH 039/595] pylint valid names such as ax --- python/ctsm/.pylintrc | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/python/ctsm/.pylintrc b/python/ctsm/.pylintrc index bc7ae54dd2..472bcb662a 100644 --- a/python/ctsm/.pylintrc +++ b/python/ctsm/.pylintrc @@ -427,6 +427,11 @@ good-names=i, k, ex, Run, + m, + n, + l1, + l2, + ax, _, # --- default list is above here, our own list is below here --- # Allow logger as a global name in each module, because this seems to follow general recommended convention: From 358834295a2377bcd8edf03200f10a91ebff63fd Mon Sep 17 00:00:00 2001 From: Negin Sobhani Date: Mon, 7 Nov 2022 19:31:32 -0700 Subject: [PATCH 040/595] pylint complaints --- python/ctsm/site_and_regional/mesh_type.py | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/python/ctsm/site_and_regional/mesh_type.py b/python/ctsm/site_and_regional/mesh_type.py index f0bc10dfef..7acba30ed5 100644 --- a/python/ctsm/site_and_regional/mesh_type.py +++ b/python/ctsm/site_and_regional/mesh_type.py @@ -3,22 +3,20 @@ This enables creating ESMF mesh file (unstructured grid file)for valid 1D or 2D lats and lons. """ import os -import sys import logging import argparse import datetime import numpy as np import xarray as xr -import pandas as pd import dask.array as da import dask.dataframe as dd from dask.diagnostics import ProgressBar # -- libraries for plotting mesh (make_mesh_plot) import cartopy.crs as ccrs -import matplotlib.pyplot as plt import cartopy.feature as cfeature +import matplotlib.pyplot as plt import matplotlib.patches as mpatches logger = logging.getLogger(__name__) From 38352a85b56343c0c601f4d9141c98a96be9ee1b Mon Sep 17 00:00:00 2001 From: Negin Sobhani Date: Tue, 8 Nov 2022 03:30:16 -0700 Subject: [PATCH 041/595] formatting --- python/ctsm/site_and_regional/mesh_type.py | 25 ++--- .../ctsm/site_and_regional/regional_case.py | 43 ++------ python/ctsm/subset_data.py | 99 ++++++++++--------- 3 files changed, 70 insertions(+), 97 deletions(-) diff --git a/python/ctsm/site_and_regional/mesh_type.py b/python/ctsm/site_and_regional/mesh_type.py index 7acba30ed5..3b2dda1dc8 100644 --- a/python/ctsm/site_and_regional/mesh_type.py +++ b/python/ctsm/site_and_regional/mesh_type.py @@ -176,16 +176,12 @@ def calculate_corners(self, unit="degrees"): # -- for the edge rows/columns. padded_lat2d = da.from_array( - np.pad( - self.center_lat2d.compute(), (1, 1), mode="reflect", reflect_type="odd" - ) + np.pad(self.center_lat2d.compute(), (1, 1), mode="reflect", reflect_type="odd") ) # -- pad center_lons for calculating edge grids padded_lon2d = da.from_array( - np.pad( - self.center_lon2d.compute(), (1, 1), mode="reflect", reflect_type="odd" - ) + np.pad(self.center_lon2d.compute(), (1, 1), mode="reflect", reflect_type="odd") ) # -- calculate corner lats for each grid @@ -406,9 +402,7 @@ def create_esmf(self, mesh_fname, area=None): ds_out.attrs["gridType"] = "unstructured mesh" ds_out.attrs["version"] = "0.9" ds_out.attrs["conventions"] = "ESMFMESH" - ds_out.attrs["date_created"] = datetime.datetime.now().strftime( - "%Y-%m-%d %H:%M:%S" - ) + ds_out.attrs["date_created"] = datetime.datetime.now().strftime("%Y-%m-%d %H:%M:%S") # -- write Xarray dataset to file if mesh_fname is not None: @@ -494,19 +488,14 @@ def make_mesh_plot(self, plot_regional, plot_global): "Center Coordinates": "tomato", # value=1 } labels, handles = zip( - *[ - (k, mpatches.Rectangle((0, 0), 1, 1, facecolor=v)) - for k, v in lc_colors.items() - ] + *[(k, mpatches.Rectangle((0, 0), 1, 1, facecolor=v)) for k, v in lc_colors.items()] ) ax.legend(handles, labels) plt.savefig(plot_regional, bbox_inches="tight") - logger.info( - "Successfully created regional plots for ESMF Mesh file : %s", plot_regional - ) + logger.info("Successfully created regional plots for ESMF Mesh file : %s", plot_regional) # -- global plot fig = plt.figure(num=None, figsize=(15, 10), facecolor="w", edgecolor="k") @@ -568,6 +557,4 @@ def make_mesh_plot(self, plot_regional, plot_global): plt.savefig(plot_global, bbox_inches="tight") - logger.info( - "Successfully created regional plots for ESMF Mesh file : %s", plot_global - ) + logger.info("Successfully created regional plots for ESMF Mesh file : %s", plot_global) diff --git a/python/ctsm/site_and_regional/regional_case.py b/python/ctsm/site_and_regional/regional_case.py index d26d964332..dca050b217 100644 --- a/python/ctsm/site_and_regional/regional_case.py +++ b/python/ctsm/site_and_regional/regional_case.py @@ -325,9 +325,7 @@ def create_mesh_at_reg(self, mesh_dir, mesh_surf): """ Create a mesh subsetted for the RegionalCase class. """ - logger.info( - "----------------------------------------------------------------------" - ) + logger.info("----------------------------------------------------------------------") logger.info("Subsetting mesh file for region: %s", self.tag) today = datetime.today() @@ -336,12 +334,7 @@ def create_mesh_at_reg(self, mesh_dir, mesh_surf): 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 - + "_c" - + today_string - + ".nc", + os.path.splitext(mesh_surf)[0] + "_" + self.tag + "_c" + today_string + ".nc", ) logger.info("mesh_in : %s", mesh_in) @@ -349,14 +342,10 @@ def create_mesh_at_reg(self, mesh_dir, mesh_surf): self.mesh = mesh_out - node_coords, subset_element, subset_node, conn_dict = self.subset_mesh_at_reg( - 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, mesh_out - ) + self.write_mesh(f_in, node_coords, subset_element, subset_node, conn_dict, mesh_out) def subset_mesh_at_reg(self, mesh_in): """ @@ -421,9 +410,7 @@ def subset_mesh_at_reg(self, mesh_in): return node_coords, subset_element, subset_node, conn_dict @staticmethod - def write_mesh( - f_in, node_coords, subset_element, subset_node, conn_dict, mesh_out - ): + def write_mesh(f_in, node_coords, subset_element, subset_node, conn_dict, mesh_out): # pylint: disable=unused-argument """ This function writes out the subsetted mesh file. @@ -545,18 +532,8 @@ def write_shell_commands(self, namelist): """ # 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 LND_DOMAIN_MESH={}".format(str(self.mesh)), nl_file - ) - self.write_to_file( - "./xmlchange MASK_MESH={}".format(str(str(self.mesh))), 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 LND_DOMAIN_MESH={}".format(str(self.mesh)), nl_file) + self.write_to_file("./xmlchange MASK_MESH={}".format(str(str(self.mesh))), nl_file) diff --git a/python/ctsm/subset_data.py b/python/ctsm/subset_data.py index 52e27b649e..107b88df93 100644 --- a/python/ctsm/subset_data.py +++ b/python/ctsm/subset_data.py @@ -394,33 +394,39 @@ def setup_files(args, defaults, cesmroot): # if the crop flag is on - we need to use a different land use and surface data file num_pft = determine_num_pft(args.crop_flag) - fsurf_in = defaults.get("surfdat", "surfdat_"+num_pft+"pft") - fluse_in = defaults.get("landuse", "landuse_"+num_pft+"pft") - - file_dict = {'main_dir': defaults.get("main", "clmforcingindir"), - '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"), - defaults.get(datm_type, 'solardir'), - defaults.get(datm_type, 'precdir'), - defaults.get(datm_type, 'tpqwdir'), - defaults.get(datm_type, 'solartag'), - defaults.get(datm_type, 'prectag'), - defaults.get(datm_type, 'tpqwtag'), - defaults.get(datm_type, 'solarname'), - defaults.get(datm_type, 'precname'), - defaults.get(datm_type, 'tpqwname')) - } + fsurf_in = defaults.get("surfdat", "surfdat_" + num_pft + "pft") + fluse_in = defaults.get("landuse", "landuse_" + num_pft + "pft") + + file_dict = { + "main_dir": defaults.get("main", "clmforcingindir"), + "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"), + defaults.get(datm_type, "solardir"), + defaults.get(datm_type, "precdir"), + defaults.get(datm_type, "tpqwdir"), + defaults.get(datm_type, "solartag"), + defaults.get(datm_type, "prectag"), + defaults.get(datm_type, "tpqwtag"), + defaults.get(datm_type, "solarname"), + defaults.get(datm_type, "precname"), + defaults.get(datm_type, "tpqwname"), + ), + } return file_dict @@ -501,19 +507,19 @@ def subset_region(args, file_dict: dict): # -- Create Region Object region = RegionalCase( - lat1 = args.lat1, - lat2 = args.lat2, - lon1 = args.lon1, - lon2 = args.lon2, - reg_name = args.reg_name, - create_domain = args.create_domain, - create_surfdata = args.create_surfdata, - 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, + lat1=args.lat1, + lat2=args.lat2, + lon1=args.lon1, + lon2=args.lon2, + reg_name=args.reg_name, + create_domain=args.create_domain, + create_surfdata=args.create_surfdata, + 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, ) logger.debug(region) @@ -528,7 +534,7 @@ def subset_region(args, file_dict: dict): file_dict["fsurf_dir"], file_dict["fsurf_in"], args.user_mods_dir ) - #if region.create_mesh: + # if region.create_mesh: # region.create_mesh_at_reg (file_dict["mesh_dir"], file_dict["mesh_surf"]) # -- Create CTSM transient landuse data file @@ -537,7 +543,6 @@ def subset_region(args, file_dict: dict): 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: @@ -552,9 +557,13 @@ def subset_region(args, file_dict: dict): region.write_shell_commands(os.path.join(args.user_mods_dir, "shell_commands")) - print ('\nFor running this regional case with the created user_mods : ') - print ('./create_newcase --case case --res CLM_USRDAT --compset I2000Clm51BgcCrop', - '--run-unsupported --user-mods-dirs ', args.user_mods_dir, '\n\n') + print("\nFor running this regional case with the created user_mods : ") + print( + "./create_newcase --case case --res CLM_USRDAT --compset I2000Clm51BgcCrop", + "--run-unsupported --user-mods-dirs ", + args.user_mods_dir, + "\n\n", + ) logger.info("Successfully ran script for a regional case.") From 543e4243a7bb981c3c3b00bffd19228dd8fbd4c8 Mon Sep 17 00:00:00 2001 From: Negin Sobhani Date: Fri, 11 Nov 2022 18:32:16 -0700 Subject: [PATCH 042/595] fix bug with CLM_USRDAT_DIR --- python/ctsm/site_and_regional/regional_case.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/python/ctsm/site_and_regional/regional_case.py b/python/ctsm/site_and_regional/regional_case.py index dca050b217..5e4b08e746 100644 --- a/python/ctsm/site_and_regional/regional_case.py +++ b/python/ctsm/site_and_regional/regional_case.py @@ -226,8 +226,6 @@ def create_surfdata_at_reg(self, indir, file, user_mods_dir): logger.info("fsurf_in: %s", fsurf_in) logger.info("fsurf_out: %s", os.path.join(self.out_dir, fsurf_out)) - fsurf_out = os.path.join(self.out_dir, fsurf_out) - if self.create_mesh: mesh_out = os.path.join( self.out_dir, From 260ecb6e6104713b2660442a837309c07f2e8fb7 Mon Sep 17 00:00:00 2001 From: Erik Kluzek Date: Tue, 6 Dec 2022 16:41:44 -0700 Subject: [PATCH 043/595] Add dask, cartopy, and matplotlib to the list of packages for conda, note that the CGD conda env is failing when I try it on cheyenne --- python/conda_env_ctsm_py.txt | 3 +++ python/conda_env_ctsm_py_cgd.txt | 3 +++ python/conda_env_ctsm_py_latest.txt | 3 +++ 3 files changed, 9 insertions(+) diff --git a/python/conda_env_ctsm_py.txt b/python/conda_env_ctsm_py.txt index d757ae1782..726a3be1cd 100644 --- a/python/conda_env_ctsm_py.txt +++ b/python/conda_env_ctsm_py.txt @@ -15,6 +15,9 @@ pandas tqdm scipy netcdf4 +dask +cartopy +matplotlib requests numpy=1.18.5 xarray=0.16.2 diff --git a/python/conda_env_ctsm_py_cgd.txt b/python/conda_env_ctsm_py_cgd.txt index e7ee4af8ab..6b339d0df2 100644 --- a/python/conda_env_ctsm_py_cgd.txt +++ b/python/conda_env_ctsm_py_cgd.txt @@ -17,6 +17,9 @@ pandas tqdm scipy netcdf4 +dask +cartopy +matplotlib requests numpy=1.18.5 xarray=0.16.2 diff --git a/python/conda_env_ctsm_py_latest.txt b/python/conda_env_ctsm_py_latest.txt index 2dc2ed518d..e174e17825 100644 --- a/python/conda_env_ctsm_py_latest.txt +++ b/python/conda_env_ctsm_py_latest.txt @@ -4,6 +4,9 @@ pandas>=1.5.1 tqdm>=4.64.1 scipy netcdf4 +dask +cartopy +matplotlib requests numpy>=1.23.0 xarray>=2022.3.0 From 304c83e32043f4f8df4c444463e5ee15563145ec Mon Sep 17 00:00:00 2001 From: Erik Kluzek Date: Thu, 26 Jan 2023 15:12:44 -0700 Subject: [PATCH 044/595] Update needed to get unit tests to work due to update --- python/ctsm/test/testinputs/default_data.cfg | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/python/ctsm/test/testinputs/default_data.cfg b/python/ctsm/test/testinputs/default_data.cfg index 7e841dca54..19b77bc47d 100644 --- a/python/ctsm/test/testinputs/default_data.cfg +++ b/python/ctsm/test/testinputs/default_data.cfg @@ -18,7 +18,9 @@ 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 landuse_16pft = landuse.timeseries_0.9x1.25_hist_16pfts_Irrig_CMIP6_simyr1850-2015_c190214.nc From 8a510c22cc940a1d612181119012e960f96bbdb3 Mon Sep 17 00:00:00 2001 From: Erik Kluzek Date: Thu, 26 Jan 2023 15:14:46 -0700 Subject: [PATCH 045/595] Change needed for update --- tools/mksurfdata_map/default_data_1850.cfg | 2 ++ 1 file changed, 2 insertions(+) diff --git a/tools/mksurfdata_map/default_data_1850.cfg b/tools/mksurfdata_map/default_data_1850.cfg index 311aeef13d..a11da9aa61 100644 --- a/tools/mksurfdata_map/default_data_1850.cfg +++ b/tools/mksurfdata_map/default_data_1850.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_simyr1850_c190214.nc surfdat_78pft = surfdata_0.9x1.25_hist_78pfts_CMIP6_simyr1850_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 13d1c99c941e2297aef08f974da8739e5c64391b Mon Sep 17 00:00:00 2001 From: Erik Kluzek Date: Thu, 26 Jan 2023 16:30:09 -0700 Subject: [PATCH 046/595] Add a test that create-user-mods requires create_mesh for regional cases --- python/ctsm/test/test_unit_subset_data.py | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/python/ctsm/test/test_unit_subset_data.py b/python/ctsm/test/test_unit_subset_data.py index fd9aef631d..95862a8ec7 100755 --- a/python/ctsm/test/test_unit_subset_data.py +++ b/python/ctsm/test/test_unit_subset_data.py @@ -155,6 +155,14 @@ def test_inputdata_setup_files_bad_inputdata_arg(self): with self.assertRaisesRegex(SystemExit, "inputdata directory does not exist"): setup_files(self.args, self.defaults, self.cesmroot) + def test_create_user_mods_without_create_mesh(self): + """ + Test that you can't run create user mods without also doing create_mesh + """ + sys.argv = ["subset_data", "regional", "--create-user-mods"] + with self.assertRaisesRegex(argparse.ArgumentError, "For regional cases, you can not create user_mods"): + check_args(self.args) + if __name__ == "__main__": unit_testing.setup_for_tests() From fc262d454b868ace3065462abe1a70dc57203f8e Mon Sep 17 00:00:00 2001 From: Erik Kluzek Date: Thu, 26 Jan 2023 17:09:37 -0700 Subject: [PATCH 047/595] Get test for create_user_mods without create_mesh working correctly, and now passes --- python/ctsm/subset_data.py | 23 +++++++++++++---------- python/ctsm/test/test_unit_subset_data.py | 11 ++++++++--- 2 files changed, 21 insertions(+), 13 deletions(-) diff --git a/python/ctsm/subset_data.py b/python/ctsm/subset_data.py index be73bb34d3..1bf5eb3463 100644 --- a/python/ctsm/subset_data.py +++ b/python/ctsm/subset_data.py @@ -412,6 +412,19 @@ def check_args(args): ) raise argparse.ArgumentError(None, err_msg) + if args.run_type == "region" and args.create_user_mods: + if not args.create_mesh: + err_msg = textwrap.dedent( + """\ + \n ------------------------------------ + \nERROR: For regional cases, you can not create user_mods + \nwithout creating the mesh file. + + \nPlease rerun the script adding --create-mesh to subset the mesh file." + """ + ) + raise argparse.ArgumentError(None, err_msg) + def setup_user_mods(user_mods_dir, cesmroot): """ @@ -655,16 +668,6 @@ def subset_region(args, file_dict: dict): # -- 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")) print("\nFor running this regional case with the created user_mods : ") diff --git a/python/ctsm/test/test_unit_subset_data.py b/python/ctsm/test/test_unit_subset_data.py index 95862a8ec7..8b1dc3c27b 100755 --- a/python/ctsm/test/test_unit_subset_data.py +++ b/python/ctsm/test/test_unit_subset_data.py @@ -159,9 +159,14 @@ def test_create_user_mods_without_create_mesh(self): """ Test that you can't run create user mods without also doing create_mesh """ - sys.argv = ["subset_data", "regional", "--create-user-mods"] - with self.assertRaisesRegex(argparse.ArgumentError, "For regional cases, you can not create user_mods"): - check_args(self.args) + sys.argv = ["subset_data", "region", "--create-user-mods", "--create-surface"] + self.args = self.parser.parse_args() + print(self.args) + with self.assertRaisesRegex( + argparse.ArgumentError, "For regional cases, you can not create user_mods" + ): + print("run check_args") + check_args(self.args) if __name__ == "__main__": From 8918db2a744b6dd07275681dd4fcde6966d6ce45 Mon Sep 17 00:00:00 2001 From: Erik Kluzek Date: Fri, 27 Jan 2023 12:21:56 -0700 Subject: [PATCH 048/595] for region test turn mesh on and user-mods --- test/tools/nl_files/subset_data_region1 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/tools/nl_files/subset_data_region1 b/test/tools/nl_files/subset_data_region1 index fce83f0e2e..19cdaf822b 100644 --- a/test/tools/nl_files/subset_data_region1 +++ b/test/tools/nl_files/subset_data_region1 @@ -1 +1 @@ -region --lat1 -40 --lat2 15 --lon1 275 --lon2 330 --create-domain --create-surface --create-landuse --verbose --overwrite --reg test1 --inputdata-dir CSMDATA +region --lat1 -40 --lat2 15 --lon1 275 --lon2 330 --create-mesh --create-surface --create-user-mods --create-landuse --verbose --overwrite --reg test1 --inputdata-dir CSMDATA From 3ee9f05ffcee68e200ffbdc8abd1d00ac5196b6e Mon Sep 17 00:00:00 2001 From: Erik Kluzek Date: Sat, 28 Jan 2023 17:49:23 -0700 Subject: [PATCH 049/595] Add a single point test that is working for mesh_maker.py --- test/tools/README.testnames | 4 ++-- test/tools/input_tests_master | 3 +++ test/tools/nl_files/mesh_maker_1x1_numaIA | 1 + 3 files changed, 6 insertions(+), 2 deletions(-) create mode 100644 test/tools/nl_files/mesh_maker_1x1_numaIA diff --git a/test/tools/README.testnames b/test/tools/README.testnames index 11d9e23d4c..7f416ea5d0 100644 --- a/test/tools/README.testnames +++ b/test/tools/README.testnames @@ -20,7 +20,7 @@ n is the configuration type: 6 -- unused 7 -- unused 8 -- unused -9 -- unused +9 -- mesh_maker.py 0 -- run_neon a -- modify_data b -- subset_data @@ -28,7 +28,7 @@ c -- mkprocdata_map d -- mkmapgrids e -- unused f -- unused -g -- mksurfdata_map +g -- unused h -- unused i -- tools scripts diff --git a/test/tools/input_tests_master b/test/tools/input_tests_master index 63e4e5173f..b4956050e9 100644 --- a/test/tools/input_tests_master +++ b/test/tools/input_tests_master @@ -38,6 +38,9 @@ blbc1 TBLscript_tools.sh site_and_regional subset_data subset_data_f09_US_pt smbd1 TSMscript_tools.sh site_and_regional subset_data subset_data_region1 blbd1 TBLscript_tools.sh site_and_regional subset_data subset_data_region1 +sm9T1 TSMscript_tools.sh site_and_regional mesh_maker.py mesh_maker_1x1_numaIA +bl9T1 TBLscript_tools.sh site_and_regional mesh_maker.py mesh_maker_1x1_numaIA + smaa2 TSMscript_tools.sh site_and_regional modify_singlept_site_neon.py modify_data_YELL blaa2 TBLscript_tools.sh site_and_regional modify_singlept_site_neon.py modify_data_YELL diff --git a/test/tools/nl_files/mesh_maker_1x1_numaIA b/test/tools/nl_files/mesh_maker_1x1_numaIA new file mode 100644 index 0000000000..9cc10df947 --- /dev/null +++ b/test/tools/nl_files/mesh_maker_1x1_numaIA @@ -0,0 +1 @@ + --input CSMDATA/lnd/clm2/surfdata_map/ctsm5.1.dev116/surfdata_1x1_numaIA_hist_78pfts_CMIP6_simyr2000_c230123.nc --output 1x1_numaIA_mesh.nc --lat lsmlat --lon lsmlon --overwrite From 3a31fd2bbb69543c489e799d628075be56ac86be Mon Sep 17 00:00:00 2001 From: Erik Kluzek Date: Sat, 28 Jan 2023 22:31:02 -0700 Subject: [PATCH 050/595] Correct area --- tools/site_and_regional/mesh_maker.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tools/site_and_regional/mesh_maker.py b/tools/site_and_regional/mesh_maker.py index 4ce8642f96..7c1c9e1960 100755 --- a/tools/site_and_regional/mesh_maker.py +++ b/tools/site_and_regional/mesh_maker.py @@ -244,7 +244,7 @@ def main(): mask = ds[mask_name].values() if area_name is not None: - area = ds[mask_name].values() + area = ds[area_name].values() this_mesh = MeshType(lats, lons, mask=None) this_mesh.calculate_corners() From 26080973aa2d4a6c64ca5a9de1e53a0fb69e552b Mon Sep 17 00:00:00 2001 From: Erik Kluzek Date: Sun, 29 Jan 2023 17:17:22 -0700 Subject: [PATCH 051/595] Move the guts of mesh_maker to the python directory --- .../ctsm}/mesh_maker.py | 4 +- test/tools/README.testnames | 2 +- test/tools/input_tests_master | 4 +- tools/site_and_regional/README | 2 +- tools/site_and_regional/mesh_maker | 40 +++++++++++++++++++ 5 files changed, 46 insertions(+), 6 deletions(-) rename {tools/site_and_regional => python/ctsm}/mesh_maker.py (98%) mode change 100755 => 100644 create mode 100755 tools/site_and_regional/mesh_maker diff --git a/tools/site_and_regional/mesh_maker.py b/python/ctsm/mesh_maker.py old mode 100755 new mode 100644 similarity index 98% rename from tools/site_and_regional/mesh_maker.py rename to python/ctsm/mesh_maker.py index 7c1c9e1960..91019044fe --- a/tools/site_and_regional/mesh_maker.py +++ b/python/ctsm/mesh_maker.py @@ -8,7 +8,7 @@ For example for running WRF-CTSM cases, the user can create a mesh file for their domain : - ./mesh_maker.py --input wrfinput_d01 --output my_region + ./mesh_maker --input wrfinput_d01 --output my_region \ --lat XLAT --lon XLONG --verbose """ @@ -31,7 +31,7 @@ def get_parser(): """ - Get the parser object for mesh_maker.py script. + Get the parser object for mesh_maker script. Returns: parser (ArgumentParser): diff --git a/test/tools/README.testnames b/test/tools/README.testnames index 7f416ea5d0..364c884b33 100644 --- a/test/tools/README.testnames +++ b/test/tools/README.testnames @@ -20,7 +20,7 @@ n is the configuration type: 6 -- unused 7 -- unused 8 -- unused -9 -- mesh_maker.py +9 -- mesh_maker 0 -- run_neon a -- modify_data b -- subset_data diff --git a/test/tools/input_tests_master b/test/tools/input_tests_master index b4956050e9..c0c567ba29 100644 --- a/test/tools/input_tests_master +++ b/test/tools/input_tests_master @@ -38,8 +38,8 @@ blbc1 TBLscript_tools.sh site_and_regional subset_data subset_data_f09_US_pt smbd1 TSMscript_tools.sh site_and_regional subset_data subset_data_region1 blbd1 TBLscript_tools.sh site_and_regional subset_data subset_data_region1 -sm9T1 TSMscript_tools.sh site_and_regional mesh_maker.py mesh_maker_1x1_numaIA -bl9T1 TBLscript_tools.sh site_and_regional mesh_maker.py mesh_maker_1x1_numaIA +sm9T1 TSMscript_tools.sh site_and_regional mesh_maker mesh_maker_1x1_numaIA +bl9T1 TBLscript_tools.sh site_and_regional mesh_maker mesh_maker_1x1_numaIA smaa2 TSMscript_tools.sh site_and_regional modify_singlept_site_neon.py modify_data_YELL blaa2 TBLscript_tools.sh site_and_regional modify_singlept_site_neon.py modify_data_YELL diff --git a/tools/site_and_regional/README b/tools/site_and_regional/README index 8349a9724a..0ceb465ca5 100644 --- a/tools/site_and_regional/README +++ b/tools/site_and_regional/README @@ -20,7 +20,7 @@ subset_data For extracting domain files, surface dataset, and DATM files at a region, use: ./subset_data region -mesh_maker.py +mesh_maker This script creates a mesh file from a netcdf file with valid lats and lons. modify_singlept_site_neon.py diff --git a/tools/site_and_regional/mesh_maker b/tools/site_and_regional/mesh_maker new file mode 100755 index 0000000000..24b36de995 --- /dev/null +++ b/tools/site_and_regional/mesh_maker @@ -0,0 +1,40 @@ +#!/usr/bin/env python3 +""" +This is a just top-level skeleton script that calls +mesh_maker.py. +The original code (mesh_maker.py) is located under the +python/ctsm folder. + +For full instructions on how to run the code and different options, +please check the python/ctsm/mesh_maker.py file. + +---------------------------------------------------------------- +Instructions for running using conda python environments: + +../../py_env_create +conda activate ctsm_py +|------------------------------------------------------------------| +|--------------------- Instructions -----------------------------| +|------------------------------------------------------------------| +This script creates ESMF unstructured GRID (mesh file) from a netcdf +file with valid lats and lons. Provided lats and lons can be 1D or 2D. + +For example for running WRF-CTSM cases, the user can create a mesh +file for their domain : + ./mesh_maker --input wrfinput_d01 --output my_region + --lat XLAT --lon XLONG --verbose + +""" +import os +import sys + +# -- add python/ctsm to path +_CTSM_PYTHON = os.path.join( + os.path.dirname(os.path.realpath(__file__)), os.pardir, os.pardir, "python" +) +sys.path.insert(1, _CTSM_PYTHON) + +from ctsm.mesh_maker import main + +if __name__ == "__main__": + main() From c10a689042830ee5aef70bd3c0cda38faf36fda8 Mon Sep 17 00:00:00 2001 From: Erik Kluzek Date: Sun, 29 Jan 2023 17:23:40 -0700 Subject: [PATCH 052/595] Run through black --- python/ctsm/mesh_maker.py | 22 +++++----------------- 1 file changed, 5 insertions(+), 17 deletions(-) diff --git a/python/ctsm/mesh_maker.py b/python/ctsm/mesh_maker.py index 91019044fe..019b194e12 100644 --- a/python/ctsm/mesh_maker.py +++ b/python/ctsm/mesh_maker.py @@ -123,13 +123,9 @@ def main(): args = parser.parse_args() if args.verbose: - logging.basicConfig( - level=logging.DEBUG, format=" %(levelname)-8s :: %(message)s" - ) + logging.basicConfig(level=logging.DEBUG, format=" %(levelname)-8s :: %(message)s") else: - logging.basicConfig( - level=logging.INFO, format=" %(levelname)-8s :: %(message)s" - ) + logging.basicConfig(level=logging.INFO, format=" %(levelname)-8s :: %(message)s") nc_file = args.input lat_name = args.lat_name @@ -141,9 +137,7 @@ def main(): area_name = args.area_name if os.path.isfile(nc_file): - ds = xr.open_dataset( - nc_file, mask_and_scale=False, decode_times=False - ).transpose() + ds = xr.open_dataset(nc_file, mask_and_scale=False, decode_times=False).transpose() else: err_msg = textwrap.dedent( """\ @@ -195,9 +189,7 @@ def main(): today_string = today.strftime("%y%m%d") if mesh_out and out_dir: - logging.error( - " Both --outdir and --output cannot be provided at the same time." - ) + logging.error(" Both --outdir and --output cannot be provided at the same time.") err_msg = textwrap.dedent( """ \n ------------------------------------ @@ -220,11 +212,7 @@ def main(): mesh_out = os.path.join( out_dir, - os.path.splitext(nc_file)[0] - + "_ESMF_UNSTRUCTURED_MESH" - + "_c" - + today_string - + ".nc", + os.path.splitext(nc_file)[0] + "_ESMF_UNSTRUCTURED_MESH" + "_c" + today_string + ".nc", ) logging.info("Creating mesh file from : %s", nc_file) From 7f69e008fbd53f3dfdb922e95a094e1d376fac7a Mon Sep 17 00:00:00 2001 From: Erik Kluzek Date: Sun, 29 Jan 2023 18:23:10 -0700 Subject: [PATCH 053/595] Add a basic unit test for mesh_maker --- python/ctsm/mesh_maker.py | 5 +-- python/ctsm/test/test_unit_mesh_maker.py | 41 ++++++++++++++++++++++++ 2 files changed, 44 insertions(+), 2 deletions(-) create mode 100755 python/ctsm/test/test_unit_mesh_maker.py diff --git a/python/ctsm/mesh_maker.py b/python/ctsm/mesh_maker.py index 019b194e12..41d879330b 100644 --- a/python/ctsm/mesh_maker.py +++ b/python/ctsm/mesh_maker.py @@ -27,6 +27,7 @@ sys.path.insert(1, _CTSM_PYTHON) from ctsm.site_and_regional.mesh_type import MeshType +from ctsm.utils import abort def get_parser(): @@ -146,7 +147,7 @@ def main(): \n ------------------------------------ """ ) - raise parser.error(err_msg) + abort(err_msg) if lat_name not in ds.coords and lat_name not in ds.variables: logging.error("Input file does not have variable named %s", lat_name) @@ -199,7 +200,7 @@ def main(): \n --output : Absolute or relative path of the ESMF mesh file created.\n """ ) - raise parser.error(err_msg) + raise argparse.ArgumentError(None, err_msg) # -- no file name and output path: if not mesh_out and not out_dir: diff --git a/python/ctsm/test/test_unit_mesh_maker.py b/python/ctsm/test/test_unit_mesh_maker.py new file mode 100755 index 0000000000..0635496565 --- /dev/null +++ b/python/ctsm/test/test_unit_mesh_maker.py @@ -0,0 +1,41 @@ +#!/usr/bin/env python3 +""" +Unit tests for mesh_maker + +You can run this by: + python -m unittest test_unit_mesh_maker.py +""" + +import unittest +import argparse +import os +import sys + +# -- add python/ctsm to path (needed if we want to run the test stand-alone) +_CTSM_PYTHON = os.path.join(os.path.dirname(os.path.realpath(__file__)), os.pardir, os.pardir) +sys.path.insert(1, _CTSM_PYTHON) + +# pylint: disable=wrong-import-position +from ctsm import unit_testing +from ctsm.mesh_maker import get_parser, main + +# pylint: disable=invalid-name + + +class TestMeshMaker(unittest.TestCase): + """ + Basic class for testing mesh_maker.py. + """ + + def test_input_file_dne(self): + """ + Test that exits if input file does not exist + """ + sys.argv = ["mesh_maker", "--input", "zztop.nc", "--lat", "lsmlat", "--lon", "lsmlon"] + with self.assertRaisesRegex(SystemExit, ".*Input file not found..*"): + main() + + +if __name__ == "__main__": + unit_testing.setup_for_tests() + unittest.main() From 050a87ba3c34db73c0d3e57ccc3a262a0dffb477 Mon Sep 17 00:00:00 2001 From: Erik Kluzek Date: Sun, 29 Jan 2023 20:07:29 -0700 Subject: [PATCH 054/595] Add a failing test for outdir and outfile at the same time --- python/ctsm/test/test_unit_mesh_maker.py | 11 ++++++++++- 1 file changed, 10 insertions(+), 1 deletion(-) diff --git a/python/ctsm/test/test_unit_mesh_maker.py b/python/ctsm/test/test_unit_mesh_maker.py index 0635496565..53b1eb5066 100755 --- a/python/ctsm/test/test_unit_mesh_maker.py +++ b/python/ctsm/test/test_unit_mesh_maker.py @@ -32,7 +32,16 @@ def test_input_file_dne(self): Test that exits if input file does not exist """ sys.argv = ["mesh_maker", "--input", "zztop.nc", "--lat", "lsmlat", "--lon", "lsmlon"] - with self.assertRaisesRegex(SystemExit, ".*Input file not found..*"): + with self.assertRaisesRegex(SystemExit, "Input file not found."): + main() + + def test_outfile_and_outdir(self): + """ + Test that exits if both outfile and outdir are provided + """ + infile = "ctsm/test/testinputs/default_data.cfg" + sys.argv = ["mesh_maker", "--input", infile, "--lat", "lsmlat", "--lon", "lsmlon", "--outdir", ".", "--output", "outthing.nc"] + with self.assertRaisesRegex(SystemExit, "You have provided both --outdir and --output."): main() From e5e7c66d8517f8b4b5636f931dc20c8298aae83a Mon Sep 17 00:00:00 2001 From: Erik Kluzek Date: Sun, 29 Jan 2023 20:28:15 -0700 Subject: [PATCH 055/595] Get test to work, but moving it up earlier --- python/ctsm/mesh_maker.py | 26 ++++++++++++------------ python/ctsm/test/test_unit_mesh_maker.py | 2 +- 2 files changed, 14 insertions(+), 14 deletions(-) diff --git a/python/ctsm/mesh_maker.py b/python/ctsm/mesh_maker.py index 41d879330b..08208ab972 100644 --- a/python/ctsm/mesh_maker.py +++ b/python/ctsm/mesh_maker.py @@ -137,6 +137,19 @@ def main(): mask_name = args.mask_name area_name = args.area_name + if mesh_out and out_dir: + logging.error(" Both --outdir and --output cannot be provided at the same time.") + err_msg = textwrap.dedent( + """ + \n ------------------------------------ + \n You have provided both --outdir and --output. + \n Please provide only one of these options to proceed: + \n --outdir : directory to save mesh file. mesh file name automatically created. + \n --output : Absolute or relative path of the ESMF mesh file created.\n + """ + ) + raise argparse.ArgumentError(None, err_msg) + if os.path.isfile(nc_file): ds = xr.open_dataset(nc_file, mask_and_scale=False, decode_times=False).transpose() else: @@ -189,19 +202,6 @@ def main(): today = datetime.today() today_string = today.strftime("%y%m%d") - if mesh_out and out_dir: - logging.error(" Both --outdir and --output cannot be provided at the same time.") - err_msg = textwrap.dedent( - """ - \n ------------------------------------ - \n You have provided both --outdir and --output. - \n Please provide only one of these options to proceed: - \n --outdir : directory to save mesh file. mesh file name automatically created. - \n --output : Absolute or relative path of the ESMF mesh file created.\n - """ - ) - raise argparse.ArgumentError(None, err_msg) - # -- no file name and output path: if not mesh_out and not out_dir: out_dir = os.path.join(os.getcwd(), "meshes") diff --git a/python/ctsm/test/test_unit_mesh_maker.py b/python/ctsm/test/test_unit_mesh_maker.py index 53b1eb5066..d52835e31f 100755 --- a/python/ctsm/test/test_unit_mesh_maker.py +++ b/python/ctsm/test/test_unit_mesh_maker.py @@ -41,7 +41,7 @@ def test_outfile_and_outdir(self): """ infile = "ctsm/test/testinputs/default_data.cfg" sys.argv = ["mesh_maker", "--input", infile, "--lat", "lsmlat", "--lon", "lsmlon", "--outdir", ".", "--output", "outthing.nc"] - with self.assertRaisesRegex(SystemExit, "You have provided both --outdir and --output."): + with self.assertRaisesRegex(argparse.ArgumentError, "You have provided both --outdir and --output."): main() From 3b41bdb7743acdb811070b16be0ce4cdead1d1ea Mon Sep 17 00:00:00 2001 From: Erik Kluzek Date: Sun, 29 Jan 2023 20:32:34 -0700 Subject: [PATCH 056/595] Add test for overwrite that fails --- python/ctsm/test/test_unit_mesh_maker.py | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/python/ctsm/test/test_unit_mesh_maker.py b/python/ctsm/test/test_unit_mesh_maker.py index d52835e31f..535232ce46 100755 --- a/python/ctsm/test/test_unit_mesh_maker.py +++ b/python/ctsm/test/test_unit_mesh_maker.py @@ -44,6 +44,15 @@ def test_outfile_and_outdir(self): with self.assertRaisesRegex(argparse.ArgumentError, "You have provided both --outdir and --output."): main() + def test_outfile_exists_and_no_overwrite(self): + """ + Test that exits if outfile already exists and overwrite was not used + """ + infile = "ctsm/test/testinputs/default_data.cfg" + sys.argv = ["mesh_maker", "--input", infile, "--lat", "lsmlat", "--lon", "lsmlon", "--outdir", ".", "--output", infile] + with self.assertRaisesRegex(SystemExit, "output meshfile exists, please choose --overwrite to overwrite the mesh file."): + main() + if __name__ == "__main__": unit_testing.setup_for_tests() From 82df5850c50757061e6cd29d73e5e1a60dbfa9f3 Mon Sep 17 00:00:00 2001 From: Erik Kluzek Date: Mon, 30 Jan 2023 00:31:13 -0700 Subject: [PATCH 057/595] Fix the overwrite test --- python/ctsm/mesh_maker.py | 49 ++++++++++++------------ python/ctsm/test/test_unit_mesh_maker.py | 6 ++- 2 files changed, 29 insertions(+), 26 deletions(-) diff --git a/python/ctsm/mesh_maker.py b/python/ctsm/mesh_maker.py index 08208ab972..112120f580 100644 --- a/python/ctsm/mesh_maker.py +++ b/python/ctsm/mesh_maker.py @@ -150,6 +150,31 @@ def main(): ) raise argparse.ArgumentError(None, err_msg) + # -- no file name and output path: + if not mesh_out and not out_dir: + out_dir = os.path.join(os.getcwd(), "meshes") + + if not mesh_out: + # -- make output path if does not exist. + if not os.path.isdir(out_dir): + os.mkdir(out_dir) + + today = datetime.today() + today_string = today.strftime("%y%m%d") + mesh_out = os.path.join( + out_dir, + os.path.splitext(nc_file)[0] + "_ESMF_UNSTRUCTURED_MESH" + "_c" + today_string + ".nc", + ) + + # -- exit if mesh_out exists and --overwrite is not specified. + print( mesh_out ) + if os.path.exists(mesh_out): + if overwrite: + os.remove(mesh_out) + else: + err_msg = "output meshfile exists, please choose --overwrite to overwrite the mesh file." + abort(err_msg) + if os.path.isfile(nc_file): ds = xr.open_dataset(nc_file, mask_and_scale=False, decode_times=False).transpose() else: @@ -199,35 +224,11 @@ def main(): ) sys.exit() - today = datetime.today() - today_string = today.strftime("%y%m%d") - # -- no file name and output path: - if not mesh_out and not out_dir: - out_dir = os.path.join(os.getcwd(), "meshes") - - if not mesh_out: - # -- make output path if does not exist. - if not os.path.isdir(out_dir): - os.mkdir(out_dir) - - mesh_out = os.path.join( - out_dir, - os.path.splitext(nc_file)[0] + "_ESMF_UNSTRUCTURED_MESH" + "_c" + today_string + ".nc", - ) logging.info("Creating mesh file from : %s", nc_file) logging.info("Writing mesh file to : %s", mesh_out) - # -- exit if mesh_out exists and --overwrite is not specified. - if os.path.exists(mesh_out): - if overwrite: - os.remove(mesh_out) - else: - logging.error( - "output meshfile exists, please choose --overwrite to overwrite the mesh file." - ) - sys.exit() if mask_name is not None: mask = ds[mask_name].values() diff --git a/python/ctsm/test/test_unit_mesh_maker.py b/python/ctsm/test/test_unit_mesh_maker.py index 535232ce46..b26f19b66c 100755 --- a/python/ctsm/test/test_unit_mesh_maker.py +++ b/python/ctsm/test/test_unit_mesh_maker.py @@ -31,7 +31,7 @@ def test_input_file_dne(self): """ Test that exits if input file does not exist """ - sys.argv = ["mesh_maker", "--input", "zztop.nc", "--lat", "lsmlat", "--lon", "lsmlon"] + sys.argv = ["mesh_maker", "--input", "zztop.nc", "--lat", "lsmlat", "--lon", "lsmlon", "--outdir", "."] with self.assertRaisesRegex(SystemExit, "Input file not found."): main() @@ -41,7 +41,8 @@ def test_outfile_and_outdir(self): """ infile = "ctsm/test/testinputs/default_data.cfg" sys.argv = ["mesh_maker", "--input", infile, "--lat", "lsmlat", "--lon", "lsmlon", "--outdir", ".", "--output", "outthing.nc"] - with self.assertRaisesRegex(argparse.ArgumentError, "You have provided both --outdir and --output."): + with self.assertRaisesRegex(argparse.ArgumentError, ".*You have provided both --outdir and --output..*"): + print( sys.argv ) main() def test_outfile_exists_and_no_overwrite(self): @@ -51,6 +52,7 @@ def test_outfile_exists_and_no_overwrite(self): infile = "ctsm/test/testinputs/default_data.cfg" sys.argv = ["mesh_maker", "--input", infile, "--lat", "lsmlat", "--lon", "lsmlon", "--outdir", ".", "--output", infile] with self.assertRaisesRegex(SystemExit, "output meshfile exists, please choose --overwrite to overwrite the mesh file."): + print( sys.argv ) main() From 2023c1421316bbc6f0a5752e773676e9675e254f Mon Sep 17 00:00:00 2001 From: Erik Kluzek Date: Mon, 30 Jan 2023 00:38:35 -0700 Subject: [PATCH 058/595] Run through black and get unit tests working again --- python/ctsm/mesh_maker.py | 10 ++--- python/ctsm/test/test_unit_mesh_maker.py | 47 ++++++++++++++++++++---- 2 files changed, 44 insertions(+), 13 deletions(-) diff --git a/python/ctsm/mesh_maker.py b/python/ctsm/mesh_maker.py index 112120f580..ecbaf73eaf 100644 --- a/python/ctsm/mesh_maker.py +++ b/python/ctsm/mesh_maker.py @@ -148,7 +148,7 @@ def main(): \n --output : Absolute or relative path of the ESMF mesh file created.\n """ ) - raise argparse.ArgumentError(None, err_msg) + abort(err_msg) # -- no file name and output path: if not mesh_out and not out_dir: @@ -167,12 +167,13 @@ def main(): ) # -- exit if mesh_out exists and --overwrite is not specified. - print( mesh_out ) if os.path.exists(mesh_out): if overwrite: os.remove(mesh_out) else: - err_msg = "output meshfile exists, please choose --overwrite to overwrite the mesh file." + err_msg = ( + "output meshfile exists, please choose --overwrite to overwrite the mesh file." + ) abort(err_msg) if os.path.isfile(nc_file): @@ -224,12 +225,9 @@ def main(): ) sys.exit() - - logging.info("Creating mesh file from : %s", nc_file) logging.info("Writing mesh file to : %s", mesh_out) - if mask_name is not None: mask = ds[mask_name].values() diff --git a/python/ctsm/test/test_unit_mesh_maker.py b/python/ctsm/test/test_unit_mesh_maker.py index b26f19b66c..54f5f689e8 100755 --- a/python/ctsm/test/test_unit_mesh_maker.py +++ b/python/ctsm/test/test_unit_mesh_maker.py @@ -31,7 +31,17 @@ def test_input_file_dne(self): """ Test that exits if input file does not exist """ - sys.argv = ["mesh_maker", "--input", "zztop.nc", "--lat", "lsmlat", "--lon", "lsmlon", "--outdir", "."] + sys.argv = [ + "mesh_maker", + "--input", + "zztop.nc", + "--lat", + "lsmlat", + "--lon", + "lsmlon", + "--outdir", + ".", + ] with self.assertRaisesRegex(SystemExit, "Input file not found."): main() @@ -40,9 +50,20 @@ def test_outfile_and_outdir(self): Test that exits if both outfile and outdir are provided """ infile = "ctsm/test/testinputs/default_data.cfg" - sys.argv = ["mesh_maker", "--input", infile, "--lat", "lsmlat", "--lon", "lsmlon", "--outdir", ".", "--output", "outthing.nc"] - with self.assertRaisesRegex(argparse.ArgumentError, ".*You have provided both --outdir and --output..*"): - print( sys.argv ) + sys.argv = [ + "mesh_maker", + "--input", + infile, + "--lat", + "lsmlat", + "--lon", + "lsmlon", + "--outdir", + ".", + "--output", + "outthing.nc", + ] + with self.assertRaisesRegex(SystemExit, "You have provided both --outdir and --output."): main() def test_outfile_exists_and_no_overwrite(self): @@ -50,9 +71,21 @@ def test_outfile_exists_and_no_overwrite(self): Test that exits if outfile already exists and overwrite was not used """ infile = "ctsm/test/testinputs/default_data.cfg" - sys.argv = ["mesh_maker", "--input", infile, "--lat", "lsmlat", "--lon", "lsmlon", "--outdir", ".", "--output", infile] - with self.assertRaisesRegex(SystemExit, "output meshfile exists, please choose --overwrite to overwrite the mesh file."): - print( sys.argv ) + sys.argv = [ + "mesh_maker", + "--input", + infile, + "--lat", + "lsmlat", + "--lon", + "lsmlon", + "--output", + infile, + ] + with self.assertRaisesRegex( + SystemExit, + "output meshfile exists, please choose --overwrite to overwrite the mesh file.", + ): main() From 0647e760eb99a44564499aa9256a3d0b43edd38d Mon Sep 17 00:00:00 2001 From: Erik Kluzek Date: Mon, 30 Jan 2023 12:44:20 -0700 Subject: [PATCH 059/595] Some important changes moved over from subset_data that allow mesh_maker to work for regional surface datasets created by subset_data, before this didn't work --- python/ctsm/mesh_maker.py | 15 ++++++++++----- 1 file changed, 10 insertions(+), 5 deletions(-) diff --git a/python/ctsm/mesh_maker.py b/python/ctsm/mesh_maker.py index ecbaf73eaf..50927cf48d 100644 --- a/python/ctsm/mesh_maker.py +++ b/python/ctsm/mesh_maker.py @@ -19,6 +19,7 @@ import textwrap from datetime import datetime import xarray as xr +import numpy as np # -- add python/ctsm to path _CTSM_PYTHON = os.path.join( @@ -209,8 +210,8 @@ def main(): len(ds[lon_name].dims).__str__(), ) - lats = ds[lat_name] - lons = ds[lon_name] + lats = ds[lat_name].astype(np.float32) + lons = ds[lon_name].astype(np.float32) if (len(lats.dims) > 2) or (len(lons.dims) > 2): time_dims = [dim for dim in lats.dims if "time" in dim.lower()] @@ -230,13 +231,17 @@ def main(): if mask_name is not None: mask = ds[mask_name].values() - + else: + mask = None + if area_name is not None: area = ds[area_name].values() + else: + area = None - this_mesh = MeshType(lats, lons, mask=None) + this_mesh = MeshType(lats, lons, mask=mask) this_mesh.calculate_corners() - this_mesh.create_esmf(mesh_out, area=None) + this_mesh.create_esmf(mesh_out, area=area) if __name__ == "__main__": From ad8ce53bfd3cbde59178531f929c58fa9fde1f5c Mon Sep 17 00:00:00 2001 From: Erik Kluzek Date: Wed, 1 Feb 2023 16:32:43 -0700 Subject: [PATCH 060/595] Update subset_data regions and tests --- test/tools/README.testnames | 6 ++++-- test/tools/input_tests_master | 14 ++++++++++---- test/tools/nl_files/subset_data_KONA | 2 +- test/tools/nl_files/subset_data_YELL | 2 +- .../nl_files/subset_data_f09_1x1pt_townshipSD | 1 + .../nl_files/subset_data_f09_58x45pt_SouthAmerica | 1 + test/tools/nl_files/subset_data_f09_US_pt | 1 - test/tools/nl_files/subset_data_region1 | 1 - 8 files changed, 18 insertions(+), 10 deletions(-) create mode 100644 test/tools/nl_files/subset_data_f09_1x1pt_townshipSD create mode 100644 test/tools/nl_files/subset_data_f09_58x45pt_SouthAmerica delete mode 100644 test/tools/nl_files/subset_data_f09_US_pt delete mode 100644 test/tools/nl_files/subset_data_region1 diff --git a/test/tools/README.testnames b/test/tools/README.testnames index 364c884b33..ad560cb852 100644 --- a/test/tools/README.testnames +++ b/test/tools/README.testnames @@ -43,9 +43,11 @@ m is the resolution 9 -- 4x5 a -- NEON YELL b -- NEON KONA -d -- region1 c -- single point from the 0.9x1.25 grid -g -- unused +d -- SouthAmerica +e -- 1850PanTropics +f -- PanBoreal +g -- AlaskaTananaValley y -- 1.9x2.5 with transient 1850-2100 for rcp=2.6 and glacier-MEC on T -- 1x1_numaIA Z -- 10x15 with crop on diff --git a/test/tools/input_tests_master b/test/tools/input_tests_master index c0c567ba29..494c0fb559 100644 --- a/test/tools/input_tests_master +++ b/test/tools/input_tests_master @@ -33,10 +33,16 @@ smbb1 TSMscript_tools.sh site_and_regional subset_data subset_data_KONA blbb1 TBLscript_tools.sh site_and_regional subset_data subset_data_KONA smb81 TSMscript_tools.sh site_and_regional subset_data subset_data_US-UMB blb81 TBLscript_tools.sh site_and_regional subset_data subset_data_US-UMB -smbc1 TSMscript_tools.sh site_and_regional subset_data subset_data_f09_US_pt -blbc1 TBLscript_tools.sh site_and_regional subset_data subset_data_f09_US_pt -smbd1 TSMscript_tools.sh site_and_regional subset_data subset_data_region1 -blbd1 TBLscript_tools.sh site_and_regional subset_data subset_data_region1 +smbc1 TSMscript_tools.sh site_and_regional subset_data subset_data_f09_1x1pt_townshipSD +blbc1 TBLscript_tools.sh site_and_regional subset_data subset_data_f09_1x1pt_townshipSD +smbd1 TSMscript_tools.sh site_and_regional subset_data subset_data_f09_58x45pt_SouthAmerica +blbd1 TBLscript_tools.sh site_and_regional subset_data subset_data_f09_58x45pt_SouthAmerica +smbe1 TSMscript_tools.sh site_and_regional subset_data subset_data_1850PanTropics +blbe1 TBLscript_tools.sh site_and_regional subset_data subset_data_1850PanTropics +smbf1 TSMscript_tools.sh site_and_regional subset_data subset_data_f09_38x288pt_PanBoreal +blbf1 TBLscript_tools.sh site_and_regional subset_data subset_data_f09_38x288pt_PanBoreal +smbg1 TSMscript_tools.sh site_and_regional subset_data subset_data_f09_4x9pt_AlaskaTananaValley +blbg1 TBLscript_tools.sh site_and_regional subset_data subset_data_f09_4x9pt_AlaskaTananaValley sm9T1 TSMscript_tools.sh site_and_regional mesh_maker mesh_maker_1x1_numaIA bl9T1 TBLscript_tools.sh site_and_regional mesh_maker mesh_maker_1x1_numaIA diff --git a/test/tools/nl_files/subset_data_KONA b/test/tools/nl_files/subset_data_KONA index c3be007869..0df59b1b17 100644 --- a/test/tools/nl_files/subset_data_KONA +++ b/test/tools/nl_files/subset_data_KONA @@ -1 +1 @@ -point --lon 263.38956 --lat 39.1082 --site KONA --dompft 17 19 23 45 --pctpft 28 12 32 28 --crop --create-domain --create-surface --outdir EXEDIR/KONA_user-mod_and_data --user-mods-dir EXEDIR/KONA_user-mod_and_data --verbose --inputdata-dir CSMDATA +point --lon 263.38956 --lat 39.1082 --site KONA --dompft 17 19 23 45 --pctpft 28 12 32 28 --crop --create-surface --outdir EXEDIR/KONA_user-mod_and_data --user-mods-dir EXEDIR/KONA_user-mod_and_data --verbose --inputdata-dir CSMDATA diff --git a/test/tools/nl_files/subset_data_YELL b/test/tools/nl_files/subset_data_YELL index 8295830c25..0d6960e7f5 100644 --- a/test/tools/nl_files/subset_data_YELL +++ b/test/tools/nl_files/subset_data_YELL @@ -1 +1 @@ -point --lon 250.45804 --lat 44.95597 --site YELL --dompft 1 --crop --create-domain --create-surface --outdir EXEDIR/YELL_user-mod_and_data --user-mods-dir EXEDIR/YELL_user-mod_and_data --verbose --inputdata-dir CSMDATA +point --lon 250.45804 --lat 44.95597 --site YELL --dompft 1 --crop --uniform-snowpack --cap-saturation --create-surface --outdir EXEDIR/YELL_user-mod_and_data --user-mods-dir EXEDIR/YELL_user-mod_and_data --silent --inputdata-dir CSMDATA diff --git a/test/tools/nl_files/subset_data_f09_1x1pt_townshipSD b/test/tools/nl_files/subset_data_f09_1x1pt_townshipSD new file mode 100644 index 0000000000..aa25c07d1e --- /dev/null +++ b/test/tools/nl_files/subset_data_f09_1x1pt_townshipSD @@ -0,0 +1 @@ +point --lon 257.5 --lat 43.822 --site f09_1x1pt_townshipSD --include-nonveg --crop --create-datm --create-user-mods --datm-syr 2000 --datm-eyr 2000 --create-surface --outdir EXEDIR/f09_US_pt_user-mod_and_data --user-mods-dir EXEDIR/f09_US_pt_user-mod_and_data --verbose --inputdata-dir CSMDATA diff --git a/test/tools/nl_files/subset_data_f09_58x45pt_SouthAmerica b/test/tools/nl_files/subset_data_f09_58x45pt_SouthAmerica new file mode 100644 index 0000000000..201dd2c76c --- /dev/null +++ b/test/tools/nl_files/subset_data_f09_58x45pt_SouthAmerica @@ -0,0 +1 @@ +region --lat1 -40 --lat2 15 --lon1 275 --lon2 330 --create-mesh --create-surface --create-user-mods --create-domain --create-landuse --verbose --overwrite --reg f09_58x45_SouthAmerica --inputdata-dir CSMDATA diff --git a/test/tools/nl_files/subset_data_f09_US_pt b/test/tools/nl_files/subset_data_f09_US_pt deleted file mode 100644 index bf6d5e2861..0000000000 --- a/test/tools/nl_files/subset_data_f09_US_pt +++ /dev/null @@ -1 +0,0 @@ -point --lon 257.5 --lat 43.822 --site 1x1_ --include-nonveg --crop --create-landuse --create-datm --create-user-mods --datm-syr 2000 --datm-eyr 2000 --create-surface --outdir EXEDIR/f09_US_pt_user-mod_and_data --user-mods-dir EXEDIR/f09_US_pt_user-mod_and_data --verbose --inputdata-dir CSMDATA diff --git a/test/tools/nl_files/subset_data_region1 b/test/tools/nl_files/subset_data_region1 deleted file mode 100644 index 19cdaf822b..0000000000 --- a/test/tools/nl_files/subset_data_region1 +++ /dev/null @@ -1 +0,0 @@ -region --lat1 -40 --lat2 15 --lon1 275 --lon2 330 --create-mesh --create-surface --create-user-mods --create-landuse --verbose --overwrite --reg test1 --inputdata-dir CSMDATA From 6592810f1027f302745eb9cd0cb06153d5cb4b1f Mon Sep 17 00:00:00 2001 From: Erik Kluzek Date: Wed, 1 Feb 2023 16:49:50 -0700 Subject: [PATCH 061/595] Change name of one of the regional tests --- test/tools/input_tests_master | 4 ++-- test/tools/nl_files/subset_data_f09_38x288pt_PanBoreal | 1 + test/tools/nl_files/subset_data_f09_4x9pt_AlaskaTananaValley | 1 + test/tools/nl_files/subset_data_f09_90x288pt_1850PanTropics | 1 + 4 files changed, 5 insertions(+), 2 deletions(-) create mode 100644 test/tools/nl_files/subset_data_f09_38x288pt_PanBoreal create mode 100644 test/tools/nl_files/subset_data_f09_4x9pt_AlaskaTananaValley create mode 100644 test/tools/nl_files/subset_data_f09_90x288pt_1850PanTropics diff --git a/test/tools/input_tests_master b/test/tools/input_tests_master index 494c0fb559..61f8e4b1ad 100644 --- a/test/tools/input_tests_master +++ b/test/tools/input_tests_master @@ -37,8 +37,8 @@ smbc1 TSMscript_tools.sh site_and_regional subset_data subset_data_f09_1x1pt_tow blbc1 TBLscript_tools.sh site_and_regional subset_data subset_data_f09_1x1pt_townshipSD smbd1 TSMscript_tools.sh site_and_regional subset_data subset_data_f09_58x45pt_SouthAmerica blbd1 TBLscript_tools.sh site_and_regional subset_data subset_data_f09_58x45pt_SouthAmerica -smbe1 TSMscript_tools.sh site_and_regional subset_data subset_data_1850PanTropics -blbe1 TBLscript_tools.sh site_and_regional subset_data subset_data_1850PanTropics +smbe1 TSMscript_tools.sh site_and_regional subset_data subset_data_f09_90x288pt_1850PanTropics +blbe1 TBLscript_tools.sh site_and_regional subset_data subset_data_f09_90x288pt_1850PanTropics smbf1 TSMscript_tools.sh site_and_regional subset_data subset_data_f09_38x288pt_PanBoreal blbf1 TBLscript_tools.sh site_and_regional subset_data subset_data_f09_38x288pt_PanBoreal smbg1 TSMscript_tools.sh site_and_regional subset_data subset_data_f09_4x9pt_AlaskaTananaValley diff --git a/test/tools/nl_files/subset_data_f09_38x288pt_PanBoreal b/test/tools/nl_files/subset_data_f09_38x288pt_PanBoreal new file mode 100644 index 0000000000..9ddf17438a --- /dev/null +++ b/test/tools/nl_files/subset_data_f09_38x288pt_PanBoreal @@ -0,0 +1 @@ +region --lat1 55 --lat2 90 --lon1 0 --lon2 360 --create-mesh --create-surface --create-domain --create-user-mods --verbose --overwrite --reg f09_38x288pt_PanBoreal --inputdata-dir CSMDATA diff --git a/test/tools/nl_files/subset_data_f09_4x9pt_AlaskaTananaValley b/test/tools/nl_files/subset_data_f09_4x9pt_AlaskaTananaValley new file mode 100644 index 0000000000..29997d1327 --- /dev/null +++ b/test/tools/nl_files/subset_data_f09_4x9pt_AlaskaTananaValley @@ -0,0 +1 @@ +region --lat1 62 --lat2 66 --lon1 -152 --lon2 -141 --create-mesh --create-surface --create-user-mods --verbose --overwrite --reg f09_4x9pt_AlaskaTananaValley --inputdata-dir CSMDATA diff --git a/test/tools/nl_files/subset_data_f09_90x288pt_1850PanTropics b/test/tools/nl_files/subset_data_f09_90x288pt_1850PanTropics new file mode 100644 index 0000000000..8a3dd76937 --- /dev/null +++ b/test/tools/nl_files/subset_data_f09_90x288pt_1850PanTropics @@ -0,0 +1 @@ +region --lat1 -55 --lat2 30 --lon1 0 --lon2 360 --crop --create-surface --overwrite --reg f09_90x288pt_1850PanTropics --inputdata-dir CSMDATA --cfg-file CTSM_ROOT/tools/mksurfdata_map/default_data_1850.cfg From 4ff8e7dcc4bee5807b128d6f4f740cb02946a9b4 Mon Sep 17 00:00:00 2001 From: Erik Kluzek Date: Wed, 1 Feb 2023 23:23:02 -0700 Subject: [PATCH 062/595] Fix some mesh_maker pylint issues --- python/ctsm/mesh_maker.py | 20 ++++++++++++-------- 1 file changed, 12 insertions(+), 8 deletions(-) diff --git a/python/ctsm/mesh_maker.py b/python/ctsm/mesh_maker.py index 50927cf48d..aba1f6c70f 100644 --- a/python/ctsm/mesh_maker.py +++ b/python/ctsm/mesh_maker.py @@ -86,7 +86,8 @@ def get_parser(): ) parser.add_argument( "--mask", - help="Name of mask varibale on netcdf input file. If none given, create a fake mask with values of 1.", + help="Name of mask varibale on netcdf input file." + + " If none given, create a fake mask with values of 1.", action="store", dest="mask_name", type=str, @@ -94,7 +95,8 @@ def get_parser(): ) parser.add_argument( "--area", - help="Name of area variable on netcdf input file. If none given, ESMF calculates element areas automatically. ", + help="Name of area variable on netcdf input file." + + " If none given, ESMF calculates element areas automatically. ", action="store", dest="area_name", type=str, @@ -183,7 +185,8 @@ def main(): err_msg = textwrap.dedent( """\ \n ------------------------------------ - \n Input file not found. Please make sure to provide the full path of Netcdf input file for making mesh. + \n Input file not found. Please make sure to provide the + \n full path of Netcdf input file for making the mesh. \n ------------------------------------ """ ) @@ -216,15 +219,16 @@ def main(): if (len(lats.dims) > 2) or (len(lons.dims) > 2): time_dims = [dim for dim in lats.dims if "time" in dim.lower()] if time_dims: - logging.debug("- time dimension found on lat {}".format(time_dims)) + logging.debug("- time dimension found on lat %s", str(time_dims)) logging.debug("- removing time dimensions from lats and lons. ") lats = lats[:, :, 0] lons = lons[:, :, 0] else: - logging.error( - "latitude or longitude has more than 2 dimensions and the third dimension cannot be detected as time." + err_msg = ( + "latitude or longitude has more than 2 dimensions and " + + "the third dimension cannot be detected as time." ) - sys.exit() + abort(err_msg) logging.info("Creating mesh file from : %s", nc_file) logging.info("Writing mesh file to : %s", mesh_out) @@ -233,7 +237,7 @@ def main(): mask = ds[mask_name].values() else: mask = None - + if area_name is not None: area = ds[area_name].values() else: From 98a69944394531b9ec1d13cc7f3b152b6106d22f Mon Sep 17 00:00:00 2001 From: Erik Kluzek Date: Wed, 1 Feb 2023 23:45:11 -0700 Subject: [PATCH 063/595] Use ctsm_logging and add a non working option to turn off plots --- python/ctsm/mesh_maker.py | 35 +++++++++++++++++------------------ 1 file changed, 17 insertions(+), 18 deletions(-) diff --git a/python/ctsm/mesh_maker.py b/python/ctsm/mesh_maker.py index aba1f6c70f..f5ea4dcd31 100644 --- a/python/ctsm/mesh_maker.py +++ b/python/ctsm/mesh_maker.py @@ -21,14 +21,13 @@ import xarray as xr import numpy as np -# -- add python/ctsm to path -_CTSM_PYTHON = os.path.join( - os.path.dirname(os.path.realpath(__file__)), os.pardir, os.pardir, "python" -) -sys.path.insert(1, _CTSM_PYTHON) - from ctsm.site_and_regional.mesh_type import MeshType from ctsm.utils import abort +from ctsm.ctsm_logging import ( + setup_logging_pre_config, + add_logging_args, + process_logging_args, +) def get_parser(): @@ -93,6 +92,13 @@ def get_parser(): type=str, required=False, ) + parser.add_argument( + "--no-plot", + help="Do not do the plots of the mesh", + action="store_true", + dest="noplot", + required=False, + ) parser.add_argument( "--area", help="Name of area variable on netcdf input file." @@ -109,27 +115,20 @@ def get_parser(): dest="overwrite", required=False, ) - parser.add_argument( - "-v", - "--verbose", - help="Increase output verbosity", - action="store_true", - dest="verbose", - required=False, - ) + add_logging_args(parser) return parser def main(): + setup_logging_pre_config() parser = get_parser() args = parser.parse_args() - if args.verbose: - logging.basicConfig(level=logging.DEBUG, format=" %(levelname)-8s :: %(message)s") - else: - logging.basicConfig(level=logging.INFO, format=" %(levelname)-8s :: %(message)s") + # --------------------------------- # + # process logging args (i.e. debug and verbose) + process_logging_args(args) nc_file = args.input lat_name = args.lat_name From 96ef2fdd7f11ff48cfd55c3ae45f278af8430806 Mon Sep 17 00:00:00 2001 From: Erik Kluzek Date: Wed, 1 Feb 2023 23:59:56 -0700 Subject: [PATCH 064/595] Add check_args for mesh_maker to shorten length of main and fix more pylint issues --- python/ctsm/.pylintrc | 1 + python/ctsm/mesh_maker.py | 70 ++++++++++++++++++++++----------------- 2 files changed, 40 insertions(+), 31 deletions(-) diff --git a/python/ctsm/.pylintrc b/python/ctsm/.pylintrc index 6fe02a09da..429d154a40 100644 --- a/python/ctsm/.pylintrc +++ b/python/ctsm/.pylintrc @@ -433,6 +433,7 @@ good-names=i, k, ex, Run, + ds, m, n, l1, diff --git a/python/ctsm/mesh_maker.py b/python/ctsm/mesh_maker.py index f5ea4dcd31..948688eca7 100644 --- a/python/ctsm/mesh_maker.py +++ b/python/ctsm/mesh_maker.py @@ -120,26 +120,9 @@ def get_parser(): return parser -def main(): - - setup_logging_pre_config() - parser = get_parser() - args = parser.parse_args() - - # --------------------------------- # - # process logging args (i.e. debug and verbose) - process_logging_args(args) - - nc_file = args.input - lat_name = args.lat_name - lon_name = args.lon_name - mesh_out = args.output - out_dir = args.out_dir - overwrite = args.overwrite - mask_name = args.mask_name - area_name = args.area_name - - if mesh_out and out_dir: +def check_args(args): + """Check the arguments""" + if args.output and args.out_dir: logging.error(" Both --outdir and --output cannot be provided at the same time.") err_msg = textwrap.dedent( """ @@ -153,31 +136,56 @@ def main(): abort(err_msg) # -- no file name and output path: - if not mesh_out and not out_dir: - out_dir = os.path.join(os.getcwd(), "meshes") + if not args.output and not args.out_dir: + args.out_dir = os.path.join(os.getcwd(), "meshes") - if not mesh_out: + if not args.output: # -- make output path if does not exist. - if not os.path.isdir(out_dir): - os.mkdir(out_dir) + if not os.path.isdir(args.out_dir): + os.mkdir(args.out_dir) today = datetime.today() today_string = today.strftime("%y%m%d") - mesh_out = os.path.join( - out_dir, - os.path.splitext(nc_file)[0] + "_ESMF_UNSTRUCTURED_MESH" + "_c" + today_string + ".nc", + args.output = os.path.join( + args.out_dir, + os.path.splitext(args.input)[0] + + "_ESMF_UNSTRUCTURED_MESH" + + "_c" + + today_string + + ".nc", ) # -- exit if mesh_out exists and --overwrite is not specified. - if os.path.exists(mesh_out): - if overwrite: - os.remove(mesh_out) + if os.path.exists(args.output): + if args.overwrite: + os.remove(args.output) else: err_msg = ( "output meshfile exists, please choose --overwrite to overwrite the mesh file." ) abort(err_msg) + +def main(): + """Main function to create a mesh file from another file""" + + setup_logging_pre_config() + parser = get_parser() + args = parser.parse_args() + + # --------------------------------- # + # process logging args (i.e. debug and verbose) + process_logging_args(args) + + check_args(args) + + nc_file = args.input + lat_name = args.lat_name + lon_name = args.lon_name + mesh_out = args.output + mask_name = args.mask_name + area_name = args.area_name + if os.path.isfile(nc_file): ds = xr.open_dataset(nc_file, mask_and_scale=False, decode_times=False).transpose() else: From 7fbdb0ca1094314f6180d59f57d7ddbafd0649d2 Mon Sep 17 00:00:00 2001 From: Erik Kluzek Date: Thu, 2 Feb 2023 00:10:32 -0700 Subject: [PATCH 065/595] Add a check that arguments get changed as expected --- python/ctsm/mesh_maker.py | 4 +++- python/ctsm/test/test_unit_mesh_maker.py | 22 ++++++++++++++++++++-- 2 files changed, 23 insertions(+), 3 deletions(-) diff --git a/python/ctsm/mesh_maker.py b/python/ctsm/mesh_maker.py index 948688eca7..799b2d514f 100644 --- a/python/ctsm/mesh_maker.py +++ b/python/ctsm/mesh_maker.py @@ -165,6 +165,8 @@ def check_args(args): ) abort(err_msg) + return args + def main(): """Main function to create a mesh file from another file""" @@ -177,7 +179,7 @@ def main(): # process logging args (i.e. debug and verbose) process_logging_args(args) - check_args(args) + args = check_args(args) nc_file = args.input lat_name = args.lat_name diff --git a/python/ctsm/test/test_unit_mesh_maker.py b/python/ctsm/test/test_unit_mesh_maker.py index 54f5f689e8..8903d1e2b6 100755 --- a/python/ctsm/test/test_unit_mesh_maker.py +++ b/python/ctsm/test/test_unit_mesh_maker.py @@ -7,7 +7,6 @@ """ import unittest -import argparse import os import sys @@ -17,7 +16,7 @@ # pylint: disable=wrong-import-position from ctsm import unit_testing -from ctsm.mesh_maker import get_parser, main +from ctsm.mesh_maker import get_parser, check_args, main # pylint: disable=invalid-name @@ -88,6 +87,25 @@ def test_outfile_exists_and_no_overwrite(self): ): main() + def test_default_outfile_as_expected(self): + """ + Test that the default outfile is as expected + """ + infile = "inputfile.nc" + sys.argv = [ + "mesh_maker", + "--input", + infile, + "--lat", + "lsmlat", + "--lon", + "lsmlon", + ] + parser = get_parser() + args = parser.parse_args() + args = check_args(args) + expected_outdir = os.path.join(os.getcwd(), "meshes") + self.assertEqual( args.out_dir, expected_outdir, "Default out_dir is not as expected") if __name__ == "__main__": unit_testing.setup_for_tests() From 013954bcd4d657ea67d7e0d0cdbe6062652b653b Mon Sep 17 00:00:00 2001 From: Erik Kluzek Date: Thu, 2 Feb 2023 00:14:57 -0700 Subject: [PATCH 066/595] Change name of subroutine --- python/ctsm/mesh_maker.py | 6 +++--- python/ctsm/test/test_unit_mesh_maker.py | 4 ++-- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/python/ctsm/mesh_maker.py b/python/ctsm/mesh_maker.py index 799b2d514f..a6adb5583e 100644 --- a/python/ctsm/mesh_maker.py +++ b/python/ctsm/mesh_maker.py @@ -120,8 +120,8 @@ def get_parser(): return parser -def check_args(args): - """Check the arguments""" +def process_and_check_args(args): + """Process and check the arguments""" if args.output and args.out_dir: logging.error(" Both --outdir and --output cannot be provided at the same time.") err_msg = textwrap.dedent( @@ -179,7 +179,7 @@ def main(): # process logging args (i.e. debug and verbose) process_logging_args(args) - args = check_args(args) + args = process_and_check_args(args) nc_file = args.input lat_name = args.lat_name diff --git a/python/ctsm/test/test_unit_mesh_maker.py b/python/ctsm/test/test_unit_mesh_maker.py index 8903d1e2b6..163493d601 100755 --- a/python/ctsm/test/test_unit_mesh_maker.py +++ b/python/ctsm/test/test_unit_mesh_maker.py @@ -16,7 +16,7 @@ # pylint: disable=wrong-import-position from ctsm import unit_testing -from ctsm.mesh_maker import get_parser, check_args, main +from ctsm.mesh_maker import get_parser, process_and_check_args, main # pylint: disable=invalid-name @@ -103,7 +103,7 @@ def test_default_outfile_as_expected(self): ] parser = get_parser() args = parser.parse_args() - args = check_args(args) + args = process_and_check_args(args) expected_outdir = os.path.join(os.getcwd(), "meshes") self.assertEqual( args.out_dir, expected_outdir, "Default out_dir is not as expected") From 1baf7c5db49c6c78c04f556fb7ce6487be9cdf8c Mon Sep 17 00:00:00 2001 From: Erik Kluzek Date: Thu, 2 Feb 2023 00:37:54 -0700 Subject: [PATCH 067/595] Start adding a system test for mesh_maker to check file data --- python/ctsm/test/test_sys_mesh_maker.py | 61 ++++++++++++++++++++++++ python/ctsm/test/test_unit_mesh_maker.py | 3 +- 2 files changed, 63 insertions(+), 1 deletion(-) create mode 100755 python/ctsm/test/test_sys_mesh_maker.py diff --git a/python/ctsm/test/test_sys_mesh_maker.py b/python/ctsm/test/test_sys_mesh_maker.py new file mode 100755 index 0000000000..9fbcbf7e9d --- /dev/null +++ b/python/ctsm/test/test_sys_mesh_maker.py @@ -0,0 +1,61 @@ +#!/usr/bin/env python3 +""" +System tests for mesh_maker + +""" + +import unittest +import os +import sys +import tempfile +import shutil + +# pylint: disable=wrong-import-position +from ctsm.path_utils import path_to_ctsm_root +from ctsm import unit_testing +from ctsm.mesh_maker import main + +# pylint: disable=invalid-name + + +class SysTestMeshMaker(unittest.TestCase): + """ + Basic class for testing mesh_maker.py. + """ + + def setUp(self): + """Setup for all tests""" + testinputs_path = os.path.join(path_to_ctsm_root(), "python/ctsm/test/testinputs") + self._testinputs_path = testinputs_path + self._fsurdat_in = os.path.join( + testinputs_path, + "surfdata_1x1_mexicocityMEX_hist_16pfts_Irrig_CMIP6_simyr2000_c221206_modified.nc", + ) + self._tempdir = tempfile.mkdtemp() + self.mesh_out = os.path.join(self._tempdir, "mesh_out.nc") + + def tearDown(self): + """ + Remove temporary directory + """ + shutil.rmtree(self._tempdir, ignore_errors=True) + + def test_basic(self): + """Do a simple basic test""" + sys.argv = [ + "mesh_maker", + "--input", + self._fsurdat_in, + "--lat", + "lsmlat", + "--lon", + "lsmlon", + "--output", + self.mesh_out, + ] + main() + + +if __name__ == "__main__": + unit_testing.setup_for_tests() + unittest.main() diff --git a/python/ctsm/test/test_unit_mesh_maker.py b/python/ctsm/test/test_unit_mesh_maker.py index 163493d601..5668e1484e 100755 --- a/python/ctsm/test/test_unit_mesh_maker.py +++ b/python/ctsm/test/test_unit_mesh_maker.py @@ -105,7 +105,8 @@ def test_default_outfile_as_expected(self): args = parser.parse_args() args = process_and_check_args(args) expected_outdir = os.path.join(os.getcwd(), "meshes") - self.assertEqual( args.out_dir, expected_outdir, "Default out_dir is not as expected") + self.assertEqual(args.out_dir, expected_outdir, "Default out_dir is not as expected") + if __name__ == "__main__": unit_testing.setup_for_tests() From 9dc675993cd8ff3a2d1cc4adf6c6b209ba91c64f Mon Sep 17 00:00:00 2001 From: Erik Kluzek Date: Wed, 15 Feb 2023 11:50:59 -0700 Subject: [PATCH 068/595] Bunch of work at starting to use conda-lock to save actual running environments rather than a minimalist conda requirements file --- py_env_create | 133 +- python/conda_env_condalock.txt | 7 + python/conda_env_ctsm_py.txt | 1 + python/conda_lock_ctsm_py_linux-64.yml | 2845 ++++++++++++++++++++++++ 4 files changed, 2975 insertions(+), 11 deletions(-) create mode 100644 python/conda_env_condalock.txt create mode 100644 python/conda_lock_ctsm_py_linux-64.yml diff --git a/py_env_create b/py_env_create index c323a374df..09619a0c82 100755 --- a/py_env_create +++ b/py_env_create @@ -20,19 +20,25 @@ if [ $error != 0 ]; then echo "For notes on installing on a user system see: https://docs.conda.io/projects/conda/en/latest/user-guide/install/index.html" echo "Error code was $error" cat condahelp.txt + rm condahelp.txt exit -1 fi rm condahelp.txt ctsm_python=ctsm_pylib +condalockenv=condalock condadir="$dir/python" domain=`domainname` +platform="linux-64" +condalockenvfile="conda_env_condalock.txt" if [[ $domain =~ cgd.* ]]; then condafile="conda_env_ctsm_py_cgd.txt" + lockfile="conda_lock_ctsm_py_cgd_${platform}.yml" else condafile="conda_env_ctsm_py.txt" + lockfile="conda_lock_ctsm_py_${platform}.yml" fi #---------------------------------------------------------------------- # Usage subroutine @@ -48,15 +54,24 @@ usage() { echo "[-v|--verbose] " echo " Run with verbose mode for the install so you see the progress bar" echo "[-f|--file ] " - echo " Conda environment file to use (can be a text format or YAML format)" + echo " Conda environment requirements text file to use (text format)" echo " Assumed to be under the directory: $condadir" echo " Default is: $condafile" + echo "[-l|--read-lock] " + echo " Read the conda-lock file for this platform rather than the requirements file" + echo " The lock file will be faster and is less likely to run into conflicts." + echo " You must have run this script previously to get an environemnt with conda-lock installed." + echo " Will use: $condadir/$lockfile if it exists" + echo "[--write-lock] " + echo " Write the conda-lock file for this platform" echo "[--option