I have a set of spatial distributed measurements. For each point p1 = [x1,y1,z1] there is a measurement v1 which is a scalar. (e.g. Temperature measurements under water.) Lets assume these measurements are on a regular grid. I would like to find out where is the most variation in this distribution. That means in what positions is the most change of temperature. I think this corresponds to the spatial derivation of temperature. Can somebody give me an advice how to proceed? What are methodologies to archive this? I tried to implement it with np.gradient() but i fail at interpreting the result...

CodePudding user response：

This is absolutely not an optimized code, but here is what I came up with, at least to explain how it works.

grid = [[[1, 2], [2, 3]], [[8, 5], [4, 1000]]]

def get_greatest_diff(g, x, y, z):
    value = g[x][y][z]
    try:
        diff_x = abs(value-g[x 1][y][z])
    except IndexError:
        diff_x = -1
    try:
        diff_y= abs(value-g[x][y 1][z])
    except IndexError:
        diff_y = -1
    try:
        diff_z = abs(value-g[x][y][z 1])
    except IndexError:
        diff_z = -1
    if diff_x>=diff_y and diff_x>=diff_z:
        return diff_x, [x 1, y, z]
    if diff_y>diff_x and diff_y>=diff_z:
        return diff_y, [x, y 1, z]
    return diff_z, [x, y, z 1]

greatest_diff = 0
greatest_diff_pos0 = []
greatest_diff_pos1 = []
for x in range(len(grid)):
    for y in range(len(grid[x])):
        for z in range(len(grid[x][y])):
            diff, coords = get_greatest_diff(grid, x, y, z)
            if diff > greatest_diff:
                greatest_diff = diff
                greatest_diff_pos0 = [x, y, z]
                greatest_diff_pos1 = coords
print(greatest_diff, greatest_diff_pos0, greatest_diff_pos1)

The try:...except:... are here to handle the edge conditions. (That's dirty but that's quick!)

For each cell, you will look at the three neighbours x 1 or y 1 or z 1 and you compute the difference with their values. You keep the largest difference in the neighborhood and you return it. (That is the explanation of get_greatest_diff)

In the main loop, you check if the difference in this neighborhood is the greatest of all, if so, store the difference, and the two cells in question.

Finally, return the greatest difference and the cells in question.

CodePudding user response：

Here is a numpy solution that returns the indices in an ndarray that has the biggest total differences with its neighbors. Say the input array is X and it is 2D. I will create D where D[i,j] = |X[i, j]-X[i-1, j]| |X[i,j]-X[i, j-1]|. And return the indices of D which give the largest value in D.

def greatest_diff(X):
    
    ndim = X.ndim
    Ds = [np.abs(np.diff(X, axis = i, prepend=0)) for i in range(ndim)]
    D = sum(Ds)
    return np.unravel_index(D.argmax(), D.shape)

X = np.zeros((5,5))
X[2,2] = 1
greatest_diff(X)
# returns (2, 2)

X = np.zeros((5,10,9))
X[2,2,7] = -1
greatest_diff(X)
# returns (2, 2, 7)

Another solution might be calculating the difference between X[i, j] and sum(X[k, l]) where k,l are the neighbors of i, j. You can achieve this by applying a gaussian filter to the X say gX then taking the squared differences: (X-gX)^2.

def greatest_diff_gaussian(X, sigma = 1):
    from scipy.ndimage import gaussian_filter
    gX = gaussian_filter(X, sigma)
    dgX = np.power(X - gX, 2)
    return np.unravel_index(dgX.argmax(), dgX.shape)