Q:

Create a matrix of neighbours

I have a sparse matrix and I need to create a new neighbor matrix of each index. 
Below I leave a representation of the data in the NxM matrix.  For each of the elements of the matrix I need to obtain the neighbors in a section of KxK. With this information, it would generate a NMxKK matrix that contains in each row the indices of the neighboring KKs of the element. 

I asked a similar question a while ago but the difference is that now the data is structured, so I can do without KdTree. 
This new matrix is ​​used to calculate the distance of non-zero neighbors, and with these distances associate a weight to each neighbor, to finally estimate the desired value as a weighted average of the neighbors. 
Thanks in advance! 
UPDATE
I have data like the ones in the image (generated with the function generate_data) and I need to perform the following operation. 

Given a filter / kernel / NxN matrix, with N being the kernel size defined by me, calculate for nonzero values the distances with respect to the central pixel.  Take as an example the value 20 that is in the position (1, 8) of the image.  Taking a matrix of 5x5, the nonzero values of interest are 40 (in (0, 6)), 37 (in (1, 6)) and 25 (in (3, 10)), with distances 2.23606798, 2 and 2.82842712 respectively (obtained making the Euclidean norm between the indices). 
What I need to get in this step is the matrix res:
[[0.          2.23606798 2.          0.          0.         ]
 [0.          0.          0.          0.          0.         ]
 [0.          0.          1.          0.          0.         ]
 [0.          0.          0.          0.          0.         ]
 [0.          0.          0.          0.          2.82842712]]

I need to obtain the 1.  in the center of the matrix too to also take into account the value where I am standing (whose distance to itself is 0.). 
With these values, I get the mask with non-zero values and calculate the weights based on a Gaussian distribution:
import scipy.stats as st 
mask = 0 < res
gauss = st.norm.pdf(res) # or st.norm.pdf(mask * kernel(5))

[[0.         , 0.03274718, 0.05399097, 0.         , 0.         ],
 [0.         , 0.         , 0.         , 0.         , 0.         ],
 [0.         , 0.         , 0.39894228, 0.         , 0.         ],
 [0.         , 0.         , 0.         , 0.         , 0.         ],
 [0.         , 0.         , 0.         , 0.         , 0.00730688]])

total = gauss.sum() # 0.4929873057962355

Finally, with these weights, I calculate the weights and the final value of the pixel by interpolating the values. 
val[1, 8] = 0.03274718 * 40 / total + 0.05399097 * 37 / total + 0.39894228 * 20 / total + 0.00730688 * 25 / total

The same thing I must do for each pixel (I guess I have to add a kernel_size padding // 2 to be able to use the whole array). 
Here is my script
import matplotlib.pylab as plt
import numpy as np
import scipy.stats as st

from scipy import sparse

def generate_data(m, n, density):
    s = 64 * sparse.random(m, n, density=density).A
    return s.astype(np.int8)

def plot_matrix(matrix):
    for (j, i), label in np.ndenumerate(s):
        plt.text(i, j, label, ha='center', va='center')

    plt.imshow(matrix)
    plt.show()

def kernel(n):
    n = n if n % 2 ! = 0 else n + 1
    mid = n // 2

    m = np.ndarray((n, n, 2))
    for i in range(n):
        for j in range(n):
            m[i, j] = np.array([i, j])

    return np.linalg.norm(m - [mid, mid], axis=2)

s = generate_data(10, 14, 0.25)
plot_matrix(s)

A:

This was really simple, although maybe not very efficient.  What I had to do was two convolutions:

In the first, it was convolucionar the Gaussian kernel with the matrix
conv_1 = convolve2d(m * mask_clean, k_gauss)
In the second, the Gaussian kernel with the mask
conv_2 = convolve2d(mask_clean, k_gauss)

In each position, conv_1 would have the sum of each value weighed by the corresponding factor of the Gaussian kernel.  conv_2 would have in each position the sum of all nonzero values.  The only thing left to do was divide them to get the final result
# m have the data
mask_clean = (0 < m) & (m_mean - 3*m_std < m) & (m < m_mean + 3*m_std)

# Custom function to create a gaussian kernel
k = gkern(kernlen=5, std=5//2)
k_gauss = st.norm.pdf(k)

conv_1 = convolve2d(m * mask_clean, k_gauss)
conv_2 = convolve2d(mask_clean, k_gauss)
final = conv_1 / conv_2

