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?

ANSWER: Ed Hyer points out in an IDL Newsgroup article that one way to do this is to use JD Smith's amazing program HIST_ND.

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=[1] & y=[2]
   IDL> print, transpose([[x],[y]])
        1   2
   IDL> help, transpose([[x],[y]])
        <Expression>    INT       = Array[2]

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.

Google
 
Web Coyote's Guide to IDL Programming