5.0

# Speeding up linear gridding of irregular points with multiple values (GRIDDATA)

Atle Borsholm

Gridding or interpolating large amounts of data is a common task for some IDL and ENVI users. Here, I am showing a trick that can speed up bi-linear interpolation using a triangulated collection of irregularly gridded points in 2-D. The assumption here is that there are multiple values for each distinct point (x,y), and instead of using GRIDDATA repeatedly for several hundred values at the same locations, the code is pre-computing weights for the triangle corners. This saves computations in the final step and thus achieves a nice speed improvement.

The speed improvement on my computer was going from 31.4 seconds to 8.4 seconds. Here is the output produced by the example code:

IDL> grid_speed

% Compiled module: GRID_SPEED.

% Time elapsed: 31.422000 seconds.

% Time elapsed: 8.4380002 seconds.

Mean, Variance, Skewness, Kurtosis

0.426667    0.0376666      1.26416     0.492909

0.426667    0.0376666      1.26416     0.492909

min, mean, max difference

-1.19209e-007 1.24474e-011 1.19209e-007

Here is the example code:

pro grid_speed

compile_opt idl2,logical_predicate

;Set up random data points

;Let's say 200,000 spatial (X,Y) points with 400 measurements each

npts = 200000

nbands = 400

im = randomu(seed, npts, nbands)

x = randomu(seed, npts)

y = randomu(seed, npts)

;Set up an output gridded space for desired locations

nx = 768

ny = 768

start = [0,0]

delta = 1d / [nx, ny]

dim = [nx, ny]

gridIm1 = fltarr(nx, ny, nbands)

gridIm2 = fltarr(nx, ny, nbands)

tic

triangulate, x, y, tr

for i=0, nbands-1 do begin

gridIm1[0,0,i] = griddata(x, y, im[*,i], triangles=tr, /linear, \$

start=start, delta=delta, dimension=dim)

endfor

toc

tic

triangulate, x, y, tr

;compute triangle numbers for each input point

;multiply by 3 so that triangles are numbered

;by the starting index 0, 3, 6, 9, ...

index = lindgen(n_elements(tr))/3*3

xt = x[tr[*]]

yt = y[tr[*]]

linTr = lindgen(size(tr, /dimensions))

tr_num = round( \$

griddata(xt, yt, float(index),triangles=linTr, /linear, \$

start=start, delta=delta, dimension=dim))

;Compute weights for each of the 3 points in the enclosing triangle

wts = ptrarr(3)

for i=0, 2 do begin

w = griddata(xt, yt, lindgen(n_elements(xt)) mod 3 eq i, \$

triangles=linTr, /linear, \$

start=start, delta=delta, dimension=dim)

wts[i] = ptr_new(w, /no_copy)

endfor

;Compute interpolation for all bands using weights

for i=0, nbands-1 do begin

gridIm2[0,0,i] = im[tr[tr_num] + i*n_elements(x)] * (*wts)

gridIm2[*,*,i] += im[tr[tr_num+1] + i*n_elements(x)] * (*wts)

gridIm2[*,*,i] += im[tr[tr_num+2] + i*n_elements(x)] * (*wts)

endfor

toc

;Verify that the results are the same for both

;methods.

print, 'Mean, Variance,Skewness, Kurtosis'

print, moment(gridIm1)

print, moment(gridIm2)

print

diff = gridIm2 - gridIm1

print, 'min, mean, maxdifference'

print, min(diff, max=maxDiff), mean(diff), maxDiff

end

MOST RECENT

### Performing Flood Mapping After Cyclone Yaas Using ENVI SARscape

9/24/2021

In May 2021, Cyclone Yaas hit the east coast of India with wind gusts up to 140 kilometers per hour (87... more »

### Cost-effectively Converting Legacy IDL Desktop Applications to Enterprise Applications

9/17/2021

Porting your IDL® application to an enterprise or cloud-based system can seem... more »