emerald-geomodelling / EmeraldTriangles

Transformations for triangle meshes build on top of scipy.spatial.Delaunay
MIT License
1 stars 0 forks source link

compute distances using cKDTree #12

Open Duke-of-Lizard opened 2 years ago

Duke-of-Lizard commented 2 years ago

right now, we're doing quite an inefficient computation of minimum distance to nearest data point in distance.py. I'm familiar with scipy.sparse and scipy.spatial to pull it off. Here's what I tried, but the method .min on a csr.array was useless because I can't exlude the blank entries of the matrix that are assumed to be 0.

import scipy.spatial.distance
from scipy.sparse import csr_array
from scipy.spatial import cKDTree
import numpy as np

def distance_to_data(col, x_col="X", y_col="Y",max_search_dist = None, **tri):
    """Calculate spatial distance (cartesian distance in the current
    projection) to vertices with non-NaN values in the column col and
    store the distance in the column col_dist.
    """
    filt = ~tri["vertices"][col].isna()
    XA = tri["vertices"].loc[~filt][[x_col,y_col]].values
    XB = tri["vertices"].loc[filt][[x_col,y_col]].values
    if len(XB):
        if max_search_dist is None:
            tri["vertices"].loc[filt,'%s_dist' % col] = 0.0
            tri["vertices"].loc[~filt,'%s_dist' % col] = np.min(scipy.spatial.distance.cdist(XA, XB),axis=1)

        else:
            A = cKDTree(XA[~filt])
            B = cKDTree(XB)
            dist_mat = A.sparse_distance_matrix(B, max_search_dist)
            dist_nearest = csr_array(dist_mat).min(axis=1).toarray()
            tri["vertices"]['%s_dist' % col] = dist_nearest
    return tri