Use PDEF For Displaying Pre-Projected Data With GrADS

Display pre-projected data with PDEF
PDEF Syntax
How grid interpolation works
How wind rotation works
PDEF BILIN option
PDEF FILE option

Display Pre-Projected Data with PDEF

Gridded data that are mapped onto a particular map projection are called 'pre-projected.' An example of pre-projected data is the output from a weather forecast model that is mapped onto a north polar stereographic grid.

In order to display pre-projected data on a map in GrADS, the descriptor file must contain a PDEF entry. A descriptor file that contains a PDEF record describes two different grids. The first grid is described by the PDEF record itself and is the "native" grid for the pre-projected data in the file. The PDEF record describes the size of this native grid, and then describes how to convert from lat/lon to the i/j of the grid (this description takes the form of stating the projection of the grid along with the projection constants or providing the mapping to the native grid in a supplementary data file). The second grid described in the desctiptor file is a rectilinear lat/lon grid, which is defined by the XDEF and YDEF records. The rectilinear grid is used by GrADS internally and can be any size or resolution -- it is completely independent of the pre-projected or native grid. GrADS uses the information about the two grids to interpolate from the PDEF-described native grid to the XDEF/ YDEF-described rectilinear grid. All displays and analyses are done using the interpolated data on the rectilinear grid. The virtue of this approach is that all built in GrADS analytic functions (e.g., aave, hcurl...) will work even though the data were not originally on a lon/lat grid. The downside is that you are looking at interpolated data.

It is possible to view the pre-projected data on its native grid. To do this, you omit the PDEF entry from the descriptor file, and use the XDEF and YDEF entries to describe the shape of the native grid. In this case, your displays must be drawn in i/j space without a map projection (set mpdraw off).

When you do a display of a pre-projected vector field, you must know whether the original vector components are defined relative to the data grid or relative to the Earth. If the data are grid-relative, they must be rotated to Earth-relative coordinates in order for the interpolation to work properly. See the "Notes" under each particular projection type for further information.

PDEF Syntax

PDEF isize jsize NPS ipole jpole lonref gridinc PDEF isize jsize SPS ipole jpole lonref gridinc
Example: | PDEF 53 45 nps 27 49 -105 190.5

Args:

isize | The size of the native grid in the x direction
jsize | The size of the native grid in the y direction
ipole | the i coordinate of the pole referenced from the lower left corner, assumed to be at (1,1)
jpole | the j coordinate of the pole referenced from the lower left corner, assumed to be at (1,1)
lonref | reference longitude
gridinc | distance between gripoints in km

Notes: | Polar stereographic projections (N and S) are defined as at NCEP. Wind rotation has also been added so that vector data will be properly displayed.
PDEF isize jsize LCCR latref lonref iref jref Struelat Ntruelat slon dx dy PDEF isize jsize LCC latref lonref iref jref Struelat Ntruelat slon dx dy
Example: | PDEF 103 69 lccr 30 -88 51.5 34.5 20 40 -88 90000 90000

Args:

isize | The size of the native grid in the x direction
jsize | The size of the native grid in the y direction
latref | reference latitude
lonref | reference longitude (in degrees, E is positive, W is negative)
iref | i of ref point
jref | j of ref point
Struelat | S true lat
Ntruelat | N true lat
slon | standard longitude
dx | grid X increment in meters
dy | grid Y increment in meters

Notes: | Starting with version 1.9b4, the LCCR option supplements the use of PDEF with data on the Lambert Conformal projection. With LCCR, wind rotation has been implemented for data files with grid-relative winds instead of Earth-relative winds. Use LCC if your vector components are already Earth-relative.
PDEF isize jsize ETA.U lonref latref dlon dlat
Example: | PDEF 181 136 eta.u -97.0 41.0 0.38888888 0.37

Args:

isize | The size of the native grid in the x direction
jsize | The size of the native grid in the y direction
lonref | reference longitude (in degrees, E is positive, W is negative)
latref | reference latitude
dlon | grid longitude increment in degrees
dlat | grid latitude increment in degrees

Notes: | The eta model native grid is awkward to work with because the variables are on staggered and non-rectangular grids. NCEP created "unstaggered" eta model fields, in which the variables are placed on a common rectangular grid. Wind rotation has also been added so that vector data will be properly displayed.
PDEF isize jsize PSE slat slon ipole jpole dx dy sign
Example: |

Args:

isize | The size of the native grid in the x direction
jsize | The size of the native grid in the y direction
slat | absolute value of the standard latitude
slon | absolute value of the standard longitude
ipole | the i coordinate of the pole referenced from the lower left corner, assumed to be at (0,0)
jpole | the j coordinate of the pole referenced from the lower left corner, assumed to be at (0,0)
dx | grid X increment in km
dy | grid Y increment in km
sign | 1 for NH; -1 for SH

