Skip to article frontmatterSkip to article content

Regrid

Tools for regridding MPAS-A output to regular lat-lon grids

Options

pyremap

pyremap simplifies the process of regridding MPAS data to other grids. It helps you compute weights (using ESMF_RegridWeightGen or the MOAB CLI[2]) and apply them, and works with global and limited-area meshes.

pyremap is available on conda-forge.

In the below example, we use the data we downloaded in Viz.

pyremap-example.py
from pathlib import Path

import xarray as xr
from pyremap import Remapper, get_lat_lon_descriptor

HERE = Path(__file__).parent

mesh_id = "x1.2562"
src_path = HERE / f"{mesh_id}.grid.nc"
out_dir = HERE
in_path = HERE / f"{mesh_id}.static.nc"
variables = ["ter"]

ds_in = xr.open_dataset(in_path, decode_times=False)

dst = get_lat_lon_descriptor(
    dlon=0.25,
    dlat=0.25,
    lon_min=-125,
    lon_max=-66,
    lat_min=25,
    lat_max=50,
)
# Note that the min/max lon/lat are cell edges, not centers.

for method in ["bilinear", "conserve"]:
    # Create weights with ESMF_RegridWeightGen
    # (It is also possible to use MOAB)
    map_path = out_dir / f"{mesh_id}_map_{method}.nc"
    remapper = Remapper(map_filename=map_path.as_posix(), method=method, ntasks=1, use_tmp=True)
    remapper.src_from_mpas(filename=src_path.as_posix(), mesh_name=mesh_id)
    remapper.dst_descriptor = dst
    remapper.build_map()

    # Apply weights with NCO ncremap, producing a file
    remapper.ncremap(
        in_filename=(HERE / f"{mesh_id}.static.nc").as_posix(),
        out_filename=(out_dir / f"out_{method[:3]}.nc").as_posix(),
        overwrite=True,
        variable_list=variables,
    )

    # Apply weights with numpy, returning an xarray Dataset
    ds = remapper.remap_numpy(ds_in[variables])
    print(ds)

NCO

NCO ncremap provides many facilities for regridding.

Examples

On Casper/Derecho:

module load nco

Let’s use the grid data we downloaded in Viz.

ln -s x1.2562.grid.nc grid.nc
ln -s x1.2562.static.nc static.nc

MPAS-Tools provides multiple ways to convert grid specs in the MPAS format to SCRIP format.

scrip_from_mpas is available when you install the mpas_tools conda-forge package [5].

# scrip_from_mpas requires [0, 2π) longitudes
# but our grid file has [-π, π) longitudes
PI=3.14159265359
ncap2 -s "where(lonCell < 0) lonCell = lonCell + 2*$PI;
where(lonVertex < 0) lonVertex = lonVertex + 2*$PI;
where(lonEdge < 0) lonEdge = lonEdge + 2*$PI" grid.nc grid_0to2pi.nc

# Create the SCRIP file (scrip.nc by default, use -s to change)
scrip_from_mpas -m grid_0to2pi.nc

mpas2esmf is not included in the conda-forge package (you must compile it from within this directory). You must supply a grid file with sphere_radius 1, but unlike scrip_from_mpas, it’s fine with negative longitudes.

# Create SCRIP file (mpas_scrip.nc) and ESMF grid file (mpas_esmf.nc)
mpas2esmf grid.nc "480-km" "$(date -I)"

Note that the results may be a bit different. For example, mpas2esmf seems to set grid_corners to the max nEdgesOnCell (generally 6), while the scrip_from_mpas SCRIP has grid_corners consistent with the original maxEdges dim, and these differences cause NCO to interpret them slightly differently. The mpas2esmf result also includes rrfac (regional refinement factor).

# Automatically generate a 1-degree target grid
# Generate weights with TempestRemap conservative monotone algorithm
ncremap -m map_con.nc -s scrip.nc -g target_grid.nc -G latlon=180,360 -a traave

