From d3945b39dbbeca3c31c46c27d7a252b583ce5e6d Mon Sep 17 00:00:00 2001 From: Lizzie Lundgren Date: Tue, 25 Jul 2023 12:12:57 -0400 Subject: [PATCH 1/5] Add cubed sphere tools from gcgridobj repository The new tools file, cstools.py, includes an example of how to extract cubed-sphere grid indexes (nf, Ydim, Xdim) given a target latitude and longitude, and how to use them to extract data for a specific lat-lon box. Using this tool requires python package pyproj not previously required. GCPy will still function without this package. However, you will not be able to call function find_index. The code added was originally developed by Liam Bindle and Sebastian Eastham and was taken directly from the gcgridobj repository. I made a minor change to remove np.int (deprecated) and add example and header. Signed-off-by: Lizzie Lundgren --- gcpy/__init__.py | 1 + gcpy/cstools.py | 342 +++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 343 insertions(+) create mode 100644 gcpy/cstools.py diff --git a/gcpy/__init__.py b/gcpy/__init__.py index 8d4d6e0d..9d9efbb1 100644 --- a/gcpy/__init__.py +++ b/gcpy/__init__.py @@ -26,3 +26,4 @@ from .benchmark import * from .file_regrid import * from .grid_stretching_transforms import * +from .cstools import * diff --git a/gcpy/cstools.py b/gcpy/cstools.py new file mode 100644 index 00000000..6c1d57f9 --- /dev/null +++ b/gcpy/cstools.py @@ -0,0 +1,342 @@ +""" +Contains tools for working with cubed-sphere data. + +Originally developed by Liam Bindle and Sebastian Eastham. + +Example: + import gcpy + c24_grid=gcpy.gen_grid(24) # Generate C24 grid + lat=40.0 # Specify target latitude + lon=150.0 # Specify target longitude + idx=gcpy.find_index(lat,lon,c24_grid) # Returns numpy ndarray [[3],[7],[1]] which can be + used to index data (nf=3, Ydim=7, Xdim=1) + + # Check and use to get data + datafile='/n/home/GeosChem.SpeciesConc.20190701_0000z.nc4' + import xarray ax xr + ds=xr.open_dataset(datafile) + nf=idx[0,0] + Ydim=idx[1,0] + Xdim=idx[2,0] + ds['lats'].isel(nf=nf,Ydim=Ydim,Xdim=Xdim].item() # 38.711082458496094 + ds['lons'].isel(nf=nf,Ydim=Ydim,Xdim=Xdim].item() # 151.61871337890625 + ds['SpeciesConcVV_O3'].isel(time=0,lev=0,nf=nf,Ydim=Ydim,Xdim=Xdim).item() # 2.7790051149167994e-08 + +""" + + +import numpy as np +import xarray as xr + +import gcpy + +try: + import pyproj + import shapely.ops + import shapely.geometry + find_point_OK = True +except: + find_point_OK = False + +# Must have: +# 1. extract_grid (returns an xarray Dataset) +# 2. grid_area (returns a 6xNxN array) +# 3. gen_grid (returns an xarray Dataset) + +def extract_grid(ds,src_var='Xdim'): + # Extract grid from xarray dataset but return a cubed-sphere grid + n_cs = ds[src_var].shape[-1] + return gen_grid(n_cs) + +def read_gridspec(gs_obj): + # Reads a gridspec object and returns an xarray dataset + n_cs = gs_obj._tiles[0].area.shape[0] + lon = np.zeros((6,n_cs,n_cs)) + lon_b = np.zeros((6,n_cs+1,n_cs+1)) + lat = np.zeros((6,n_cs,n_cs)) + lat_b = np.zeros((6,n_cs+1,n_cs+1)) + area = np.zeros((6,n_cs,n_cs)) + for i in range(6): + tile = gs_obj._tiles[i] + lon_b[i,...] = tile.supergrid_lons[::2,::2] # Identical to original definition + lat_b[i,...] = tile.supergrid_lats[::2,::2] # Identical to original definition + lon[i,...] = tile.supergrid_lons[1::2,1::2] # NOT identical to original definition + lat[i,...] = tile.supergrid_lats[1::2,1::2] # NOT identical to original definition + area[i,...] = tile.area[...] + ds = xr.Dataset( + data_vars=dict( + area=(['nf','Ydim','Xdim'],area), + lon=(['nf','Ydim','Xdim'],lon), + lat=(['nf','Ydim','Xdim'],lat), + lon_b=(['nf','Ydim_b','Xdim_b'],lon_b), + lat_b=(['nf','Ydim_b','Xdim_b'],lat_b), + ), + coords=dict( + nf=(['nf'],list(range(6))), + Ydim=(['Ydim'],list(range(n_cs))), + Xdim=(['Xdim'],list(range(n_cs))), + Ydim_b=(['Ydim_b'],list(range(n_cs+1))), + Xdim_b=(['Xdim_b'],list(range(n_cs+1))), + ), + attrs=dict(description=f'c{n_cs:d} grid data'), + ) + return ds + +def face_area(lon_b, lat_b, r_sphere = 6.375e6): + """Calculate area of cubed-sphere grid cells on one face + Inputs must be in degrees. Edge arrays must be + shaped [N+1 x N+1] + """ + + # Convert inputs to radians + lon_b_rad = lon_b * np.pi / 180.0 + lat_b_rad = lat_b * np.pi / 180.0 + + r_sq = r_sphere * r_sphere + n_cs = lon_b.shape[1] - 1 + + # Allocate output array + cs_area = np.zeros((n_cs,n_cs)) + + # Ordering + valid_combo = np.array([[1,2,4],[2,3,1],[3,2,4],[4,1,3]]) - 1 + + for i_lon in range(n_cs): + for i_lat in range(n_cs): + lon_corner = np.zeros(4) + lat_corner = np.zeros(4) + xyz_corner = np.zeros((4,3)) + for i_vert in range(4): + x_lon = i_lon + (i_vert > 1) + x_lat = i_lat + (i_vert == 0 or i_vert == 3) + lon_corner[i_vert] = lon_b_rad[x_lon,x_lat] + lat_corner[i_vert] = lat_b_rad[x_lon,x_lat] + for i_vert in range(4): + xyz_corner[i_vert,:] = ll2xyz(lon_corner[i_vert],lat_corner[i_vert]) + tot_ang = 0.0 + for i_corner in range(4): + curr_combo = valid_combo[i_corner,:] + xyz_mini = np.zeros((3,3)) + for i_mini in range(3): + xyz_mini[i_mini,:] = xyz_corner[curr_combo[i_mini],:] + curr_ang = sphere_angle(xyz_mini[0,:],xyz_mini[1,:],xyz_mini[2,:]) + tot_ang += curr_ang + cs_area[i_lon,i_lat] = r_sq * (tot_ang - (2.0*np.pi)) + + return cs_area + +def ll2xyz(lon_pt,lat_pt): + """Converts a lon/lat pair (in radians) to cartesian co-ordinates + Vector should point to the surface of the unit sphere""" + + xPt = np.cos(lat_pt) * np.cos(lon_pt) + yPt = np.cos(lat_pt) * np.sin(lon_pt) + zPt = np.sin(lat_pt) + return [xPt,yPt,zPt] + +def sphere_angle(e1,e2,e3): + # e1: Mid-point + # e2 and e3 to either side + pVec = np.ones(3) + qVec = np.ones(3) + pVec[0] = e1[1]*e2[2] - e1[2]*e2[1] + pVec[1] = e1[2]*e2[0] - e1[0]*e2[2] + pVec[2] = e1[0]*e2[1] - e1[1]*e2[0] + + qVec[0] = e1[1]*e3[2] - e1[2]*e3[1] + qVec[1] = e1[2]*e3[0] - e1[0]*e3[2] + qVec[2] = e1[0]*e3[1] - e1[1]*e3[0] + ddd = np.sum(pVec*pVec) * np.sum(qVec*qVec) + if ddd <= 0.0: + angle = 0.0; + else: + ddd = np.sum(pVec*qVec)/np.sqrt(ddd); + if (np.abs(ddd)>1.0): + angle = np.pi/2.0; + else: + angle = np.arccos(ddd); + + return angle + +def grid_area(cs_grid=None,cs_res=None): + """Return area in m2 for each cell in a cubed-sphere grid + Uses GMAO indexing convention (6xNxN) + """ + # Calculate area on a cubed sphere + if cs_res is None: + cs_res = cs_grid['lon_b'].shape[-1] - 1 + elif cs_grid is None: + cs_grid = gcpy.csgrid_GMAO(cs_res) + elif cs_grid is not None and cs_res is not None: + assert cs_res == cs_grid['lon_b'].shape[-1], 'Routine grid_area received inconsistent inputs' + cs_area = np.zeros((6,cs_res,cs_res)) + cs_area[0,:,:] = face_area(cs_grid['lon_b'][0,:,:],cs_grid['lat_b'][0,:,:]) + for i_face in range(1,6): + cs_area[i_face,:,:] = cs_area[0,:,:].copy() + return cs_area + +def gen_grid(n_cs, stretch_factor=None, target_lon=None, target_lat=None): + if stretch_factor is not None: + cs_temp, ignore = gcpy.make_grid_SG(n_cs,stretch_factor,target_lon,target_lat) + else: + cs_temp = gcpy.csgrid_GMAO(n_cs) + return xr.Dataset({'nf': (['nf'],np.array(range(6))), + 'Ydim': (['Ydim'],np.array(range(n_cs))), + 'Xdim': (['Xdim'],np.array(range(n_cs))), + 'Ydim_b': (['Ydim_b'],np.array(range(n_cs+1))), + 'Xdim_b': (['Xdim_b'],np.array(range(n_cs+1))), + 'lat': (['nf','Ydim','Xdim'], cs_temp['lat']), + 'lon': (['nf','Ydim','Xdim'], cs_temp['lon']), + 'lat_b': (['nf','Ydim_b','Xdim_b'], cs_temp['lat_b']), + 'lon_b': (['nf','Ydim_b','Xdim_b'], cs_temp['lon_b']), + 'area': (['nf','Ydim','Xdim'], grid_area(cs_temp))}) + +def corners_to_xy(xc, yc): + """ Creates xy coordinates for each grid-box. The shape is (n, n, 5) where n is the cubed-sphere size. + Developed, tested, and supplied by Liam Bindle. + + :param xc: grid-box corner longitudes; shape (n+1, n+1) + :param yc: grid-box corner latitudes; shape (n+1, n+1) + :return: grid-box xy coordinates + """ + p0 = slice(0, -1) + p1 = slice(1, None) + boxes_x = np.moveaxis(np.array([xc[p0, p0], xc[p1, p0], xc[p1, p1], xc[p0, p1], xc[p0, p0]]), 0, -1) + boxes_y = np.moveaxis(np.array([yc[p0, p0], yc[p1, p0], yc[p1, p1], yc[p0, p1], yc[p0, p0]]), 0, -1) + return np.moveaxis(np.array([boxes_x, boxes_y]), 0, -1) + + +def central_angle(x0, y0, x1, y1): + """ Returns the distance (central angle) between coordinates (x0, y0) and (x1, y1). This is vectorizable. + Developed, tested, and supplied by Liam Bindle. + + :param x0: pt0's longitude (degrees) + :param y0: pt0's latitude (degrees) + :param x1: pt1's longitude (degrees) + :param y1: pt1's latitude (degrees) + :return: Distance (degrees) + """ + RAD2DEG = 180 / np.pi + DEG2RAD = np.pi / 180 + x0 = x0 * DEG2RAD + x1 = x1 * DEG2RAD + y0 = y0 * DEG2RAD + y1 = y1 * DEG2RAD + return np.arccos(np.sin(y0) * np.sin(y1) + np.cos(y0) * np.cos(y1) * np.cos(np.abs(x0-x1))) * RAD2DEG + +def find_index_single(lat,lon,x_centers_flat,y_centers_flat,xy,cs_size,latlon_crs,jitter_size=0.0): + #import pyproj + #import shapely.ops + #import shapely.geometry + # Center on x_find, y_find + x_find = lon + y_find = lat + gnomonic_crs = pyproj.Proj(f'+proj=gnom +lat_0={y_find} +lon_0={x_find}') + + # Generate all distances + distances = central_angle(x_find, y_find, x_centers_flat, y_centers_flat) + four_nearest_indexes = np.argpartition(distances, 4)[:4] + + # Unravel 4 smallest indexes + four_nearest_indexes = np.unravel_index(four_nearest_indexes, (6, cs_size, cs_size)) + four_nearest_xy = xy[four_nearest_indexes] + four_nearest_polygons = [shapely.geometry.Polygon(polygon_xy) for polygon_xy in four_nearest_xy] + + # Transform to gnomonic projection + gno_transform = pyproj.Transformer.from_proj(latlon_crs, gnomonic_crs, always_xy=True).transform + four_nearest_polygons_gno = [shapely.ops.transform(gno_transform, polygon) for polygon in four_nearest_polygons] + + # Figure out which polygon contains the point + Xy_find = shapely.geometry.Point(x_find, y_find) + Xy_find_GNO = shapely.ops.transform(gno_transform, Xy_find) + polygon_contains_point = [polygon.contains(Xy_find_GNO) for polygon in four_nearest_polygons_gno] + + if np.count_nonzero(polygon_contains_point) == 0: + if jitter_size>0.0: + # Move longitude by ~1 m + nf, YDim, XDim = find_index_single(y_find,x_find+jitter_size,x_centers_flat,y_centers_flat,xy,cs_size,latlon_crs,jitter_size=0.0) + else: + raise ValueError(f'Point at {x_find:8.2f} E, {y_find:8.2f} N could not be matched') + # The first will be selected, if more than one + polygon_with_point = np.argmax(polygon_contains_point) + + # Get original index + nf = four_nearest_indexes[0][polygon_with_point] + YDim = four_nearest_indexes[1][polygon_with_point] + XDim = four_nearest_indexes[2][polygon_with_point] + + return nf, YDim, XDim + +def find_index(lat,lon,grid,jitter_size=0.0): + # For point-finding + #import pyproj + #import shapely.ops + #import shapely.geometry + assert find_point_OK, "Cannot perform index finding - need pyproj and shapely" + + # Based on a routine developed, tested, and supplied by Liam Bindle. + lon_vec = np.asarray(lon) + lat_vec = np.asarray(lat) + n_find = lon_vec.size + + # Get the corners + x_corners = grid['lon_b'].values + y_corners = grid['lat_b'].values + x_centers = grid['lon'].values + y_centers = grid['lat'].values + x_centers_flat = x_centers.flatten() + y_centers_flat = y_centers.flatten() + + cs_size = x_centers.shape[-1] + + # Generate everything that will be reused + # Get XY polygon definitions for grid boxes + xy = np.zeros((6, cs_size, cs_size, 5, 2)) # 5 (x,y) points defining polygon corners (first and last are same) + for nf in range(6): + xy[nf, ...] = corners_to_xy(xc=x_corners[nf, :, :], yc=y_corners[nf, :, :]) + latlon_crs = pyproj.Proj("+proj=latlon") + + # Find 4 shortest distances to (x_find, y_find) + #idx = np.full((3,n_find),np.int(0)) # np.int no longer in numpy (ewl) + idx = np.full((3,n_find),0) + for x_find, y_find, i_find in zip(np.nditer(lon_vec),np.nditer(lat_vec),list(range(n_find))): + ## Center on x_find, y_find + #gnomonic_crs = pyproj.Proj(f'+proj=gnom +lat_0={y_find} +lon_0={x_find}') + + ## Generate all distances + #distances = central_angle(x_find, y_find, x_centers_flat, y_centers_flat) + #four_nearest_indexes = np.argpartition(distances, 4)[:4] + + ## Unravel 4 smallest indexes + #four_nearest_indexes = np.unravel_index(four_nearest_indexes, (6, cs_size, cs_size)) + #four_nearest_xy = xy[four_nearest_indexes] + #four_nearest_polygons = [shapely.geometry.Polygon(polygon_xy) for polygon_xy in four_nearest_xy] + + ## Transform to gnomonic projection + #gno_transform = pyproj.Transformer.from_proj(latlon_crs, gnomonic_crs, always_xy=True).transform + #four_nearest_polygons_gno = [shapely.ops.transform(gno_transform, polygon) for polygon in four_nearest_polygons] + + ## Figure out which polygon contains the point + #Xy_find = shapely.geometry.Point(x_find, y_find) + #Xy_find_GNO = shapely.ops.transform(gno_transform, Xy_find) + #polygon_contains_point = [polygon.contains(Xy_find_GNO) for polygon in four_nearest_polygons_gno] + + ##assert np.count_nonzero(polygon_contains_point) == 1 + ##assert np.count_nonzero(polygon_contains_point) > 0 + #if np.count_nonzero(polygon_contains_point) == 0: + # if allow_jitter: + # idx[:,i_find] = + # raise ValueError(f'Point at {x_find:8.2f} E, {y_find:8.2f} N could not be matched') + ## The first will be selected, if more than one + #polygon_with_point = np.argmax(polygon_contains_point) + + ## Get original index + #nf = four_nearest_indexes[0][polygon_with_point] + #YDim= four_nearest_indexes[1][polygon_with_point] + #XDim= four_nearest_indexes[2][polygon_with_point] + nf, YDim, XDim = find_index_single(y_find,x_find,x_centers_flat,y_centers_flat, + xy,cs_size,latlon_crs,jitter_size=jitter_size) + idx[:,i_find] = [nf,YDim,XDim] + return idx + + From 925fd88ea28bd8f784caadec87899104deae99ac Mon Sep 17 00:00:00 2001 From: Bob Yantosca Date: Tue, 25 Jul 2023 12:30:14 -0400 Subject: [PATCH 2/5] Add pyproj==3.6.0 to the gcpy_env environment.yml ifle docs/environment_files/environment.yml - The updates to the cubed-sphere tools require the pyproj library. We have now added pyproj==3.6.0 as a PyPi dependency. Signed-off-by: Bob Yantosca --- docs/environment_files/environment.yml | 1 + 1 file changed, 1 insertion(+) diff --git a/docs/environment_files/environment.yml b/docs/environment_files/environment.yml index a9497bf1..f77acee8 100644 --- a/docs/environment_files/environment.yml +++ b/docs/environment_files/environment.yml @@ -27,6 +27,7 @@ dependencies: - numpy==1.21.1 # Optimized mathematical functions - pandas==1.3.1 # Tables/timeseries manipulation - pycodestyle==2.9.1 # Python style checker (formerly PEP8) + - pyproj==3.6.0 # Python map projections library - pylint==2.15.3 # Python linter - pypdf2==1.26.0 # PDF utilities (bookmarks, etc.) - recommonmark==0.7.1 # Dependency for Sphinx From 3699f98e3a24e21af48d0fe4587f007c65dccf15 Mon Sep 17 00:00:00 2001 From: Bob Yantosca Date: Tue, 25 Jul 2023 16:19:38 -0400 Subject: [PATCH 3/5] Cleaned up cstools.py and adopted style suggestions from Pylint gcpy/cstools.py - Ran through the Pylint linter and adopted several style suggestions: - Updated Pydoc strings / Added missing Pydoc strings - Redefined variable names to adhere to snake_case style - Now define radians to degrees & vice versa constants at the top of the module - Make sure most lines are within 80 characters --- gcpy/cstools.py | 678 ++++++++++++++++++++++++++++++++++-------------- 1 file changed, 488 insertions(+), 190 deletions(-) diff --git a/gcpy/cstools.py b/gcpy/cstools.py index 6c1d57f9..0ef94124 100644 --- a/gcpy/cstools.py +++ b/gcpy/cstools.py @@ -2,14 +2,19 @@ Contains tools for working with cubed-sphere data. Originally developed by Liam Bindle and Sebastian Eastham. +Included into GCPy by Lizzie Lundgren and Bob Yantosca. + +Style updates suggested by the "Pylint" python linter have been adopted. Example: import gcpy c24_grid=gcpy.gen_grid(24) # Generate C24 grid lat=40.0 # Specify target latitude lon=150.0 # Specify target longitude - idx=gcpy.find_index(lat,lon,c24_grid) # Returns numpy ndarray [[3],[7],[1]] which can be - used to index data (nf=3, Ydim=7, Xdim=1) + idx=gcpy.find_index(lat,lon,c24_grid) # Returns numpy ndarray + # [[3],[7],[1]] which can be + # used to index data + # (nf=3, Ydim=7, Xdim=1) # Check and use to get data datafile='/n/home/GeosChem.SpeciesConc.20190701_0000z.nc4' @@ -18,52 +23,85 @@ nf=idx[0,0] Ydim=idx[1,0] Xdim=idx[2,0] - ds['lats'].isel(nf=nf,Ydim=Ydim,Xdim=Xdim].item() # 38.711082458496094 - ds['lons'].isel(nf=nf,Ydim=Ydim,Xdim=Xdim].item() # 151.61871337890625 - ds['SpeciesConcVV_O3'].isel(time=0,lev=0,nf=nf,Ydim=Ydim,Xdim=Xdim).item() # 2.7790051149167994e-08 - + ds['lats'].isel(nf=nf,Ydim=Ydim,Xdim=Xdim].item() + # prints 38.711082458496094 + ds['lons'].isel(nf=nf,Ydim=Ydim,Xdim=Xdim].item() + # prints 151.61871337890625 + ds['SpeciesConcVV_O3'].isel(time=0,lev=0,nf=nf,Ydim=Ydim,Xdim=Xdim).item() + # prints 2.7790051149167994e-08 """ - - import numpy as np import xarray as xr - import gcpy - try: import pyproj import shapely.ops import shapely.geometry - find_point_OK = True -except: - find_point_OK = False - -# Must have: -# 1. extract_grid (returns an xarray Dataset) -# 2. grid_area (returns a 6xNxN array) -# 3. gen_grid (returns an xarray Dataset) - -def extract_grid(ds,src_var='Xdim'): - # Extract grid from xarray dataset but return a cubed-sphere grid - n_cs = ds[src_var].shape[-1] +except ImportError as exc: + raise ImportError( + "gcpy.cstools needs packages 'pyproj' and 'shapely'!" + ) from exc + +# Constants +RAD_TO_DEG = 180.0 / np.pi +DEG_TO_RAD = np.pi / 180.0 + + +def extract_grid( + data, + src_var='Xdim' +): + """ + Extracts the grid information from an xarray.Dataset object and + returns a new xarray.Dataset object on a cubed-sphere grid. + + Args: + ----- + data : xarray.Dataset + The input dataset + + data_cs: xarray.Dataset + Same data as in argument "ds", but on a cubed-sphere grid + """ + if not isinstance(data, xr.Dataset): + return TypeError("Argument 'ds' is not of type 'xarray.Dataset'!") + n_cs = data[src_var].shape[-1] return gen_grid(n_cs) + def read_gridspec(gs_obj): - # Reads a gridspec object and returns an xarray dataset + """ + Reads a GridSpec object and returns an xarray.Dataset object. + + Args: + ----- + gs_obj : GridSpec + The GridSpec object as input + + Returns: + -------- + ds : xarray.Dataset + The same data as an xarray.Dataset object. + """ n_cs = gs_obj._tiles[0].area.shape[0] - lon = np.zeros((6,n_cs,n_cs)) - lon_b = np.zeros((6,n_cs+1,n_cs+1)) - lat = np.zeros((6,n_cs,n_cs)) - lat_b = np.zeros((6,n_cs+1,n_cs+1)) - area = np.zeros((6,n_cs,n_cs)) + lon = np.zeros((6, n_cs, n_cs)) + lon_b = np.zeros((6, n_cs+1, n_cs+1)) + lat = np.zeros((6, n_cs, n_cs)) + lat_b = np.zeros((6, n_cs+1, n_cs+1)) + area = np.zeros((6, n_cs, n_cs)) for i in range(6): tile = gs_obj._tiles[i] - lon_b[i,...] = tile.supergrid_lons[::2,::2] # Identical to original definition - lat_b[i,...] = tile.supergrid_lats[::2,::2] # Identical to original definition - lon[i,...] = tile.supergrid_lons[1::2,1::2] # NOT identical to original definition - lat[i,...] = tile.supergrid_lats[1::2,1::2] # NOT identical to original definition + # lon_b is identical to original definition + lon_b[i,...] = tile.supergrid_lons[::2,::2] + # lat_b is dentical to original definition + lat_b[i,...] = tile.supergrid_lats[::2,::2] + # lon is NOT identical to original definition + lon[i,...] = tile.supergrid_lons[1::2,1::2] + # lat is NOT identical to original definition + lat[i,...] = tile.supergrid_lats[1::2,1::2] area[i,...] = tile.area[...] - ds = xr.Dataset( + + data = xr.Dataset( data_vars=dict( area=(['nf','Ydim','Xdim'],area), lon=(['nf','Ydim','Xdim'],lon), @@ -80,27 +118,48 @@ def read_gridspec(gs_obj): ), attrs=dict(description=f'c{n_cs:d} grid data'), ) - return ds + return data -def face_area(lon_b, lat_b, r_sphere = 6.375e6): - """Calculate area of cubed-sphere grid cells on one face - Inputs must be in degrees. Edge arrays must be - shaped [N+1 x N+1] + +def face_area( + lon_b, + lat_b, + r_sphere=6.375e6 +): + """ + Calculates area of cubed-sphere grid cells on one face. + Inputs must be in degrees. Edge arrays must be shaped [N+1 x N+1]. + + Args: + ----- + lon_b, lat_b : list of float + Longitude and latitude bounds (degrees) + + Keyword Args (optional): + ------------------------ + r_sphere : float + Radius of Earth (meters). Default value: 6.375e6 + + Returns: + -------- + cs_area : numpy.ndarray + Array of surface area (m2) in each grid box of a + cubed-sphere grid face. """ - + # Convert inputs to radians - lon_b_rad = lon_b * np.pi / 180.0 - lat_b_rad = lat_b * np.pi / 180.0 - + lon_b_rad = lon_b * DEG_TO_RAD + lat_b_rad = lat_b * DEG_TO_RAD + r_sq = r_sphere * r_sphere n_cs = lon_b.shape[1] - 1 - + # Allocate output array cs_area = np.zeros((n_cs,n_cs)) - + # Ordering valid_combo = np.array([[1,2,4],[2,3,1],[3,2,4],[4,1,3]]) - 1 - + for i_lon in range(n_cs): for i_lat in range(n_cs): lon_corner = np.zeros(4) @@ -112,55 +171,120 @@ def face_area(lon_b, lat_b, r_sphere = 6.375e6): lon_corner[i_vert] = lon_b_rad[x_lon,x_lat] lat_corner[i_vert] = lat_b_rad[x_lon,x_lat] for i_vert in range(4): - xyz_corner[i_vert,:] = ll2xyz(lon_corner[i_vert],lat_corner[i_vert]) + xyz_corner[i_vert,:] = ll2xyz( + lon_corner[i_vert], + lat_corner[i_vert] + ) tot_ang = 0.0 for i_corner in range(4): curr_combo = valid_combo[i_corner,:] xyz_mini = np.zeros((3,3)) for i_mini in range(3): xyz_mini[i_mini,:] = xyz_corner[curr_combo[i_mini],:] - curr_ang = sphere_angle(xyz_mini[0,:],xyz_mini[1,:],xyz_mini[2,:]) + curr_ang = sphere_angle( + xyz_mini[0,:], + xyz_mini[1,:], + xyz_mini[2,:] + ) tot_ang += curr_ang cs_area[i_lon,i_lat] = r_sq * (tot_ang - (2.0*np.pi)) - + return cs_area -def ll2xyz(lon_pt,lat_pt): - """Converts a lon/lat pair (in radians) to cartesian co-ordinates - Vector should point to the surface of the unit sphere""" - - xPt = np.cos(lat_pt) * np.cos(lon_pt) - yPt = np.cos(lat_pt) * np.sin(lon_pt) - zPt = np.sin(lat_pt) - return [xPt,yPt,zPt] - -def sphere_angle(e1,e2,e3): - # e1: Mid-point - # e2 and e3 to either side - pVec = np.ones(3) - qVec = np.ones(3) - pVec[0] = e1[1]*e2[2] - e1[2]*e2[1] - pVec[1] = e1[2]*e2[0] - e1[0]*e2[2] - pVec[2] = e1[0]*e2[1] - e1[1]*e2[0] - - qVec[0] = e1[1]*e3[2] - e1[2]*e3[1] - qVec[1] = e1[2]*e3[0] - e1[0]*e3[2] - qVec[2] = e1[0]*e3[1] - e1[1]*e3[0] - ddd = np.sum(pVec*pVec) * np.sum(qVec*qVec) + +def ll2xyz( + lon_pt, + lat_pt +): + """ + Converts a lon/lat pair (in radians) to Cartesian co-ordinates. + This is vectorizable. + + Args: + ----- + lon_pt, lat_pt : float + Longitude & latitude in radians. + + Returns: + -------- + [x_pt, y_pt, z_pt] : list of numpy.float64 + Cartesian vector coordinates (X,Y.Z) normalized to + the unit sphere. + """ + x_pt = np.cos(lat_pt) * np.cos(lon_pt) + y_pt = np.cos(lat_pt) * np.sin(lon_pt) + z_pt = np.sin(lat_pt) + + return [x_pt, y_pt, z_pt] + + +def sphere_angle( + e_1, + e_2, + e_3 +): + """ + Computes the angle between 3 points on a sphere. + + Args: + ----- + e_1 : list of float + (x, y) coordinates at mid point + + e_2, e_3: : list of float + (x, y) coordinates at points on either side of midpoint + + Returns: + -------- + angle : float + The spherical angle at point e_1. + """ + p_vec = np.ones(3) + q_vec = np.ones(3) + p_vec[0] = e_1[1]*e_2[2] - e_1[2]*e_2[1] + p_vec[1] = e_1[2]*e_2[0] - e_1[0]*e_2[2] + p_vec[2] = e_1[0]*e_2[1] - e_1[1]*e_2[0] + + q_vec[0] = e_1[1]*e_3[2] - e_1[2]*e_3[1] + q_vec[1] = e_1[2]*e_3[0] - e_1[0]*e_3[2] + q_vec[2] = e_1[0]*e_3[1] - e_1[1]*e_3[0] + + ddd = np.sum(p_vec * p_vec) * np.sum(q_vec * q_vec) if ddd <= 0.0: - angle = 0.0; + angle = 0.0 else: - ddd = np.sum(pVec*qVec)/np.sqrt(ddd); - if (np.abs(ddd)>1.0): - angle = np.pi/2.0; + ddd = np.sum(p_vec * q_vec) / np.sqrt(ddd) + if np.abs(ddd) > 1.0: + angle = np.pi / 2.0 else: - angle = np.arccos(ddd); + angle = np.arccos(ddd) return angle -def grid_area(cs_grid=None,cs_res=None): - """Return area in m2 for each cell in a cubed-sphere grid - Uses GMAO indexing convention (6xNxN) + +def grid_area( + cs_grid=None, + cs_res=None +): + """ + Return area (m2) for each cell in a cubed-sphere grid + + Args: + ----- + cs_grid : dict + Cubed-sphere grid definition as a dict of: + {'lat' : lat midpoints, + 'lon' : lon midpoints, + 'lat_b' : lat edges, + 'lon_b' : lon edges} + where each value has an extra face dimension of length 6. + + Returns: + -------- + grid_area : numpy.ndarray + Surface area (m2) for each cell in the cubed-sphere grid. + NOTE: Uses GMAO convention, array shape = (6, n, n), + where n is the number of cells along a face edge. """ # Calculate area on a cubed sphere if cs_res is None: @@ -168,113 +292,309 @@ def grid_area(cs_grid=None,cs_res=None): elif cs_grid is None: cs_grid = gcpy.csgrid_GMAO(cs_res) elif cs_grid is not None and cs_res is not None: - assert cs_res == cs_grid['lon_b'].shape[-1], 'Routine grid_area received inconsistent inputs' + assert cs_res == cs_grid['lon_b'].shape[-1], \ + 'Routine grid_area received inconsistent inputs' cs_area = np.zeros((6,cs_res,cs_res)) - cs_area[0,:,:] = face_area(cs_grid['lon_b'][0,:,:],cs_grid['lat_b'][0,:,:]) + cs_area[0,:,:] = face_area( + cs_grid['lon_b'][0,:,:], + cs_grid['lat_b'][0,:,:] + ) for i_face in range(1,6): cs_area[i_face,:,:] = cs_area[0,:,:].copy() + return cs_area -def gen_grid(n_cs, stretch_factor=None, target_lon=None, target_lat=None): + +def gen_grid( + n_cs, + stretch_factor=None, + target_lon=None, + target_lat=None +): + """ + Returns an xarray.Dataset object specifying a cubed-sphere + stretched grid. + + Args: + n_cs : int + Number of grid boxes along a single face of the cubed-sphere. + + stretch_factor : int + Specifies the stretching factor. Default value: None + + target_lon, target_lat : float + Specifies the longitude and latitude at the center of the + cubed-sphere grid face that will be stretched. + Default values: None, None + """ if stretch_factor is not None: - cs_temp, ignore = gcpy.make_grid_SG(n_cs,stretch_factor,target_lon,target_lat) + cs_temp, ignore = gcpy.make_grid_SG( + n_cs, + stretch_factor, + target_lon, + target_lat + ) else: cs_temp = gcpy.csgrid_GMAO(n_cs) - return xr.Dataset({'nf': (['nf'],np.array(range(6))), - 'Ydim': (['Ydim'],np.array(range(n_cs))), - 'Xdim': (['Xdim'],np.array(range(n_cs))), - 'Ydim_b': (['Ydim_b'],np.array(range(n_cs+1))), - 'Xdim_b': (['Xdim_b'],np.array(range(n_cs+1))), - 'lat': (['nf','Ydim','Xdim'], cs_temp['lat']), - 'lon': (['nf','Ydim','Xdim'], cs_temp['lon']), - 'lat_b': (['nf','Ydim_b','Xdim_b'], cs_temp['lat_b']), - 'lon_b': (['nf','Ydim_b','Xdim_b'], cs_temp['lon_b']), - 'area': (['nf','Ydim','Xdim'], grid_area(cs_temp))}) - -def corners_to_xy(xc, yc): - """ Creates xy coordinates for each grid-box. The shape is (n, n, 5) where n is the cubed-sphere size. + + return xr.Dataset( + {'nf': (['nf'], np.array(range(6))), + 'Ydim': (['Ydim'], np.array(range(n_cs))), + 'Xdim': (['Xdim'], np.array(range(n_cs))), + 'Ydim_b': (['Ydim_b'], np.array(range(n_cs+1))), + 'Xdim_b': (['Xdim_b'], np.array(range(n_cs+1))), + 'lat': (['nf','Ydim','Xdim'], cs_temp['lat']), + 'lon': (['nf','Ydim','Xdim'], cs_temp['lon']), + 'lat_b': (['nf','Ydim_b','Xdim_b'], cs_temp['lat_b']), + 'lon_b': (['nf','Ydim_b','Xdim_b'], cs_temp['lon_b']), + 'area': (['nf','Ydim','Xdim'], grid_area(cs_temp)) + } + ) + + +def corners_to_xy( + x_c, + y_c +): + """ + Creates xy coordinates for each grid-box Developed, tested, and supplied by Liam Bindle. - :param xc: grid-box corner longitudes; shape (n+1, n+1) - :param yc: grid-box corner latitudes; shape (n+1, n+1) - :return: grid-box xy coordinates + Args: + ----- + x_c : numpy.ndarray + Grid-box corner longitudes; array shape = (n+1, n+1), + where n is the cubed-sphere grid size. + + y_c : numpy.ndarray + Grid-box corner longitudes; array shape = (n+1, n+1), + where n is the cubed-sphere grid size. + + Returns: + -------- + x_y : numpy.ndarray + Grid-box cartesian coordinates; array shape = (n, n, 5), + where n is the cubed-sphere grid size. """ - p0 = slice(0, -1) - p1 = slice(1, None) - boxes_x = np.moveaxis(np.array([xc[p0, p0], xc[p1, p0], xc[p1, p1], xc[p0, p1], xc[p0, p0]]), 0, -1) - boxes_y = np.moveaxis(np.array([yc[p0, p0], yc[p1, p0], yc[p1, p1], yc[p0, p1], yc[p0, p0]]), 0, -1) - return np.moveaxis(np.array([boxes_x, boxes_y]), 0, -1) + p_0 = slice(0, -1) + p_1 = slice(1, None) + boxes_x = np.moveaxis( + np.array( + [ + x_c[p_0, p_0], + x_c[p_1, p_0], + x_c[p_1, p_1], + x_c[p_0, p_1], + x_c[p_0, p_0] + ] + ), + 0, -1 + ) + boxes_y = np.moveaxis( + np.array( + [ + y_c[p_0, p_0], + y_c[p_1, p_0], + y_c[p_1, p_1], + y_c[p_0, p_1], + y_c[p_0, p_0] + ] + ), + 0, -1 + ) + return np.moveaxis( + np.array( + [boxes_x, boxes_y] + ), 0, -1 + ) -def central_angle(x0, y0, x1, y1): - """ Returns the distance (central angle) between coordinates (x0, y0) and (x1, y1). This is vectorizable. +def central_angle( + x_0, + y_0, + x_1, + y_1): + """ + Returns the distance (central angle) between cartesian + coordinates (x_0, y_0) and (x_1, y_1). This is vectorizable. Developed, tested, and supplied by Liam Bindle. - :param x0: pt0's longitude (degrees) - :param y0: pt0's latitude (degrees) - :param x1: pt1's longitude (degrees) - :param y1: pt1's latitude (degrees) - :return: Distance (degrees) + Args: + ----- + x_0, y_0 : float + Longitude and latitude (degrees) of coordinates (x_0, y_0). + + x_1, y_1: float + Longitude and latitude (degrees) of coordinates (x_1, y_1). + + Returns: + -------- + distance : float + Distance (degrees) between (x_0, y_0) and (x_1, y_1). + """ + x_0 = x_0 * DEG_TO_RAD + x_1 = x_1 * DEG_TO_RAD + y_0 = y_0 * DEG_TO_RAD + y_1 = y_1 * DEG_TO_RAD + + return np.arccos( + np.sin(y_0) * np.sin(y_1) + \ + np.cos(y_0) * np.cos(y_1) * \ + np.cos(np.abs(x_0 - x_1)) + ) * RAD_TO_DEG + + +def find_index_single( + lat, + lon, + x_centers_flat, + y_centers_flat, + xy_polygon_defs, + cs_size, + latlon_crs, + jitter_size=0.0 +): + """ + Returns the cubed-sphere grid box corresponding to a given + latitude and longitude. Called by routine find_index. + + Args: + ----- + lat, lon : float or list(float) + Latitude and longitude (degrees) of the point for which + cubed-sphere grid indices are desired. + + x_centers_flat, y_centers_flat : float or list(float) + Flattened (i.e. in Fortran column-major notation) arrays + of cubed-sphere xDim and yDim values (degrees). + + xy_polygon_defs : float or list(float) + XY polygon definitions for cubed-sphere grid boxes + (i.e. the output of function corners_to_xy). + + cs_size : int or list(int) + Cubed-sphere grid size (i.e. the number of points along + a face edge). + + latlon_crs : ? + Ouptut of pyproj.Proj("+proj=latlon") + + jitter_size : float + ?? + + Returns: + -------- + nf_cs, xdim_cs, ydim_cs : int, list(float), list(float) + nf_cs is the number of cube-sphere face + xdim_cs (aka XDim) and ydim_cs (aka YDim) are the longitude + and latitude arrays for each cell of the cubed-sphere grid. """ - RAD2DEG = 180 / np.pi - DEG2RAD = np.pi / 180 - x0 = x0 * DEG2RAD - x1 = x1 * DEG2RAD - y0 = y0 * DEG2RAD - y1 = y1 * DEG2RAD - return np.arccos(np.sin(y0) * np.sin(y1) + np.cos(y0) * np.cos(y1) * np.cos(np.abs(x0-x1))) * RAD2DEG - -def find_index_single(lat,lon,x_centers_flat,y_centers_flat,xy,cs_size,latlon_crs,jitter_size=0.0): - #import pyproj - #import shapely.ops - #import shapely.geometry # Center on x_find, y_find x_find = lon y_find = lat gnomonic_crs = pyproj.Proj(f'+proj=gnom +lat_0={y_find} +lon_0={x_find}') # Generate all distances - distances = central_angle(x_find, y_find, x_centers_flat, y_centers_flat) + distances = central_angle( + x_find, + y_find, + x_centers_flat, + y_centers_flat + ) four_nearest_indexes = np.argpartition(distances, 4)[:4] # Unravel 4 smallest indexes - four_nearest_indexes = np.unravel_index(four_nearest_indexes, (6, cs_size, cs_size)) - four_nearest_xy = xy[four_nearest_indexes] - four_nearest_polygons = [shapely.geometry.Polygon(polygon_xy) for polygon_xy in four_nearest_xy] + four_nearest_indexes = np.unravel_index( + four_nearest_indexes, + (6, cs_size, cs_size) + ) + four_nearest_xy = xy_polygon_defs[four_nearest_indexes] + four_nearest_polygons = [ + shapely.geometry.Polygon(polygon_xy) for polygon_xy in four_nearest_xy + ] # Transform to gnomonic projection - gno_transform = pyproj.Transformer.from_proj(latlon_crs, gnomonic_crs, always_xy=True).transform - four_nearest_polygons_gno = [shapely.ops.transform(gno_transform, polygon) for polygon in four_nearest_polygons] + gno_transform = pyproj.Transformer.from_proj( + latlon_crs, + gnomonic_crs, + always_xy=True + ).transform + four_nearest_polygons_gno = [ + shapely.ops.transform(gno_transform, polygon) \ + for polygon in four_nearest_polygons + ] # Figure out which polygon contains the point - Xy_find = shapely.geometry.Point(x_find, y_find) - Xy_find_GNO = shapely.ops.transform(gno_transform, Xy_find) - polygon_contains_point = [polygon.contains(Xy_find_GNO) for polygon in four_nearest_polygons_gno] + xy_find = shapely.geometry.Point(x_find, y_find) + xy_find_gno = shapely.ops.transform(gno_transform, xy_find) + polygon_contains_point = [ + polygon.contains(xy_find_gno) for polygon in four_nearest_polygons_gno + ] if np.count_nonzero(polygon_contains_point) == 0: - if jitter_size>0.0: + if jitter_size > 0.0: # Move longitude by ~1 m - nf, YDim, XDim = find_index_single(y_find,x_find+jitter_size,x_centers_flat,y_centers_flat,xy,cs_size,latlon_crs,jitter_size=0.0) + nf_cs, ydim_cs, xdim_cs = find_index_single( + y_find, + x_find+jitter_size, + x_centers_flat, + y_centers_flat, + xy_polygon_defs, + cs_size, + latlon_crs, + jitter_size=0.0 + ) else: - raise ValueError(f'Point at {x_find:8.2f} E, {y_find:8.2f} N could not be matched') + msg = f'Point at {x_find:8.2f} E, {y_find:8.2f} N ' + msg+= 'could not be matched' + raise ValueError(msg) + # The first will be selected, if more than one polygon_with_point = np.argmax(polygon_contains_point) # Get original index - nf = four_nearest_indexes[0][polygon_with_point] - YDim = four_nearest_indexes[1][polygon_with_point] - XDim = four_nearest_indexes[2][polygon_with_point] + nf_cs = four_nearest_indexes[0][polygon_with_point] + ydim_cs = four_nearest_indexes[1][polygon_with_point] + xdim_cs = four_nearest_indexes[2][polygon_with_point] - return nf, YDim, XDim + return nf_cs, ydim_cs, xdim_cs -def find_index(lat,lon,grid,jitter_size=0.0): - # For point-finding - #import pyproj - #import shapely.ops - #import shapely.geometry - assert find_point_OK, "Cannot perform index finding - need pyproj and shapely" - # Based on a routine developed, tested, and supplied by Liam Bindle. +def find_index( + lat, + lon, + grid, + jitter_size=0.0 +): + """ + Returns the cubed-sphere grid box indices corresponding to + given latitude and longitude coordinates. + + Based on a routine developed, tested, and supplied by Liam Bindle. + + Args: + ----- + lat, lon : float + Latitude and longitude (degrees) of the point for which + cubed-sphere indices are desired. + + grid : dict + Cubed-sphere grid definition as a dict of: + {'lat' : lat midpoints, + 'lon' : lon midpoints, + 'lat_b' : lat edges, + 'lon_b' : lon edges} + where each value has an extra face dimension of length 6. + + + Returns: + -------- + ind : numpy.ndarray + Array containing (nf, YDim, XDim), where: + nf is the number of cubed-sphere faces (= 6) + YDim is the cubed-sphere longitude index at (lat, lon) + XDim is the cubed-sphere latitude index at (lat, lon) + """ lon_vec = np.asarray(lon) lat_vec = np.asarray(lat) n_find = lon_vec.size @@ -291,52 +611,30 @@ def find_index(lat,lon,grid,jitter_size=0.0): # Generate everything that will be reused # Get XY polygon definitions for grid boxes - xy = np.zeros((6, cs_size, cs_size, 5, 2)) # 5 (x,y) points defining polygon corners (first and last are same) - for nf in range(6): - xy[nf, ...] = corners_to_xy(xc=x_corners[nf, :, :], yc=y_corners[nf, :, :]) + # 5 (x,y) points defining polygon corners (first and last are same) + xy_polygon_defs = np.zeros((6, cs_size, cs_size, 5, 2)) + for nf_cs in range(6): + xy_polygon_defs[nf_cs, ...] = corners_to_xy( + x_c=x_corners[nf_cs, :, :], + y_c=y_corners[nf_cs, :, :] + ) latlon_crs = pyproj.Proj("+proj=latlon") # Find 4 shortest distances to (x_find, y_find) - #idx = np.full((3,n_find),np.int(0)) # np.int no longer in numpy (ewl) - idx = np.full((3,n_find),0) - for x_find, y_find, i_find in zip(np.nditer(lon_vec),np.nditer(lat_vec),list(range(n_find))): - ## Center on x_find, y_find - #gnomonic_crs = pyproj.Proj(f'+proj=gnom +lat_0={y_find} +lon_0={x_find}') - - ## Generate all distances - #distances = central_angle(x_find, y_find, x_centers_flat, y_centers_flat) - #four_nearest_indexes = np.argpartition(distances, 4)[:4] - - ## Unravel 4 smallest indexes - #four_nearest_indexes = np.unravel_index(four_nearest_indexes, (6, cs_size, cs_size)) - #four_nearest_xy = xy[four_nearest_indexes] - #four_nearest_polygons = [shapely.geometry.Polygon(polygon_xy) for polygon_xy in four_nearest_xy] - - ## Transform to gnomonic projection - #gno_transform = pyproj.Transformer.from_proj(latlon_crs, gnomonic_crs, always_xy=True).transform - #four_nearest_polygons_gno = [shapely.ops.transform(gno_transform, polygon) for polygon in four_nearest_polygons] - - ## Figure out which polygon contains the point - #Xy_find = shapely.geometry.Point(x_find, y_find) - #Xy_find_GNO = shapely.ops.transform(gno_transform, Xy_find) - #polygon_contains_point = [polygon.contains(Xy_find_GNO) for polygon in four_nearest_polygons_gno] - - ##assert np.count_nonzero(polygon_contains_point) == 1 - ##assert np.count_nonzero(polygon_contains_point) > 0 - #if np.count_nonzero(polygon_contains_point) == 0: - # if allow_jitter: - # idx[:,i_find] = - # raise ValueError(f'Point at {x_find:8.2f} E, {y_find:8.2f} N could not be matched') - ## The first will be selected, if more than one - #polygon_with_point = np.argmax(polygon_contains_point) - - ## Get original index - #nf = four_nearest_indexes[0][polygon_with_point] - #YDim= four_nearest_indexes[1][polygon_with_point] - #XDim= four_nearest_indexes[2][polygon_with_point] - nf, YDim, XDim = find_index_single(y_find,x_find,x_centers_flat,y_centers_flat, - xy,cs_size,latlon_crs,jitter_size=jitter_size) - idx[:,i_find] = [nf,YDim,XDim] - return idx - + idx = np.full((3,n_find), 0) + for x_find, y_find, i_find in \ + zip(np.nditer(lon_vec), np.nditer(lat_vec), list(range(n_find))): + + nf_cs, ydim_cs, xdim_cs = find_index_single( + y_find, + x_find, + x_centers_flat, + y_centers_flat, + xy_polygon_defs, + cs_size, + latlon_crs, + jitter_size=jitter_size + ) + idx[:,i_find] = [nf_cs, ydim_cs, xdim_cs] + return idx From 0fbe4842dcca63c15863b1d40ad8f05337801b0a Mon Sep 17 00:00:00 2001 From: Bob Yantosca Date: Tue, 25 Jul 2023 16:26:51 -0400 Subject: [PATCH 4/5] Updated CHANGELOG w/ info about cubed-sphere tools updates CHANGELOG.py - Added sentence about adding gcpy/cstools.py Signed-off-by: Bob Yantosca --- CHANGELOG.md | 1 + 1 file changed, 1 insertion(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index 9ce3805a..f6ccedb8 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -21,6 +21,7 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/), - Added `multi_index_lat` keyword to `reshape_MAPL_CS` function in `gcpy/util.py` - Added FURA to `emission_species.yml` and `benchmark_categories.yml` - Added new routine `format_number_for_table` in `util.py` +- Added module `gcpy/cstools.py` with utility functions for cubed-sphere grids ### Changed - Simplified the Github issues templates into two options: `new-feature-or-discussion.md` and `question-issue.md` From 36196440e4ee43407a89446ed3d7c33f72ecffe3 Mon Sep 17 00:00:00 2001 From: Bob Yantosca Date: Tue, 25 Jul 2023 16:31:26 -0400 Subject: [PATCH 5/5] Update CHANGELOG.md with info about pyproj install CHANGELOG.md - Add sentence about the environment.yml file being updated to install pyproj version 3.6.0 via pip. Signed-off-by: Bob Yantosca --- CHANGELOG.md | 1 + 1 file changed, 1 insertion(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index f6ccedb8..3b907240 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -36,6 +36,7 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/), - Routine `print_totals` now prints small and/or large numbers in scientific notation - Truncate names in benchmark & emissions tables to improve readability - Add TransportTracers species names to `gcpy/emissions_*.yml` files +- Updated `docs/environment_files/environment.yml` to install `pyproj==3.6.0` via pip ### Fixed - Generalized test for GCHP or GCClassic restart file in `regrid_restart_file.py`