Notes: | The polar stereo projection used by the original NMC models is not very precise because it assumes the earth is round (eccentricity = 0). While this approximation was reasonable for coarse resolution NWP models, it is inadequate to work with higher resolution data such as SSM/I. Wind rotation has not been implemented!!! Use only for scalar fields.
PDEF isize jsize OPS latref lonref xoff yoff iref jref dx dy
Example: | PDEF 26 16 ops 40.0 -100.0 90000.0 90000.0 14.0 9.0 180000.0 180000.0

Args:

isize | The size of the native grid in the x direction
jsize | The size of the native grid in the y direction
latref | reference latitude
lonref | reference longitude (in degrees, E is positive, W is negative)
xoff | lonref offset in meters
yoff | latref offset in meters
iref | the i coordinate of the reference point
jref | the j coordinate of the reference point
dx | grid X increment in km
dy | grid Y increment in km
dy | grid Y increment in km

Notes: | The CSU RAMS model uses an oblique polar stereo projection. Wind rotation has not been implemented!!! Use only for scalar fields.
PDEF isize jsize BILIN format byteorder fname
Example: | PDEF 100 100 BILIN sequential binary-big ^mygrid.interp.values

Args:

isize | The size of the native grid in the x direction
jsize | The size of the native grid in the y direction
format | Must be either STREAM (direct access) or SEQENTIAL (fortran formatted)
byteorder | If set to BINARY, byte odering is assumed to be same as local machine If set to BINARY-BIG, byte ordering is assumed to be big-endian If set to BINARY-LITTLE, byte ordering is assumed to be little-endian
fname | The name of the supplementary file

Notes: | The supplementary file contains three lat-lon floating-point grids: i values, j values, and wind rotation values. The size of these grids must match the XDEF and YDEF entries in the descriptor file.
PDEF size 1 FILE num format byteorder fname
Example: | PDEF 5000 1 file 4 sequential binary-big ^mygrid.interp.values PDEF 15238 1 file 1 stream binary ^gtd.filepdef

Args:

size | size of the native grid, expressed as a vector (eg, all gridpoints in an x-y grid)
num | number of sets of interpolation grids supplied
format | Must be either STREAM (direct access) or SEQENTIAL (fortran formatted)
byteorder | If set to BINARY, byte odering is assumed to be same as local machine If set to BINARY-BIG, byte ordering is assumed to be big-endian If set to BINARY-LITTLE, byte ordering is assumed to be little-endian
fname | The name of the supplementary file; it may be mixed case.

Notes: | See below

How PDEF Grid Interpolation Works
To illustrate how the data is interpolated from the native grid to the rectilinear grid, let's consider an example. Here are a set of relevant records from a descriptor file:
PDEF 100 100 nps ...
XDEF 181 linear -180 1
YDEF 90 linear 0 1
These three entries describe data on a native 100x100 North Polar stereographic projection and a rectilinear lat/lon grid that is 181 by 90 and has an interval of 1 degree in both lat and lon. Consider one point within the rectilinear grid, the point -90,40. GrADS calls an internal routine to calculate the i and j values in the native grid that correspond to this lat/lon point. Let's say we get i,j values of 31.24 and 67.88. To do the interpolation to the lat/lon point -90,40, GrADS uses the data values from the following four native grid points: 31,67 - 31,68 - 32,67 - 32,68. Bi-linear interpolation is used within this grid box to get down to the position 31.24,67.88. The interpolation is linear within the i,j grid. When a descriptor file is opened that contains a PDEF record, GrADS calculates the i/j values in the native grid that correspond to the lat/lon pair for each gridpoint in the rectilinear grid.

How PDEF Wind Rotation Works
There is a third value calculated for every lat/lon point, and that is the wind rotation value. With some "pre-projected" or native grids, the winds are given relative to the i/j space and not the lat/lon space. To insure correct interpolation, the winds must be rotated to lat/lon space. This is done by determining a rotation amount for each lat/lon point. When u or v wind components are displayed, the values are not just interpolated but also rotated.

To do the wind rotation properly, GrADS requires both the u and v components. Even if you are just displaying u, GrADS has to retrieve (internally) both the u and v component in order to do the rotation calculation. GrADS determines how to match u and v variables by checking the units field of the variable record in the descriptor file. The u variable must have a units value of 33, and the v variable must have a units value of 34. (This is the GRIB convention). If there are more than one u/v pairs, secondary units values are used. For example:

u | 18 | 33,100 | U-Wind Components on Pressure Levels
v | 18 | 34,100 | V-Wind Components on Pressure Levels
u10 | 0 | 33,105 | 10 Meter U Wind
v10 | 0 | 34,105 | 10 Meter V Wind

