Grid tools

NEDAS.grid.grid module

This module contain classes to handle grid with various geometries and spatial structures.

Classes

• Grid: A factory class for 2D grids.

• GridBase: Base class for handling 2D grids.

• RegularGrid: Class for handling 2D regular grid.

• IrregularGrid: Class for handling 2D irregular grid (triangular mesh).

• Grid1D: Class for handling 1D grid.

Map projections

The projection is typically handled by a pyproj.Proj object, see list of supported projections. However, if you are working with some unsupported grid geometries, you can also define a custom proj function. What Grid needs from proj() is the conversion between longitude,latitude and the x,y coordinates, and the inverse, i.e.

>>> x, y = proj(lon, lat)
>>> lon, lat = proj(x, y, inverse=True)

Coordinates

The coordinates x,y has the same dimensions as the 2D grid.

For a regular grid, when we have a scalar field indexed by fld[j,i], its coordinate is x[j,i], y[j,i]. Note that for Python the rightmost dimension is the quickest (row index first), if you have data that is column index first fld[i,j] then you will need to transpose it before using Grid.

So far we only support regular grid x,y where all the x[j,:] and all the y[:,i] are identical, namely x, y = np.meshgrid(xc, yc) where xc,yc are the 1D coordinates arrays. If you have a rotated x,y plane, you need to let the projection take care of the rotation first.

Cyclic boundary condition is supported, by setting the cyclic_dim parameter you can tell Grid which dimension(s) should wrap around when the coordinates fall outside of the given range (done by wrap_cyclic_xy())

The longitude,latitude (geodetic) grid is special since it is the only projection not on a flat surface. For this +proj=longlat projection we support longitude conventions both 0 ~ 360 and -180 ~ 180 as long as its unit is “degrees east”, you shall set cyclic_dim=’x’ if x is your longitude dimension. The latitude has units “degrees north”. We also support special treatment of the two poles (lat=-90, 90) at which rotation angle is not well defined. You shall specify which dimension contains the pole and which index in that dimension is the pole. For example, if your data is [lat=40~90, long=0~360], you can set pole_dim=’y’ and pole_index=(-1).

For an unstructured mesh, the x,y coordinates are for the nodes (vertices of the triangles), the 2D field can be either defined on nodes or elements (triangle faces). If the later is true, x_elem,y_elem will be used in handling the field. Note that x_elem has different size with x.

Examples

Define a proj function, e.g. stereographic:

>>> import pyproj
>>> proj = pyproj.Proj('+proj=stere +lat_0=90 +lon_0=0')

Create a regular mesh x, y coordinates, with 10 km resolution, in meter units:

>>> import numpy as np
>>> x, y = np.meshgrid(np.arange(-3e6, 3e6, 1e4), np.arange(-3e6, 3e6, 1e4))

Create a grid instance:

>>> from NEDAS.grid import Grid
>>> grid = Grid(proj, x, y)

self.nx, self.ny are the grid dimension for regular grids.

self.xmin, self.xmax, self.ymin, self.ymax define the grid extent (bounding box) for regular grids.

self.Lx, self.Ly are the grid length in x and y directions.

self.dx is the grid spacing (resolution), for unstructured mesh it is the averaged edge length of the triangles. self.dy is also provided for regular grid, although it is typical that we use dy`=`dx.

