Regridding E3SM Data with ncremap

This page is devoted to instruction in NCO’s regridding operator, ncremap. It describes steps necessary to create grids, and to regrid datasets between different grids with ncremap. Some of the simpler regridding options supported by ncclimo are also described at Generate, Regrid, and Split Climatologies (climo files) with ncclimo. This page describes those features in more detail, and other, more boutique features often useful for custom regridding solutions.

The Zen of Regridding

Most modern climate/weather-related research requires a regridding step in its workflow. The plethora of geometric and spectral grids on which model and observational data are stored ensures that regridding is usually necessary to scientific insight, especially the focused and variable resolution studies that E3SM models conduct. Why does such a common procedure seem so complex? Because a mind-boggling number of options are required to support advanced regridding features that many users never need. To defer that complexity, this HOWTO begins with solutions to the prototypical regridding problem, without mentioning any other options. It demonstrates how to solve that problem simply, including the minimal software installation required. Once the basic regridding vocabulary has been introduced, we solve the prototype problem when one or more inputs are "missing", or need to be created. The HOWTO ends with descriptions of different regridding modes and workflows that use features customized to particular models, observational datasets, and formats. The overall organization, including TBD sections (suggest others, or vote for prioritizing, below), is:

Software Requirements:

At a minimum, install a recent version NCO on your executable $PATH with the corresponding library on your $LD_LIBRARY_PATH. NCO installation instructions are here. Unless you have reason to do otherwise, we recommend installing NCO through the Conda package. The Conda NCO package automagically installs two other important regridding tools, the ESMF_RegridWeightGen (aka ERWG) executable and the TempestRemap (aka TR) executables. Execute 'ncremap --config' to verify you have a working installation:

zender@aerosol:~$ ncremap --config
ncremap, the NCO regridder and grid, map, and weight-generator, version 4.9.3-alpha02 "Fuji Rolls"
...[Legal Stuff]...
Config: ncremap script located in directory /Users/zender/bin
Config: NCO binaries located in directory /Users/zender/bin, linked to netCDF library version 4.7.3
Config: No hardcoded machine-dependent path/module overrides. (If desired, turn-on NCO hardcoded paths at supported national labs with "export NCO_PATH_OVERRIDE=Yes").
Config: External (non-NCO) program availability:
Config: ESMF weight-generation command ESMF_RegridWeightGen version 7.1.0r found as /opt/local/bin/ESMF_RegridWeightGen
Config: MOAB-Tempest weight-generation command mbtempest not found
Config: MPAS depth coordinate addition command add_depth.py found as /Users/zender/bin/add_depth.py
Config: TempestRemap weight-generation command GenerateOfflineMap found as /usr/local/bin/GenerateOfflineMap

Only NCO is required for many operations including applying existing regridding weights (aka, regridding) and/or generating grids, maps, or conservative weights with the NCO algorithms. Generating new weights (and map-files) with ERWG or TR requires that you install those packages (both of which come with the NCO Conda package). MOAB-Tempest (MBTR) is not yet available in a pre-packaged format and must currently be built from scratch (MBTR is expected to come as a Conda package sometime in 2020). MBTR is only required for power-users. Make sure ncremap reports a sufficiently working status as above before proceeding further.

Prototypical Regridding I: Use Existing Map-file

The regridding problem most commonly faced is converting output from a standard resolution model simulation to equivalent data on a different grid for visualization or intercomparison with other data. The EAM v1 model low-resolution simulations are performed and output on the ne30np4 SE (spectral element) grid, aka the "source grid". The corresponding EAM v2 simulations were conducted on the ne30pg2 FV (finite volume) grid. E3SM source grids like these are called “unstructured” because they have only one horizontal dimension (i.e., 1D) which makes them difficult to visualize. The recommended 2D (latitude-longitude) grids for analysis (aka the "destination grid") of v1 simulations was the 129x256 Cap grid (the gridcells at the poles look like yarmulke caps), and since v2 is the 180x360 equi-angular grid used by CMIP. The single most important capability of a regridder is the intelligent application of weights that transform data on the input grid to the desired output grid. These weights are stored in a "map-file", a single file which contains all the necessary weights and grid information necessary. While most regridding problems revolve around creating the appropriate map-file, this prototype problem assumes that has already been done, so the appropriate map-file (map.nc) already exists and ncremap can immediately transform the input dataset (dat_src.nc) to the output (regridded) dataset (dat_rgr.nc):

ncremap -m map.nc dat_src.nc dat_rgr.nc

This solution is deceptively simple because it conceals the choices, paths, and options required to create the appropriate map.nc for all situations. We will discuss creating map.nc later after showing more powerful and parallel ways to solve the prototype problem. The solution above only works for users savvy enough to know how to find appropriate pre-built map-files. E3SM map-files are available at https://web.lcrc.anl.gov/public/e3sm/inputdata/cpl/gridmaps/ . Many commonly used maps and grids can also be found in my (@czender's) directories as ~zender/data/[maps,grids] at most DOE High Performance Computing (HPC) centers. Take a minute now to look at these locations.

Pre-built map-files use the (nearly) standardized naming convention map_srcgrd_to_dstgrd_algtyp.YYYYMMDD.nc, where srcgrd and dstgrd are the source and destination grid names, algtyp is a shorthand for the numerical regridding algorithm, and YYYYMMDD is the date assigned to the map (i.e., the date the map was created). The source grid in the example above is called ne30np4, the destination is called fv129x256. A pre-built map for the v1 combination is map.nc = map_ne30np4_to_fv129x256_aave.20150901.nc, and for v2 is map_ne30pg2_to_cmip6_180x360_aave.20200201.nc. What is aave? Weight generators can use about a dozen interpolation algorithms for regridding, and each has a shorthand name. For now, it is enough to know that the two most common algorithms are (first-order) conservative area-average regridding (aave) and bilinear interpolation (bilin or blin). Hence this conservatively regrids dat_src.nc to dat_rgr.nc with first order accuracy:

ncremap -m map_ne30np4_to_fv129x256_aave.20150901.nc dat_src.nc dat_rgr.nc # EAMv1 ncremap -m map_ne30pg2_to_cmip6_180x360_aave.20200201.nc dat_src.nc dat_rgr.nc # EAMv2 ncremap -m map_ne30pg2_to_cmip6_180x360_traave.20231201.nc dat_src.nc dat_rgr.nc # EAMv3

Before looking into map-file generation in the next section, try a few ncremap features. For speed's sake, regrid only selected variables:

ncremap -v FSNT,AODVIS -m map.nc dat_src.nc dat_rgr.nc

To regrid multiple files with a single command, supply ncremap with the source and regridded directory names (drc_src, drc_rgr). It regrids every file in drc_src and places the output in drc_rgr:

Or supply specific input filenames on the command-line, piped through standard input, or redirected from a file:

When an output directory is not specified, ncremap writes to the current working directory. When the output and input directories are the same, ncremap appends a string (based on the destination grid resolution) to each input filename to avoid name collisions. Finally, be aware that multiple-file invocations of ncremap execute in parallel by default. Power users will want to tune this as described in the section on "Intermediate Regridding".

Prototypical Regridding II: Create Map-file from Known Grid-files

The simplest regridding procedure applies an existing map-file to your data, as in the above example. E3SM map-files are publicly available at https://web.lcrc.anl.gov/public/e3sm/inputdata/cpl/gridmaps/. At most DOE High Performance Computing (HPC) centers these can also be found in my (@czender's) directory, ~zender/data/maps. If the desired map-file cannot be found, then you must create it. Creating a map-file requires a complete specification of both source and destination grids (meshes). The files that contain these grid specifications are called "grid-files". Many E3SM grid-files are publicly available within model-specific directories of the previous location, e.g., https://web.lcrc.anl.gov/public/e3sm/inputdata/ocn/mpas-o/oEC60to30v3/ . At most DOE High Performance Computing (HPC) centers these can also be found in my (@czender's) directory, ~zender/data/grids. Take a minute now to look there for the prototype problem grid-files, i.e., for FV 129x256 and ne30np4 grid-files.

You might find multiple grid-files that contain the string 129x256. Grid-file names are often ambiguous. The grid-file global metadata (ncks -M grid.nc) often displays the source of the grid. These metadata, and sometimes the actual data (fxm: link), are usually more complete and/or accurate in files with a YYYYMMDD-format date-stamp. For example, the metadata in file 129x256_SCRIP.20150901.nc clearly state it is an FV grid and not some other type of grid with 129x256 resolution. The metadata in 129x256_SCRIP.130510 tell the user nothing about the grid boundaries, and some of the data are flawed. When grids seem identical except for their date-stamp, use the grid with the later date-stamp. The curious can examine a grid-file (ncks -M -m grid.nc) and easily see it looks completely different from a typical model or observational data file. Grid-files and data-files are not interchangeable.

