Convention for reduced Gaussian grids

[LoreaGSM]

Convention for reduced Gaussian grids

Moderator

Moderator Status Review [last updated: YYYY-MM-DD]

Requirement Summary

The reduced Gaussian grid is a widely used format in atmospheric and climate data, yet it is not currently represented by any of the grid types defined in the CF Conventions. Unlike regular grids, the reduced Gaussian grid has a varying number of longitude points per latitude, decreasing toward the poles. This irregularity makes it both more efficient for spectral transforms and more complex to describe with existing CF structures. Because no existing CF grid mapping can capture these characteristics, a new convention is proposed to ensure consistent representation, efficient storage, and interoperability across software and datasets that rely on this grid type.

Technical Proposal Summary

Following discussions between NCAS, University of Reading (@JonathanGregory, @davidhassell and @sadielbartholomew) and ECMWF (@sebvi, @LoreaGSM), we propose the following approach for incorporating reduced Gaussian grids into the CF Conventions:
As the grid is not a regular but complex 2D structure, storing full latitude-longitude vectors is inefficient. Instead, the grid can be defined by:

• The latitudes, which correspond to the zeros of the Legendre polynomial of order 2N, where N is the number of latitude lines from pole to equator. These grids exclude points at the poles and equator.
• The number of longitude points per latitude as defined above, stored in a points-per-latitude vector. Longitudes can then be reconstructed using a known formula. The first longitude is always 0 for each of the latitude lines.
• The grid subtype (e.g., Octahedral, Normal) and the number of latitudes from pole to equator (the “resolution” of the grid).

The reduced Gaussian grids should also include an indexing vector (reduced_gaussian_index) in global or non-global domains.This indexing vector can be stored with a reduced size by the compression by coordinate subsampling (CF chapter 8.3). In the particular case of a subdomain or a masked domain, data is given only for the non missing points. These indexes of these grid-cells are referenced to the global index (reduced_gaussian_index) and compression by gathering (CF chapter 8.2) can be used. The pl and lat vectors are always global.

Examples

Example for the global most general case:

dimensions:
   reduced_gaussian_index = 6599680 ;
   n_lats = 2560 ;
variables:
   char reduced_gaussian ;
        reduced_gaussian:grid_mapping_name = “reduced_gaussian” ;
        reduced_gaussian:grid_subtype = "octahedral” ;
        reduced_gaussian:grid_resolution = 1280 ;
        reduced_gaussian:points_per_latitude = “pl” ;
        reduced_gaussian:latitudes = “lat” ;
   float lat(n_lats) ;
        lat:units = "degrees_north" ;
        lat:long_name = "latitude" ; 
   int pl(n_lats) ;
        pl:long_name = "number of points per latitude" ;
        pl:units = “1” ;
   int reduced_gaussian_index(reduced_gaussian_index) ;
        reduced_gaussian_index:long_name = “Reduced Gaussian grid point index”;
        reduced_gaussian_index:units = “1”;

   float data(reduced_gaussian_index) ;
        data:grid_mapping = “reduced_gaussian” ;
        data:coordinates = “reduced_gaussian_index” ;
data:
   latitude = …. ;                     #  vector of 2560 values of latitude
   points_per_latitude = … ;           #  vector of 2560 values of number of points per each latitude
   data = … ;                          #  variable with  6599680 values
   reduced_gaussian_index = ... ;      #  vector with values [0, …, 6599679]

Either points_per_latitude (pl) or the accumulated_points_per_latitude (accum_pl) vectors can be provided indistinctively, the latter being the cumulative sum of the first one. Both vectors define the grid, but the accumulated vector may ease user computations of latitudes and longitudes. The points-per-latitude vectors and the latitude both have a length of twice the grid_resolution (2560 for grid_resolution=1280).
The reduced_gaussian_index can be stored with a reduced size using the compression by coordinate subsampling (CF chapter 8.3)

Example for the regional case or a masked global (a particular case of the global example):

dimensions:
   cells = 3 ;
   n_lats = 2560 ;
