4. Regridding Data¶
4.1. Overview¶
CDMS provides several methods for interpolating gridded data:
From one rectangular, latlon grid to another (CDMS regridder)
Between any two latlon grids (SCRIP regridder)
From one set of pressure levels to another
From one vertical (lat/level) crosssection to another vertical crosssection.
4.2. CDMS Horizontal Regrider¶
The simplest method to regrid a variable from one rectangular, lat/lon grid to another is to use the regrid function defined for variables. This function takes the target grid as an argument, and returns the variable regridded to the target grid:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16  >>> # wget "http://cdat.llnl.gov/cdat/sample_data/clt.nc"
>>> # wget "http://cdat.llnl.gov/cdat/sample_data/geos5sample.nc"
>>> import cdms2
>>> import cdat_info
>>> f1=cdms2.open("clt.nc")
>>> f2=cdms2.open("geos5sample.nc")
>>> clt=f1('clt') # Read the data
>>> clt.shape
(120, 46, 72)
>>> ozone=f2['ozone'] # Get the file variable (no data read)
>>> outgrid = ozone.getGrid() # Get the target grid
>>> cltnew = clt.regrid(outgrid)
>>> cltnew.shape
(120, 181, 360)
>>> outgrid.shape
(181, 360)

A somewhat more efficient method is to create a regridder function. This has the advantage that the mapping is created only once and can be used for multiple arrays. Also, this method can be used with data in the form of an MV2.MaskedArray. The steps in this process are:
Given an input grid and output grid, generate a regridder function.
Call the regridder function on a Numpy array, resulting in an array defined on the output grid. The regridder function can be called with any array or variable defined on the input grid.
The following example illustrates this process. The regridder function is generated at line 9, and the regridding is performed at line 10:
1 2 3 4 5 6 7 8 9 10 11 12 13  >>> # wget "http://cdat.llnl.gov/cdat/sample_data/clt.nc"
>>> # wget "http://cdat.llnl.gov/cdat/sample_data/geos5sample.nc"
>>> import cdms2
>>> from regrid2 import Regridder
>>> f = cdms2.open("clt.nc")
>>> cltf = f['clt']
>>> ingrid = cltf.getGrid()
>>> g = cdms2.open('geos5sample.nc')
>>> outgrid = g['ozone'].getGrid()
>>> regridfunc = Regridder(ingrid, outgrid)
>>> cltnew = regridfunc(cltf)
>>> f.close()
>>> g.close()