Multiple grid-files also contain the string ne30. These are either slightly different grids, or the same grids store in different formats meant for different post-processing tools. The different SE (and FV) grid types are described with figures and described here (https://acme-climate.atlassian.net/wiki/spaces/Docs/pages/34113147/Atmosphere+Grids). As explained there, many people will want the "dual-grid" with pentagons. The correct grid-file for this is ne30np4_pentagons.091226.nc. Do not be tempted by SE grid-files named with latlon.

All grid-files discussed so far are in SCRIP-format, named for the Spherical Coordinate Remapping and Interpolation Package (authored by @pjones). Other formats exist and are increasingly important, especially for SE grids. For now just know that these other formats are also usually stored as netCDF, and that some tools allow non-SCRIP formats to be used interchangeably with SCRIP.

Once armed with source and destination grid-files, one can generate their map-file with

Regrid a datafile at the same time as generating the map-file for archival:

Regrid a datafile without archiving the (internally generated) map-file:

For the prototype problem, the map-file generation procedure becomes

The map-files above are named with alg_typ=nco because the ncremap default interpolation algorithm is the first-order conservative NCO algorithm (NB: before NCO 4.9.1 the default algorithm was ESMF bilin). To re-create the aave map in the first example, invoke ncremap with -a esmfaave (the newest v3 naming convention) or -a conserve (same algorithm, different name in v1, v2):

This takes a few minutes, so save custom-generated map-files for future use. Computing weights to generate map-files is much more computationally expensive and time-consuming than regridding, i.e., than applying the weights in the map-file to the data. We will gloss over most options that weight-generators can take into consideration, because their default values often work well. One option worth knowing now is -a. The invocation synonyms for -a are --alg_typ, --algorithm, and --regrid_algorithm. These are listed in the square brackets in the self-help message that ncremap prints when it is invoked without argument, or with --help:

At least one valid option argument for each supported interpolation type is shown separated by vertical bars. The arguments shown have multiple synonyms, separated by commas, that are equivalent. For example, -a esmfaave is equivalent to -a aave and --alg_typ=conserve. Use the longer option form for clarity and precision, and the shorter form for conciseness. The full list of synonyms, and the complete documentation, is at http://nco.sf.net/nco.html#alg_typ. The NCO algorithm ncoaave is the default. Commonly-used algorithms that invoke ERWG are esmfbilin and esmfaave. TR options are discussed below. Peruse the list of options now, though defer a thorough investigation until you reach the "Intermediate Regridding" section.

Prototypical Regridding III: Infer Grid-file from Data-file

Thus far we have explained how to apply a map-file to data, and how, if necessary, to generate a map-file from known grids. What if there is no map-file and the source or the destination grid-files (or both) are unavailable? Often, one knows the basic information about a grid (e.g., resolution) but lacks the grid-file that contains the complete information for that grid geometry. In such cases, one must create the grid-file via one of two methods. First, one can let ncremap attempt to infer the grid-file from a data file known to be on the desired grid. This procedure is called "inferral" and is fairly painless. Second, one can feed ncremap all the required parameters (for rectangular grids only) and it will generate a grid-file. This requires a precise specification of the grid geometry, and will be covered the sub-section on "Manual Grid-file Generation".

Before we describe what the inferral procedure does, here is an example that demonstrates how easy it is. You can regrid an SE (or FV) dataset from our prototype example to the same grid as an evaluation dataset. Pick any 2D (e.g., MxN latitude-by-longitude) dataset to compare the E3SM simulations to. Inferral uses the grid information in the evaluation dataset, which is already on the desired destination grid, to create the (internally generated) destination grid-file. Supply ncremap with any dataset on the desired destination grid (dat_dst.nc) with -d (for "dataset") instead of -g (for "grid"):

This tells ncremap to infer the destination grid-file and to use it to generate the desired map-file, named with the supposed destination resolution (here, MxN degrees or gridcells). To archive the inferred destination grid for later use, supply a name for the grid with -g:

One can infer a grid without having to regrid anything. Supply ncremap with a data-template file (dat.nc) and a grid-file name (grd.nc). Since there are no input files to regrid, ncremap exits after inferring the grid-file:

Grid-inferral is easier to try than manual grid-generation, and will work if the data file contains the necessary information. The only data needed to construct a SCRIP grid-file are the vertices of each gridcell. The gridcell vertices define the gridcell edges and these in turn define the gridcell area which is equivalent to the all-important weight parameter necessary to regrid data. Of course the gridcell vertices must be stored with recognizable names and/or metadata indicators. The Climate & Forecast (CF) metadata convention calls the gridcell vertices the "cell-bounds". Coordinates (like latitude and longitude) usually store cell-center values, and should, according to CF, have bounds attributes whose values point to variables (e.g., lat_bounds or lon_vertices) that contain the actual vertices. Relatively few datasets "in the wild" contain gridcell vertices, though the practice is, happily, on the rise. Formally SE models have nodal points with weights without any fixed perimeter or vertices assigned to the nodes, so the lack of vertices in SE model output is a feature, not an oversight. The dual-grid (referenced above) addresses this by defining "pretend" gridcell vertices for each nodal point so that an SE dataset can be treated like an FV dataset. However, dual-grids are difficult to generate, and may not exist or be accurate for many SE grids. In that case, ncremap cannot infer the grid (because the vertices are unknown) and one needs to use a different package (such as TempestRemap, below) to construct the grid-file and the mapping weights.

Inferral works well on important categories of grids for which ncremap can guess the missing grid information. In the absence of gridcell vertice information, ncremap examines the location of and spacing between gridcell centers and can often determine from these what type of grid a data-file (not a grid-file!) is stored on. A data-file simply means the filetype preferred for creation/distribution of modeled/observed data. Hence ncremap has the (original and unique, so far as we know) ability to infer all useful rectangular grid-types from data-files that employ the grid. The key constraint here is "rectangular", meaning the grid must be orthogonal (though not necessarily regularly spaced) in latitude and longitude. This includes all uniform angle grids, FV grids, and Gaussian grids. For curvilinear grids (including most swath data), ncremap infers the underlying grid to be the set of curves that bisect the curves created by joining the gridcell centers. This often works well for curvilinear grids that do not cross a pole. Inferral works for unstructured (i.e., 1D) grids only when the cell-bounds are stored in the datafile as described above. Hence inferral will not work on raw output from SE models.

A few more examples will flesh-out how inferral can be used. First, ncremap can infer both source and destination grids in one command:

Here the user provides only data files (no grid- or map-files) yet still obtains regridded output! The first positional (i.e., not immediately preceded by an option indicator) argument (dat_src.nc) is interpreted as the data to regrid, and the second positional argument (dat_rgr.nc) is the output name. The -d argument is the name of any dataset (dat_dst.nc) on the desired destination grid. ncremap infers the source grid from dat_src.nc, then infers the destination grid from dat_dst.nc, then generates weights (with the default algorithm since none is specified) and creates a map-file that it uses to regrid dat_src.nc to dat_rgr.nc. No grid-file or map-file names were specified (with -g or -m) so both grid-files and the map-file are generated internally in temporary locations and erased after use.

Second, this procedure, like most ncremap features, works for multiple input files:

Unless a map-file or source grid-file is explicitly provided (with -m or -s, respectively), ncremap infers a separate source grid-file (and computes a corresponding map-file) for each input file. This allows it to regrid lists of uniquely gridded data (e.g., satellite swaths each on its own grid) to a common destination grid. When all source files are on the same grid (as is typical with models), then turn-off the expensive multiple inferral and map-generation procedures with the -M switch to save time:

Prototypical Regridding IV: Manual Grid-file Generation

If a desired grid-file is unavailable, and no dataset on that grid is available (so inferral cannot be used), then one must manually create a new grid. Users create new grids for many reasons including dataset intercomparisons, regional studies, and fine-tuned graphics. NCO and ncremap support manual generation of the most common rectangular grids as SCRIP-format grid-files. Create a grid by supplying ncremap with a grid-file name and "grid-formula" (grd_sng) that contains, at a minimum, the grid-resolution. The grid-formula is a hash-separated string of name-value pairs each representing a grid parameter. All parameters except grid resolution have reasonable defaults, so a grid-formula can be as simple as "latlon=180,360":

Congratulations! The new grid-file grd.nc is a valid source or destination grid for ncremap commands.

Grid-file generation documentation in the NCO Users Guide at http://nco.sf.net/nco.html#grid describes all the grid parameters and contains many examples. Note that the examples in this section use grid generation API for ncremap version 4.7.6 (August, 2018) and later. Earlier versions can use the ncks API explained in the Users Guide.

The most useful grid parameters (besides resolution) are latitude type (lat_typ), longitude type (lon_typ), title (ttl), and, for regional grids, the SNWE bounding box (snwe). The three supported varieties of global rectangular grids are Uniform/equiangular (lat_typ=uni), Cap/FV (lat_typ=cap), and Gaussian (lat_typ=gss). The four supported varieties of longitude types are the first (westernmost) gridcell centered at Greenwich (lon_typ=grn_ctr), western edge at Greenwish (grn_wst), or at the Dateline (lon_typ=180_ctr and lon_typ=180_wst, respectively). Grids are global, uniform, store latitudes from south-to-north, and have their first longitude centered at Greenwich by default. The grid-formula for this is 'lat_typ=uni#lon_typ=grn_ctr#lat_drc=s2n'. Some examples (remember, this API requires NCO 4.7.6+):

Regional grids are a powerful tool in regional process analyses, and can be much smaller in size than global datasets. Regional grids are always uniform. Specify the rectangular bounding box, i.e., the outside edges of the region, in SNWE order:

Intermediate Regridding:

The sections on Prototypical Regridding were intended to be read sequentially and introduced the most frequently required ncremap features. The Intermediate and Advanced regridding sections are an a la carte description of features most useful to particular component models, workflows, and data formats. Peruse these sections in any order.

Intermediate Regridding I: Treatment of Missing Data

The terminology of this section can be confusing. The most important point to clarify is that the treatments of missing data described here are independent of the regridding algorithm (and weight-generator) which is specified by alg_typ. alg_typ only defines the algorithm used to generate the weights contained in a map-file, not how to apply those weights. ncremap applies the map-file weights (i.e., regrids) to data-files that may contain fields with spatio-temporally missing values, such as AODVIS in EAM, all  fields with ocean gridpoints in ELM output, most fields in MPAS-Seaice, and many fields in satellite-retrieved datasets. Unless the locations of a field's (like AODVIS) missing values were supplied as a mask to the weight-generator algorithm at map-creation time, the map-file weights cannot and will not automatically take into account that some of the source gridcells contain invalid (aka missing) values. This section describes three options for how the weight-applicator (not weight-generator) is to reconstruct values in destination gridcells affected by missing source gridcells. One option preserves the integral value of the input gridcells and is called "integral-preserving" or "conservative" because it is first-order conservative in the affected gridcells locally, as well as globally conservative. At the other extreme is the option to preserve (locally) the gridcell mean of the input values in the affected output gridcells. This option is called "mean-preserving" or "renormalized" and is not conservative. The third option allows the user to specify a sliding scale between the integral- and mean-preserving options, which is to say it is a more general form of renormalization. Once again, these weight-application approaches are independent of the weight-generation algorithm defined by alg_typ.

Conservative weight-generation is, for first-order accurate algorithms, a straightforward procedure of identifying gridcell overlap and apportioning values correctly from source to destination. The absence of valid values (or presence of missing values) forces a choice on how to construct destination gridcell values where some but not all contributing source cells are valid. ncremap supports three distinct approaches: "Integral-preserving", "Mean-preserving", and "Sliding Scale". The integral-preserving approach uses all valid data from the input grid on the output grid once and only once. Destination cells receive the weighted valid values of the source cells. This is a first-order conservative procedure, and will, in conjunction with weights from conservative regridding algorithm, conserve and preserve the local and global integrals of the source and destination fields. The mean-preserving approach divides the destination value by the sum of the valid weights. This preserves and extends the mean of the valid input values throughout the entire destination gridcell. In other words, it extrapolates valid data to missing regions while preserving the mean valid data. Input and output integrals are unequal and mean-preserving regridding is not conservative and does not preserve the integrals when missing data are involved. Furthermore, mean-preserving regridding only preserves the local (i.e., gridcell-level) mean of the input data, not the global mean. Both approaches, integral-preserving and mean-preserving (aka conservative and re-normalized) produce identical answers when no missing data maps to the destination gridcell. Before explaining the nuance of the third approach ("sliding-scale"), we demonstrate the symmetric ways of invoking integral- or mean-preserving regridding: 

By default, ncremap implements the integral-preserving, conservative approach because it has useful properties, is simpler to understand, and requires no additional parameters. However, this approach will often produce unrealistic gridpoint values (e.g., ocean temperatures < 100 K and that is not a typo) values along coastlines or near data gaps where state variables are regridded to/from small fractions of a gridcell. The mean-preserving approach solves this problem, yet incurs others, like non-conservation. The mean-preserving approach (aka re-normalization) ensures that the output values are physically consistent and do not stick-out like sore thumbs in a plot, although the integral of their value times area is not conserved.

The third option, the "sliding scale" approach, is slightly more complex and much more flexible than the all-or-nothing integral- or mean-preserving examples above. The sliding-scale refers to the fraction of the destination gridcell that must be overlapped by valid source gridcells. Users can specify the renormalization threshold weight rnr_thr which is the required valid fraction of destination gridcells with the renormalization threshold option "--rnr_thr=rnr_thr". The weight-application algorithm then ensures that valid values overlap at least the fraction rnr_thr of each destination gridcell for that gridcell to meet the threshold for a non-missing destination value. When rnr_thr is exceeded, the mean valid value in the valid area is placed in the destination gridcell. If the valid area covers less than rnr_thr, then the destination gridcell is assigned the missing value. Valid values of rnr_thr range from zero to one. A threshold weight of rnr_thr = 0.0 indicates that any amount (no matter how small) of valid data should be represented by its mean value on the output grid. Keep in mind though, that the actual valid area divides a sum (the area-weighted integral of the valid values), and values of zero or very near to zero in the divisor can lead to floating-point underflow and divide-by-zero errors. A threshold weight of rnr_thr = 1.0 indicates that the entire destination gridcell (to within machine precision) must be overlapped by valid data in order for that destination gridcell to be assigned a valid (non-missing) value. Threshold weights 0.0 < rnr_thr < 1.0 invoke the real-power of the sliding scale. Remote sensing classification schemes applied to L2 data may require, for example, valid retrievals over at least 50% of source pixels before assigning a valid value to a destination gridcell. This is equivalent to rnr_thr = 0.5. And so, 

These sliding scale examples specify that valid values must cover at least 25% and 50% of the destination gridcell to meet the threshold for a non-missing destination value. With actual valid destination areas of 25% or 50%, this approach would produce destination values greater than the conservative algorithm by factors of four and two, respectively. Careful readers may already have observed that the mean-preserving approach (--preserve=mean) is exactly equivalent to the sliding scale approach with a renormalization threshold weight rnr_thr = 0.0. The latter approach is just a more numerical way of expressing the former approach. On the other hand, no numerical threshold weight produces answers equivalent to the integral-preserving approach. However, it is convenient for power users to be able invoke all three missing data approaches via the rnr_thr option alone. Hence, we made setting the threshold weight to the string "none" (--rnr_thr=none) to be exactly equivalent to specifying the integral-preserving approach. In summary, the --preserve option with its two valid arguments integral and mean that invoke the integral-preserving and mean-preserving algorithms will suffice for most situations. The sliding-scale algorithm must be invoked via the --rnr_thr option with a numerical argument, although the string "none" will default to the integral-preserving algorithm.

In practice, it may make sense to use the default integral-preserving treatment with conservative weights, and the mean-preserving or sliding-scale treatment with other (non-conservative) weights such as those produced by bilinear interpolation or nearest-neighbor. Another consideration in selecting the weight-application algorithm for missing values is whether the fields to regrid are fluxes or state variables. For example, temperature (unlike heat) and concentrations (amount per unit volume) are not physically conserved quantities under areal-regridding so it often makes sense to interpolate them in a non-conservative fashion, to preserve their fine-scale structures. Few researchers can digest the unphysical values of temperature that the integral-preserving treatment produces in regions rife with missing values. On the other hand, mass and energy fluxes should be physically conserved under areal-regridding. Hence, one must consider both the type of field and its conservation properties when choosing a missing value treatment.

Intermediate Regridding II: TempestRemap

TempestRemap (TR) is the chief alternative to ERWG and NCO for regridding weight-generation. TR is replacing ERWG as the default on-line weight generator in E3SMv2. Tempest algorithms, written by @paulullrich, have many numerical advantages described in papers and at Transition to TempestRemap for Atmosphere grids . Verify that ncremap can access your Tempest installation as described in the above section on "Software Requirements" before trying the examples below. Once installed, TR can be as easy to use as ERWG or NCO with FV grid-files, e.g.,

This command, the same as shown in the "Create Map-file from Known Grid-files" section above, except using alg_typ='fv2fv_flx', is a jumping-off point for understanding Tempest features and quirks. First, simply note that the ncremap interfaces for ERWG, NCO, and TR weight-generators are the same even though the underlying ERWG, NCO, and TR applications have different APIs. Second, Tempest accepts SCRIP-format grids (as shown) and Exodus-format grid-files, also stored in netCDF though typically with a '.g' suffix, e.g., ne30.g as described at Transition to TempestRemap for Atmosphere grids Exodus grid-files contain grid "connectivity" and other information required to optimally treat advanced grids like SE. Hence this also works

This produces subtly different weights because ne30.g encodes the SE ne30np4 grid-specification, not its dual-grid FV representation. TR generates superior weights when instructed to use an algorithm optimized for the type of remapped variable and the grid representation as described below. The above example employs the recommended algorithm to remap fluxes on the SE grid to the FV destination grid, i.e., se2fv_flx.

The exact options that NCO invokes for a specific TR algorithm like se2fv_flx can be discovered in multiple ways: First, all TR options that NCO employs are recommended on the Transition to TempestRemap for Atmosphere grids; Second, the NCO Users Guide documents TR options at http://nco.sf.net/nco.html#tr; Third, the options --dbg_lvl=1 or --dbg_lvl=2 cause ncremap to print the sub-commands it executes, including the TR commands with options. Experiment, intercompare, and find the algorithm that works best for your purposes. Advanced users may be interested in the quantitative evaluations of the quality of the weights in the map-file provided by the map-checker (ncks --chk_map map.nc) described below.

As mentioned at  the Tempest overlap-mesh generator expects to be sent the two grid-files in the order smaller then larger (ERWG and NCO have no corresponding restriction). For example, Tempest considers the global ocean to be a smaller domain than the global atmosphere since it covers less area (due to masked points). Hence ncremap must be told when the source grid-file covers a larger domain than the destination. Do this with the "--a2o" or "--l2s" switch (for "atmosphere-to-ocean" or "large-to-small", respectively):

Transition to TempestRemap for Atmosphere grids describes eight specific global FV<->SE mappings that optimize different combinations of accuracy, conservation, and monotonicity desirable for remapping flux-variables (flx), state-variables (stt), and an alternative (alt) mapping. A plethora of boutique options and switches control the Tempest weight-generation algorithms for these six cases. To simplify their invocation, ncremap names these eight algorithms fv2se_stt, fv2se_flx, fv2se_alt, fv2fv_flx, fv2fv_stt, se2fv_stt, se2fv_flx, se2fv_alt, and se2se. E3SM maps with these algorithms have adopted the suffixes "mono" (for fv2fv_flx, se2fv_flx, and fv2se_alt), "highorder" (fv2fv_stt, se2fv_stt, fv2se_stt), "intbilin" (se2fv_alt), and "monotr" (fv2se_flx). The relevant Tempest map-files can be generated with

These maps, applied to appropriate flux and state variables, should exactly reproduce the online remapping in the E3SM v2/v3 coupler. However, explicitly generating all standard maps this way is not recommended because ncremap includes an MWF-mode (for "Make All Weight Files") described below. MWF-mode generates and names, with one command and in a self-consistent manner, all combinations of E3SM global atmosphere<->ocean maps for ERWG and Tempest.

Intermediate Regridding III: MOAB/mbtempest support

Weight-generator computations (e.g., finding grid intersections) increase non-linearly with the grid size, so the largest grids are most efficiently computed with parallel algorithms. ERWG has long supported distributed computation in an MPI-environment, and NCO has always supported multi-threaded weight computation via OpenMP. A growing subset of TempestRemap algorithms have now been ported to the parallel MOAB (Mesh Oriented datABase) tool mbtempest (aka MBTR). ncremap can invoke the MOAB regridding toolchain as of NCO version 5.0.2 from September, 2021. The "spoiler" answer to "How do I get ncremap to invoke mbtempest?" is simple: Add the --mpi_nbr=mpi_nbr option to an ncremap command that already calls a TempestRemap algorithm. However, MOAB involves complex, MPI-enabled software, and support for it is continually changing and subject to important limitations. Read on to understand what to currently expect.

First, MOAB requires an MPI environment to perform well. Invoking MOAB (i.e., using --mpi_nbr=mpi_nbr with a TR algorithm name) in a non-MPI environment will result in an error. One can easily obtain MPI-enabled MOAB environment with Conda. For example, install Conda MPI versions of the MOAB package with

conda install -c conda-forge moab=5.3.0=mpich_tempest esmf

Please ensure you have the latest version of ERWG, MOAB, and/or TempestRemap before reporting any related problems to NCO.

This section is intentionally placed after the TR section because mbtempest re-uses the TR algorithm names described in the previous section. For example, ncremap invokes TR to generate weights to remap fluxes from FV to FV grids when invoked with

ncremap -a fv2fv_flx --src_grd=src.g --dst_grd=dst.nc -m map.nc # Invoke TR

Add the --mpi_nbr option and ncremap will instead invoke the MOAB toolchain to compute weights for any TempestRemap algorithm (otherwise the TR toolchain would be used):

ncremap --mpi_nbr=8 -a fv2fv_flx --src_grd=src.g --dst_grd=dst.nc -m map.nc # Invoke MOAB

Although mapping weights generated by MOAB and TempestRemap use the same numerical algorithms, they are likely to produce slightly different weights due to round-off differences. MOAB is heavily parallelized and computes and adds terms together in an unpredictable order compared to the serial TempestRemap.

Transparently to the user, ncremap supplies the selected weight generator with the recommended options, which can be quite complex. For the preceding example, ncremap either invokes TempestRemap's GenerateOverlapWeights with the boutique options --in_type fv --in_np 1 --out_type fv --out_np 1 --correct_areas that E3SM recommends for conservative and monotone remapping of fluxes, or it invokes multiple components of the MOAB toolchain:

mbconvert -B -o PARALLEL=WRITE_PART -O PARALLEL=BCAST_DELETE -O PARTITION=TRIVIAL -O PARALLEL_RESOLVE_SHARED_ENTS "src.g" "src.h5m"
mbconvert -B -o PARALLEL=WRITE_PART -O PARALLEL=BCAST_DELETE -O PARTITION=TRIVIAL -O PARALLEL_RESOLVE_SHARED_ENTS "dst.nc" "dst.h5m"
mbpart 8 --zoltan RCB "src.h5m" "src_8p.h5m"
mbpart 8 --zoltan RCB --recompute_rcb_box --scale_sphere --project_on_sphere 2 "dst.h5m" "dst_8p.h5m"
mpirun -n 8 mbtempest --type 5 --weights --load "src_8p.h5m" --load "dst_8p.h5m" --method fv --order 1 --method fv --order 1 --file "map.nc"

The purpose of the ncremap front-end to MOAB is to hide this complexity from the user while preserving the familiar look and feel of invoking other weight-generators. Once again, the MOAB toolchain should produce a map-file identical (to rounding precision) to one produced by TR. When speed matters (i.e., large grids), and the algorithm is supported (i.e., fv2fv_flx), invoke MOAB, otherwise invoke TR.

This section comes before the parallelism section because mbtempest supports MPI-enabled parallelism that is distinct from the ncremap workflow parallelism described in the next section.

Caveat lector: As of September, 2021 Weights generated by MOAB (version 5.3.0 and earlier) are only trustworthy for the fv2fv_flx algorithm. The options for all other algorithms are implemented as indicated though they should be invoked for testing purposes only. High order and spectral element map algorithms are not fully implemented and will produce unreliable results. MOAB anticipates supporting more TempestRemap algorithms in the future. Always use the map-checker to test maps before use, e.g., with ncks --chk_map map.nc.

Another limitation is that mbtempest (version 5.3.0) currently only generates map-files that regrid to rank 1 (unstructured) grids. mbtempest "unrolls" any rank 2 (e.g., latxlon) destination grid into its rank 1 equivalent. Users seeking to regrid to rank 2 grids can manually alter the MOAB-generated mapfile with something like:

ncks -O -x -v dst_grid_dims ~/map.nc ~/map_fix.nc
ncap2 -O -s 'defdim("dst_grid_dims",2);dst_grid_dims($dst_grid_dims)={256,129};' ~/map_tmp.nc ~/map_fix.nc

Replace 256,129 above by the lon,lat (Fortran order!) of your grid. The resulting map-file, map_fix.nc, will almost regrid as intended. The regridded output will have curvilinear coordinates (e.g, lat(lat,lon)) instead of rectangular coordinates (e.g., lat(lat)), and the coordinate bounds variables (e.g., lat_bnds) will not be correct and may cause plotting issues at poles and Greenwich. Nevertheless, the weights are correct and of unsurpassed accuracy. MOAB anticipates supporting rank 2 mapfiles and TR high-order and spectral element algorithms in the future.

Intermediate Regridding IV: Parallelism

ncremap can exploit three types of parallelism: multiples nodes in a cluster, multiple simultaneous file regriddings on a single node, and multiple simultaneous variables regridded within a single file. Each level of parallelism reduces the wallclock time to complete the regridding workflow at the expense of increase resource requirements.

File-level parallelism accelerates throughput when regridding multiple files in one ncremap invocation, and has no effect when only one file is to be regridded. Note that the ncclimo and ncremap semantics for selecting file-level parallelism are identical, though their defaults differ (Background mode for ncclimo and Serial mode for ncremap). Select the desired mode with the argument to --par_typ=par_typ. Explicitly select Background mode with par_typ values of bck, background, or Background. The values mpi or MPI select MPI mode, and the srl, serial, Serial, nil,  and none will all select Serial mode (which disables file-level parallelism, though still allows intra-file OpenMP parallelism).

The default file-level parallelism for ncremap is Serial mode (i.e., no file-level parallelism), in which ncremap processes one input file at a time. Background and MPI modes implement true file-level parallelism. Typically both these parallel modes scale well with sufficent memory unless and until I/O contention becomes the bottleneck. In Background mode ncremap issues all commands to regrid the input file list as UNIX background processes on the local node. Nodes with mutiple cores and sufficient RAM take advantage of this to simultaneously regrid multiple files. In MPI mode ncremap issues commands to regrid the input file list in round-robin fashion to all available compute nodes. Prior to NCO version 4.9.0 (released December, 2019), Background and MPI parallelism modes both regridded all the input files at one time and there was no way to limit the number of files being simultaneously regridded. Subsequent versions allow finer grained parallelism by introducing the ability to limit the number of discrete workflow elements or ``jobs'' (i.e., file regriddings) to perform simultaneously within an ncremap invocation or ``workflow''.

As of NCO version 4.9.0 (released December, 2019), the --job_nbr=job_nbr option specifies the maximum number of files to regrid simultaneously on all nodes being harnessed by the workflow. Thus job_nbr is an additional parameter to fine-tune file level parallelism (it has no effect in Serial mode). In both parallel modes ncremap spawns processes in batches of job_nbr jobs, then waits for those processes to complete. Once a batch finishes, ncremap spawns the next batch. In Background mode, all jobs are spawned to the local node. In MPI mode, all jobs are spawned in round-robin fashion to all available nodes until job_nbr jobs are running.

If regridding consumes so much RAM (e.g., because variables are large and/or the number of threads is large) that a single node can perform only one regridding job at a time, then a reasonable value for job_nbr is the number of nodes, node_nbr. Often, however, nodes can regrid multiple files simultaneously. It can be more efficient to spawn multiple jobs per node than to increase the threading per job because I/O contention for write access to a single file prevents threading from scaling indefinitely.

By default job_nbr = 2 in Background mode, and job_nbr = node_nbr in MPI mode. This helps prevent users from overloading nodes with too many jobs. Subject to the availability of adequate RAM, expand the number of jobs per node by increasing job_nbr until, ideally, each core on the node is used. Remember that processes and threading are multiplicative in core use. Four jobs each with four threads each consumes sixteen cores.

We have thus far demonstrated how to control file-level parallelism (with par_typ) and workflow level parallelism (with job_nbr). The third level of parallelism mentioned above is that ncremap can use OpenMP shared-memory techniques to simultaneosly regrid multiple variables within a single file. This shared memory parallelism is efficient because it requires only a single copy of the regridding weights in physical memory to regrid multiple variable simultaneously. Even so, regridding multiple variables at high resolution may become memory-limited, meaning that the insufficient RAM can often limit the number of variables that the system can simultaneously regrid.

By convention all variables to be regridded share the same regridding weights stored in a map-file, so that only one copy of the weights needs to be in memory, just as in Serial mode. However, the per-thread (i.e., per-variable) OpenMP memory demands are considerable, with the memory required to regrid variables amounting to no less than about 5-7 times (for type NC_FLOAT) and 2.5-3.5 times (for type NC_DOUBLE) the size of the uncompressed variable, respectively. Memory requirements are so high because the regridder performs all arithmetic in double precision to retain the highest accuracy, and must allocate separate buffers to hold the input and output (regridded) variable, a tally array to count the number of missing values and an array to sum the of the weights contributing to each output gridcell (the last two arrays are only necessary for variables with a _FillValue attribute). The input, output, and weight-sum arrays are always double precision, and the tally array is composed of four-byte integers. Given the high memory demands, one strategy to optimize thr_nbr for repetitious workflows is to increase it to keep doubling it (1, 2, 4, ...) until throughput stops improving. With sufficient memory, the NCO regridder scales well up to 8-16 threads.

As an example, consider regridding 100 files with a single map. Say you have a five-node cluster, and each node has 16 cores and can simultaneously regrid two files using eight threads each. (One needs to test a bit to optimize these parameters.) Then an optimal (in terms of wallclock time) invocation would request five nodes with 10 simultaneous jobs of eight threads. On PBS or SLURM batch systems this would involve a scheduler command like qsub -l nodes=5 ... or sbatch --nodes=5 ..., respectively, followed by ncremap --par_typ=mpi --job_nbr=10 --thr_nbr=8 .... This job will likely complete between five and ten-times faster than a serial-mode invocation of ncremap to regrid the same files. The uncertainty range is due to unforeseeable, system-dependent load and I/O charateristics. Nodes that can simultaneously write to more than one file fare better with multiple jobs per node. Nodes with only one I/O channel to disk may be better exploited by utilizing more threads per process.

Advanced Regridding:

Advanced procedures have in common that they activate non-standard processing modes for ncremap. These modes do something different, or in addition to, the standard weight-generation and application. Generally these modes were created in order to automate frequently recurring workflows that can leverage the ncremap infrastructure so long as various bells and whistles are introduced along the way. Please let us know if you have ideas for new or improved processing modes.

Advanced Regridding I: MPAS-mode

MPAS models produce output in their own format distinct from CESM-heritage models. MPAS-mode invokes three pre-processing steps to massage MPAS datasets until they are amenable to regridding. These steps include missing value annotation, missing value treatment, and dimension permutation. We will shortly describe these steps in order. First, though, MPAS-mode, like most other ncremap modes, is explicitly invoked with the -P md_nm option where the md_nm is some variation of the MPAS component model name:

Many model and observational datasets use missing values that are not annotated in the standard manner. For example, the MPAS ocean and ice models use -9.99999979021476795361e+33 as the missing value, yet (at least from 2015-2018) do not store this in a _FillValue attribute with any variables. To prevent arithmetic from treating these values as valid, MPAS-mode automatically puts this value in a _FillValue attribute for all floating-point variables via

Oddly, the MPAS land-ice model uses -1.0e36 for missing values, so currently MPAS-LI users must explicitly supply this missing value, or invoke ncremap with the -P mali option

Next, MPAS datasets usually have masked regions (e.g., non-ocean cells) yet MPAS users like to visualize regridded data with realistic (not blocky) boundaries along those cells and so they decided to, by default, treat missing values with the renormalization approach described above in the section on Treatment of Missing Values. Hence MPAS-mode automatically invokes ncremap with maximum renormalization, equivalent to

Finally, ncremap requires the horizontal spatial dimension(s), whether latitude and longitude or some unstructured dimension, to be the final (most-rapidly-varying) dimension(s) of input datasets. MPAS datasets natively place their horizontal spatial dimension (typically nCells) closer to the least-rapidly-varying position. While this makes perfect sense from an I/O-efficiency point-of-view for unstructured models, it does not play well with regridders. Hence all MPAS-modes permute the input dimensions to a regridder-friendly order (i.e., ending with nCells) with a command of the form

The combination of these three data manipulations defines MPAS-mode. It can be difficult to learn the ocean mesh (i.e., grid) names and thus to find the appropriate pre-made map-files. The standard low resolution ocean map-files for post-processing each version of MPAS are:

ncremap -P mpasocean --map=map_oEC60to30v3_to_cmip6_180x360_aave.20181001.nc mpov1.nc out.nc # MPAS Ocean v1
ncremap -P mpasocean --map=map_EC30to60E2r2_to_cmip6_180x360_aave.20220301.nc mpov2.nc out.nc # MPAS Ocean v2
ncremap -P mpasocean --map=map_IcoswISC30E3r5_to_cmip6_180x360_traave.20231201.nc mpov3.nc out.nc # MPAS Ocean v3

ncremap -P mpasseaice --map=map_oEC60to30v3_to_cmip6_180x360_aave.20181001.nc msiv1.nc out.nc # MPAS Seaice v1
ncremap -P mpasseaice --map=map_EC30to60E2r2_to_cmip6_180x360_aave.20220301.nc msiv2.nc out.nc # MPAS Seaice v2
ncremap -P mpasseaice --map=map_IcoswISC30E3r5_to_cmip6_180x360_traave.20231201.nc msiv3.nc out.nc # MPAS Seaice v3

Advanced Regridding II: EAMXX-mode

EAMXX storage conventions and files differ from EAM (and CAM) in only two ways. ncremap introduced a -P eamxx mode to support these gotchas in 2022. Dimension permutation is the primary pre-processing step necessary to massage EAMXX datasets so that they are amenable to regridding.

First, ncremap requires the horizontal spatial dimension(s), whether latitude and longitude or some unstructured dimension, to be the final (most-rapidly-varying) dimension(s) of input datasets. EAMXX datasets natively place their horizontal spatial dimension (typically ncol) closer to the least-rapidly-varying position. While this makes perfect sense from an I/O-efficiency point-of-view for unstructured models, it does not play well with horizontal regridding. Hence EAMXX-mode permutes the input dimensions to a regridder-friendly order (i.e., ending with ncol) with a command of the form

Second, EAMXX names the surface pressure variable ps by default, not PS as in EAM. The distinction is important whenever vertical interpolation is invoked. Hence EAMXX mode automatically tells the vertical interpolaion routine to look for ps not PS. The combination of these two pre-processing steps defines EAMXX-mode.

Advanced Regridding III: Sub-Gridscale Regridding (SGS-mode)

ncremap has a sub-gridscale (SGS) mode that performs the special pre-processing and weighting necessary to conserve fields that represent fractional spatial portions of a gridcell, and/or fractional temporal periods of the analysis. Spatial fields output by most geophysical models are intensive, and so by default the regridder attempts to conserve the integral of the area times the field value such that the integral is equal on source and destination grids. However some models (like ELM, CLM, CICE, and MPAS-Seaice) output gridcell values intended to apply to only a fraction sgs_frc (for "sub-gridscale fraction'') of the gridcell. The sub-gridscale fraction usually changes spatially with the distribution of land and ocean, and spatiotemporally with the distribution of sea ice and possibly vegetation. For concreteness consider a sub-grid field that represents the land fraction. Land fraction is less than one in gridcells that resolve coastlines or islands. ELM and CLM happily output temperature values valid only for a small (i.e., sgs_frc << 1) island within the larger gridcell. Model architecture dictates this behavior and savvy researchers expect it. The goal of the NCO weight-application algorithm is to treat SGS fields as seamlessly as possible so that those less familiar with sub-gridscale models can easily regrid them correctly.

Fortunately, models like ELM and CLM that run on the same horizontal grid as the overlying atmosphere can use the same mapping-file as the atmosphere, so long as the SGS weight-application procedure is invoked. Not invoking an SGS-aware weight application algorithm is equivalent to assuming sgs_frc = 1 everywhere. Regridding sub-grid values correctly versus incorrectly (e.g., with and without SGS-mode) alters global-mean answers for land-based quantities by about 1% for horizontal grid resolutions of about one degree. The resulting biases are in intricately shaped regions (coasts, lakes, sea-ice floes) and so are easy to overlook.

To invoke SGS mode and correctly regrid sub-gridscale data, specify the names of the fractional area sgs_frc and, if applicable, the mask variable sgs_msk (strictly, this is only necessary if these names differ from their respective defaults landfrac and landmask). Trouble will ensue if sgs_frc is a percentage or an absolute area rather than a fractional area (between zero and one). ncremap must know the normalization factor sgs_nrm by which sgs_frc must be divided (not multiplied) to obtain a true, normalized fraction. Datasets (such as those from CICE) that store sgs_frc in percent should specify the option --sgs_nrm=100 to instruct ncremap to normalize the sub-grid area appropriately before regridding. ncremap will re-derive sgs_msk based on the regridded values of sgs_frc: sgs_msk = 1 is assigned to destination gridcells with sgs_frc > 0.0, and all others sgs_msk = 0. As of NCO version 4.6.8 (released June, 2017), invoking any of the options --sgs_frc, --sgs_msk, or --sgs_nrm, automatically triggers SGS-mode, so that also invoking -P sgs is redundant though legal. As of NCO version 4.9.0 (released December, 2019), the values of the sgs_frc and sgs_msk variables should be explicitly specified. In previous versions they defaulted to landfrac and landmask, respectively, when -P sgs was selected. This behavior still exists but will likely be deprecated in a future version.

The area and sgs_frc fields in the regridded file will be in units of sterradians and fraction, respectively. However, ncremap offers custom options to reproduce the idiosyncratic data and metadata format of two particular models, ELM and CICE. When invoked with -P elm (or -P clm), a final step converts the output area from sterradians to square kilometers. When invoked with -P cice, the final step converts the output area from sterradians to square meters, and the output sgs_frc from a fraction to a percent.

It is sometimes convenient to store the sgs_frc field in an external file from the field(s) to be regridded. For example, CMIP-style timeseries are often written with only one variable per file. NCO supports this organization by accepting sgs_frc arguments in the form of a filename followed by a slash and then a variable name:

Files regridded using explicitly specified SGS options will differ slightly from those regridded using the -P elm or -P cice options. The former will have an area field in sterradians, the generic units used internally by the regridder. The latter produces model-specific area fields in square kilometers (for ELM) or square meters (for CICE), as expected in the raw output from these two models. To convert from angular to areal values, NCO assumes a spherical Earth with radius 6,371,220 m or 6,371,229 m, for ELM and CICE, respectively. The ouput sgs_frc field is expressed as a decimal fraction in all cases except for -P cice which stores the fraction in percent. Thus the generic SGS and model-specific convenience options produce equivalent results, and the latter is intended to be indistinguishable (in terms of metadata and units) to raw model output. This makes it more interoperable with many existing analysis scripts.

Advanced Regridding IV: Regional Unstructured Output (RRG-mode)

EAM (and CAM-SE) will produce regional output if requested to with the finclNlonlat namelist parameter. Output for a single region can be higher temporal resolution than the host global simulation. This facilitates detailed yet economical regional process studies. Regional output files are in a special format that we call RRG (for "regional regridding"). An RRG file may contain any number of rectangular regions. The coordinates and variables for one region do not interfere with other (possibly overlapping) regions because all variables and dimensions are named with a per-region suffix string, e.g., lat_128e_to_134e_9s_to_16s. ncremap can easily regrid RRG 2D logically rectangular output from an FV-dycore because ncremap can infer (as discussed above) the grid from any well-annotated regional FV data file. Regridding unstructured regional grid data, however, is more complex because unstructured grids without cell vertices and unstructured grid weight-generators are not yet flexible enough to to output only regional (as opposed to global) grids with weights. To summarize, regridding RRG data leads to three difficulties (#1-3 below) and two difficulties (#4-5) shared with FV RRG files:

  1. RRG files contain only regional gridcell center locations, not weights

  2. Global SE grids have well-defined weights not vertices for each gridpoint

  3. Grid generation software (ESMF and TempestRemap) only create global not regional SE grid files

  4. Non-standard variable names and dimension names

  5. Regional files can contain multiple regions

ncremap's RRG mode resolves these issues to allow trouble-free regridding of SE RRG files. The user must provide two additional input arguments, --dat_glb=dat_glb and --grd_glb=grd_glb to point to a global SE dataset and grid, respectively, of the same resolution as the model that generated the RRG datasets. Hence a typical RRG regridding invocation is:

Here grd_rgn is a regional destination grid-file, dat_rgn is the RRG file to regrid, and dat_rgr is the regridded output. Typically grd_rgn is a uniform rectangular grid covering the same region as the RRG file. Generate this as described in the last example in the section above on "Manual Grid-file Generation". grd_glb is the standard dual-grid grid-file for the SE resolution of the simulation, e.g., ne30np4_pentagons.091226.nc. ncremap regrids the global data file dat_glb to the global dual-grid in order to produce a intermediate global file annotated with gridcell vertices. Then it hyperslabs the lat/lon coordinates (and vertices) from the regional domain to use with regridding the RRG file. A grd_glb file with only one 2D field suffices (and is fastest) for producing the information needed by the RRG procedure. One can prepare an optimal dat_glb file by subsetting any 2D variable (e.g., ncks -v FSNT in.nc dat_glb.nc) from a full global SE output dataset.

ncremap RRG mode supports two additional options to override parameters set internally. First, the per-region suffix string may be set with --rnm_sng=rnm_sng. RRG mode will, by default, regrid the first region it finds in an RRG file. Explicitly set the desired region with rnm_sng for files with multiple regions, e.g., --rnm_sng= . Second, the bounding-box of the region may be explicitly set with --bb_wesn=lon_wst,lon_est,lat_sth,lat_nrt. The normal parsing of the bounding-box string from the suffix string may fail in (as yet undiscovered) corner cases, and the --bb_wesn option provides a workaround. The bounding-box string must include the entire RRG region, specified in WESN order. The two override options may be used independently or together, as in:

RRG-mode supports most normal ncremap options, including input and output methods and regridding algorithms.

Advanced Regridding V: Make All Weight Files (MWF-mode)

As mentioned above in the TempestRemap section, ncremap includes an MWF-mode (for "Make All Weight Files") that generates and names, with one command and in a self-consistent manner, all combinations of E3SM global atmosphere<->ocean maps with both ERWG and Tempest. MWF-mode automates the laborious and error-prone process of generating numerous map-files with various switches. Its chief use occurs when developing and testing new global grid-pairs for the E3SM atmosphere and ocean components. Invoke MWF-mode with a number of specialized options to control the naming of the output map-files:

where grd_ocn is the "global" ocean grid, grd_atm, is the global atmosphere grid, nm_src sets the shortened name for the source (ocean) grid as it will appear in the output map-files, nm_dst sets, similarly, the shortend named for the destination (atmosphere) grid, and dt_sng sets the date-stamp in the output map-file name map_${nm_src}_to_${nm_dst}_${alg_typ}.${dt_sng}.nc. Setting nm_src, nm_dst, and dt_sng is optional though highly recommended. For example,

produces the 10 ERWG map-files:

The ordering of source and destination grids is immaterial for ERWG maps since MWF-mode produces all map combinations. However, as described above in the TempestRemap section, the Tempest overlap-mesh generator must be called with the smaller grid preceding the larger grid. For this reason, always invoke MWF-mode with the smaller grid (i.e., the ocean) as the source, otherwise some Tempest map-file will fail to generate. The six optimized SE<->FV Tempest maps described above in the TempestRemap section will be generated when the destination grid has a ".g" suffix which ncremap interprets as indicating an Exodus-format SE grid (NB: this assumption is an implementation convenience that can be modified if necessary). For example,

produces the 6 TempestRemap map-files:

MWF-mode takes significant time to complete (~20 minutes on my MacBookPro) for the above grids. To accelerate this, consider installing the MPI-enabled instead of the serial version of ERWG. Then use the --wgt_cmd option to tell ncremap the MPI configuration to invoke ERWG with, for example:

Background and distributed node parallelism (as described above in the the Parallelism section) of MWF-mode are possible though not yet implemented. Please let us know if this feature is desired.

Advanced Regridding VI: CMIP6 Timeseries

This section describes the recommended procedures to construct and regrid E3SM timeseries data to CMIP6 specifications. Most models provide data to CMIP6 in timeseries format, meaning one variable-per-file with multiple years per file. These timeseries must be regridded to at least one of the CMIP6 standard grids. The E3SM project chose to supply its v1 experiments to CMIP6 archived on rectangular, uniform (i.e., equiangular in latitude and longitude), one-degree (for standard-resolution) and quarter-degree (for high-resolution) grids. Generating these timeseries from experiments as lengthy as 500 model years, formatted to CMIP6 specifications, requires many non-standard options to both ncclimo (to construct the timeseries) and to ncremap (to regrid timeseries), and is a natural capstone exercise in using both together. This section is arranged in reverse order where first we present the final actual commands, followed by the descriptions, meanings, and reasons for particular options.

The recommended procedures for generating EAM and MPAS-O timeseries of the 500-yr DECK pre-industrial simulations for CMIP6 are:

Take a moment to compare the methods for EAM and for MPAS. They are nearly identical except for the variable names, experiment names and directories, map-files (so far nothing surprising or important) AND the additional MPAS options in ${mpas_opt}. We will discuss those soon. Each command-set begins with setting experiment-dependent I/O directories and a map-files. Other experiments will require changing these to the appropriate I/O directories, yet the map-file remains the same unless the native or destination grid changes. The next three or four lines in each command-set configure the splitter and regridder with options that many ncclimo/ncremap users have never before tried. Finally the list of input files and all the configuration options are sent to ncclimo. The entire procedure for the user boils down to creating then executing this one splitter command for each desired variable.

Regridding is performed only if the splitter (i.e., ncclimo) is invoked with the --map option that supplies a suitable mapfile from the native to the desired destination grid. CMIP6 will only distribute data on 2D structured grids, yet E3SM will itself distribute the timeseries on native (unstructured grids). Hence the commands above construct both the native timeseries (stored in ${drc_out}) and the regridded timeseries for CMIP6 (stored in ${drc_rgr}). Internally, the splitter constructs the native grid timeseries for the same time segement for all requested variables in parallel, waits for completion, and then calls the regridder (ncremap) in parallel with all timeseries for that segment, waits, then continues to the next segment.

The first configuration options to discuss are the MPAS-specific options. In order to automatically trigger a number of MPAS-specific behaviors, the regridder must first know that the model type is MPAS. When invoked with '-m mpas' the splitter will pass the MPAS-flag to the regridder. The splitter itself simply creates timeseries and does nothing different for MPAS files other than pass options through to the regridder. MPAS outputs its native grid data in double precision, not single precision like EAM/ELM/CAM. Thus raw MPAS datasets are twice the size needed by most analyses. The --d2f flag tells the regridder to demote doubles (unless they are coordinate variables) to floats in an additional pre-processing step. Otherwise, regridded MPAS output would be twice the size with no appreciable benefits for analysis.

The CMIP6-specific options (${cmip6_opt}) collectively ensure that the timeseries are compliant, compact, and concise. CMIP6 requires datasets be in netCDF4-Classic format, i.e., netCDF4 storage constrained to the netCDF3 API. This is achieved with the -7 switch (mnemonic: 7=4+3). Additionally, CMIP6 requires datasets use netCDF4's internal compression, the DEFLATE algorithm (same as in gzip). We recommend deflation level 1 (i.e., --dfl_lvl=1) since higher levels compress only marginally better yet require significantly more wallclock time.

The next three CMIP6 options trim the timeseries to exclude variables that could otherwise be included. The --no_cll_msr (no-cell-measures) switch excludes variables typically listed in the CF cell_measures attribute such as gridcell area and volume. The --no_frm_trm (no-formula-terms) switch excludes variables that appear in the CF formula_terms attribute, notably the 2D surface pressure (PS) for EAM. The --no_stg_grd (no-staggered-grid) switch excludes the offset (aka staggered) grid that ncremap normally adds to rectangular output grids. The specific variables excluded are slat, slon, and w_stag. There is no downside to this option for MPAS data, although it can cause problems for older versions of AMWG diagnostics. Thus timeseries processed with these options include no "extras" that might inflate their size or, alas, their convenience.

The splitter options (${spl_opt}) configure the timeseries length and number of segments. The splitter expects the number of (monthly) input files to equal the number of years (between ${yr_srt} and ${yr_end}, inclusive) times twelve. This sanity check helps prevent inadvertent omissions/inclusions of unwanted months. Each timeseries is split, if necessary, into a number of segments of equal length and possibly one shorter length tail segment. The --ypf option specifies the number of years per file (i.e., segment). CMIP6 recommends file sizes be no greater than a few gigabytes. Factors that influence the segment filesize include the segment length, the variable rank and number of layers if 3D, and, for regridded timeseries, the grid resolution and the presence of missing values (e.g., due to ocean bathymetry). A compromise that meets these criteria is segment lengths up to 500-years for 2D variables, and 25-years for 3D variables. For consistency, these same segments length limits are used for all E3SM v1 models and experiments in CMIP6. Note that because the segment lengths differ for 2D and 3D variables, it is necessary to call the splitter at least twice per experiment, once with 2D variables (supplied with --var) and segment size, and likewise for 3D variables and segment size.

These 500-year and 25-year segment lengths yield native-grid files of sizes ~800 MB and 2.3 GB for v1 EAM 2D and 3D variables, respectively, that have regridded (to 1x1 degree) sizes of ~1.0 GB and 3.0 GB. For v1 MPAS-Ocean data, these segment lengths yield native-grid files of sizes ~9.7 GB and 22 GB for MPAS-Ocean 2D and 3D variables, respectively, that have regridded (to 1x1 degree) sizes of ~900 MB and 1.5 GB. Hence all the regridded E3SM data distributed by CMIP6 will be in files of sizes between 1-3 GB. Note MPAS regridded data consumes ~90% less space than native grid data. This is due to two factors: 1. Raw MPAS data do not utilize the netCDF _FillValue attribute (which would substantially improve compression), and 2. Raw MPAS data are double precision not single precision.

The RAM overhead of timeseries generation can also be a factor on small nodes. Splitting does most of its work on disk and so requires only as much RAM as required to store a single timestep of a single variable. Regridding is a different kettle of fish, a bird of another feather, and potentially a can of worms. The maximum RAM usage is about three times the uncompressed size of the entire timeseries. For the 500-yr 2D and 25-yr 3D segments considered here, expect peak RAM usage of 20 GB and 64 GB, respectively, for MPAS data. If the regridder exhausts available memory when called with multiple variables, then reduce the parallelization over variables using the --job=${job_nbr} option (not shown). This is unlikely to occur on beefy nodes because job_nbr defaults to 2 (i.e., variables are split and regridded in groups of two). The splitter parallelizes well for typical timeseries of 2D variables, and can be invoked with ${job_nbr} exceeding 100 when no regridding (which consumes much more memory than splitting) is performed.

Now that the content of the rather lengthy CMIP6 splitter/regridder commands has been explained, it is worthwhile describing the method of invocation. The splitter accepts filenames supplied in numerous ways (command-line arguments, pipes to stdin, directory contents, redirection operators) as described above. For large numbers of input files typical of CMIP6 experiments, piping filenames as output by ls from the input file directory into the splitter is preferred for two reasons. First, ls automatically sorts files into alphanumeric order. This is equivalent to timeseries order because of the filename conventions employed by E3SM. Thus ls ensures that timeseries monotonically increase. Moreover, ls understands command-line globbing to simplify culling only required time periods from directories with longer simulations. Second, issuing ls from the input file directory removes the lengthy path component of each filename received by the splitter. For a 500-year pre-industrial DECK simulation, this removes 500*12=6000 copies of the same ~100-character directory path from the provenance metadata maintained in the history attribute of each downstream file.

Finally, note that we explicitly set ${TMPDIR} to a capacious writable directory prior to execution. The regridder writes all intermediate files to this directory, and removes them only upon successful completion. (The splitter itself never writes to ${TMPDIR}). However, for MPAS files, the regridder may write as many as three or four intermediate files per output file to ${TMPDIR}. Since some 3D MPAS-Ocean DECK PI files are 10's of GB in size, it is best to ensure the intermediate files are written to volatile storage. They will be automatically deleted upon successful completion of regridding. Should the splitter or regridder fail for any reason, the files will remain in ${TMPDIR} to assist in debugging. Thus it is best if ${TMPDIR} is automatically scrubbed every so often, e.g., on re-boots as with most Linux and MacOS workstations.

This discussion of splitting and regridding has focused on "one-off" experiments such as the DECK 500-yr pre-industrial simulation. The above methods with minor modifications also apply to ensemble experiments such as those with historical forcing since 1850. For example, this generates CMIP6 timeseries to analyze historical cloud radiative effects in the ensemble of five E3SM v1 simulations designated H1-H5:

The main difference between generating the timeseries for the Historical ensemble and the DECK PI experiment is the need to loop over the ensemble. Here the splitter command is not backgrounded so that one member experiment is processed at a time (to avoid overwhelming nodes with I/O and RAM demands). Set the input directory in the ensemble loop, and ensure the globbing pattern for filenames matches the naming convention used for all five members. Consider whether to output to member-specific directories or to a single, ensemble-wide directory. If the former, then nothing special need be done. If the latter, use the --fml_nm (family-name) option as above to avoid identical timeseries names (that will overwrite one another) and to create instead member-specific timeseries names like

CLDLOW_H1_185001_201412.nc
CLDLOW_H2_185001_201412.nc
...

Advanced Regridding VII: Initial Condition Files

First, use the right tool. ncremap can regrid an initial conditions (IC) file, both vertically and horizontally. However, generating scientifically validated IC files for new model resolutions is best done with a lengthy workflow (https://acme-climate.atlassian.net/wiki/spaces/ED/pages/872579110/Adding+support+for+new+grids) in which regridding plays a relatively small role. That said, ncremap is a good tool to place an atmospheric state onto a new grid where it can then be nudged into a valid IC file, or to place the contents of an IC file on a rectangular grid (as shown below) where it is easier to plot. Regridding atmosphere IC files was straightforward until E3SM v2 when the atmosphere separated the dynamical and physical grids. This resulted in IC files containing two sets of grid variables, so now we must regrid v2 IC files with two successive invocations of ncremap:

ncremap --map=map_ne30pg2_to_cmip6_180x360_nco.20200901.nc foo.eam.i.2001-01-01-00000.nc foo.nc
ncremap --map=map_ne30np4_to_cmip6_180x360_nco.20200901.nc foo.nc foo.eam.i.2001-01-01-00000.rgr.nc

The first command remaps the variables on the PG2 physics grid to the desired output grid. This is straightforward since the physics grid variables (SICTHK, ICEFRAC, TS1...) use ncol (the ncremap default) as the horizontal coordinate. The second command remaps the remaining variables, stored on the NP4 dynamics grid with horizontal dimension ncol_d, to the same output grid. The second invocation automatcally regrids all the variables on the dynamics grid (with the ncol_d dimension) because the intermediate file foo.nc no longer has an ncol dimension.

To regrid only the variables with a non-default (e.g., not ncol) horizontal dimension in a file with multiple horizontal dimensions, one would explicitly select the non-default dimension using the ncremap -R option, whose argument is passed directly to the underlying regridder (ncks):

ncremap -R '--rgr col_nm=ncol_d' --map=map_ne30np4_to_cmip6_180x360_nco.20200901.nc foo.eam.i.2001-01-01-00000.nc foo.rgr.nc

Advanced Regridding VIII: Fixing Grid Files to Work as Intended

Gridfiles store a wealth of highly precise information using loosely standardized rules that are open to interpretation. One may encounter gridfiles that conform to one regridder's expectations though not another's. The section provides guidance on how to adjust or repair some of the most frequently encountered problems with gridfiles. The problems currently addressed are floating-point masks, _FillValue, and imprecise RLL grids.

Floating-point mask variables (grid_imask) in SCRIP files---they are non-standard and may break some software. TempestRemap's GenerateOverlapMesh program, for example, breaks (as of this writing, 20210428) when asked to ingest a SCRIP grid-file with a floating point mask. This problem can occur when using grid-files from TR's old ConvertExodusToSCRIP program which output floating point masks. The new TR program, ConvertMeshToSCRIP, appears to have fixed this problem (as of this writing in 20240110). Also, users often create masks from floating-point variables (as described in the next section) and inadvertently leave the mask as a floating point variable. The solution to the problem of floating-point masks can be as simple as a straightforward conversion of the mask to an integer:

ncap2 -s 'grid_imask=int(grid_imask)' grd.nc # Integerize mask

However, ncap2 uses implicit type-conversion rules that truncate floating-point variables, i.e., round-down towards zero, so an input floating point mask value of 0.99999 will be converted to an integer mask value of zero in the output file. If this is not the desired behavior, consider

ncap2 -s 'grid_imask=int(round(grid_imask))' grd.nc # Integerize mask

This converts input floating point mask values between 0.5 and 1.5 to integer mask values of one in the output file.

The attribute to identify data whose values are missing, impossible, or not-yet defined in a netCDF dataset has long been _FillValue. Many people are used to adding _FillValue attributes to all variables in data files. That is usually fine and it is understandable how that practice sometimes gets carried into grid-files. However, there is no reason for a SCRIP grid-file to include any _FillValue attributes. The SCRIP format allows the incorporation of a mask variable (grid_imask, discussed further in the next section) that performs some of the roles of a _FillValue attribute. Masks designate cells intended to be ignored by the weight generation algorithm; this does not require that data values be ill-defined in masked-out cells, just that masked cells not contribute to (for source-grid masks) or receive (for destination-grid masks) mapping weights.

In any case, the interpretation of values identified as missing due to their equality with the _FillValue attribute is not-standardized and will likely vary among weight-generation software. To avoid this ambiguity, we recommend that the _FillValue attribute be removed from all variables in the grid-file:

ncatted -a _FillValue,,d,, grd.nc # Remove all _FillValue attributes

The last issue we address here is imprecise RLL grids. Within and before the last decade, climate modelers often generated rectangular (including Gaussian) grids with only single-precision arithmetic. That often works well enough, especially for low-resolution grids. However, we design weight generation software to be as precise and repeatable as possible, and this requires double-precision accuracy. Moreover, for reasons not completely understood (though perhaps related to the density of RLL meridians near the poles), high resolution single-precision RLL grids break some weight generation software. For example, single precision globally uniform 3 minute (0.05x0.05 degree) RLL grids tend to either break NCO or cause it to generate bad weights. ERWG usually completes weight generation with such grids, though the weights are often bad.

The first step in avoiding imprecise RLL grids is recognizing these symptoms. If ERWG or NCO produce bad weights, e.g., max(frac_a,frac_b) > 1.01, then this might be due to an imprecise grid. The next step is verify whether the RLL gridcell corners exhibit random behavior after the seventh digit or so. For example, examine the first (or last) few gridcells and verify that the increments between corners (and/or centers) is as expected:

zender@cori06:~$ grid=/global/cfs/cdirs/e3sm/inputdata/lnd/clm2/mappingdata/grids/SCRIPgrid_3minx3min_GLOBE-Gardner_c120922.nc
zender@cori06:~$ ncks --trd -H -d grid_size,,2 -v grid_corner_lat,grid_corner_lon,grid_center_lat,grid_center_lon ${grid}
grid_size[0] grid_center_lat[0]=-89.9749984741
grid_size[1] grid_center_lat[1]=-89.9749984741
grid_size[2] grid_center_lat[2]=-89.9749984741

grid_size[0] grid_center_lon[0]=-179.975006104
grid_size[1] grid_center_lon[1]=-179.925003052
grid_size[2] grid_center_lon[2]=-179.875
...

These gridcell centers random behavior after the seventh digit. If we regenerate this grid with full double-precision accuracy, the grid vertices and centers will be exactly 0.05 degrees apart as desired:

zender@cori06:~$ ncremap -G latlon=3600,7200#lat_typ=uni#lon_typ=180_wst -g grd.nc
zender@cori06:~$ ncks --trd -H -d grid_size,,2 -v grid_corner_lat,grid_corner_lon,grid_center_lat,grid_center_lon grd.nc
grid_size[0] grid_center_lat[0]=-89.975
grid_size[1] grid_center_lat[1]=-89.975
grid_size[2] grid_center_lat[2]=-89.975

grid_size[0] grid_center_lon[0]=-179.975
grid_size[1] grid_center_lon[1]=-179.925
grid_size[2] grid_center_lon[2]=-179.875
...

Once the grid contains coordinates accurate to double-precision, other grid-file fields such as the mask (grid_imask) can be appended at will, including from the original, imprecise file:

ncks -C -m -A -v grid_imask ${grid} grid.nc # Replace newly generated default mask with GLOBE/Gardner landmask
ncks -O -7 -L 1 grid.nc grid.nc # Optional step to compress large grids

Advanced Regridding IX: Creating and Using Land Surface Grids with Masks

Regridding software is commonly used to convert continuous fields into binary (True/False) masks for use within numerical models. Land surface models, for example, use a "land mask" to distinguish gridcells containing land (mask is True = 1) from non-land gridcells (mask is False = 0). Similarly, glacier models models use an "ice mask" to distinguish gridcells available to the ice sheet from non-ice gridcells. This section describes the steps to follow and potential crevasses to avoid when manipulating fields with masks into (primarily SCRIP-format) grid-files that can then be used as described above in the creation of weights to map between grids.

Masks intended for use in regridding must satisfy constraints imposed by most or all weight-generation software. Such masks must be integer-valued and time-invariant (no time dimension). Regridding masks are thus one-dimensional for unstructured grids, or two-dimensional for rectangular or curvilinear grids. Some weight-generation software (TempestRemap, for example) will crash if it encounters a floating-point valued mask. No known software will create, or apply, time-varying map-files, so time dimensions in grid files must be avoided.

Three generic tasks often occur during mask creation: 1) Using rules to combine continuous field(s) into a mask, 2) Merging a mask into a grid-file, 3) Inferring a mask from an exising continuous field.

As an example of defining a mask from a set of rules applied to continuous input data, consider the problem of defining an ice mask (e.g., for ELM/E3SM) from fields like bathymetry/topography and ice thickness, which MALI produces in the output variables bedTopography and thickness, respectively. Say one wishes the mask to include all gridcells that have multi-year land-based ice as well as gridcells that might form multi-year ice in a simulation. NCO can process these (and other) conditions into an integer-valued output field using the ncap2 where() operator:

ncap2 -s 'icemask[lat,lon]=0;where(thickness > 0 || bedTopography > 0) icemask=1' in.nc mask.nc

This command ensures an integer-valued mask by initializing the icemask variable to 0 (note the lack of a decimal point), an integer. Initializing to 0.0f or 0.0 would result in icemask being defined as a single or double-precision floating-point variable, respectively. If the input variables include a leading time dimension, then above command could be modified to define a time-invariant mask in terms of the variables at the initial timestep:

ncap2 -s 'icemask[lat,lon]=0;where(thickness(0,:,:) > 0 || bedTopography(0,:,:) > 0) icemask=1' in.nc mask.nc

A mask defined from a data file in this way must then be added to gridfile in order to be utilized when generating remapping weights. This is straightforward so long as the name, type, rank, and dimensionality of the mask variable match the expectations of the gridfile:

ncrename -v icemask,grid_imask mask.nc
ncks -A -C -v grid_imask mask.nc grd.nc

Next, consider inferring a mask from a single continuous input field. Suppose the input field has non-zero floating point values where the mask is to be true, and _FillValue or zero-values where the mask is to be false. ncremap will automatically infer the correct mask from a data file containing that input field. For example, RACMO ice sheet model data contains a field named Icemask_GR that is true only over multi-year Greenland ice. One can create a gridfile that contains a mask variable (grid_imask) derived from the Icemask_GR field like this:

ncremap --msk_dst=Icemask_GR --dst_fl=racmo_data.nc --grd_dst=grid_icemask.nc

In this case the grid-file is complete and maps may be created from it, e.g.,

ncremap --grd_src=grid_icemask.nc --grd_dst=1x1.nc --map=map.nc

Such a map applied to a RACMO dataset would allow only points inside the Greenland Ice Sheet to contribute to the destination (analysis) grid.

Epilogue: User-Suggested Examples

E3SM employs some the most advanced grid formulations used in ESM modeling today, and since regridding is a complex subject full of details, it is unlikely that the foregoing documentation already answers everyone's questions. Moreover, since some of the greatest E3SM innovations and optimization rely on even newer grid treatments, we are unlikely to ever stop needing to learn newer regridding techniques. If you read this far you know about ncremap's main capabilities, yet either you and U2 still haven't found what you're looking for, or you may have a hunch that one of those relevant-sounding-though-undocumented (at least in Confluence) options that appears with ncremap --help, is for.In that spirit, we welcome user requests for annotated examples for their real-world regridding issues. Fire away!