# Apply weights, selecting a specific variable (terrain height from the static file)
ncremap -P mpasa -m map_con.nc -v ter static.nc out_con.nc
# Use the same target grid and regrid using ESMF bilinear interpolation
ncremap -m map_bil.nc -s scrip.nc -g target_grid.nc -a bilinear
ncremap -P mpasa -m map_bil.nc -v ter static.nc out_bil.nc

Notes

In practice, it may make sense to use the default “conservative” algorithm when performing conservative regridding, and the “renormalized” algorithm when performing other regridding such as bilinear interpolation or nearest-neighbor. Another consideration is whether the fields being regridded 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 structure. Few researchers can digest the unphysical values of temperature that the “conservative” option will produce in regions rife with missing values. A counter-example is fluxes, which should be physically conserved under areal-regridding. One should consider both the type of field and its conservation properties when choosing a regridding strategy.

NCO doc, Section 3.26 Regridding

Renormalization is specified using --rnr=<threshold> where the threshold ranges from 0 to 1 and indicates the fraction of given destination cell covered by valid (data not missing) source cells. Cell-times where the threshold is met preserve the mean, while others are set to missing. pyremap supports renormalization as well (example uses 0.01).

It may be possible to use NCO RRG mode to generate weights for limited-area meshes, but it is more involved than using pyremap. However, you can use NCO to apply pyremap-generated weights.

It is also possible to use TempestRemap or MOAB directly instead of through NCO.

CDO

On Casper/Derecho:

module load cdo

Global regridding with CDO is straightforward. As in the NCO example, we can use the data from Viz, aliased to grid.nc and static.nc.

# Generate weights for a 0.25-degree regular lat-lon grid (SCRIP format)
# (lon centers 0--359.75; lat centers -89.875--89.875)
cdo gencon,r1440x720 -setgrid,mpas:grid.nc -selgrid,1 grid.nc map_con.nc

# Apply weights
cdo remap,r1440x720,map_con.nc -selvar,ter -setgrid,mpas:grid.nc static.nc out_con.nc

genbil (bilinear) doesn’t support unstructured grids. But gennn (nearest neighbor) does.

convert_mpas

convert_mpas[6] provides a simple way to convert to a 0.5°x0.5° regular lat-lon grid (or other rectangular lat-lon grids).

To obtain it:

git clone https://github.com/mgduda/convert_mpas
cd convert_mpas
make
export PATH="$PWD:$PATH"

On Casper/Derecho this should just work. Now we can use the convert_mpas command[7].

Examples

Here we use the data from our global run in Run.

# Default 0.5-degree global, one file
convert_mpas x1.10242.init.nc diag.2017-09-20_12.00.00.nc
# 2-degree global, all diag files
convert_mpas x1.10242.init.nc diag.*.nc nlat=90 nlon=180
# CONUS rectangle, all diag files
echo "startlat=25
endlat=50
nlat=25
startlon=-125
endlon=-66
nlon=59
" > target_domain

convert_mpas x1.10242.init.nc diag.*.nc

Note that the start/end points are cell edges, not centers.

Recommendations

convert_mpas for a quick remap, e.g. for checking your data with tools like ncview

Otherwise:

Footnotes

  1. Or Derecho, but Derecho is usually overkill for this kind of task.

    ↩

  2. mbconvert (to the MOAB native format), mbpart (partition into pieces for parallel processing), mbtempest (generate weights using the TempestRemap library)

    ↩↩

  3. GenerateOverlapMesh (overlay the source and destination meshes and compute intersections), GenerateOfflineMap (generate weights)

    ↩

  4. -a nco_con (first-order conservative), -a nco_idw (inverse distance weighting, can extrapolate)

    ↩

  5. On Casper/Derecho:

    module load conda/latest
    mamba create -n mpas_tools -c conda-forge mpas_tools
conda activate mpas_tools
    ↩

  6. Provided by the MPAS-A lead developer, and used in the official virtual tutorial.

    ↩

  7. For the current shell session. Add an export line to your .bashrc to make it permanent.

    ↩
MPAS Tutorial
Viz