4.2.1. Notes¶
Line #3 Makes the CDMS module available.
Line #4 Makes the Regridder class available from the regrid module.
Line #5 Opens the input dataset.
Line #6 Gets the variable object named ‘clt’. No data is read.
Line #7 Gets the input grid.
Line #8 Opens a dataset to retrieve the output grid.
Line #9 The output grid is the grid associated with the variable named ‘ozone’ in dataset g. Just the grid is retrieved, not the data.
Line #10 Generates a regridder function regridfunc.
Line #11 Reads all data for variable cltf, and calls the regridder function on that data, resulting in a transient variable cltnew.
4.3. SCRIP Horizontal Regridder¶
To interpolate between grids where one or both grids is nonrectangular, CDMS provides an interface to the SCRIP regridder package developed at Los Alamos National Laboratory (https://oceans11.lanl.gov/trac/SCRIP).
Figure 3 illustrates the process:
Obtain or generate the source and target grids in SCRIP netCDF format. A CDMS grid can be written to a netCDF file, in SCRIP format, using the writeScripGrid method.
Edit the input namelist file scrip_in to reference the grids and select the method of interpolation, either conservative, bilinear, bicubic, or distanceweighted. See the SCRIP documentation for detailed instructions.
Run the scrip executable to generate a remapping file containing the transformation coefficients.
CDMS, open the remapping file and create a regridder function with the readRegridder method.
Call the regridder function on the input variable, defined on the source grid. The return value is the variable interpolated to the new grid. Note that the variable may have more than two dimensions. Also note that the input arguments to the regridder function depend on the type of regridder. For example, the bicubic interpolation has additional arguments for the gradients of the variable.
4.4. Regridding Data with SCRIP¶
Example:
Regrid data from a T42 to POP4/3 grid, using the firstorder, conservative interpolator.
In this example:
The input grid is defined in remap_grid_T42.nc.
The output grid is defined in remap_grid_POP43.nc.
The input data is variable src_array in file sampleT42Grid.nc.
The file scrip_in has contents:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16  >>> &remap_inputs
>>> num_maps = 1
>>>
>>> grid1_file = 'remap_grid_T42.nc'
>>> grid2_file = 'remap_grid_POP43.nc'
>>> interp_file1 = 'rmp_T42_to_POP43_conserv.nc'
>>> interp_file2 = 'rmp_POP43_to_T42_conserv.nc'
>>> map1_name = 'T42 to POP43 Conservative Mapping'
>>> map2_name = 'POP43 to T42 Conservative Mapping'
>>> map_method = 'conservative'
>>> normalize_opt = 'frac'
>>> output_opt = 'scrip'
>>> restrict_type = 'latitude'
>>> num_srch_bins = 90
>>> luse_grid1_area = .false.
>>> luse_grid2_area = .false.

num_maps
specifies the number of mappings generated, either 1 or 2.
For a single mapping, grid1_file
and grid2_file
are the source
and target grid definitions, respectively. The map_method
specifies
the type of interpolation, either ‘conservative’, ‘bilinear’, ‘bicubic’,
or ‘distwgt’ (distanceweighted). The remaining parameters are described
in the SCRIP documentation.
Once the grids and input file are defined, run the scrip executable to generate the remapping file ‘rmp_T42_to_POP43_conserv.nc’
1 2 3 4 5 6 7 8 9  >>> % scrip
>>> Using latitude bins to restrict search.
>>> Computing remappings between:
>>> T42 Gaussian Grid
>>> and
>>> POP 4/3 DisplacedPole T grid
>>> grid1 sweep
>>> grid2 sweep
>>> Total number of links = 63112

Next, run CDAT and create the regridder:
1 2 3 4 5 6 7 8 9 10 11  >>> # wget "http://cdat.llnl.gov/cdat/sample_data/remap_grid_POP43.nc"
>>> # wget "http://cdat.llnl.gov/cdat/sample_data/remap_grid_T42.nc"
>>> # wget "http://cdat.llnl.gov/cdat/sample_data/rmp_POP43_to_T42_conserv.nc"
>>> # wget "http://cdat.llnl.gov/cdat/sample_data/rmp_T42_to_POP43_conserv.nc"
>>> # wget "http://cdat.llnl.gov/cdat/sample_data/xieArkinT42.nc"
>>> # Import regrid package for regridder functions
>>> import regrid2, cdms2
>>> # Read the regridder from the remapper file
>>> remapf = cdms2.open('rmp_T42_to_POP43_conserv.nc')
>>> regridf = regrid2.readRegridder(remapf)
>>> remapf.close()

Then read the input data and regrid:
1 2 3 4 5 6  >>> # Get the source variable
>>> f = cdms2.open('xieArkinT42.nc')
>>> t42prc = f('prc')
>>> f.close()
>>> # Regrid the source variable
>>> popdat = regridf(t42prc)

Note that t42dat
can have rank greater than 2. The trailing
dimensions must match the input grid shape. For example, if t42dat
has shape (12, 64, 128), then the input grid must have shape (64,128).
Similarly if the variable had a generic grid with shape (8092,), the
last dimension of the variable would have length 8092.
4.5. PressureLevel Regridder¶
To regrid a variable which is a function of latitude, longitude,
pressure level, and (optionally) time to a new set of pressure levels,
use the pressureRegrid
function defined for variables. This function
takes an axis representing the target set of pressure levels, and
returns a new variable d
regridded to that dimension.
1 2 3 4 5 6 7 8 9 10  >>> # wget "http://cdat.llnl.gov/cdat/sample_data/ta_ncep_876884.nc"
>>> f=cdms2.open("ta_ncep_876884.nc")
>>> ta=f('ta')
>>> ta.shape
(11, 17, 73, 144)
>>> ta.getAxisIds()
['time', 'level', 'latitude', 'longitude']
>>> result = ta.pressureRegrid(cdms2.createAxis([1000.0]))
>>> result.shape
(11, 1, 73, 144)

4.6. CrossSection Regridder¶
To regrid a variable which is a function of latitude, height, and
(optionally) time to a new latitude/height crosssection, use the
crossSectionRegridder
defined for variables. This function takes as
arguments the new latitudes and heights, and returns the variable
regridded to those axes.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15  >>> # wget "http://cdat.llnl.gov/cdat/sample_data/ta_ncep_876884.nc"
>>> f=cdms2.open("ta_ncep_876884.nc")
>>> ta=f('ta')
>>> ta.shape
(11, 17, 73, 144)
>>> levOut=cdms2.createAxis([1000.0,950.])
>>> levOut.designateLevel()
>>> latOut=cdms2.createAxis(ta.getLatitude()[10:20])
>>> latOut.designateLatitude()
>>> ta0 = ta[0,:]
>>> ta0.getAxisIds()
['level', 'latitude', 'longitude']
>>> taout = ta0.crossSectionRegrid(levOut, latOut)
>>> taout.shape
(2, 10, 144)

4.7. Regrid Module¶
The regrid
module implements the CDMS regridding functionality as
well as the SCRIP interface. Although this module is not strictly a part
of CDMS, it is designed to work with CDMS objects.
4.8. CDMS Horizontal Regridder¶
from regrid2 import Regridder
Makes the CDMS Regridder class available within a Python program. An instance of Regridder is a function which regrids data from rectangular input to output grids.
4.8.1. CDMS Regridder Constructor¶
Constructor 
Description 



4.9. SCRIP Regridder¶
SCRIP regridder functions are created with the regrid.readRegridder
function:
4.9.1. SCRIP Regridder Constructor¶
Constructor 
Description 



4.10. Regridder Functions¶
It is only necessary to specify the map method if it is not defined in the file.
If checkGrid
is 1 (default), the grid cells are checked for
convexity, and ‘repaired’ if necessary. Grid cells may appear to be
nonconvex if they cross a 0 / 2pi
boundary. The repair consists of
shifting the cell vertices to the same side modulo 360 degrees.
4.11. CDMS Regridder Functions¶
A CDMS regridder function is an instance of the CDMS Regridder
class. The function is associated with rectangular input and output
grids. Typically its use is straightforward:
The function is passed an input array and returns the regridded array. However, when the array has missing data, or the input and/or output grids are masked, the logic becomes more complicated.
4.11.1. Step 1¶
The regridder function first forms an input mask. This mask is either twodimensional or ndimensional, depending on the rank of the usersupplied mask. If no mask or missing value is specified, the mask is obtained from the data array mask if present.
Twodimensional case:
Let mask_1 be the twodimensional user mask supplied via the mask argument, or the mask of the input grid if no user mask is specified.
If a missingdata value is specified via the missing argument, let the implicit_mask be the twodimensional mask defined as 0 where the first horizontal slice of the input array is missing, 1 elsewhere.
The input mask is the logical AND(mask_1, implicit_mask)
Ndimensional case:
If the user mask is 3 or 4dimensional with the same shape as the input array, it is used as the input mask.
4.11.2. Step 2¶
The data is then regridded. In the twodimensional case, the input mask is ‘broadcast’ across the other dimensions of the array. In other words, it assumes that all horizontal slices of the array have the same mask. The result is a new array, defined on the output grid. Optionally, the regridder function can also return an array having the same shape as the output array, defining the fractional area of the output array which overlaps a nonmissing input grid cell. This is useful for calculating areaweighted means of masked data.
4.11.3. Step 3¶
Finally, if the output grid has a mask, it is applied to the result array. Where the output mask is 0, data values are set to the missing data value, or 1.0e20 if undefined. The result array or transient variable will have a mask value of 1 (invalid value) for those output grid cells which completely overlap input grid cells with missing values
4.11.4. CDMS Regridder Function¶
Type 
Function 
Description 

Array or TransientVariable 


4.11.5. DMS Regridder Function(cont’d)¶
Type 
Function 
Description 

Array or TransientVariable 


Array, Array 

If called with the optional

4.12. SCRIP Regridder Functions¶
A SCRIP regridder function is an instance of the ScripRegridder class. Such a function is created by calling the regrid.readRegridder method. Typical usage is straightforward:
1 2 3 4 5 6 7 8 9  >>> import cdms2
>>> import regrid2
>>> remapf = cdms2.open('rmp_T42_to_POP43_conserv.nc')
>>> regridf = regrid2.readRegridder(remapf)
>>> f = cdms2.open('xieArkinT42.nc')
>>> t42prc = f('prc')
>>> f.close()
>>> # Regrid the source variable
>>> popdat = regridf(t42prc)

The bicubic regridder takes four arguments:
>>> # outdat = regridf(t42prc, gradlat, gradlon, gradlatlon)
A regridder function also has associated methods to retrieve the following fields:
Input grid
Output grid
Source fraction: the fraction of each source (input) grid cell participating in the interpolation.
Destination fraction: the fraction of each destination (output) grid cell participating in the interpolation.
In addition, a conservative regridder has the associated grid cell areas for source and target grids.
4.12.1. SCRIP Regridder Functions¶
Return Type 
Method 
Description 

Array or TransientVariable 
[conservative, bilinear, and distanceweighted regridders] 

Array or TransientVariable 
[bicubic regridders] 

4.12.2. SCRIP Regridder Functions(con’td)¶
Return Type 
Method 
Description 

Numpy array 


Numpy array 


CurveGrid or GenericGrid 

Return the input grid, or None if no input grid is associated with the regridder. 
CurveGrid or GenericGrid 

Return the output grid. 
Numpy array 


4.13. Examples¶
4.13.1. CDMS Regridder¶
Example:
Regrid data to a uniform output grid.
1 2 3 4 5 6 7 8 9  >>> import cdms2
>>> from regrid2 import Regridder
>>> f = cdms2.open('clt.nc')
>>> cltf = f.variables['clt']
>>> ingrid = cltf.getGrid()
>>> outgrid = cdms2.createUniformGrid(90.0, 46, 4.0, 0.0, 72, 5.0)
>>> regridFunc = Regridder(ingrid, outgrid)
>>> newrls = regridFunc(cltf)
>>> f.close()

4.13.2. Regridder Constructure¶
Line 
Notes 

3 
Open a netCDF file for input. 
6 
Create a 4 x 5 degree output grid. Note that this grid is not associated with a file or dataset. 
7 
Create the regridder function. 
8 
Read all data and regrid. The missing data value is obtained from variable rlsf 
Return the area fraction of the source (input) grid cell that participates in the regridding. The array is 1D, with length equal to the number of cells in the input grid.
Example:
Get a mask from a separate file, and set as the input grid mask.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17  >>> # wget http://cdat.llnl.gov/cdat/sample_data/clt.nc
>>> # wget http://cdat.llnl.gov/cdat/sample_data/geos5sample.nc
>>> import cdms2
>>> from regrid2 import Regridder
>>> #
>>> f = cdms2.open('clt.nc')
>>> cltf = f.variables['clt']
>>> outgrid = cltf.getGrid()
>>> g = cdms2.open('geos5sample.nc')
>>> ozoneg = g.variables['ozone']
>>> ingrid = ozoneg.getGrid()
>>> regridFunc = Regridder(ingrid,outgrid)
>>> uwmaskvar = g.variables['uwnd']
>>> uwmask = uwmaskvar[:]<0
>>> outArray = regridFunc(ozoneg.subSlice(time=0),mask=uwmask)
>>> f.close()
>>> g.close()

Line 
Notes 

7 
Get the input grid. 
10 
Get the output grid. 
11 
Create the regridder function. 
14 
Get the mask. 
15 
Regrid with a user mask. The subslice call returns a transient variable corresponding to variables of at time 0. 
Note: Although it cannot be determined from the code, both mask and the input arrays of are fourdimensional. This is the ndimensional case.
Example:
Generate an array of zonal mean values.
1 2 3 4 5 6 7  >>> f = cdms2.open(‘rls_ccc_per.nc’)
>>> rlsf = f.variables[‘rls’]
>>> ingrid = rlsf.getGrid()
>>> outgrid = cdms2.createZonalGrid(ingrid)
>>> regridFunc = Regridder(ingrid,outgrid)
>>> mean = regridFunc(rlsf)
>>> f.close()

Line 
Notes 

3 
Open a netCDF file for inputGet the input grid. Return the area fraction of the source (input) grid cell that participates in the regridding. The array is 1D, with length equal to the number of cells in the input grid. 
4 
Create a zonal grid. outgrid has the same latitudes as ingrid, and a singleton longitude dimension. createGlobalMeanGrid could be used here to generate a global mean array. 
5 
Generate the regridder function. 
6 
Generate the zonal mean array. 
Example:
Regrid an array with missing data, and calculate the areaweighted mean of the result.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17  >>> import cdms2
>>> from cdms2.MV2 import *
>>> from regrid2 import Regridder
>>> f = cdms2.open("ta_ncep_876884.nc")
>>> var = f('ta')
>>> outgrid = cdms2.createUniformGrid(90.0, 46, 4.0, 0.0, 72, 5.0)
>>> outlatw, outlonw = outgrid.getWeights()
>>> outweights = outerproduct(outlatw, outlonw)
>>> grid = var.getGrid()
>>> sample = var[0,0]
>>> latw, lonw = grid.getWeights()
>>> weights = outerproduct(latw, lonw)
>>> inmask = where(greater(absolute(sample),1.e15),0,1)
>>> mean = add.reduce(ravel(inmask*weights*sample))/add.reduce(ravel(inmask*weights))
>>> regridFunc = Regridder(grid, outgrid)
>>> outsample, outmask = regridFunc(sample, mask=inmask, returnTuple=1)
>>> outmean = add.reduce(ravel(outmask*outweights*outsample)) / add.reduce(ravel(outmask*outweights))

Line 
Notes 

2 
Create a uniform target grid. 
3 
Get the latitude and longitude weights. 
4 
Generate a 2D weights array. 
5 
Get the input grid. 
6 
Get the first horizontal slice from 
78 
Get the input weights, and generate a 2D weights array. 
9 
Set the 2D input mask. 
10 
Calculate the input array areaweighted mean. 
11 
Create the regridder function. 
12 
Regrid. Because returnTuple is set to 1, the result is a tuple (dataArray, maskArray). 
13 
Calculate the areaweighted mean of the regridded data. mean and outmean should be approximately equal. 
4.13.3. SCRIP Regridder¶
Example:
Regrid from a curvilinear to a generic grid, using a conservative remapping. Compute the areaweighted means on input and output for comparison.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28  >>> # wget "http://cdat.llnl.gov/cdat/sample_data/remap_grid_T42.nc"
>>> # wget http://cdat.llnl.gov/cdat/sample_data/rmp_T42_to_C02562_conserv.nc
>>> # wget "http://cdat.llnl.gov/cdat/sample_data/xieArkinT42.nc"
>>> import cdms2, regrid2, MV2
>>> # Open the SCRIP remapping file and data file
>>> fremap = cdms2.open('rmp_T42_to_C02562_conserv.nc')
>>> fdat = cdms2.open('xieArkinT42.nc')
>>> # Input data array
>>> dat = fdat('prc')[0,:]
>>> # Read the SCRIP regridder
>>> regridf = regrid2.readRegridder(fremap)
>>> # Regrid the variable
>>> outdat = regridf(dat)
>>> # Get the cell area and fraction arrays. Areas are computed only
>>> # for conservative regridding.
>>> srcfrac = regridf.getSourceFraction()
>>> srcarea = regridf.getSourceArea()
>>> dstfrac = regridf.getDestinationFraction()
>>> dstarea = regridf.getDestinationArea()
>>> # calculate areaweighted means
>>> inmean = MV2.sum(srcfrac*srcarea*MV2.ravel(dat)) / MV2.sum(srcfrac*srcarea)
>>> outmean = MV2.sum(dstfrac*dstarea*MV2.ravel(outdat)) / MV2.sum(dstfrac*dstarea)
>>> print 'Input mean:', inmean
Input mean: 2.60376502339
>>> print 'Output mean:', outmean
Output mean: 2.60376502339
>>> fremap.close()
>>> fdat.close()
