python - Efficient adaptative (row dependant) dilation of a numpy array

Question

Let's say I have a 2D array reprezenting equally-spaced (in °) lat/lon points, the "top" of the array being the North, the left being the West, etc. For instance with North = 55.4°, South = 37.5°, West = -12°, East = 16°, and resolution = 0.025°, the shape of the 2D array will be (717, 1121).

Now, considering some points of this 2D array, I would like to expend/dilate their value through a circular neighborhood in km, but since the distance in km between two lat/lon points depends on the latitude, the "array reprezentation" of the neighborhood will more be an ellipse than a disk, and beside will depend on the row of the 2D array. Indeed, the more we get closer to latitudes +/-90°, the more we'll need points to cross an East - West distance in km.

To perform that, I start by "binarizing" my 2D array, affecting the value "1" to the points where I want to apply the neighborhood. Once done I loop over the different rows (eg the latitudes), compute the "array reprezentation" of the circular neighborhood, then use it with function binary_dilation() from scipy.ndimage through argument 'structure' to perform the dilation. Here is my code:

import numpy as np
from geopy.distance import distance
from scipy import ndimage

def circular_mask(r):
    mask = np.zeros((2*r+1, 2*r+1))
    y,x = np.ogrid[-r:r+1, -r:r+1]
    disk = x**2 + y**2 <= r**2
    mask[disk] = 1
    return(mask)
def elliptic_mask(a,b):
    #NB: 'a' is the semi-minor axis, 'b' the semi-major axis
    mask = np.zeros((2*a+1, 2*b+1))
    y,x = np.ogrid[-a:a+1, -b:b+1]
    ellipse = (x/b)**2 + (y/a)**2 <= 1
    mask[ellipse] = 1
    return(mask)

North = 55.4
South = 37.5
West = -12.
East = 16.
Resol = 0.025
lon = West  #random longitude (the distance in km between tow lat/lon points depends only on the lat)
radius = 30.  #(radius in km for circular neighborhood)

latitudes = np.arange(South, North + Resol/2., Resol)  #list of different latitudes
array = np.random.randint(20, size = (717,1121))  #((North-South)/Resol)+1, ((East-West)/Resol)+1
array_01 = np.where(array >= 15, 1, 0)  #[binarization] dilate points >= 15
array_01_dilated = array_01.copy()  #will be filled during the loop

idx_row = array.shape[0]
for lat in latitudes:
    idx_row = idx_row + 1  #first latitude corresponds to South = the bottom of the array

    EW_resol_to_km = distance((lat, lon), (lat, lon + lon_resol)).km #number of km between two points side by side on latitude 'lat'
    EW_nb_resol_to_reach_radius = radius/EW_resol_to_km #number of side by side points on latitude 'lat' needed to reach distance 'radius'
    EW_nb_points = int(round(EW_nb_resol_to_reach_radius)) #convert into integer
    NS_resol_to_km = distance((lat, lon), (lat + lat_resol, lon)).km #number of km between two points side by side on longitude 'lon'
    NS_nb_resol_to_reach_radius = radius/NS_resol_to_km #number of side by side points on longitude 'lon' needed to reach distance 'radius'
    NS_nb_points = int(round(NS_nb_resol_to_reach_radius))

    if (NS_nb_points != EW_nb_points):
        ellipse = elliptic_mask(a = min(NS_nb_points,EW_nb_points), b = max(NS_nb_points,EW_nb_points)) #array reprezentation of circular neighborhood (perfectible)
    else:
        ellipse = circular_mask(radius = NS_nb_points)

    neighborhood_band = np.zeros((ellipse.shape[0], array.shape[1]))
    neighborhood_band[NS_nb_points + 1, :] = array_01[idx_row, :]
    neighborhood_band = ndimage.binary_dilation(neighborhood_band, structure = ellipse) #corresponds to the latitude 'lat' and the nearest ones modified due to dilation

    # fill 'array_01_dilated', taking into account possible borderies overtakings #
    idx_row_top_band = idx_row - NS_nb_points
    idx_row_bottom_band = idx_row + NS_nb_points
    if (idx_row_top_band < 0):
        neighborhood_band = neighborhood_band[(-idx_row_top_band):,:]
        idx_row_top_band = 0
    if (idx_row_bottom_band > array.shape[0]-1):
        neighborhood_band = neighborhood_band[:-(idx_row_bottom_band - (array.shape[0]-1)),:]
        idx_row_bottom_band = array.shape[0]-1
    array_01_dilated[idx_row_top_band:(idx_row_bottom_band + 1),:] = array_01_dilated[idx_row_top_band:(idx_row_bottom_band + 1),:] + neighborhood_band

array_binarized_neighborhood[array_binarized_neighborhood >= 1] = 1.  #for possible "dilation coverings"

This method being not very fast, I would like to know if there is another way to solve this problem (using scipy.ndimage.binary_dilation() with a row depending structure? a vectorized approach? ...).

Thank you very much !

question from:https://stackoverflow.com/questions/66045372/efficient-adaptative-row-dependant-dilation-of-a-numpy-array

Answer 1

Waitting for answers