;---------------------------------------------------------------------------- ; SCCS/s.nn_interp.pro 1.14 01/06/20 17:57:04 ; ; function nn_interp ; ; Purpose: ; Interpolates data from one grid (can be regular or irregular) to ; another grid (which can also be regular or irregular) using a near- ; neighbor scheme. ; ; Algorithm: ; For each point defined by XPout and YPout, procedure finds all the ; points in in_z that are within Search_Radius away from the X/YPout ; point, and that are not NaN. For all those values, a distance- ; weighted arithmetic mean is done, and that value is entered into ; out_z_interp at that location. If no in_z values are found that ; are within Search_Radius, that location of out_z_interp is given ; the value NaN. The weight for the distance-weighting is the ratio ; of the inverse of the distance over the total of all the inverses. ; For out_z_interp points that are directly at the same location as ; z, the interp. value is just the mean of the z at the same location ; (which is usually the mean of one value). ; ; File I/O: ; None. ; ; Entry and Exit State: ; N/A. ; ; Input Parameters: ; in_z Array of z values at 2-D (x,y) locations on the surf. ; of the Earth. Required. Unchanged. Array can be any ; number of dimensions. Source values which will be used ; to interpolate from. ; in_x Array of x (longitude) locations. In decimal degrees. ; Same dims and type as in_z. Required. Unchanged. ; in_y Array of y (latitude) locations. In decimal degrees. ; Same dims and type as in_z. Required. Unchanged. ; ; Output Parameter: ; out_z_interp Array of z values interpolated onto the (x,y) grid ; specified by XPout and YPout. Type DOUBLE. Created. ; Has same dimensions as X/YPout. ; ; Keywords: ; ALL_DIST_FROM_I Array dimensioned (N_src, N_interp). The vector ; in_all_dist_from_i[*,i] gives the value of local ; variable dist_from_i[*] for the location defined ; by XPout[i] and YPout[i]. Sometimes, you're ; interested in using this function on many datasets ; that have the same source and interpolated ; locations. Thus, precalculating dist_from_i ahead ; of time and just reading it in can save substantial ; time. NB: Since this array is often very large, ; this function does NOT convert its copy into DOUBLE ; for local use, which is the default practice for ; other input variables. Though I think the ; keywords isn't changed, it is not protected by ; duplication. ; SEARCH_RADIUS Radius of near-neighbor searching. In meters. Not ; changed. ; SUPPRESS_EXACT_LOC If keyword is set, procedure does not include ; the values of z that are identically at the ; location specified by XPout and YPout when ; calculating the interpolation. This is useful in ; cases when the point is question has NaN at that ; location in z, and the procedure is being run to ; try and fill in that value. If no points inside ; the search_radius exists besides the exact loca- ; tion, the procedure returns NaN. ; VERBOSE If keyword is set, prints diagnostic messages. Set ; this especially if the array being operated on is huge, ; since this helps you keep track of what's happening. ; XPOUT Array of x locations that we're interpolating in_z to, ; i.e. the x locations of out_z_interp. Required ; keyword. Longitude in decimal degrees. Not changed. ; Can be any number of dimensions. ; YPOUT Array of y locations that we're interpolating in_z to, ; i.e. the y locations of out_z_interp. Required ; keyword. Latitude in decimal degrees. Same dimensions ; and type as XPout. Not changed. ; ; Revision History: ; - 8 Feb 2001: Orig. ver. by Johnny Lin, CIRES/Univ. of Colo. Email: ; air_jlin@yahoo.com. Passed moderately reasonable tests. ; - 9 Feb 2001: Changed so input can any dimension of array. Fixed ; error that gave incorrect value for zero_dist_pts cases. Passed ; moderately reasonable tests. ; - 12 Feb 2001: Added Suppress_Exact_Loc keyword. Passed moderately ; reasonable tests. ; - 15 Feb 2001: Slight change to MEAN enable IDL 5.2.1 to run. ; - 20 Jun 2001: Add All_Dist_From_i keyword. Passed passably reasonable ; tests. ; ; Notes: ; - Written for IDL 5.2.1. ; - Copyright (c) 2001 Johnny Lin. For licensing and contact information ; see http://www.johnny-lin.com/lib.html. ; - Calculations done in DOUBLE. ; - No procedures called with _EXTRA keyword invoked. ; - User-written procedures called: WTD_MEAN, DISTANCE ; - No common blocks are used in this program. ;---------------------------------------------------------------------------- FUNCTION NN_INTERP, in_z, in_x, in_y $ , ALL_DIST_FROM_I = in_all_dist_from_i $ , SEARCH_RADIUS = in_search_radius $ , SUPPRESS_EXACT_LOC = suppress_exact_loc $ , VERBOSE = verbose $ , XPOUT = in_xpout $ , YPOUT = in_ypout $ , _EXTRA = extra ; -------------------- Error Check and Parameter Setting -------------------- ON_ERROR, 0 if (N_PARAMS() ne 3) then MESSAGE, 'error-bad param list' if (N_ELEMENTS(in_search_radius) eq 0) then MESSAGE, 'error-no SEARCH_RADIUS' if (N_ELEMENTS(in_xpout) eq 0) then MESSAGE, 'error-no XPOUT' if (N_ELEMENTS(in_ypout) eq 0) then MESSAGE, 'error-no YPOUT' z = DOUBLE(in_z) ;- protect input parameters x = DOUBLE(in_x) y = DOUBLE(in_y) search_radius = DOUBLE(in_search_radius) xpout = DOUBLE(in_xpout) ypout = DOUBLE(in_ypout) sd_z = SIZE(z) sd_x = SIZE(x) sd_y = SIZE(y) if (WHERE(sd_z ne sd_x) ne -1) then MESSAGE, 'error-input different' if (WHERE(sd_z ne sd_y) ne -1) then MESSAGE, 'error-input different' sd_xp = SIZE(xpout) sd_yp = SIZE(ypout) if (WHERE(sd_xp ne sd_yp) ne -1) then MESSAGE, 'error-input different' z_interp_dims = SIZE(xpout, /Dimensions) ;- dims for future z_interp array N_src = N_ELEMENTS(z) ;- source array number elements N_interp = N_ELEMENTS(xpout) ;- interp. array number elements suppress_exact_loc_state = KEYWORD_SET(suppress_exact_loc) ;- suppress_exact_ ; loc keyword ; state if (WHERE(FINITE(x) eq 0) ne -1) then MESSAGE, 'error-loc have NaN values' if (WHERE(FINITE(y) eq 0) ne -1) then MESSAGE, 'error-loc have NaN values' if (WHERE(FINITE(xpout) eq 0) ne -1) then MESSAGE, 'error-loc have NaN values' if (WHERE(FINITE(ypout) eq 0) ne -1) then MESSAGE, 'error-loc have NaN values' ; ----------- Calculate Weighted Average Based on Near Neighbors ------------ z_interp = DBLARR(N_interp) ;- initialize eventual output array for i=0L,N_interp-1L do begin ;- cycle each pt. on interpol. (output) grid if (KEYWORD_SET(verbose) eq 1) then begin ;- print out PRINT, ' Processing output grid (i) = (' $ ; diagnostic if , STRING(i, Format='(g0.0)'), ')' ; VERBOSE is set endif if (N_ELEMENTS(in_all_dist_from_i) gt 0) then begin if (N_ELEMENTS(in_all_dist_from_i) ne N_src*N_interp) then $ MESSAGE, 'error--bad in_all_dist_from_i' dist_from_i = REFORM(in_all_dist_from_i[*,i]) endif else begin lon1 = REPLICATE(xpout[i], N_src) ;+ the location of the point on lat1 = REPLICATE(ypout[i], N_src) ; the interpol. grid lon2 = REFORM(x, N_src) ;+ location of the src. data grid lat2 = REFORM(y, N_src) dist_from_i = DISTANCE( lon1, lat1, lon2, lat2 ) endelse near_neighbors = WHERE( dist_from_i le search_radius, count ) if (count gt 0) then begin ;+ wtd. mean of near neighbors values = z[near_neighbors] weights = dist_from_i[near_neighbors] case suppress_exact_loc_state of 0: begin zero_dist_pts = WHERE(weights eq 0.0, zp_count) case zp_count of 0: begin weights = 1.d0 / weights z_interp[i] = WTD_MEAN(values, weights) end 1: z_interp[i] = values[zero_dist_pts] else: z_interp[i] = MEAN(values[zero_dist_pts]) endcase end 1: begin not_zero_dist_pts = WHERE(weights ne 0.0, np_count) if (np_count gt 0) then begin values = values[not_zero_dist_pts] weights = weights[not_zero_dist_pts] weights = 1.d0 / weights z_interp[i] = WTD_MEAN(values, weights) endif else $ z_interp[i] = !VALUES.F_NAN end else: MESSAGE, 'error-bad setting' endcase endif else begin PRINT, 'NN_INTERP: Warning - NaN set at (i) = (' $ , STRING(i, Format='(g0.0)'), ')' z_interp[i] = !VALUES.F_NAN endelse endfor ; -------------------------------- Clean-Up --------------------------------- out_z_interp = REFORM(z_interp, z_interp_dims) RETURN, out_z_interp END ; ===== end of procedure ===== ; ========== end file ==========