# Applying a filter to a large array with elements that are not regularly spaced

I need to apply a Gaussian filter to a 2D numpy array where the distance between adjacent array elements depends on the row of the array. (Specifically, the data are evenly spaced in latitude and longitude but are not evenly spaced in terms of distance on the surface of the sphere.)

This means that I need a different filtering array for each row of data. My approach has been to pre-compute the filtering array for each row of pixels, then pass these arrays to SciPy's scipy.ndimage.generic_filter, along with a function that selects the appropriate filtering array.

import gc
import numpy as np
from scipy.ndimage import generic_filter

degs = ['deg', 'degs', 'degree', 'degrees']
coord_array = np.pi/2 - coord_array
elif angle in degs:
coord_array = 90 - coord_array
else:
print('angle must be radians or degrees!')
quit()

return coord_array

'''
Calculate the haversine-based distance between two points on the surface of a sphere. Should be more accurate than the arc cosine strategy. See, for example: http://en.wikipedia.org/wiki/Haversine_formula
'''

degs = ['deg', 'degs', 'degree', 'degrees']
degrees = False
elif angle in degs:
degrees = True
else:
print('angle must be radians or degrees!')
quit()

spherical_distance = 2.0 * sphere_radius * np.arcsin(np.sqrt( ((1 - np.cos(theta_2-theta_1))/2.) + np.sin(theta_1) * np.sin(theta_2) * ( (1 - np.cos(phi_2-phi_1))/2.)  ))

return spherical_distance

class Surface(object):
'''
the surface class defines a surface on a segment of a spherical shell. it
contains a numpy meshgrid of coordinates in theta and phi, and a
corresponding array of the surface values at each coordinate pair.

should passed as [rs, thetas, phis] (or equivalent ndarray)
'''

def __init__(self, data, coords='spherical', r=6371000, angle='radians'):

degs = ['deg', 'degs', 'degree', 'degrees']
pass
elif angle in degs:
else:
print("angle must be 'radians' or 'degrees'!")
quit()

self.surface = np.array(data, dtype="float16")
self.thetas = np.array(data, dtype="float16")
self.phis = np.array(data, dtype="float16")

self.r = r

self.dtheta = self.thetas[1, 0] - self.thetas[0, 0]
self.dphi = self.phis[0, 1] - self.phis[0, 0]

def smooth(self, sigma=50000, n_sigma=4):

# here we construct the 'private' filtering function
## here 'middle_theta' is (2*m)+1 due to the way the flattening works
def spherical_gaussian_filter(image_slice, filter_arrays, middle_theta):
n_middle_theta = int(image_slice[middle_theta])
return (image_slice * filter_arrays[n_middle_theta]).sum()

# we can get the required height of the Gaussian array immediately,
# since degrees of latitude are evenly spaced (on a perfect sphere)
## first get the required change in theta to cover n*sigma on the sphere
Delta_theta = n_sigma*sigma/self.r

## next, get the required size of the array's outer dimension (its
## vertical size) AWAY from the centre pixel (i.e. the total size of the
## outer dimension is 2m+1)
## applying ceil guarantees we are at least as good as n*sigma
m = np.ceil(Delta_theta/np.abs(self.dtheta)).astype('int')

# now get the maximum required width of the Gaussian array
## the maximum width will be at the most extremal value of theta
## convert to latitude so finding the angle farthest from the equator is
## easy
latitudes = convert_between_latitude_and_polar_angle(self.thetas)
max_theta = convert_between_latitude_and_polar_angle(np.max(np.abs(latitudes[..., 0])))
del latitudes
gc.collect()

## now get the required change in phi
Delta_phi = n_sigma*sigma/(self.r*np.sin(max_theta))

print(np.degrees(Delta_theta))
print(np.degrees(Delta_phi))

## and then the size of the inner dimension of the array (away from the
## centre pixel)
n = np.ceil(Delta_phi/self.dphi).astype('int')

# next create all the arrays of coordinates around the central element
# note that we only need an array for each row of pixels, since they
# are independant of longitude
## first make arrays corresponding to the coordinates of the 'central'
## pixel for each row
## the None suffix in the array slicing functions like numpy.expand_dim
central_thetas = np.broadcast_to(self.thetas[..., 0, None, None], (len(self.thetas), 2*m+1, 2*n+1))
central_phis = np.zeros_like(central_thetas, dtype="float16")

