Gridding Irregular Data
QUESTION: I am interpolating satellite data (latitude, longitude, and temperature) from an irregular grid to a regular lat/lon grid using GRIDDATA, like this:
pro get_gridded_temp_image, ....
.
.
.
grid_lat=findgen(31)*2.0+30.0 ;30 to 90 deg @2 deg
grid_lon=findgen(181)*2.0 ;0 to 360 deg @2 deg
qhull, lon, lat, tri, /delaunay, /sphere
temp_grid = griddata( lon, lat, temperature, /sphere, /degree, $
/grid, xout = grid_lon, yout = grid_lat, $
method = 'NaturalNeighbor', $
triangles=tri, $
missing = !values.f_nan)
Which is patterned after Example 5 in the IDL help for GRIDDATA.
Usually this procedure works fine, over thousands of different cases. However, today I tried it and got this error message:
GRIDDATA: Triangle 652 not in counterclockwise order. Execution halted at: GET_GRIDDED_TEMP_IMAGE 184
Any ideas on what could be wrong?
ANSWER: This is, apparently, an extremely common problem with GRIDDATA. It appears to be caused by overlapping, or nearly overlapping, data points in the triangulation part of the routine. There are, generally, two things you can try to eliminate this problem. First, you can try to eliminate these overlapping or duplicate points by preproccessing your data with GRID_INPUT, whose purpose is to do just this. The following code snippet might fix your problem.
pro get_gridded_temp_image, ....
.
.
.
grid_lat=findgen(31)*2.0+30.0 ;30 to 90 deg @2 deg
grid_lon=findgen(181)*2.0 ;0 to 360 deg @2 deg
grid_input, lon, lat, temperature, xyz, newtemp, /degree, /sphere, epsilon = 0.5
qhull, xyz[0,*], xyz[1,*], tri, /delaunay, /sphere
temp_grid = griddata( xyz[0,*], xyz[1,*], newtemp, /degree, /sphere, $
/grid, xout = grid_lon, yout = grid_lat, $
method = 'NaturalNeighbor', $
triangles=tri, $
missing = !values.f_nan)
Indeed, the person who reported the problem indicated that there were duplicate points in his data he had not noticed before, and that this solved the problem for him.
Sometimes, using GRID_INPUT isn't enough to solve this problem, and you can try a second solution, which is to choose a different method, other than QHULL, to create the triangulation that goes into GRIDDATA. Ben Tupper reports such a problem gridding a dataset that is regular in longitude and irregular in latitude. You can download the program, named GridData_Example, Ben used to illustrate the problem. The program, as Ben original wrote it, threw these errors.
IDL> gridded = GridData_Example(/Sphere) % GRIDDATA: Triangle 1 not in counterclockwise order. % GRIDDATA: Triangle 2 not in counterclockwise order. % GRIDDATA: Triangle 4 not in counterclockwise order. % GRIDDATA: Triangle 5 not in counterclockwise order. % GRIDDATA: Triangle 7 not in counterclockwise order. % GRIDDATA: Triangle 11 not in counterclockwise order. % GRIDDATA: Triangle 19 not in counterclockwise order. % GRIDDATA: Triangle 20 not in counterclockwise order. % GRIDDATA: Triangle 61 not in counterclockwise order. % GRIDDATA: Triangle 65 not in counterclockwise order. % Execution halted at: $MAIN$ 137 C:\griddata_example.pro
In this case, the solution was to use the TRIANGULATE method of producing the triangles required for GRIDDATA.
IDL> gridded = GridData_Example(/Sphere, TRIMETHOD='triangulate')
Although there is not suppose to be a difference in the Delaunay triangulations produced by the two routines, in practice, there is occasionally. For example, in Ben's code we find the two methods produce a different number of triangles.
QHULL Triangles: Array[3, 744557] TRIANGULATE Triangles: Array[3, 742122]
Sometimes, it is just impossible to get spherical gridding to work (i.e, using the SPHERE keyword). No one seems to know why. In those cases, you have to convince yourself that linear scaling is close enough for your purposes. A combination of these methods usually solves the problem.
An Alternative Gridding Method
Ken Bowman points out that there is another common way to grid latitude and longitude data of the sort Ben uses in his program.
If your grid is rectangular and separable (in the sense that all the longitudes in each "column" of data are the same and all of the latitudes in each "row" of data are same), even if the coordinates are not regularly spaced, then it is actually quite easy to interpolate to any set of points (regular or irregular) using INTERPOLATE. This should be much faster than triangulating.
Note that the requirement Ken imposes for the grid is not the case in the first example above, although it is the case with Ben's data. Ken suggests the following code outline for this purpose. Assuming original variables lon, lat, and temperature, the first step is to create vectors for constructing a regular grid.
nx = 360 ny = 181 slon = FINDGEN(nx) slat = -90.0 + FINDGEN(ny)
Next, we compute the fractional dimension “interpolation coordinates” from the original data.
x = Interpol(Findgen(N_Elements(lon)), lon, slon) y = Interpol(Findgen(N_Elements(lat)), lat, slat)
The next step is to expand the new x and y vectors into 2D arrays.
xx = Rebin(x, nx, ny, /SAMPLE) yy = Rebin(Reform(y, 1, ny), nx, ny, /SAMPLE)
And, finally, we interpolate the new temperature data at the new regular grid locations.
newTemperature = INTERPOLATE(temperature, xx, yy)
I used Ken's method to interpolate output from a CCCMA climate model, which is in a 96 by 48 array, where the latitudes are spaced according to a Gaussian grid. For my project, I needed a regularly spaced grid, so I wanted to interpolate to a regular grid. My approach looks like this.
DATADIR = 'G:\data\cccma'
orig = Filepath(ROOT_DIR=DATADIR, $
'pcmdi.ipcc4.cccma_cgcm3_1.1_t47_1850_2000.nc')
; Read the latitude and longitude data required for regridding.
cdf = Obj_New('NCDF_DATA', orig)
lat = cdf -> ReadVariable('lat')
lon = cdf -> ReadVariable('lon')
ice = cdf -> ReadVariable('sit')
ice = ice[*,*,1000] ; For a particular year in the analysis.
; Interpolate to a regular grid.
nx = N_Elements(lon)
ny = N_Elements(lat)
slon = Scale_Vector(Findgen(nx), Min(lon), Max(lon))
slat = Scale_Vector(Findgen(ny), Min(lat), Max(lat))
x = Interpol(Findgen(N_Elements(lon)), lon, slon)
y = Interpol(Findgen(N_Elements(lat)), lat, slat)
xx = Rebin(x, nx, ny, /SAMPLE)
yy = Rebin(Reform(y, 1, ny), nx, ny, /SAMPLE)
newIce = Interpolate(ice, xx, yy)
Version of IDL used to prepare this article: IDL 7.0.1.
Copyright © 2008 David W. Fanning
Last Updated 8 March 2008