Reduced Gaussian Grids

Reduced Gaussian Grids#

Interpolation of the reduced Gaussian grids of the ECMWF analyses and forecasts (reduced_gg in GRIB, the N320, N640, O1280 grids): a fixed set of parallels, the Gaussian latitudes, each holding equally spaced points whose number decreases towards the poles. The grids are interpolated in 2D (a value per point), in 3D (a column of model levels per point) and in time, from numeric inputs only: the latitudes of the parallels, their number of points and the values. Reading the files, GRIB or otherwise, is left to the caller; Reduced Gaussian Grids gives a complete recipe with ecCodes.

Mesh, zone, grid and series#

  • A Mesh describes the geometry and never holds values. Its points are numbered parallel after parallel from the north, and on each parallel eastwards from longitude 0, as a global GRIB field stores them.

  • A Zone is the subset of the points a regional grid holds: sorted, disjoint ranges [start, stop) of the flat indices of the mesh. It is computed from the mesh alone, before any value is read, from a box (from_box()) or from the points to interpolate (from_points(), exactly the corners they read). A parallel may hold several ranges: a satellite pass crosses most parallels twice, on opposite sides of the globe. Read the values of the zone with its indices() (fancy indexing of an array, Zarr or netCDF orthogonal selection), or with its ranges() for readers that seek; coalesce() merges close ranges for fewer, larger reads without changing the zone.

  • A Grid holds the fields of one time step on the mesh or on a zone: float32 or float64 arrays of shape (points,) (2D) or (points, levels) (3D), C-contiguous, read without copy, and optionally the hybrid coefficients a and b of the half levels.

  • A Series holds grids in time on one mesh. Each step may hold its own zone: time_zones() computes them, each from the points whose date reads the step.

A point is located once on the mesh. Its cell, the four points around it on the two parallels that bracket its latitude, then serves every field, every level and both time steps it reads. The cell is a trapezoid, its two parallels holding different numbers of points, and the bilinear weights account for it. A NaN value is left out and the weights of the other corners renormalized, level by level, so that a level under the orography of one corner does not void the whole column.

Time steps#

By default a date reads the two steps that bracket it, blended linearly, or the nearest one (time_method="nearest"). interval reads instead the two ends of the analysis interval the date falls in, on a lattice whose phase is the first step (00, 06, 12, 18 h for analyses; 03, 09, 15, 21 h for forecasts): an interval that lost one of its ends is read off the other alone, one that lost both is NaN, and map_availability() counts the missing ends, the flag the altimetry products publish. max_step_gap treats a wider pair as a missing step. The dates are numpy.datetime64 (handled at nanosecond resolution) or integers in a unit of your choice; the durations are then numpy.timedelta64 or integers, respectively.

Computing from the columns#

interpolate() returns the interpolated fields: a value per point for a 2D field, a column per point for a 3D field. Anything computed point by point from them with numpy, such as a dry tropospheric correction from the surface pressure, needs nothing more. apply() computes from the columns without holding them all in memory, in two ways:

  • a Python callable, called on blocks of points with the columns, the coordinates of the points and their extra arguments, written with numpy operations over the block;

  • a native kernel, a function compiled by numba with column_kernel() (pip install numba), called for each point from several threads without the GIL: for computations that loop over the levels and do not vectorize over points, such as a vertical integration. A kernel cannot raise; it returns 0, or a non-zero code that makes the point NaN.

For a series, blend="columns" (default) blends the columns of the two steps before calling the function, blend="results" calls it on each step and blends its results, which costs twice as much: it interpolates in time the result rather than the fields, which differs when the function is nonlinear in them.

NumPy functions#

Mesh(self, latitudes[, shape, order, shape, ...])

A reduced Gaussian mesh: parallels of equally spaced points, north to south, starting at longitude 0.

Zone(self, mesh, ranges[, shape, writable])

The points of a mesh a regional grid holds: sorted, disjoint ranges [start, stop) of flat indices.

Grid(self, mesh, fields, zone, a[, shape, ...])

The fields of one time step on a mesh or a zone of it.

Series(mesh)

A time series of grids on one mesh.

interpolate(source, lon, lat[, time, ...])

Interpolate the fields of a grid or a series at points.

apply(source, callback, lon, lat[, time, ...])

Compute values per point from the interpolated columns of the fields.

column_kernel([n_out])

Compile a Python function into a native column kernel with numba.

ColumnKernel(address, n_out, cfunc)

A native column kernel, called per point without the GIL.

time_zones(mesh, step_times, time, lon, lat, *)

Compute the zone each time step of a series needs.

map_availability(step_times, time, interval)

Count the steps missing from the analysis interval of each date.

coalesce(ranges, io_gap)

Merge the ranges of a zone for reading, without changing the zone.

gaussian_latitudes(n)

Compute the latitudes of a Gaussian grid.

octahedral_points(n)

Compute the number of points per parallel of an octahedral grid.

regular_points(n)

Compute the number of points per parallel of a regular Gaussian grid.

Low-level functions#

interpolate(source, lon[, shape, writable, ...])

Interpolate the fields of a reduced Gaussian grid or series at points.

apply(source, callback, lon[, shape, ...])

Compute values per point from the interpolated columns of the fields.

bracket(step_times[, shape, writable, ...])

The steps each date reads, by the rule of a Series: (i0, i1, weight), with -1 for a date that reads no step and i1 == i0 for a date that reads a single one.