gcpy.cstools

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.

Examples

Generate a C24 grid and find the indices for a given lat/lon:

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 as 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()
# 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

Functions

central_angle(x_0, y_0, x_1, y_1)

Returns the central angle (great-circle distance) between two points on the sphere.

corners_to_xy(x_c, y_c)

Creates XY coordinates for each grid box from corner arrays.

extract_grid(data)

Extracts grid information from an xarray Dataset and returns it as a cubed-sphere xarray Dataset.

face_area(lon_b, lat_b[, r_sphere])

Calculates the area of cubed-sphere grid cells on one face.

find_index(lat, lon, grid[, jitter_size])

Returns the cubed-sphere grid box indices for one or more latitude/longitude coordinates.

find_index_single(lat, lon, x_centers_flat, ...)

Returns the cubed-sphere grid box indices for a single latitude/longitude point.

gen_grid(n_cs[, stretch_factor, target_lon, ...])

Returns an xarray Dataset specifying a cubed-sphere grid, optionally with stretching applied.

get_cubed_sphere_res(data)

Returns the number of grid cells along one edge of a cubed-sphere grid face.

grid_area([cs_grid, cs_res, stretch_factor, ...])

Returns the surface area in m2 for each cell of a cubed-sphere grid.

is_cubed_sphere(data)

Determines whether data is placed on a cubed-sphere grid.

is_cubed_sphere_diag_grid(data)

Determines whether a dataset has cubed-sphere History diagnostic dimensions (i.e. a dimension named "nf").

is_cubed_sphere_rst_grid(data)

Determines whether a dataset has cubed-sphere restart file dimensions (i.e. lat and lon with lat = lon * 6).

is_gchp_lev_positive_down(data)

Determines whether GCHP data levels are ordered from the top of the atmosphere downwards.

ll2xyz(lon_pt, lat_pt)

Converts a lon/lat pair in radians to Cartesian coordinates.

read_gridspec(gs_obj)

Reads a GridSpec object and returns an xarray Dataset.

sphere_angle(e_1, e_2, e_3)

Computes the angle between three points on a sphere.
gcpy.cstools.extract_grid(data)[source]

Extracts grid information from an xarray Dataset and returns it as a cubed-sphere xarray Dataset.

Parameters:

data (xarray.Dataset or xarray.DataArray) – The input dataset.

Returns:

data_cs – Cubed-sphere grid extracted from data. Returns None if the data is not on a cubed-sphere grid.

Return type:

xarray.Dataset or None

gcpy.cstools.read_gridspec(gs_obj)[source]

Reads a GridSpec object and returns an xarray Dataset.

Parameters:

gs_obj (GridSpec) – The GridSpec object to read.

Returns:

ds – The same data as an xarray Dataset object.

Return type:

xarray.Dataset

gcpy.cstools.face_area(lon_b, lat_b, r_sphere=6375000.0)[source]

Calculates the area of cubed-sphere grid cells on one face.

Inputs must be in degrees. Edge arrays must be shaped [N+1 x N+1].

Parameters:

  • lon_b (array_like) – Longitude bounds of grid cell edges in degrees.

  • lat_b (array_like) – Latitude bounds of grid cell edges in degrees.

  • r_sphere (float, optional) – Radius of the Earth in metres. Default: 6.375e6

Returns:

cs_area – Surface area in m:sup`2` for each grid box of a cubed-sphere grid face.

Return type:

numpy.ndarray

gcpy.cstools.ll2xyz(lon_pt, lat_pt)[source]

Converts a lon/lat pair in radians to Cartesian coordinates.

This function is vectorizable.

Parameters:

  • lon_pt (float) – Longitude in radians.

  • lat_pt (float) – Latitude in radians.

Returns:

coords – Cartesian vector coordinates [x_pt, y_pt, z_pt] normalised to the unit sphere.

Return type:

list of numpy.float64

gcpy.cstools.sphere_angle(e_1, e_2, e_3)[source]

Computes the angle between three points on a sphere.

Parameters:

  • e_1 (array_like) – Cartesian (x, y, z) coordinates of the midpoint.

  • e_2 (array_like) – Cartesian (x, y, z) coordinates of a point on one side of the midpoint.

  • e_3 (array_like) – Cartesian (x, y, z) coordinates of a point on the other side of the midpoint.

