Interpolation

Previously, various tools in this package used scipy for interpolation. However, the interpolation routines in scipy are not well suited to interpolation from regular grids to MPAS meshes—they are slow and very memory intensive, particularly for large meshes.

For bilinear interpolation from a tensor lon/lat grid to an MPAS mesh, it will be faster to use the function mpas_tools.mesh.interpolation.interp_bilin() Here is an example where we define cell width for an EC mesh (see Defining an Eddy-closure Mesh), read in longitude and latitude from an MPAS mesh, and interpolate the cell widths to cell centers on the MPAS mesh.

import numpy as np
import netCDF4 as nc4
from mpas_tools.mesh.interpolation import interp_bilin

dlon = 1.
dlat = dlon
earth_radius = constants['SHR_CONST_REARTH']
nlon = int(360./dlon) + 1
nlat = int(180./dlat) + 1
lon = np.linspace(-180., 180., nlon)
lat = np.linspace(-90., 90., nlat)

cellWidth = mdt.EC_CellWidthVsLat(lat)

# broadcast cellWidth to 2D
_, cellWidth = np.meshgrid(lon, cellWidthVsLat)

ds = nc4.Dataset('base_mesh.nc', 'r+')
lonCell = ds.variables['lonCell'][:]
latCell = ds.variables['latCell'][:]

lonCell = np.mod(np.rad2deg(lonCell) + 180., 360.) - 180.
latCell = np.rad2deg(latCell)

cellWidthOnMpas = interp_bilin(lon, lat, cellWidth, lonCell, latCell)