Finding nearest xy-point in numpy array and second nearest with condition

Question

My problem is like the problem in the thread Finding index of nearest point in numpy arrays of x and y coordinates, but it's extended:

For better visualization here's an image (manipulated image, original from: by 112BKS - Eigenes WerkOriginal graph/Data from [.. ? ..], CC BY-SA 3.0, link to page):

On the one hand there is a array datafield. It consists of a numpy array with elements [value x y]. That are the thin blue lines with the numbers (they are the value). On the other hand there is the array orangeline in a numpy array with elements [x y].

What I want to do is to calculate the value of any elements in orangeline.
I visualized one concrete element of orangeline with the green circle. The value for it can I interpolate with the two elements from datafield, visualized with the triangles. As result I get for the green circle a value between 225 and 230.

First step:
Find for every element in orangeline the closest element in datafield.
(In the example that is the pink triangle.)

Second step:
Find for every element in 'orangeline' the closest element in datafield but with another value than the one from the first step.
(In the example that is the brown triangle.)

Third step: Interpolate the value for every element in orangeline from those the two founded values and the distances to those elements.

First step can be solved with

mytree = scipy.spatial.cKDTree(datafield[:, 1:3])
dist1, indexes1 = mytree.query(orangeline)

But now I don't know how to filter the datafield for the second step. Is there a solution?

Answer 1

With help from @unutbu's comment I found this solution which works quite good also in those cases where the orangeline goes not through the field.

Here are the functions for the grid:

import matplotlib.mlab as mlab
import numpy as np
import scipy

def define_grid(rawdata):
    xmin, xmax = np.amin(rawdata[:, 1]), np.amax(rawdata[:,1])
    ymin, ymax = np.amin(rawdata[:, 2]), np.amax(rawdata[:,2])

    x, y, z = rawdata[:, 1], rawdata[:, 2], rawdata[:, 0]

    # Size of regular grid
    ny, nx = (ymax - ymin), (xmax - xmin)

    # Generate a regular grid to interpolate the data.
    xi = np.linspace(xmin, xmax, nx)
    yi = np.linspace(ymin, ymax, ny)
    xi, yi = np.meshgrid(xi, yi)

    # Interpolate using delaunay triangularization
    zi = mlab.griddata(x,y,z,xi,yi)
    return xi, yi, zi

def grid_as_array(xi,yi,zi):
    xi_flat, yi_flat, zi_flat = np.ravel(xi), np.ravel(yi), np.ravel(zi)

    # reduce arrays for faster calculation, take only every second element
    xi_red, yi_red, zi_red = xi_flat[1::2], yi_flat[1::2], zi_flat[1::2]

    # stack to array with elements [x y z], but there are z values that are 'nan'
    xyz_with_nan = np.hstack((xi_red[:, np.newaxis], yi_red[:, np.newaxis],
                              zi_red[:, np.newaxis]))

    # sort out those elements with 'nan'
    xyz = xyz_with_nan[~np.isnan(xyz_with_nan).any(axis=1)]
    return xyz

Another function to find the closest point from the grid for the values from orangeline:

def closest_node(points, datafield):
    mytree = scipy.spatial.cKDTree(datafield)
    dist, indexes = mytree.query(points)
    return indexes

And now the code:

# use function to create from the raw data an interpolated datafield
xi, yi, zi = define_grid(datafield)

# rearrange those values to bring them in the form of an array with [x y z]
xyz = grid_as_array(xi, yi, zi)    

# search closest values from grid for the points of the orangeline
# orangeline_xy is the array with elements [x y]
indexes = self.closest_node(orangeline_xy, xyz[:,0:2])

# take z values from the grid which we found before
orangeline_z = xyz[indexes, 2]

# add those z values to the points of the orangeline
orangeline_xyz = np.hstack((orangeline_xy,orangeline_z[:, np.newaxis]))