## create the arrays of distances (relative coordinates) in theta and
## phi from the central pixel
relative_theta = np.arange(-m, m + 1, dtype="float16")*np.abs(self.dtheta)
relative_phi = np.arange(-n, n + 1, dtype="float16")*self.dphi
relative_phis, relative_thetas = np.meshgrid(relative_phi, relative_theta)
del relative_phi, relative_theta
gc.collect()

## now we can get the actual coordinates of each element of the
## filtering array
absolute_thetas = (central_thetas + relative_thetas).astype("float16")
del relative_thetas
gc.collect()
absolute_phis = (central_phis + relative_phis).astype("float16")
del relative_phis
gc.collect()

# next we need the distance on the surface of the sphere from the
# central pixel to each element
## use the Haversine formula, slightly adapted from the formulation used
## in 2017 for my honours project
print("About to do distances between points")
print(absolute_thetas.shape)
distances_from_central_pixel = calculate_haversine_distances_between_theta_phi_points(absolute_thetas, absolute_phis, central_thetas, central_phis, self.r)
del absolute_thetas, absolute_phis, central_thetas, central_phis
gc.collect()

# now we need to calculate the gaussian weights of each distance array
gaussian_weights = ( 1 / (sigma*np.sqrt(2*np.pi)) ) * np.exp( -0.5 * (distances_from_central_pixel / sigma)**2 )
# and normalise
sums = gaussian_weights.sum(axis=(1,2))[..., None, None]
gaussian_weights /= sums
del sums
gc.collect()
# print(gaussian_weights, end='\n\n')

# we need to add in the zeros that will take care of the theta_ns in the
# surface that is passed to generic_filter
gaussian_weights = gaussian_weights.flatten()
gaussian_weights = np.insert(gaussian_weights, np.arange(1, len(gaussian_weights) + 1), np.zeros(len(gaussian_weights)))
gaussian_weights = gaussian_weights.reshape(len(self.thetas), 2*(2*m+1)*(2*n+1))

# lastly, need to create an array with the n of each theta at each point
# on the surface
surface = np.copy(self.surface).flatten()
theta_ms = np.arange(len(self.surface)).reshape(len(self.surface), 1)
theta_ms = np.repeat(theta_ms, len(self.surface)).flatten()
surface_with_theta_ms = np.insert(surface, np.arange(1, len(surface) + 1), theta_ms)
surface_with_theta_ms = surface_with_theta_ms.reshape(len(self.surface), len(self.surface), 2)
del surface, theta_ms
gc.collect()
# print(surface_with_theta_ms, end='\n\n')

# pass generic_filter an array which is pretending to be 2d but each element actually contains the height, theta
self.filtered_surface = generic_filter(surface_with_theta_ms.astype("float16"), spherical_gaussian_filter, size=(2*m+1, 2*n+1, 2), extra_arguments=(gaussian_weights, (2*m)+1))[:,:,0]
del surface_with_theta_ms
gc.collect()

if __name__ == '__main__':
m, n = 300, 400
x, y = np.meshgrid(np.linspace(-15, 15, num=n), np.linspace(120, 150, num=m))
surface = np.ones((m,n))
surface[:,int(n/2):] = surface[:,int(n/2):] * 0
test_data = np.array([surface, y, x])
print(test_data, end='\n\n')
test_surface = Surface(test_data, angle='deg')
test_surface.smooth(sigma=50000, n_sigma=4)
print(test_surface.filtered_surface)


It works fine for small arrays of data, but when I scale it up to run on my data (~ 3000 x 2000 pixels), my pre-computed filtering arrays are large, and there's around 2000 of them, resulting in using up all 22 GB of free RAM.

As you can see I've made some efforts to reduce memory usage already by using less accurate data types in my ndarrays by using gc.collect(), but the program still uses all 22 GB of RAM when it tries to calculate the distances on the surface of the sphere for each filtering array.

This is my first time posting a question to Code Review, so I appreciate that it might be inappropriate as written, but I really am at a loss for how to reduce memory usage further.

• @Graipher Thanks for the feedback! I've included the two functions as well as the test data I was running it on – Notso Sep 26 '18 at 10:02