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.
Tim Robishaw points out to me, however, that there is a problem in this code. I'll let Tim explain:
My student and I have been working hard on doing some heavy-duty interpolation of astronomical data on an irregular spherical grid. We had a heck of a time trying to figure out how to get GRIDDATA to work for us. In the end the QHULL scheme outlined in the IDL help and sketched out [in this article] turned out to be the way to go. However, as you point out, even doing it this way sometimes causes those not-in-counter-clockwise-order errors from GRIDDATA. After beating our heads against the whiteboard for a couple days, we've figured out that there's a problem in the code snippet that you posted on your page.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)The problem is that when the SPHERE keyword is set, the XYZ value returned from GRID_INPUT is a variable that contains a 3-by-n array of Cartesian coordinates representing points on a sphere.
But when the SPHERE keyword is set for QHULL, it's expecting the two input parameters V0 and V1 to be (from the documentation):
- The V0 argument contains the longitude, in degrees, and V1 contains the latitude, in degrees, of each point.
In the example on your page, you are sending QHULL the X and Y Cartesian components (xyz[0,*] and xyz[1,*]) rather than the latitude and longitude. Instead, you need to define the longitude and latitude in degrees and pass those into QHULL:
lon = !radeg * atan(xyz[1,*],xyz[0,*]) lat = !radeg * asin(xyz[2,*]) qhull, lon, lat, tri, /delaunay, SPHERE=sphThen you need to pass these lon and lat variables into GRIDDATA because that's what it's expecting when /SPHERE and /DEGREES are set and only 3 input arguments are provided:
temp_grid = griddata( lon, lat, newtemp, /degree, /sphere, $ /grid, xout = grid_lon, yout = grid_lat, $ method = 'NaturalNeighbor', $ triangles=tri, $ missing = !values.f_nan)When we do this, it works like a charm. GRID_INPUT has thus far successfully removed all duplicates such that GRIDDATA has thrown no weirdo counter-clockwise-order errors no matter how weird we choose our irregular grid. QHULL takes forever, but it does a very nice job. To be specific, we're actually using the KRIGING and MIN_POINTS=5 keywords in the place of your method = 'NaturalNeighbor'.
Also, maybe IDL has changed things since the code snippet was written, but the SPHERE keyword to QHULL has to be set to a named variable; it can't be set with a switch. If it is, you get the error:
% QHULL: Expression must be named variable in this context:The output sent back from QHULL via this SPHERE keyword seems to be of no use in the GRIDDATA call as illustrated in the IDL help for GRIDDATA, but someone with a better grasp of black magic might be able to parse of what use the SPHERE output might be.
It sounds like this might solve the counterclockwise order problem. But, if it doesn't, 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 = cgScaleVector(Findgen(nx), Min(lon), Max(lon))
slat = cgScaleVector(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)
For additional information on this topic, please see the article Gridding Map Data.
Version of IDL used to prepare this article: IDL 7.0.1.
Copyright © 2008 David W. Fanning
Written 8 March 2008
Last Updated 5 February 2013