self.mfx, self.mfy are the map factors in x and y directions. Due to map projection the actual length is distorted on the 2D grid plane (see [Tissot’s indicatrix](https://en.wikipedia.org/wiki/Tissot%27s_indicatrix) for an illustration). These map factors are useful in computing actual distances on Earth, using self.dx / self.mfx and self.dy / self.mfy

self.tri for unstructured mesh this is the triangle indices, computed by matplotlib.tri.Triangulation

self.x_elem, self.y_elem are the center coordinates for each triangular element.

Typically we use the Grid class for conversion of a 2D field from one grid to another, which involves interpolation, rotation of vectors and coarse-graining if the two grid has different resolution. Additionally, we also provide basic tools for visualizing a field on the map. See tutorials/grid_convert for some illustrated examples.

Once you created a grid1 object on which a 2D field fld1 is defined, and a grid2 object for the destination grid. You convert fld1 by first setting

grid1.set_destination_grid(grid2), then

fld2 = grid1.convert(fld1)

It is recommended that you initialize a grid object, set the destination grid by set_destination_grid, then apply the same grid.convert to many fields needing conversion. Since the interpolation and rotation matrices are computed by set_destination_grid only once, the following application of grid.convert can be very efficient.

### 2.1 Interpolation

Two methods are supported for now. method=linear does bi-linear interpolation for the regular grids and barycentric interpolation for the triangular meshes. method=nearest performs nearest neighbor interpolation.

If you just want to do interpolation, you can also call grid.interp(fld, x, y, method) directly, where fld is the original field, x,y are the target coordinates (can be just one number, or an array, multi-dimension arrays will be flattened), method is the interpolation method.

If you only want the indices of given x,y coordinates, you can also use grid.find_index(x, y). But note that we return more than just the indices, also some other information. Five things will be returned by find_index:

• inside: True if the point x,y is inside the grid

• indices: None for regular grid; triangle index which the x,y position falls in.

• vertices: the indices of rectangle(4) or triangle(3) vertices the x,y position falls in.

• in_coords: the internal coordinates of the x,y point within the rectangle (in_x, in_y) or triangle (in1, in2, in3)

• nearest: index in x.flatten(), y.flatten() for the nearest grid point to the given x,y position.

### 2.2 Rotation of Vectors

We support both scalar and vector fields. For a vector field, its u- and v- components are stacked in the first dimension, so that vec_fld has dimension [2, ny, nx]. Parameter is_vector = True will let grid.convert know it’s operating on a vector field and need to rotate vectors if projection changes.

The rotation matrix is computed by displacing a small amount in x,y directions on each grid point, and project both the original point and displaced point to the new grid, then measuring the displacement and figure out the rotation.

If you just want to rotate vectors, you can apply grid.rotate_vectors(vec_fld).

### 2.3 Coarse-graining

Setting coarse_grain=True in grid.convert will perform an additional coarse graining (by grid.coarsen) during interpolation. When the destination grid has lower resolution than the source grid, more points will fall in the same destination grid element. The coarse-graining will average them to represent the field on the low-resolution grid, instead of interpolating only the nearest source grid points. For a field with small-scale variability on the source grid, the interpolation to low-resolution grid will be sensitive to those sub-grid-scale variability, causing the interpolated value to have representation errors. Taking the average in the low-resolution grid will remedy this issue.

However, if you are dealing with a sparse observing network where you want to interpolate in the source grid to find the observed values, you don’t need coarse graining, because observation should pick up sub-grid scale variability in its measurement anyway. Remember to set coarse_grain=False if this is the case.

### 2.4 Visualization

To avoid intricate dependencies from installing the cartopy package, we instead provide some internal solutions for basic plotting.

ax is a matplotlib.pyplot.Axes object for plotting.

Scalar fields can be plotted using grid.plot_field(ax, fld, vmin, vmax, cmap), we use pyplot.pcolor for regular grids and pyplot.tripcolor for unstructured meshes. vmin, vmax are minimum and maximum values shown for the field. cmap is the color map.

Vector fields can be plotted using grid.plot_vectors(ax, vec_fld, V, L, spacing, num_steps, linecolor, linewidth, showred, ref_xy, refcolor, showhead, headwidth, headlength). This replaces the pyplot.quiver and provides more control over how a vector plot looks. vec_fld has a dimension [2, …] for the two components. V is the velocity scale that is the typical velocity value in this field (same units as vec_fld), L is the length scale that is the length (in x,y units) for vectors with velocity V. If you don’t provide V and L, a standard value will be used based on the data provided, which is often good enough. spacing controls the density of vectors in the plot, they are spaced with interval spacing * L. num_steps`=1 displays straight vectors (as in `quiver); >1 lets you display curved trajectories, at each sub-step the velocity is re-interpolated at the new position along the trajectories. linecolor, linewith controls the style of vector lines. show_ref=True displays the reference vector at position ref_xy=(x,y) and with background color refcolor. showhead=True displays the vector heads with headwidth, headlength controlling their styles.

To show the map overlay, you can use grid.plot_land(ax, color, linecolor, linewidth, showriver, rivercolor, showgrid, dlon, dlat). There are several shape files downloaded from [www.naturalearthdata.com](https://www.naturalearthdata.com) are stored in grid directory, where the plot_land reads and gather land_data, river_data and lake_data. They data are cached after first time use, so making several plots in a role will be more efficient. color controls the face color of the landmass, linecolor and linewidth controls the style of the coastlines. showriver=True will plot the rivers and lakes with rivercolor. showgrid=True turns on the latitude longitude grid at dlon, dlat intervals (degrees). The x,y axes will show the coordinate tick marks in its units.

class NEDAS.grid.grid.Grid(proj, x, y, regular=True, bounds=None, cyclic_dim=None, pole_dim=None, pole_index=None, distance_type='cartesian', triangles=None, neighbors=None, dst_grid=None)

Bases: object

Factory class for creating a grid object.

This class serves as a factory for creating either a RegularGrid or an IrregularGrid instance based provided parameters. It allows for easy instantiation of different grid types without needing to directly interact with the specific grid classes.

Examples

The factory class provides two classmethods to help creating a grid instance more easily. To create a regular grid with coordinates x from 0 to 100 and y from 0 to 100, resolution dx=1:

>>> from NEDAS.grid import Grid
>>> grid = Grid.regular_grid(None, 0, 100, 0, 100, 1)

To create a irregular grid with same range of coordinates, in total 2000 grid points randomly placed:

>>> grid = Grid.random_grid(None, 0, 100, 0, 100, 2000)

classmethod regular_grid(proj, xstart: float, xend: float, ystart: float, yend: float, dx: float, centered: bool = False, **kwargs) → RegularGrid

Create a regular grid within specified boundaries.

Parameters:

• proj (Callable) – Projection from lon,lat to x,y.

• xstart (float) – Lower bound for X coordinates.

• xend (float) – Upper bound for X coordinates.

• ystart (float) – Lower bound for Y coordinates.

• yend (float) – Upper bound for Y coordinates.

• dx (float) – Grid spacing.

• centered (bool) – Optional, toggle for grid points to be on vertices (False) or in the middle of each grid box (True). Default is False.

• **kwargs – Additional keyword arguments.

Returns:

A Grid object representing the regular grid.

Return type:

RegularGrid

classmethod random_grid(proj, xstart: float, xend: float, ystart: float, yend: float, npoints: int, min_dist: float | None = None, **kwargs) → IrregularGrid

Create a grid with randomly positioned points within specified boundaries.

Parameters:

• proj (pyproj.Proj) – Projection from lon,lat to x,y.

• xstart (float) – Lower bound for X coordinates.

• xend (float) – Upper bound for X coordinates.

• ystart (float) – Lower bound for Y coordinates.

• yend (float) – Upper bound for Y coordinates.

• npoints (int) – Number of grid points.

• min_dist (float) – Optional, minimal distance allowed between each pair of grid points.

• **kwargs – Additional keyword arguments.

Returns:

A Grid object representing the randomly positioned grid.

Return type:

IrregularGrid

NEDAS.grid.grid_base module

NEDAS.grid.grid_regular module

class NEDAS.grid.grid_regular.RegularGrid(proj, x, y, bounds=None, cyclic_dim=None, distance_type='cartesian', pole_dim=None, pole_index=None, neighbors=None, dst_grid=None)

Bases: Grid2DBase

Regular 2D grid class

Parameters:

• pole_dim (str, optional) – None (default), if one of the dimension has poles, 'x' or 'y'

• pole_index (str, optional) – None (default), tuple of the pole index(s) in pole_dim

• neighbors (np.ndarray, optional) – None (default), for regular grid with special geometry (e.g. tripolar ocean grid), neighbors stores the j,i index of 4 neighors (east, north, west, and south) on each grid point. Since specifying neighbors already take care of cyclic boundary conditions, cyclic_dim will be discarded if neighbors is set.

pole_dim: str | None

pole_index: list[int] | None

nx: int

ny: int

change_resolution_level(nlevel)

Generate a new grid with changed resolution.

Parameters:

nlevel (int) – Positive number, downsample grid abs(nlevel) times, each time doubling the grid spacing; Negative number, upsample grid abs(nlevel) times, each time halving the grid spacing

Returns:

A new grid object with changed resolution.

find_index(x_, y_)

Find indices of self.x, self.y corresponding to the given x_, y_.

Parameters:

• x (float or np.ndarray) – x-coordinates of target point(s).

• y (float or np.ndarray) – y-coordinates of target point(s).

Outputs:

inside (np.ndarray of bool): Boolean array of shape (x_.size,) indicating whether

each x_, y_ point lies inside the grid.

indices (np.ndarray of int or None): Indices of grid elements containing the input points.

• For regular grids, this is None since vertices suffice to locate the grid box.

• For unstructured meshes, these are indices into tri.triangles, from tri_finder.

vertices (np.ndarray of int): Array of shape (inside_size, n), where

n = 4 for regular grid boxes or n = 3 for mesh triangles. These are indices into self.x, self.y (flattened) for the vertices of the grid element that each point falls in.

in_coords (np.ndarray of float): Array of shape (inside_size, n) giving internal coordinates

of each point within the containing element. Used to compute interpolation weights.

nearest (np.ndarray of int): Array of shape (inside_size,) with indices of the grid nodes

closest to each point.

Notes

• This function assumes self.x, self.y define either a regular or triangular grid.

• Internal coordinates are used for interpolation and vary in dimension based on the grid type.

rotate_vectors(vec_fld)

Apply the rotate_matrix to a vector field

Parameters:

vec_fld (np.array) – The input vector field, shape (2, self.x.shape)

Returns:

The vector field rotated to the dst_grid.

get_corners(fld)

given fld defined on a regular grid, obtain its value on the 4 corners/vertices

interp(fld, x=None, y=None, method='linear')

Interpolation of 2D field data (fld) from one grid (self or given x,y) to another (dst_grid). This can be used for grid refining (low->high resolution) or grid thinning (high->low resolution). This also converts between different grid geometries.

Parameters:

• fld (np.array) – Input field defined on self, should have same shape as self.x

• x (float or np.array) – Optional; If x,y are specified, the function computes the weights and apply them to fld If x,y are None, the self.dst_grid.x,y are used. Since their interp_weights are precalculated by dst_grid.setter it will be efficient to run interp for many different input flds quickly.

• y (float or np.array) – Optional; If x,y are specified, the function computes the weights and apply them to fld If x,y are None, the self.dst_grid.x,y are used. Since their interp_weights are precalculated by dst_grid.setter it will be efficient to run interp for many different input flds quickly.

• method (str) – Interpolation method, can be ‘nearest’ or ‘linear’

Returns:

The interpolated field defined on the destination grid

coarsen(fld)

Coarse-graining is sometimes needed when the dst_grid is at lower resolution than self. Since many points of self.x,y falls in one dst_grid box/element, it is better to average them to represent the field on the low-res grid, instead of interpolating only from the nearest points that will cause representation errors.

Parameters:

fld (np.array) – Input field to perform coarse-graining on, it is defined on self.

Returns:

The coarse-grained field defined on self.dst_grid.

plot_field(ax, fld, vmin=None, vmax=None, cmap='viridis', **kwargs)

Plot a scalar field using pcolor/tripcolor

Parameters:

• ax (matplotlib.pyplot.Axes) – Axes handle for plotting

• fld (np.array) – The scalar field for plotting

• vmin (float, optional) – The minimum and maximum value range for the colormap, if not specified (None) the np.min, np.max of the input fld will be used.

• vmax (float, optional) – The minimum and maximum value range for the colormap, if not specified (None) the np.min, np.max of the input fld will be used.

• cmap (matplotlib colormap, or str, optional) – Colormap used in the plot, default is ‘viridis’

NEDAS.grid.grid_irregular module

class NEDAS.grid.grid_irregular.IrregularGrid(proj, x, y, bounds=None, cyclic_dim=None, distance_type='cartesian', triangles=None, dst_grid=None)

Bases: Grid2DBase

Class for handling irregular grids.

Parameters:

triangles (np.ndarray, optional) – None (default), if the triangle indices ind[grid.x.size, 3] for the unstructured mesh are given here, the grid will skip computing triangulation.

change_resolution_level(nlevel)

find_index(x_, y_)

Find indices of self.x, self.y corresponding to the given x_, y_.

Parameters:

• x (float or np.ndarray) – x-coordinates of target point(s).

• y (float or np.ndarray) – y-coordinates of target point(s).

Outputs:

inside (np.ndarray of bool): Boolean array of shape (x_.size,) indicating whether

each x_, y_ point lies inside the grid.

indices (np.ndarray of int or None): Indices of grid elements containing the input points.

• For regular grids, this is None since vertices suffice to locate the grid box.

• For unstructured meshes, these are indices into tri.triangles, from tri_finder.

vertices (np.ndarray of int): Array of shape (inside_size, n), where

n = 4 for regular grid boxes or n = 3 for mesh triangles. These are indices into self.x, self.y (flattened) for the vertices of the grid element that each point falls in.

in_coords (np.ndarray of float): Array of shape (inside_size, n) giving internal coordinates

of each point within the containing element. Used to compute interpolation weights.

nearest (np.ndarray of int): Array of shape (inside_size,) with indices of the grid nodes

closest to each point.

Notes

• This function assumes self.x, self.y define either a regular or triangular grid.

• Internal coordinates are used for interpolation and vary in dimension based on the grid type.

interp(fld, x=None, y=None, method='linear')

Interpolation of 2D field data (fld) from one grid (self or given x,y) to another (dst_grid). This can be used for grid refining (low->high resolution) or grid thinning (high->low resolution). This also converts between different grid geometries.

Parameters:

• fld (np.array) – Input field defined on self, should have same shape as self.x

• x (float or np.array) – Optional; If x,y are specified, the function computes the weights and apply them to fld If x,y are None, the self.dst_grid.x,y are used. Since their interp_weights are precalculated by dst_grid.setter it will be efficient to run interp for many different input flds quickly.

• y (float or np.array) – Optional; If x,y are specified, the function computes the weights and apply them to fld If x,y are None, the self.dst_grid.x,y are used. Since their interp_weights are precalculated by dst_grid.setter it will be efficient to run interp for many different input flds quickly.

• method (str) – Interpolation method, can be ‘nearest’ or ‘linear’

Returns:

The interpolated field defined on the destination grid

coarsen(fld)

Coarse-graining is sometimes needed when the dst_grid is at lower resolution than self. Since many points of self.x,y falls in one dst_grid box/element, it is better to average them to represent the field on the low-res grid, instead of interpolating only from the nearest points that will cause representation errors.

Parameters:

fld (np.array) – Input field to perform coarse-graining on, it is defined on self.

Returns:

The coarse-grained field defined on self.dst_grid.

plot_field(ax, fld, vmin=None, vmax=None, cmap='viridis', **kwargs)

Plot a scalar field using pcolor/tripcolor

Parameters:

• ax (matplotlib.pyplot.Axes) – Axes handle for plotting

• fld (np.array) – The scalar field for plotting

• vmin (float, optional) – The minimum and maximum value range for the colormap, if not specified (None) the np.min, np.max of the input fld will be used.

• vmax (float, optional) – The minimum and maximum value range for the colormap, if not specified (None) the np.min, np.max of the input fld will be used.

• cmap (matplotlib colormap, or str, optional) – Colormap used in the plot, default is ‘viridis’

NEDAS.grid.grid_1d module

class NEDAS.grid.grid_1d.Grid1D(x: ndarray, bounds=None, regular=True, cyclic=False, distance_type='cartesian', dst_grid=None)

Bases: object

Grid class to handle fields defined on a 1D grid

Methods convert, interp and distance has the same input args as their counterparts in the Grid class. This introduces y coordinates in Grid1D class, although y is not defined for 1D grid, this makes the code that calles Grid/Grid1D methods to be easier to maintain.

mask: ndarray

x: ndarray

regular: bool

cyclic: bool

nx: int

dx: float

Lx: float

classmethod regular_grid(xstart, xend, dx, centered=False, **kwargs)

classmethod random_grid(xstart, xend, npoints, **kwargs)

change_resolution_level(nlevel)

property dst_grid

set_destination_grid(grid)

wrap_cyclic(x_)

find_index(x_)

interp(fld, x=None, y=None, method='linear')

coarsen(fld)

convert(fld, is_vector=False, method='linear', coarse_grain=False, **kwargs)

distance(ref_x, x, ref_y=None, y=None, p=1)

plot_field(ax, fld, vmin, vmax)

plot_vectors(ax, vec_fld)