I have a Python code that works, but it's quite slow and I believe there has to be a way of doing this more efficiently.
The idea is to apply a filter to an image. The filter is an average of the points which fall within a specified radius. The input is a mx2
array representing x,y and mx1
array z, representing coordinates of m observation points.
The program that works is the following
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.)
The problem is that I have 6 million points per data set, and the code is horribly slow. There must be a way to do this more efficiently. I thought of using scipy.spatial.distance.cdist() which is fast, but then I'd have to reproject the data to UTM, and I'd like to avoid reprojection. Any suggestions?
Thanks, Reniel