Returns:

angle – The spherical angle at point e_1 in radians.

Return type:

float

gcpy.cstools.grid_area(cs_grid=None, cs_res=None, stretch_factor=None, target_lon=None, target_lat=None)[source]

Returns the surface area in m2 for each cell of a cubed-sphere grid.

Parameters:

  • cs_grid (dict, optional) – Cubed-sphere grid definition containing keys: 'lat', 'lon', 'lat_b', 'lon_b', where each value has an extra face dimension of length 6. Default: None

  • cs_res (int, optional) – Cubed-sphere resolution (number of grid cells along one face edge). Default: None

  • stretch_factor (float, optional) – Stretching factor for a stretched-grid simulation. Default: None

  • target_lon (float, optional) – Longitude of the centre of the stretched grid face in degrees. Default: None

  • target_lat (float, optional) – Latitude of the centre of the stretched grid face in degrees. Default: None

Returns:

cs_area – Surface area in m2 for each cell of the cubed-sphere grid. Array shape is (6, n, n) following the GMAO convention, where n is the number of cells along a face edge.

Return type:

numpy.ndarray

gcpy.cstools.gen_grid(n_cs, stretch_factor=None, target_lon=None, target_lat=None)[source]

Returns an xarray Dataset specifying a cubed-sphere grid, optionally with stretching applied.

Parameters:

  • n_cs (int) – Number of grid boxes along a single face of the cubed-sphere.

  • stretch_factor (float, optional) – Stretching factor for a stretched-grid simulation. Default: None

  • target_lon (float, optional) – Longitude of the centre of the stretched grid face in degrees. Default: None

  • target_lat (float, optional) – Latitude of the centre of the stretched grid face in degrees. Default: None

Returns:

grid – Cubed-sphere grid definition containing variables: 'lat', 'lon', 'lat_b', 'lon_b', 'area', where each has a face dimension of length 6.

Return type:

xarray.Dataset

gcpy.cstools.corners_to_xy(x_c, y_c)[source]

Creates XY coordinates for each grid box from corner arrays.

Developed, tested, and supplied by Liam Bindle.

Parameters:

  • 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 latitudes. Array shape (n+1, n+1), where n is the cubed-sphere grid size.

Returns:

x_y – Grid-box Cartesian coordinates. Array shape (n, n, 5, 2), where n is the cubed-sphere grid size and the last two dimensions are the 5 corner points (first and last identical) and the XY pair respectively.

Return type:

numpy.ndarray

gcpy.cstools.central_angle(x_0, y_0, x_1, y_1)[source]

Returns the central angle (great-circle distance) between two points on the sphere.

This function is vectorizable. Developed, tested, and supplied by Liam Bindle.

Parameters:

  • x_0 (float) – Longitude of the first point in degrees.

  • y_0 (float) – Latitude of the first point in degrees.

  • x_1 (float) – Longitude of the second point in degrees.

  • y_1 (float) – Latitude of the second point in degrees.

Returns:

distance – Great-circle distance between (x_0, y_0) and (x_1, y_1) in degrees.

Return type:

float

gcpy.cstools.find_index_single(lat, lon, x_centers_flat, y_centers_flat, xy_polygon_defs, cs_size, latlon_crs, jitter_size=0.0)[source]

Returns the cubed-sphere grid box indices for a single latitude/longitude point.

Called internally by find_index().

Parameters:

  • lat (float) – Latitude of the query point in degrees.

  • lon (float) – Longitude of the query point in degrees.

  • x_centers_flat (numpy.ndarray) – Flattened array of cubed-sphere grid cell centre longitudes in column-major (Fortran) order.

  • y_centers_flat (numpy.ndarray) – Flattened array of cubed-sphere grid cell centre latitudes in column-major (Fortran) order.

  • xy_polygon_defs (numpy.ndarray) – XY polygon definitions for cubed-sphere grid boxes, as returned by corners_to_xy().

  • cs_size (int) – Number of grid cells along one cubed-sphere face edge.

  • latlon_crs (pyproj.Proj) – Projection object from pyproj.Proj("+proj=latlon").

  • jitter_size (float, optional) – If the point cannot be matched to a grid box, the longitude is shifted by this distance in metres before retrying. Useful for points near the poles. Default: 0.0

