python - Numpy griddata interpolation up to certain radius -


i'm using griddata() interpolate (irregular) 2-dimensional depth-measurements; x,y,depth. method great job - interpolates on entire grid can find opposing points. don't want behaviour. i'd have interpolation around existing measurements, extent of radius.

is possible tell numpy/scipy: don't interpolate if you're far existing measurement? resulting in nodata-value? ideal = griddata(.., .., .., radius=5.0)

edit example: in image below; black dots measurements. shades of blue interpolated cells numpy. area marked in green in fact part of picture considered nodata numpy (because there's no points in between). now, red areas, interpolated, want rid of them. ideas? enter image description here

ok cool. don't think there built-in option griddata() want, need write yourself.

this comes down calculating distances between n input data points , m interpolation points. simple enough if have lot of points can slow @ ~o(m*n). here's example calculates distances alln data points, each interpolation point. if number of data points withing radius @ least neighbors, keeps value. otherwise writes value of nodata.

neighbors 4 because griddata() use biilinear interpolation needs points bounding interpolants in each dimension (2*2 = 4).

#invec - input points nx2 numpy array #mvec - interpolation points mx2 numpy array  #just random points example n=100 invec = 10*np.random.random([n,2])  m=50 mvec = 10*np.random.random([m,2])  # --- here put griddata() call, returning interpolated_values interpolated_values = np.zeros(m) nodata=np.nan  radius = 5.0 neighbors = 4  m in range(m):     data_in_radius = np.sqrt(np.sum( (invec - mvec[m])**2, axis=1)) <= radius      if np.sum(data_in_radius) < neighbors :         interpolated_values[m] = nodata 

edit: ok re-read , noticed input 2d. example modified.

just additional comment, accelerated if first build coarse mapping each point mvec[m] subset of relevant data points. costliest step in loop change from

np.sqrt(np.sum( (invec - mvec[m])**2, axis=1))

to like

np.sqrt(np.sum( (invec[subset[m]] - mvec[m])**2, axis=1))

there plenty of ways this, example using quadtree, hashing function, or 2d index. whether gives performance advantage depends on application, how data structured, etc.


Comments

Popular posts from this blog

angular - Ionic slides - dynamically add slides before and after -

minify - Minimizing css files -

Add a dynamic header in angular 2 http provider -