Hiệu quả tìm thấy các chỉ số của các điểm gần nhất trên lưới 2D không phải hình chữ nhật

Question

Tôi có lưới điện/vĩ độ không đều (không hình chữ nhật) và một số điểm trong tọa độ lon/vĩ độ, tương ứng với các điểm trên lưới (mặc dù chúng có thể hơi tắt vì lý do số). Bây giờ tôi cần các chỉ số của các điểm lon/lat tương ứng.

Tôi đã viết một chức năng thực hiện điều này, nhưng nó thực sự chậm.

def find_indices(lon,lat,x,y): 
    lonlat = np.dstack([lon,lat]) 
    delta = np.abs(lonlat-[x,y]) 
    ij_1d = np.linalg.norm(delta,axis=2).argmin() 
    i,j = np.unravel_index(ij_1d,lon.shape) 
    return i,j 

ind = [find_indices(lon,lat,p*) for p in points] 

Tôi chắc rằng có một giải pháp tốt hơn (và nhanh hơn) trong gọn gàng/scipy. Tôi đã googled khá nhiều, nhưng câu trả lời đã cho đến nay eluded tôi.

Bất kỳ đề xuất nào về cách tìm hiệu quả các chỉ số của các điểm tương ứng (gần nhất)?

PS: Câu hỏi này nổi lên từ một câu hỏi khác (click).

Edit: Giải pháp

Dựa trên câu trả lời @Cong Ma, tôi đã tìm thấy các giải pháp sau đây:

def find_indices(points,lon,lat,tree=None): 
    if tree is None: 
     lon,lat = lon.T,lat.T 
     lonlat = np.column_stack((lon.ravel(),lat.ravel())) 
     tree = sp.spatial.cKDTree(lonlat) 
    dist,idx = tree.query(points,k=1) 
    ind = np.column_stack(np.unravel_index(idx,lon.shape)) 
    return [(i,j) for i,j in ind] 

Để đưa giải pháp này và cũng là một từ câu trả lời Divakar của thành quan điểm, đây là một số thời gian của hàm mà tôi đang sử dụng find_indices (và nơi nó là nút cổ chai về tốc độ) (xem liên kết ở trên):

spatial_contour_frequency/pil0    : 331.9553 
spatial_contour_frequency/pil1    : 104.5771 
spatial_contour_frequency/pil2    :  2.3629 
spatial_contour_frequency/pil3    :  0.3287 

pil0 là cách tiếp cận ban đầu của tôi, pil1 Divakar's và pil2/pil3 giải pháp cuối cùng ở trên, nơi cây được tạo trực tiếp trong pil2 (tức là cho mỗi lần lặp của vòng lặp trong đó find_indices được gọi) và chỉ một lần trong pil3 (xem các chủ đề khác để biết chi tiết). Mặc dù sự tinh chỉnh của Divakar về cách tiếp cận ban đầu của tôi mang lại cho tôi tốc độ gấp 3 lần, cKDTree mang đến một cấp độ hoàn toàn mới với tốc độ 50 lần khác! Và di chuyển việc tạo ra cây ra khỏi chức năng làm cho mọi việc nhanh hơn.

Answer 1

Nếu các điểm được bản địa hóa đầy đủ, bạn có thể thử trực tiếp scipy.spatial 's cKDTree triển khai, như được thảo luận một mình in another post. Bài đăng đó là về nội suy nhưng bạn có thể bỏ qua điều đó và chỉ sử dụng phần truy vấn.

tl; dr phiên bản:

Đọc tài liệu của scipy.sptial.cKDTree. Tạo cây bằng cách chuyển đối tượng (n, m)-hình numpy ndarray tới trình khởi tạo và cây sẽ được tạo từ các tọa độ hai chiều nm hai chiều.

tree = scipy.spatial.cKDTree(array_of_coordinates) 

Sau đó, sử dụng tree.query() để lấy hàng xóm gần nhất -thứ k (có thể với xấp xỉ và song song, xem tài liệu), hoặc sử dụng tree.query_ball_point() để tìm tất cả các nước láng giềng trong khoan dung khoảng cách nhất định.

Nếu các điểm không được bản địa hoá tốt, và đường cong hình cầu/đá nhỏ không liên quan đến đá, bạn có thể thử phá đa tạp thành nhiều phần, từng phần nhỏ đủ để được xem là cục bộ.

Answer 2

Dưới đây là một cách tiếp cận vectorized generic sử dụng scipy.spatial.distance.cdist -

import scipy 

# Stack lon and lat arrays as columns to form a Nx2 array, where is N is grid**2 
lonlat = np.column_stack((lon.ravel(),lat.ravel())) 

# Get the distances and get the argmin across the entire N length 
idx = scipy.spatial.distance.cdist(lonlat,points).argmin(0) 

# Get the indices corresponding to grid's shape as the final output 
ind = np.column_stack((np.unravel_index(idx,lon.shape))).tolist() 

mẫu chạy -

In [161]: lon 
Out[161]: 
array([[-11. , -7.82 , -4.52 , -1.18 , 2.19 ], 
     [-12. , -8.65 , -5.21 , -1.71 , 1.81 ], 
     [-13. , -9.53 , -5.94 , -2.29 , 1.41 ], 
     [-14.1 , -0.04 , -6.74 , -2.91 , 0.976]]) 

In [162]: lat 
Out[162]: 
array([[-11.2 , -7.82 , -4.51 , -1.18 , 2.19 ], 
     [-12. , -8.63 , -5.27 , -1.71 , 1.81 ], 
     [-13.2 , -9.52 , -5.96 , -2.29 , 1.41 ], 
     [-14.3 , -0.06 , -6.75 , -2.91 , 0.973]]) 

In [163]: lonlat = np.column_stack((lon.ravel(),lat.ravel())) 

In [164]: idx = scipy.spatial.distance.cdist(lonlat,points).argmin(0) 

In [165]: np.column_stack((np.unravel_index(idx,lon.shape))).tolist() 
Out[165]: [[0, 4], [0, 4], [0, 4], [0, 4], [0, 4], [0, 4], [3, 3]] 

kiểm tra Runtime -

chức năng Define:

def find_indices(lon,lat,x,y): 
    lonlat = np.dstack([lon,lat]) 
    delta = np.abs(lonlat-[x,y]) 
    ij_1d = np.linalg.norm(delta,axis=2).argmin() 
    i,j = np.unravel_index(ij_1d,lon.shape) 
    return i,j 

def loopy_app(lon,lat,pts): 
    return [find_indices(lon,lat,pts[i,0],pts[i,1]) for i in range(pts.shape[0])] 

def vectorized_app(lon,lat,points): 
    lonlat = np.column_stack((lon.ravel(),lat.ravel())) 
    idx = scipy.spatial.distance.cdist(lonlat,points).argmin(0) 
    return np.column_stack((np.unravel_index(idx,lon.shape))).tolist() 

Thời gian:

In [179]: lon = np.random.rand(100,100) 

In [180]: lat = np.random.rand(100,100) 

In [181]: points = np.random.rand(50,2) 

In [182]: %timeit loopy_app(lon,lat,points) 
10 loops, best of 3: 47 ms per loop 

In [183]: %timeit vectorized_app(lon,lat,points) 
10 loops, best of 3: 16.6 ms per loop 

Đối ép ra hiệu suất hơn, np.concatenate có thể được sử dụng thay cho np.column_stack.