| re.gex - A 2D regridding function for GrADS |
re.gex - A 2D regridding function for GrADS
re(expr,dlon), while dlon=dlat
re(expr,dlon,dlat,['ig',nyig],['ba'|'bl'|'bs'|'vt',vtmax,vtmin|'ma',min]
re(expr,dlon,gYY,['ig',nyig],['ba'|'bl'|'bs'|'vt',vtmax,vtmin|'ma',min]
re(expr,nx,'linear',lon,dlon,ny,'linear',lat,dlat, ['ig',nyig],['ba'|'bl'|'bs'|'vt',vtmax,vtmin|'ma',min]
re(expr,nx,'linear',lon,dlon,ny,'gaus',gstart,njog, ['ig',nyig],['ba'|'bl'|'bs'|'vt',vtmax,vtmin|'ma',min]
The regrid function re solves a common problem of transforming
horizontal 2-D gridded fields from/to different resolutions/grid
types for quantitative intercomparison. For example, a model monthly
mean precipitation field on a T126 gaussian grid can be compared to an
observed climatology on a 2.5x2.5 grid using re. The function re
offers many transform options ranging from simple bilinear
interpolation to box averaging with voting. Additional methods can
be added to re as needed.
re transforms two-dimensional (2-D) lat/lon GrADS grids from
one grid type/resolution to another. The input is any 2-D
lat/lon grid defined by the current GrADS lat/lon dimension
environment. re handles input grids which are cyclically
continuous in longitude and excludes undefined input grid values
from participation in the transform. If a valid transform cannot
be made (i.e., insufficient defined data points), the output grid
is set to undefined. re supports two output grid types: 1)
lat/lon with fixed increments; and 2) gaussian. Four transforms
are available: 1) box averaging for regridding fine to coarse
grids; 2) box averaging with ``voting'' for noncontinuous/index
data such, as soil type; 3) bilinear interpolation; and 4)
4-point bessel interpolation.
Any valid GrADS grid expression (e.g., z or ave(z.3(t+0,t=120,1yr)), real number, 'undef' or '-u'.
The number of points in longitude (integer)
Beginning longitude (center of the lower left hand corner grid or the grid (1,1)) of the output domain (float)
Delta longitude (dlon) or number of gaussian longitudes on the GLOBE (float)
The number of points in latitude (integer)
Beginning latitude (center of the lower left hand corner grid or the grid (1,1)) of the output domain (float)
Delta latitude (dlat) or the number of gaussian latitudes on the GLOBE (float)
The first gaussian grid number. If the data span all latitudes, start would be 1, indicating the southernmost gaussian grid latitude (integer)
The number of GLOBAL gaussian latitudes on the output grid. (integer)
Linear mapping (string)
Gaussian latitide mapping (string)
All strings are CASE INSENSITIVE.
Input grid is gaussian with nyig being the number of gaussian latitudes (e.g., ig92 for the NMC T62 grid). nyig must be >= 8 and a multiple of four. This parameter is used to invoke a more precise calculation of the boundary between gaussian grid boxes.
Box averaging (the default, while regrids to coarse resoultion)
Bi-linear interpolation (the default, when regridding to a finer resolution)
3rd order Bessel interpolation
Vote interpolation or box averaging with voting. The parameters (vtmax,vtmin) (range: [0-1]) set the fraction of an output grid box that must be covered by defined input grid data for a ``winner'' to be chosen in the election. The default is vtmax=vtmin=1/2.
The parameter vtmin must be the same as vtmax except for three or more candidates. The fraction for two candidates is midway between vtmax and vtmin.
When there is only one candidate, vtmax is the minimum fraction of an output grid point hat must be covered by defined input grid data for a ``winner'' to be chosen in the election.
Specifying vtmax = vtmin = 1 would require that 100% of the output grid box must be covered by a single, unique value from the input grid whereas vtmax = vtmin = 0 would allow a winner to be chosen if ANY candidates where running. The default value of 0.5 means that a simple majority is required to have a winner.
This option applies ONLY to box averaging without voting when the input grid has undefined points. The parameter fraction (range: [0-1]) specifies the minimum area which must be covered with DEFINED data boxes to be considered a valid interpolation. The old regrid v1.0 assumed fraction was 0 or that if ANY input grid boxes contained defined data which intersected the output grid produced a box average. This was clearly too liberal and fraction is now set by default to 50% or that half the output grid box must be covered with defined data to produced a defined output grid point.
Regrid a global T62 gaussian grid (192x94) to a 2.5 deg lat/lon by box averaging,
open /reanl1/pilot20/fluxgrb8508.ctl set x 1 192 set y 1 94 define p25=re(p,144,linear,0,2.5,72,linear,-88.75,2.5,ba)
or set lon 0 360 set lat -90 90 define p25=re(p,2.5,2.5,ba)
or more simply,
define p25=re(p,2.5)
Note: The lat/lon dimension environment is set using grid coordinates (x,y) to make the input and output grids global. To minimize future graphics/analysis calculations with the regridded field, we use the GrADS define function to store the grid in memory where it can be referenced as any other GrADS grid.
Regrid a 4x5 SiB vegetation type to a R15 (48x40) gaussian grid using box averaging with ``voting.'' Require that at least 60% of the output grid box must be covered with a single candidate value from the input grid for an election to occur. Otherwise the output grid box is set to undefined. Relax the one-candidate election requirement to 20% when there are three or more candidates,
open /export/sgi18/fiorino/data/sib/sib.param.annual.ctl set lon 0 360 set lat -90 90 define i21=re(index,48,linear,0,7.5,40,gaus,1,40,vt,0.60,0.20) set gxout grfill d index d i21
Note : During candidate selection, undefined input grid points do not contribute to the fraction of the output grid box covered with input grid boxes. The best way to display index type data in GrADS is to use the ``grid fill'' display option (set gxout grfill). GrADS will draw the grid box and color it according to the grid value.
Regrid 1x1 Aviation run 500 mb z to 2.5x2.5 grid for the region (-140, 20) to (-40, 70) using bessel interpolation,
open /export/sgi39/wd22sl/grads/avn/avn.93092800.fcst.ctl set lev 500 set lat -180 180 set lon -90 90
d re(z,40,linear,-138.75,2.5,20,linear,21.25,2.5,bs)
or
set lat 20 70 set lon -140 -40
d re(z,40,linear,-138.75,2.5,20,linear,21.25,2.5,bs)
or
d re(z,2.5,2.5,bs)
Note: The above three regrid commands produce exactly the same results. Box averaging would be more appropriate when regridding to a coarser grid.
Regrid 1x1 Aviation run 500 mb z to 2.5x2.5 grid using box averaging and setting the grid to start at a specific lat/lon,
open /export/sgi39/wd22sl/grads/avn/avn.93092800.fcst.ctl set lev 500 set lat -20 70 set lon -140 -40 d re(z,40,linear,-138.75,2.5,20,linear,21.25,2.5,ba) set lat 30 50 set lon -50 50 d re(z,40,linear,-138.75,2.5,20,linear,21.25,2.5,ba)
Note: The above two regrids produce DIFFERENT results since the input domain does not cover the entire output domain. Missing values will be filled for the uncovered regions.
There is no restriction in the dimension of input/output grids while there is sufficient memory. Note that there was a restriction of input/output grids of dimension 730x380 (~T225) in version 2.0 and earlier.
Any valid GrADS grid can be regridded. However, GrADS (V1.5)
currently only supports lat/lon grids where the mapping between
grid and world coordinates is one dimensional, i.e.,
longitude(i,j)=f(i) vice longitude(i,j)=f(i,j).
Only two output grid types have been implemented: 1) lat/lon
(dlat does NOT have to equal dlon); and 2) gaussian grids.
Mercator output grids could be added as lon(i,j)=f(i) and
lat(i,j)=f(j) in this projection.
The first step in the regrid transform is to define a relationship between the input and output grids within a common frame of reference or coordinate system. regrid bases the inter-grid relationship on ``world'' coordinates, and the GrADS map between grid coordinates (i,j) and world coordinates (lon, lat). As noted above, the world-grid mapping in GrADS is one-dimensional. Thus, the world-grid map of an input GrADS grid takes the form,
lat(i,j)=f(j) and lon(i,j)=g(i).
By specifying a similar mapping for an output GrADS grid of the form
LAT(I,J)=F(J) and LON(I,J)=G(I),
as defined by the input parameters X1, X2 and X3-6, regrid calculates,
X(I)=i(G(I)) and Y(J)=j(F(J)),
where i(G(I)) is the location of the output grid with respect to
the input grid dimension i and j(F(J)) for j.
For simplicity, and greater generality, regrid assumes that the grid point is at the center of a rectangular grid box and maps the location of the boundaries of the output grid box to that of the input grid box. By default the boundaries are assumed to lie midway between grid points and while this is not strictly true for a gaussian grid near the poles, it is close nonetheless. The boundaries for gaussian grids can be calculated by specifying ig XXX in options. The reason why this cannot be automatic is that GrADS does not directly support gaussian grids (i.e., there is no ydef gauss 40 option in the data descriptor .ctl file, just linear and levels).
Given the inter-grid map X(I) and Y(J), regrid uses two basic
methods for doing the transform: 1) box averaging; or 2)
interpolation. Box averaging is simply the area-weighted
integral of all input grid boxes which intersect an output grid
box, divided by the area of the output grid box. This approach
is most appropriate: 1) for transforming from fine (e.g., dlon =
1 deg) to coarse grids (e.g., dlon = 2.5 deg); and 2) when
approximate conservation of an integral quantity (e.g., global
average) is desired.
Box averaging is also useful for regridding noncontinuous,
parametric or ``index'' data. For example, suppose you have been
given a 0.5x0.5 deg global grid of vegetation type and want to
use these data in an R43 global model. The intuitive solution is
to assign the output grid the value of the intersecting input
grid box(es) which account(s) for the greatest percentage of the
output grid box surface area. In the example of vegetation data,
if 70% of the output grid box is covered by deciduous forest,
then it might be reasonable to call the output grid deciduous
forest. However, if there were 5 distinct vegetation types or
``candidates'' available, then regrid, being an American function,
holds an ``election'' and select a ``winner'' based on the candidate
covering the greatest percentage of the total surface area in the
output grid box. Of course, coming from an imperfect democracy,
the election can be ``rigged'' for a desired effect....
This grid transform strategy is invoked using the ``vote'' option in box averaging (vt in C1). Conditions on the percentage of the output grid box (number of candidates and what it takes to get elected) can be finely controlled by the X4 and X5 parameters.
Perhaps the most conventional way of regridding meteorological data (e.g., 500 mb geopotential heights) is interpolation because weather data tend to be continuous . regrid features a 4x4 point bessel interpolation routine developed at Fleet Numerical Meteorology and Oceanography Center (courtesy of D. Hensen, FNMOC). While this routine is in wide use at FNMOC, the regrid implementation has been substantially improved to handle more general input grids.
First, bilinear interpolation is performed at all points to produce a ``first guess.'' Improvements to the bilinear ``first guess'' are made using the higher-order terms in the bessel interpolator, but only if the bessel option is set (i.e., bs in options). Second, an undefined value in the 2x2 bilinear stencil causes the output grid to be undefined. If the bilinear value is defined, but any of the points in the larger 4x4 bessel stencil are undefined, the output grid is assigned the bilinear value. The third improvement is that the routine handles grids which are cyclically continuous in longitude.
It is generally recommended that you use the bessel option when interpolating because the higher-order terms in the polynomial interpolation function produce a closer fit to the original data in regions of rapid changes (e.g., large gradients of the gradient around low pressure centers).
By default, the box averaging is used while the resolution of input grid is finer than the out grid. Otherwise, the bessel interlopation is used.
The only down side to a regridded field is that its dimension environment cannot be controlled by its ``grid'' coordinate system. The best way to understand this is by an example. Suppose you regrid a T62 global Gaussian grid (192x94) to a uniform 2.5 deg grid (144x73) using box averaging and the GrADS define capability, e.g.,
define p25=re(p,2.5,2.5,ba)
You now want to calculate the global average of the original field p and the defined regridded field p25. The only unambiguous way (using all grid boxes) of doing this calculation for p would be,
d aave(p,x=1,x=192,y=1,y=94)
and not,
d aave(p,lon=0,lon=360,lat=-90,lat=90)
This is because the cyclic continuity feature in GrADS would count grid boxes at the prime meridian twice, i.e., GrADS would really be doing,
d aave(p,x=1,x=193,y=1,y=94)
Trying to find the global average of the 2.5 deg regridded field p25 using,
d aave(p25,x=1,x=144,y=1,y=73)
would not yield a global average even though p25 IS 144x73! However,
d aave(p25,x=1,x=192,y=1,y=94)
would because GrADS converts the grid coordinate range to (x=1,x=192) to world coordinates (lon=0,lon=360-1.875) and grabs all grid boxes in p25 within that range when putting together the data for the areal averaging calculation. Despite this restriction on grid coordinates, you can access specific chunks of a regridded defined variable using world coordinates. Suppose you want to look at the latitudinal variation of the u wind component at 180 deg and 500 mb on a 1 deg grid, e.g.,
set lev 500 set lon 180 set lat -90 90 d u
if the you had,
define u5=regrid2(u,5)
you could then,
d u5(lon=175)
but not,
d u5(x=1)
regrid sends information on the transform process (starting/ending lat/lon, number of grid points and the regridding method) to the terminal window where you are running GrADS. Additionally, errors during the call to regrid (e.g., trying to regrid a two-dimensional longitude-pressure cross section) will be displayed, the process terminated, and control returned back to GrADS.
Regrid, a Grid Analysis and Display System (GrADS) function developed originally for the Development Division of the former National Meteorological Center (now NOAA's National Centers for Environmental Predection, NCEP), was substantially improved by Mike Fiorino at the Program for Climate Model Diagnosis and Intercomparison (PCMDI), Lawrence Livermore National Laboratory, from 1995 to about 2003. Fiorino introduced a simpler calling syntax and made it available through the Internet as a GrADS User Defined Function (UDF).
Starting in 2003 Ben-Jei Tsuang converted the original regrid from
Fortran to C and implemented it as an intrinsic function in GrADS
v1.8, making his source code patches available to the
community. However, this regrid patch was not adopted by COLA for
GrADS v1.9. The porting was done with the utility f2c to convert
original f77 code to C, and the resulting C code was modified to
eliminate the dependency on the f2c library. Among the improvements
Tsuang made to regrid are:
Removed the restriction on the sizes of input/output dimensions.
The function was renamed from regrid() to re(), and the
arguments were reorganized to a syntax that closely followed the
conventions used by the standard GrADS file descriptor. In this new
syntax, the size and the locations of output grids can be exactly
described.
The function is now aware of the longitudinal cyclic continuity in GrADS
(where the first and last longitudinal grid point is repeated), and does
not generate an extra-x column as its predecessor regrid.
This feature is useful when using fwrite to create binary files.
The output domain can be smaller or larger than the input domain. If the input domain is smaller than the output domain missing values will be filled for the uncovered regions.
Due to the cyclic continuity feature in GrADS, you may not see your
results correctly on the screen when using display. It is strongly
suggested to use lats4d or fwrite to write your results to a
disk file, and then visualize your results afterwards.
For convenience, the function writes out XDEF, YDEF and UNDEF records
which can be cut-and-pasted to create a ctl file when using this
function in conjunction with fwrite or lats4d.
In 2006 Arlindo da Silva implemented re() as dynamic user defined
function and it became a poster child for the new User Defined Extensions
(UDXTs) in GrADS.
Copyright (c) 1995-2003 by Mike Fiorino <mfiorino@gmail.com>
Copyright (c) 2003-2007 by Ben-Jei Tsuang <btsuang@yahoo.com>
This is free software; see the source for copying conditions. There is NO warranty; not even for MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
| re.gex - A 2D regridding function for GrADS |