Matching Point Source Data with Map Grids
QUESTION: I have a point source data set and I would like to match it with a particular map projection grid, or visa versa. The FOR loop method I am using is exceedingly slow. Is there a faster way to do this in IDL?
I use HIST_ND a great deal. One of the important things I do with it is match gridded maps to point data sets and vice versa. Like this:IDL> histogram = HIST_ND(TRANSPOSE([[point_x],[point_y]]), $ [dx_grid,dy_grid], MIN=[minx_grid,miny_grid], $ MAX=[maxx_grid,maxy_grid],REVERSE_INDICES=ri) IDL> matched_pts = point_x * 0.0 ; Initialize array to hold sampled output IDL> FOR i=0L,N_ELEMENTS(histogram)-1 DO IF (histogram[i] GT 0) THEN $ matched_points[ri[ri[i]:(ri[i+1]-1)]] = gridded_data[i]
or:IDL> grid_totals = gridded_data * 0.0 ; Initialize array to hold gridded output IDL> FOR i=0L,N_ELEMENTS(histogram)-1 DO IF(histogram[i] GT 0) THEN $ grid_totals[i] = TOTAL(point_data[ri[ri[i]:(ri[i+1]-1)]])
The only thing to watch out for is that HIST_ND will be unhappy when the inputs [point_x] and [point_y] contain only one point. JD Smith offers this explanation and solution.
HIST_ND wants an NxP input array. This is to warn against accidentally sending it: [point_x,point_y] or some such. The problem with your method of forming an NxP array with a single point is easily illustrated:IDL> x= & y= IDL> print, transpose([[x],[y]]) 1 2 IDL> help, transpose([[x],[y]]) <Expression> INT = Array
IDL has happily trimmed your final shallow dimension for you, and HIST_ND can no longer tell that you've intended a 2x1 array. To rememdy this, you'll need to add that trailing dimension back on with REFORM.IDL> help, reform(transpose([[x],[y]]), 2, 1) <Expression> INT = Array[2, 1]
I suppose I could do this for you inside of HIST_ND, but then it would fail to trap accidental usage of a pure vector input.
Copyright © 2007 David W. Fanning
Last Updated 26 April 2007