Montag, 24. November 2014

Regridding spatial data

This post is about one of the most difficult tasks in processing spatial data, in particular gridded data: regridding. Regridding is all about "warping" a given data blob into a different grid. First of all, it is tremendously helpful to really know what grid the original data is in. Here are a few examples:
Regular lon/lat grid:
Given the corner coordinate as a longitude/latitude pair and the per-data-pixel offset in degrees, this is one of the easiest format for gridded spatial data. In terms of projections, this can be expressed as a standard epsg:4326 projection (disregarding subtile differences in the meaning of longitude and latitude depending on the used "datum", i.e. the used model of the earth), meaning that all longitude/latitude coordinate pairs are simply "flattened" to the plane. Dealing with this kind of data is usually pleasant, gdal and other spatial libraries fully support this, no big deal. Except, the notion of longitude/latitude in commonly used data files in the data format "grib" (which is often used in numerical weather models) is slightly (usually neglectible) different from the standard WGS84 datum.
Rotated lon/lat grid:
The regular lon/lat grid has (at least) one major shortcoming: Considering the whole globe, it is far too detailed near the poles and too coarse at the equator. The reason is simply that the number of grid points per latitude, starting from the north pole going south at constant latitude steps, is constant. So if the region of interest is a small area it might be a considerable improvement to the accuracy if the pole is "virtually" rotated to that area, such that the highest accuracy is obtained for this ROI. Although this rotation does not pose a huge problem computationally, not all libraries can properly deal with data files that contain data in this form (again, "grib" supports this). Also, it is usually a bit more tricky to formulate a simple projection definition for this kind of data.
Gaussian grid:
The Gaussian grid basically uses non-evenly spaced latitude values in combination with evenly spaced longitude values (per one latitude the longitude values are evently spaced, but the number of longitude values may be different at different latitudes). Again, this grid type is not terribly complicated, but it is also not commonly supported (no easy projection definition). Luckily, the errors introduced if treating a Gaussian grid like a regular grid are usually neglectible (and mostly close to the poles anyway).
Triangular/hexagonal grid:
This is probably the most complicated (apart from completely irregular grids) grid. There is no projection definition, no support in most libraries, and dealing with this kind of grid usually boils down to handle it just like a complete irregular grid, meaning the need to extract the coordinates of the cells (usually: the center).

So, in order to use the data, some kind of transformation is necessary. That can be either the knowledge of the coordinates (in lon/lat) for each point, or the warping of the whole data blob into some well-known grid. For simple grids, it is as easy as using the standard gdal tool gdalwarp. For more complicated grids however, this is not enough. One solution that seems to work for transforming many kinds of grids back and forth is "cdo", which is also able to interpolate data from one grid to another. It appears that "cdo" is not very widely known, but if grid transformations is an issue, it is definitely worth a try!

Although cdo is indeed capable of dealing with/regridding triangular or unstructured grids, unfortunately it has some issues when the input is a grib file (issues = segmentation faults). So the "easy" solution is to convert the grib file to a different file format (this, for some reason, is possible with cdo), such as NetCDF. Again, there is a "but": gribs with triangular meshes cannot be translated to NetCDF directly. Instead, the grib data needs to be transferred to a different format, along with a manual definition of the grid. The result can then be converted to NetCDF and finally warped to a different grid.

Update #2:
Luckily, cdo is under active development and the bug using grib files has been fixed (will be in version 1.6.7)!