variables:
   char reduced_gaussian ;
        reduced_gaussian:grid_mapping_name = “reduced_gaussian” ;
        reduced_gaussian:grid_subtype = "octahedral” ;
        reduced_gaussian:grid_resolution = 1280 ;
        reduced_gaussian:points_per_latitude = “pl” ;
        reduced_gaussian:latitudes = “lat” ;
   float lat(n_lats) ;
        lat:units = "degrees_north" ;
        lat:long_name = "latitude" ; 
   int pl(n_lats) ;
        pl:long_name = "number of points per latitude" ;
        pl:units = “1” ;
   int cells(cells);
        cells:compress = “reduced_gaussian_index”;
        cells:long_name = “indices of reduced_gaussian_index for gathered data”
   float data(cells) ;
        data:grid_mapping = “reduced_gaussian” ;

data:
   cells= 3507, 6789, 10689 ;
   latitude = …. ;                     # vector of 2560 values of latitude
   points_per_latitude = … ;           # vector of 2560 values of number of points per each latitude
   data = 10, 20, 30 ;                 # variable with 3 non missing values

In this case, only 3 points have valid data from the 6599680 of the global case, and the cells vector represents the position in the global vector. The pl and lat vectors in the grid definition are referred to as the global case.

Calculations of the latitude and longitude for a particular grid cell from the points_per_latitude vector:

For a given grid cell point defined by index i, the latitude index k can be calculated as the minimum k-th element so that the accumulation points_per_latitude vector (accum_pl) for that element is higher than i. The longitude index is then calculated for that latitude index as the remaining number of elements in the following latitude line.

$$ {accum}\_{pl}[0]= pl[0] $$

$$ {accum}\_{pl}[n] = {accum}\_{pl}[n-1] + pl[n] , \quad n \in {1, \dots, n\_{lats}-1} $$

Then:

$$ k = \min \{ s \mid {accum}\_{pl}[s] > i \} $$

$$ m = i - {accum}\_{pl}[k-1] - 1 $$

The latitude for that grid point would be the k-th element of the latitude vector, and the longitude would be the longitude index times the increment of longitude for that latitude line.

$$ lat_i = lat[k], \quad lon_i = 0^\circ + m \cdot \frac{360}{pl[k]} $$

Longitude and latitude bounds are located in the middle point between each consecutive longitude and each consecutive latitude.

Benefits

The reduced Gaussian grid is the native format for most ECMWF data products. With the upcoming ERA6 reanalysis release (expected late 2026 or early 2027), a substantial volume of data will be disseminated using this type of grid in GRIB2 but also NetCDF. Users will widely benefit from this proposal.

Status Quo

At present, the reduced Gaussian grid is not supported by CF, and a large amount of data produced in this type of grid (ERA6 and ECMWF products) would not be able to be covered by the convention.

Associated pull request

A pull request has not yet been created.

Detailed Proposal

All covered above.

[ChrisBarker-NOAA]

I'm a bit confused:

At present, the reduced Gaussian grid is not supported by CF, and a large amount of data produced in this type of grid (ERA6 and ECMWF products) would not be able to be covered by the convention.

First: -- what are the ERA6 and ECMWF products doing now?

Second: IIUC (and I may not) -- these grids could be represented by the UGRID standard, or maybe even similarly to how curvilinear grids are represented -- essentially, that the coordinates of the grid nodes are explicitly stored.

I think that this proposal is based on the idea that the locations of the grid nodes can be computed from a small amount of information, and thus it's unnecessary and inefficient to store all the coordinates explicitly.

However, for the most part, CF has been about storing everything explicitly.

