Use scipy ndimage generic_filter to compute set intersection across pixel-neighborhoods

Question

I have raster images reporting the start day and end day of some periods across the globe. For each pixel (read numpy array element) I'd like to get a measure of how much its day-unit period overlap/intersect with that of the neighboring pixels. Formally, I want to obtain:

X_p = Σ_{k∈N_p} | T_p ∩ T_k|

where: T_p = {S_p , S_p + 1, ..., E_p} is the period in pixel p and N_p is the set identifying the neighbor of pixel p.

I've tried to compute my measure using a custom function in scipy.ndimage.generic_filter(). My first approach was to work directly with an array reporting T_p sets; this didn't work because I think that the generic_filter() function does not accept arrays with sets (I got RuntimeError: array type not supported).
I have thus tried to work with a custom function that: (i) has the two arrays for period start and end as inputs, (ii) builds the T sets and (iii) computes their intersection. However, I don't know how to pass it correctly into the generic_filter(). Below is my code:

import numpy as np
import scipy.ndimage as ndimage

# Preliminaries
def calendar(s,e):    # build the T array
    return set(range(s, e))
vect_calendar = np.vectorize(calendar)

s = np.array([[99, 99, 1], [2, 2, 4], [1, 99, 99], [99, 3, 3]]) # start day of T (99=missing day)
e = np.array([[0, 0, 3], [6, 4, 7], [3, 0, 0], [0, 4, 5]])
t = vect_calendar(s,e)
x_expected = np.array([[0, 0, 2], [6, 4, 2], [2, 0, 0], [0, 1, 1]])

# Function for computing intersection across pixels (note I exclude the central pixel from computations)
ns = 9
ray = int((ns**(1/2)-1)/2)

def neig_func(arr_s, arr_e):
    p = int((ns-1)/2)  # central pixel of the neighborhood
    t_p = set(range(int(arr_s[p]), int(arr_e[p])))
    
    agg = []
    for k in range(ns):
        t_k = set(range(int(arr_s[k]), int(arr_e[k])))
        f_pk = set.intersection(t_p, t_k)
        f_pk = len(f_pk)
        agg.append(f_pk)
    return np.sum(agg) - len(t_p)

# Run the function across all pixel neighbors except those on the top & bottom edges of the array
s_bord = np.pad(s, pad_width=((ray,ray),(0,0)), mode="constant", constant_values=0)
e_bord = np.pad(e, pad_width=((ray,ray),(0,0)), mode="constant", constant_values=0)
e_bord_flat = e_bord.flatten()

neig = np.ones((int(ns**(1/2)), int(ns**(1/2))))
x = ndimage.generic_filter(s_bord, neig_func, footprint=neig, mode="wrap", extra_arguments=(e_bord_flat,))
print("Input data is: \n", t, "\n\n****\n Expected output is: \n", x_expected, "\n\n****\n Actual output is \n", x[ns:-ns])

The output is an array with all zeros! I guess that the problem is that the array e is read by the generic_filter() always in its first ns elements.
An alternative solution could be have a 3-d array with s and e stacked one over the other. Yet, it seems that I can't run the generic_filter() on 3-dimensional arrays.

Use scipy ndimage generic_filter to compute set intersection across pixel-neighborhoods

0 Answers0