Returns:

  • nf_cs (int) – Index of the cubed-sphere face (0–5).

  • ydim_cs (int) – YDim (latitude) index of the matched grid box.

  • xdim_cs (int) – XDim (longitude) index of the matched grid box.

Raises:

ValueError – If the point cannot be matched to any grid box and jitter_size is 0.0.

gcpy.cstools.find_index(lat, lon, grid, jitter_size=0.0)[source]

Returns the cubed-sphere grid box indices for one or more latitude/longitude coordinates.

Based on a routine developed, tested, and supplied by Liam Bindle.

Parameters:

  • lat (float or array_like) – Latitude(s) of the query point(s) in degrees.

  • lon (float or array_like) – Longitude(s) of the query point(s) in degrees.

  • grid (xarray.Dataset) – Cubed-sphere grid definition containing variables: 'lat', 'lon', 'lat_b', 'lon_b', where each has a face dimension of length 6.

  • jitter_size (float, optional) – If a point cannot be matched to a grid box, the longitude is shifted by this distance in metres before retrying. Useful for points near the poles. Default: 0.0

Returns:

ind – Array of shape (3, N) containing (nf, YDim, XDim) for each of the N query points, where:

  • nf — cubed-sphere face index (0–5)

  • YDim — latitude index on the matched face

  • XDim — longitude index on the matched face

Return type:

numpy.ndarray

gcpy.cstools.is_cubed_sphere(data)[source]

Determines whether data is placed on a cubed-sphere grid.

Parameters:

data (xarray.Dataset or xarray.DataArray) – The input data to be tested.

Returns:

True if the data is on a cubed-sphere grid, False otherwise.

Return type:

bool

Notes

A cubed-sphere data file satisfies at least one of:

  1. Has a dimension named "nf" (GCHP/GEOS diagnostic files).

  2. The lat/lon size ratio is exactly 6 (GCHP/GEOS checkpoints).

gcpy.cstools.is_cubed_sphere_diag_grid(data)[source]

Determines whether a dataset has cubed-sphere History diagnostic dimensions (i.e. a dimension named "nf").

Parameters:

data (xarray.DataArray or xarray.Dataset) – The input data.

Returns:

True if the grid has History diagnostic dimensions (contains an "nf" dimension), False otherwise.

Return type:

bool

gcpy.cstools.is_cubed_sphere_rst_grid(data)[source]

Determines whether a dataset has cubed-sphere restart file dimensions (i.e. lat and lon with lat = lon * 6).

Parameters:

data (xarray.DataArray or xarray.Dataset) – The input data.

Returns:

True if the grid has restart file dimensions, False otherwise.

Return type:

bool

Notes

If the "lat" coordinate is present, the check is performed by comparing the lengths of lat and lon. Otherwise, variable names are searched for the prefix "SPC_".

gcpy.cstools.get_cubed_sphere_res(data)[source]

Returns the number of grid cells along one edge of a cubed-sphere grid face.

For example, returns 24 for a C24 grid (which has 24 × 24 grid cells per face).

Parameters:

data (xarray.DataArray or xarray.Dataset) – The input data.

Returns:

cs_res – The cubed-sphere resolution. Returns 0 if the data is not on a cubed-sphere grid.

Return type:

int

Notes

dims is a dict in Dataset objects but a tuple in DataArray objects. Returning the length of the corresponding coords array works in both cases.

gcpy.cstools.is_gchp_lev_positive_down(data)[source]

Determines whether GCHP data levels are ordered from the top of the atmosphere downwards.

The vertical ordering follows these conventions:

  • Checkpoint files: lev:positive = "down"

  • Emissions collection: lev:positive = "down"

  • All other collections: lev:positive = "up"

Parameters:

data (xarray.DataArray or xarray.Dataset) – The input data.

Returns:

True if levels are ordered top-of-atmosphere downwards, False if ordered surface upwards.

Return type:

bool