Pythonでイメージに円形フィルタを適用する/ numpy配列の各要素に関数を適用する

Question

私はPythonコードを使用していますが、それはかなり遅く、これをより効率的に行う方法が必要だと思います。 アイデアは画像にフィルタを適用することです。フィルタは、指定された半径内にある点の平均です。入力は、x、yを表すmx2配列であり、m個の観測点の座標を表すmx1配列zです。作品

プログラムは

import numpy as np 
def haversine(point, xy_list): 
    earth_radius = 6378137.0 
    dlon = np.radians(xy_list[:,0]) - np.radians(point[0]) 
    dlat = np.radians(xy_list[:,1]) - np.radians(point[1]) 
    a = np.square(np.sin(dlat/2.0)) + np.cos(np.radians(point[0])) * np.cos(np.radians(xy_list[:,0])) * np.square(np.sin(dlon/2.0)) 
    return 2 * earth_radius * np.arcsin(np.sqrt(a)) 

def circular_filter(xy, z, radius): 
    filtered = np.zeros(xy.shape[0]) 
    for q in range(xy.shape[0]): 
     dist = haversine(xy[q,:],xy) 
     masked_z = np.ma.masked_where(dist>radius, z) 
     filtered[q] = masked_z.mean() 
    return filtered 

x = np.random.uniform(low=-90, high=0, size=(1000,1)) # x represents longitude 
y = np.random.uniform(low=0, high=90, size=(1000,1)) # y represents latitude 
xy = np.hstack((x,y)) 
z = np.random.rand(1000,) 
fitered_z = circular_filter(xy, z, radius=100.) 

次のような問題は、私は、データセットごとに600万ポイントを持っているということです、そしてコードが恐ろしく遅いです。これをより効率的に行う方法が必要です。高速であるscipy.spatial.distance.cdist（）を使用することを考えましたが、データをUTMに再投影する必要がありました。再投影を避けたいと思います。助言がありますか？

おかげで、 Reniel

Answer 1

私はようやく私のコードが実行される永遠にかかった理由を発見した検索を読んだ多くの後。なぜなら、私はフィルターカーネルの概念を理解し、適用する必要があったからです。基本的に私は私の問題と、このポストの間の接続があった実現： Local Maxima with circular window

欠点：ユーザーは、適切なEPSGコードを提供する必要があるが、私はこの後の回避策を見つけることができると思います。

アップサイド：非常に高速で効率的です。

私にとっては、円形のカーネルを作成してscipyからgeneric_filterを適用できるように、lat longをUTMに変換していました。

import numpy as np 
from pyproj import Proj, transform 
from scipy.ndimage.filters import generic_filter 
def circular_filter(tile, radius): 
    x, y = np.meshgrid(tile['lon'], tile['lat']) 
    x = x.reshape(x.size) 
    y = np.flipud(y.reshape(y.size)) 
    z = tile['values'].reshape(tile['values'].size) 
    wgs84 = Proj(init='epsg:4326') 
    utm18N = Proj(init='epsg:26918') 
    x,y = transform(wgs84,utm18N,x,y) 

    dem_res = np.abs(x[1]-x[0]) # calculates the raster resolution, (original data is geoTiff read using gdal). 

    radius = int(np.ceil(radius/dem_res)) # user gives meters, but we figure out number of cells 
    print(radius) 
    kernel = np.zeros((2*radius+1, 2*radius+1)) 
    y,x = np.ogrid[-radius:radius+1, -radius:radius+1] 
    mask = x**2 + y**2 <= radius**2 
    kernel[mask] = 1 
    print('Commence circular filter.'); start = time.time() 
    tile['values'] = generic_filter(tile['values'], np.mean, footprint=kernel) 
    print('Took {:.3f} seconds'.format(time.time()-start)) 

私もここからクラスタリング技術を見ていた：http://geoffboeing.com/2014/08/clustering-to-reduce-spatial-data-set-size/

しかし、私はこれらのクラスタリング技術は完全に異なる目的を果たす実現。