For instance, in the GIS world, they often store raster data on a regular grid by specifying a corner, delta-x, delta-y and size. [https://en.wikipedia.org/wiki/Esri_grid] This completely specifies the "grid", but I don't think it's CF compliant, unless I missed something.

This proposal seems to me to be analogous.

To be fair, in the raster case, storing the full x and y vector is relatively small, so why not? whereas a full Gaussian grid could be large.

However, this is an issue with unstructured and curvilnear grids as well, and I'm not sure a huge one -- the data associated with the grid [*] is usually a lot larger than the grid specification itself anyway.

And the point of having the full grid (and coordinates) explicit is that it makes it a lot easier for client code to work with / visualize the data.

Another analogous example is for information on grids that is "on the face" or "on an edge" -- e.g. averages over a cell, or fluxes through a cell wall. In that case, the data are not at a particular point, but the CF spec still requires that coordinates be specified. And, in fact, the cells themselves need to be explicitly specified as well.

Personally, I think CF could relax some of these constraints, but they have the reasons.

[*] -- well, if you have only a variable or two, and only a few time steps, then the grid definition can be large compared to the data -- but I think the solution to that is to be able to specify the grid in a separate file, so it can be re-used rather than repeated -- which is the topic of another proposal(s). @davidhassell -- what's the status of that?

[JonathanGregory]

Dear @LoreaGSM

Thanks for your proposal. I support its inclusion as a new grid mapping.

I have some suggestions (mostly the same as when I commented on the draft):

1. You could say that the bounds in longitude for cell i are 0° + (m ± ½) × 360° / pl[k]. For latitude, since it's a one-dimensional array, I think it would be good to store it explicitly as bounds on the latitude variable.

2. I suggest that the latitude variable should have the same name as its dimension, making it a NUG coordinate variable, and that the attribute should be regarded as naming the dimension. (We can continue to debate whether thathat has any implication for the CF data model, but it's not problematic.)

3. Since the grid resolution is exactly half the latitude dimension, it is redundant, and I suggest omitting it.

If you adopted these suggestions, the regional example would look like this:

dimensions:
  cells = 3 ;
  lat = 2560 ;
  two = 2;
variables:
  char reduced_gaussian ;
    reduced_gaussian:grid_mapping_name = "reduced_gaussian" ;
    reduced_gaussian:grid_subtype = "octahedral" ;
    reduced_gaussian:points_per_latitude = "pl" ;
    reduced_gaussian:latitude_dimension = "lat" ;   //  latitude_dimension instead of latitudes
  float lat(lat) ;
    lat:units = "degrees_north" ;
    lat:long_name = "latitude" ; 
    lat:bounds="lat_bounds";
  float lat_bounds(lat,two);
  int pl(lat) ;
    pl:long_name = "number of points per latitude" ;
    pl:units = "1" ;
  int cells(cells);
    cells:compress = "reduced_gaussian_index";
    cells:long_name = "indices of reduced_gaussian_index for gathered data"
  float data(cells) ;
    data:grid_mapping = "reduced_gaussian" ;

data:
  cells= 3507, 6789, 10689 ;
  lat = ... ;                 // vector of 2560 values of latitude
  lat_bounds = ...            // array of 2560 pairs of latitude bounds
  points_per_latitude = ... ; // vector of 2560 values of number of points per each latitude
  data = 10, 20, 30 ;         // variable with 3 non missing values

Best wishes

Jonathan

[JonathanGregory]

Here are some brief responses to points of Chris's.

• This proposal is for a new grid mapping. For each grid mapping, there's a formula or algorithm for obtaining the latitude and longitude from other coordinates (often projection_x_coordinate and projection_x_coordinate) using parameters provided by the grid mapping. In this new case, there is just one coordinate, namely the reduced_gaussian_index, from which latitude and longitude can be obtained. It's similar to the HEALPix grid mapping, which was added in CF 1.13.

• It used to be mandatory to provide 2D latitude and longitude auxiliary coordinate variables if a grid mapping is in use. This was dropped because of use-cases where the 2D coordinate variables are very large, potentially making the data three times larger if a file contains a single data variable. Issue Location / Identification for grid specs outside of a file #357 (which you opened) suggests that it should be possible to store information defining the domain in another file. That issue stalled 18 months ago. It needs a summary to restart it.

• Regardless of that discussion, I think it should still be recommended to include 2D latitude and longitudes, if it's practical to do so, because they make the data more generally usable.

• UGRID is for unstructured grids, which require more explicit information than structured ones. Any horizontal rectilinear or curvilinear grid could be described by UGRID, but it would take more space, and it wouldn't provide a description of the grid construction.

Cheers

Jonathan

[larsbarring]

If I am not mistaken, 2560 grid cells at the equator corresponds to about 15 km resolution. This means that the shape of the earth is becoming a factor that needs to be considered. As we saw in the issue on HEALPIX the "concept of latitude" is becoming more complex (e.g. this comment) at high resolution. Gaussian grids are not (currently?) aiming at the same very high resolution, and I guess that the ERA6 system is using a spherical earth, which means that the full complexity discussed in relation to HEALPIX is not needed here. Still I do think that it is relevant to specify the shape of the earth, which in this case would be the radius of the sphere using attribute earth_radius (cf. Appendix F: Grid Mapping). This would benefit all users that may wish to reproject the data using some (other) map projection, e.g for publication purposes.

[LoreaGSM]

Thank you for your input @ChrisBarker-NOAA . The main idea behind this proposal is that the data that will be shared in GRIB2 format in this type of grid, can be shared with the community also in CF-NetCDF format without any need to interpolate it to a different grid (e.g. a regular one) which would mean altering the data. We also believe that a structured grid, even if not a regular one, should be encoded as such, not as an unstructured grid. This will not only mean a reduction in size, which is important when sharing large amount of files and data, but also a more optimal way of using this data by the different existing tools.

Thank you also for your comments @JonathanGregory , I think you are right and the latitude should have the same name as its dimension, this would strengthen the grid specifications by making it a NUG coordinate variable. As for the grid resolution, as you have suggested before, maybe we should include only a reference to the specification of the grid, something like reduced_gaussian: grid_specification= “O1280” and omit both, the subtype and the resolution. For an easier and faster recognition of the type of grid by the user and/or tools. I like your suggestion for the regional example, and I think latitude bounds should be strongly recommended.

As for the shape of the Earth, @larsbarring , I think it is a very relevant suggestion. I agree with you and this can be a really important attribute to include.

[ChrisBarker-NOAA]

Thanks @JonathanGregory for the clarifications -- I was a bit confused.

However, it does seem like we've got two different, but similar, concepts captured by grid_mapping (which is probably fine, but worth being explicit about).

1. Mapping of coordinates from, e.g. one projection to another in this case, you are storing the same number of coordinates, that is: x, y but specifying how they map to lon, lat, so any space savings is by not having to store both, but otherwise it's essentially the same. And, indeed, if an end-user wants to work with the native coordinate system (projection), then they don't have to do anything special.

2. Mapping an indexed system to actual coordinates -- I think this is what HealPIx is, and what I hadn't realized was a thing.

So this is a proposal for a new mapping of type (2). Makes sense to me.

I think it should still be recommended to include 2D latitude and longitudes

Agree on that.

Side note: as per my example of a GIS raster: specifying a corner, delta-x, delta-y and size -- maybe we should provide that as a grid_mapping? It would allow easier / less "lossy" specification for GIS rasters.

@LoreaGSM wrote:

The main idea behind this proposal is that the data that will be shared in GRIB2 format in this type of grid, can be shared with the community also in CF-NetCDF format without any need to interpolate it to a different grid (e.g. a regular one) which would mean altering the data.

I fully agree there, that is critically important.

We also believe that a structured grid, even if not a regular one, should be encoded as such, not as an unstructured grid.

I agree there as well, you definitely don't want to throw out structure if it's there to begin with. I should have taken the time to understand what a reduced gaussian grid actually is before making the UGRID comment.

But the question is, other than space storage issues (which, of course, are real), is there any problem with the storing the full grid in lon-lat coordinates in a structured way, similarly to how we do with curvilinear grids?

From my very brief attempt to understand the structure, it seems storing the grid points would be straightforward, but maybe not the cells -- I haven't yet found a diagram of the cells for a reduced gaussian grid. Maybe the whole concept of cells is irrelevant?

[davidhassell]

As Lars said, the shape of the Earth is indeed important, and it's worth highlighting that the attributes which describe the ellipsoid and prime meridian may be included, when applicable, with any grid mapping, and so would also be available to reduced_gaussian by definition. These are currently earth_radius, inverse_flattening, longitude_of_prime_meridian, prime_meridian_name, reference_ellipsoid_name, semi_major_axis, semi_minor_axis.

It may or may not be relevant here (I don't know), but for HEALPix, the shape of the Earth needed a special mention in its grid mapping definition ("If the healpix grid mapping is used on an ellipsoid, the latitudes returned by the [GHB05] algorithm must be interpreted as authalic latitude, instead of geodetic latitude (for a sphere, the two latitudes are the same). This is necessary to retain the equal area property"). Would the reduced Gaussian algorithm need any similar (but obviously different to the HEALPix case!) qualification?

[JonathanGregory]

Dear @ChrisBarker-NOAA

But the question is, other than space storage issues (which, of course, are real), is there any problem with the storing the full grid in lon-lat coordinates in a structured way, similarly to how we do with curvilinear grids?

The problem is that it isn't a rectangular grid in index space; the number of longitudes depends on the index of latitude. Hence it could be stored in a 2D array only if you padded most latitude rows with missing data. In addition, you would then require a 2D auxiliary coordinate variable for longitude, with the same padding. Curvilinear grids like the NEMO tripolar one are truly 2D in index space, although of course with lots of missing data on land if the data is for the ocean.

Best wishes

Jonathan

[ChrisBarker-NOAA]

Thanks all for accommodating my ignorance.

First: This does seem like an appropriate use for a grid mapping, akin to HEALPix, so no issue there.

However, I'm still a bit confused, and I think maybe there could be a way to accommodate this sort of grid in a more CF-like way as well.

If I'm really off-base, then tell me to stop :-).

Agreed -- this is not an unstructured grid, nor is is a logically rectangular grid, so we can't simply use those specs to store it.

However, the term "unstructured" grid is really a bit of a misnomer -- UGRids don't have no structure, they very much have a structure, it's just not a regular structure. And the UGRID spec is a way to specify what that structure is.

So I'm wondering if a similar principle could be applied to Reduced Gaussian Grids.

But where I'm confused is what the structure IS to a reduced Gaussian grids -- with a little googling, I've found discussions of the grid, and diagrams that show points, but not one that shows how those points relate to one another -- e.g. "cells". Is there a cell definition???

IIUC, an RGG has points (nodes?) at a defined number of (regularly spaced?) latitudes, and at each latitude there are regularly spaced points, with a reduced number of points as latitude increases. All in a way so that the the lat-lon locations can be computed from a small number of parameters.

All good, and why a grid mapping makes sense.

And you could also simply store the lat-long coordinates of each point (and I think it's been suggested that that's a recommendation) but this would "throw away" the structure. But not completely, it would be apparent from the data which points lie on which latitude line, for instance.

But what I'm not clear on is what that structure is? -- how the points along one latitude line relate to the points on the next latitude line? That is, if I want to interpolate a value to a point in between the latitude lines, how do I do that? what points do I want to interpolate from? And if this is clearly defined, maybe we can find a way to specify that other than doing the RGG math.

NOTE: this isn't idle speculation -- in my work, what I generally need is exactly that -- to interpolate values from model results to arbitrary points -- I haven't yet encountered a RGG, but if I did, I'd want to know how to do that. And ideally a CF compliant dataset would contain all the information I need.

[taylor13]

If regridding is of interest, models are often formulated such that when one model component is calculated on one grid and another on a different grid, "conservative" remapping algorithms are used so that energy, mass, and the like are conserved. These algorithms consider which original grid cells contribute to each target grid cell (and their relative importance). For general algorithms of this sort the cells are defined by the locations of their vertices. The only approximations typically made by these algorithms is that the vertices are joined along great circles which isn't always the case. Some of this is discussed in, for example, Taylor (2024); https://doi.org/10.5194/gmd-17-415-2024).

But you are probably already aware of how regridding typically works.

[ChrisBarker-NOAA]

I'm not talking about regridding, but rather, interpolation to particular points for:

• comparison to measurements or other models
• in my case: for a particle tracking model, where we need values at arbitrary points as the particles move.

But the issues are pretty much the same except that for general interpolation, conservation is less important (or not as well defined, anyway).

so:

... "which original grid cells contribute to each [point]. ... the cells are defined by the locations of their vertices."

So that should be defined in a "complete" CF grid specification.

Which is the question I have -- how are the cells defined for a reduced Gaussian grid?

Maybe that's part of the grid spec, but that isn't clear to me from (a tiny bit) of googling. All see is how "points" or "nodes" are defined.

And even if it is, I think there should be a way to specify it more directly via variables / attributes in CF, even if optional, like we do for other grid types.

[JonathanGregory]

Dear @ChrisBarker-NOAA

Lorea wrote in her initial proposal that, "Longitude and latitude bounds are located in the middle point between each consecutive longitude and each consecutive latitude" in this grid (at least, in ECMWF's use of it). The longitudes are computed rather than stored explicitly in the proposed grid-mapping, which means there has to be a stated method for computing them. In my comment I suggested stating that

The bounds in longitude for cell i are 0° + (m ± ½) × 360° / pl[k],

where m is the longitude index and k the latitude index, from which the cell index i is computed.

The latitudes are stored explicitly in a one-dimensional variable, which has the attributes of a coordinate variable. It can therefore specify the latitude bounds explicitly in a variable named by its bounds attribute, as usual. I believe that Lorea is happy with these suggestions. Does that answer the question about the definition of cells?

Best wishes

Jonathan