might be some variable records in the descriptor file. If wind rotation is called for, u and v would be paired, and u10 and v10 would be paired (since the secondary values would be checked, ie, the 105,100 values).

The PDEF BILIN Option
When a descriptor file is opened that contains a PDEF record, we have explained that GrADS internally generates three grids, each one the size of the rectilinear lat/lon grid. The first two grids contain the i and j values (respectively) from the native grid that correspond to each grid point in the rectilinear grid; the third grid contains wind rotation values. But this only works for a small set of well-defined native grids. GrADS will generate these three internal grids automatically for polar stereographic, lamber conformal, and some eta grids. If the native grid for your data is not one of the predefined projections, it is still possible for GrADS to handle the data. All you have to do is supply these three grids to GrADS with a supplementary data file and use the bilin option in your PDEF record.

The supplementary file will contains three lat-lon floating-point grids sized according to the XDEF and YDEF records in the descriptor file.The three grids contain: i values, j values, and wind rotation values. A value of -999 in the i-value grid indicates not to interpolate to that lat-lon point (will end up missing on output) and a value of -999 in the wind-rotation grid indicates not to do wind rotation for that point. If the wind-rotation grid is all -999 values, no rotation is ever done and a flag is set not to even attempt rotation.

The PDEF FILE Option
All of the PDEF examples discussed so far involve the same method for grid interpolation: a grid point value in the rectilinear grid is calculated by finding the four neighboring grid points in the native grid and averaging them, with weights applied bi-linearly according to their proximity to the to rectilinear grid point. The PDEF FILE option generalizes this method so that an arbitrary number of native grid point values and their weights are averaged to generate the interpolated rectilinear grid point values. The index values for the native grid values that are to be used and their weights are specified by the user in a supplementary data file (fname).

The num argument in the PDEF FILE entry specifies the number of native grid points that will be used to calcuate each interpolated rectilinear grid point value. For each num, the supplementary data file will contain two grids -- both will be the size of the rectilinear grid (as defined by XDEF and YDEF). The first grid contains the index or offset values that point to the native grid value that will be used in the interpolation; the second grid contains the weights for those native grid values. The first grid contains integer values, the second grid contains floating-point values. Finally, the supplementary data file must also contain one grid of floating-point wind rotation values. Thus if num equals 1, there will be 3 grids in fname. If num equals 3, there will be 7 grids in fname (3 sets of 2 grids plus the wind rotation grid).

To do the grid interpolation, GrADS will first read the data in the native grid (vector) along with the values in the supplementary grids. To calculate the interpolated value for a particular lat-lon point, GrADS will get num native grid point values, multiply each by their weight, sum over all the weighted native grid points, and divide by the sum of the weights.

Native grid offset values of -999 indicate not to use an input point for that portion of the calculation (thus you can apply less than the "num" number of interpolation points for some of the points). So in theory this can be used for such things as thinned grids, staggered grids, etc.

An Example:

The original data are set up as a vector of land points only, taken from a 1-degree lat/lon grid. There are 15238 land points in the native grid (vector). We use the PDEF FILE option to repopulate a rectilinear lat/lon grid with the land values, leaving the ocean grid points as missing. In this case, ther eis no interpolation done. The PDEF option is used simply to convert a vector of land points into a 2D grid with ocean points missing. The descriptor file looks like this:

DSET ^gswp_vegetation_parameters.nc
DTYPE netcdf
TITLE Monthly Vegetation parameters at 1x1 degree
UNDEF 1.e+20
PDEF 15238 1 file 1 stream binary ^gswp.filepdef
XDEF 360 linear -179.5 1
YDEF 150 linear -59.5 1
ZDEF 1 linear 1 1
TDEF 204 linear 00Z01jan1982 1mo
VARS 1
NDVI=>ndvi 0 t,x Monthly vegetation index
ENDVARS

The supplementary file gtd.filepdef contains three grids -- the first contains the index values that associate each location in the lat/lon grid with it's point in the vector. All of the ocean points will have a missing value of -999. The second grid will contain the weights, which will be 1 for land points, 0 for ocean points. The third grid will contain all missing values since wind rotation is not a issue in this example. Here is a descriptor file for the supplementary data file (a useful strategy for making sure you've got everything written out correctly):

DSET ^gswp.filepdef
TITLE PDEF file for GSWP Vegetation Parameters
UNDEF -999
XDEF 360 linear -179.5 1
YDEF 150 linear -59.5 1
ZDEF 1 linear 1 1
TDEF 1 linear 00z01jul1982 3hr
VARS 3
i 0 -1,40,4 Index Values
w 0 99 Weights
r 0 99 Wind Rotation Values
ENDVARS