1
\$\begingroup\$

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

def convert_between_latitude_and_polar_angle(coord_array, angle='radians'):

    rads = ['rad', 'rads', 'radian', 'radians']
    degs = ['deg', 'degs', 'degree', 'degrees']
    if angle in rads:
        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

def calculate_haversine_distances_between_theta_phi_points(theta_2, phi_2, theta_1, phi_1, sphere_radius, angle='radians'):
    '''
    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
    '''

    rads = ['rad', 'rads', 'radian', 'radians']
    degs = ['deg', 'degs', 'degree', 'degrees']
    if angle in rads:
        degrees = False
    elif angle in degs:
        degrees = True
        theta_2 = np.radians(theta_2)
        phi_2 = np.radians(phi_2)
        theta_1 = np.radians(theta_1)
        phi_1 = np.radians(phi_1)
    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'):

        rads = ['rad', 'rads', 'radian', 'radians']
        degs = ['deg', 'degs', 'degree', 'degrees']
        if angle in rads:
            pass
        elif angle in degs:
            data[1] = np.radians(data[1])
            data[2] = np.radians(data[2])
        else:
            print("angle must be 'radians' or 'degrees'!")
            quit()

        self.surface = np.array(data[0], dtype="float16")
        self.thetas = np.array(data[1], dtype="float16")
        self.phis = np.array(data[2], 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[0])).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[0]), 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[0], 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.

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.