Numpy gradient in masked array

Question

When I compute the gradient of a masked array in numpy as

import numpy as np
import numpy.ma as ma
x = np.array([100, 2, 3, 5, 5, 5, 10, 100])
mx = ma.masked_array(x, mask=[1, 0, 0, 0, 0, 0, 0, 1])

the mask of the resulting array is different from the original mask:

np.gradient(mx) 
masked_array(data = [-- -- 1.5 1.0 0.0 2.5 -- --],
         mask = [ True  True False False False False  True  True],
   fill_value = 999999)

Why is the gradient not computed at the 'new' boundaries? How can I change that?

For masked entries in the middle it gets weirder:

x = np.array([100, 2, 3, 5, 5, 6, 6, 6, 6, 6, 7, 5, 10, 100])
mx = ma.masked_array(x, mask=[1, 0, 0, 0,0, 0,1,0,0,0, 0, 0, 0, 1])
np.gradient(mx)
masked_array(data = [-- -- 1.5 1.0 0.5 -- 0.0 -- 0.0 0.5 -0.5 1.5 ----],
mask = [ True  True False False False  True False  True False False False False
True  True],
fill_value = 1e+20)

I would expect np.gradient to just treat the masked cells as boundaries.

Update:

What I want to do: I need to compute the gradient on the array and not change the mask, nor the shape of the array (I want 2d in the end) Masked cells should not contribute to the gradient. Points next to a masked cell should be considered boundaries and a one-sided difference applied.

- - - - - - - - - - -
- - - - o o - - o - -
- - - o x x o o x o -
- - - o x o - - o - -
- - - - o - - - - - -
- - - - - - - - - - -

In this sketch x represent a masked cell and o are cells where a one-sided difference should be computed (cells at the edges of the area need to be one-sided too, but I don't paint them in here for clarity).

Hi, I'm having the same issue (in 2d like you). Did you ever find a (performant) solution for this? — Kroppeb
– Kroppeb, Commented Feb 17, 2023 at 14:30

York's · Accepted Answer · 2018-04-29 19:34:18Z

2

I think the reason is that you expect masked elements of the mx array to be skipped during the computation of the gradient, so that instead of computing gradient on x = np.array([100, 2, 3, 5, 5, 5, 10, 100]) we will compute it on x = np.array([2, 3, 5, 5, 5, 10]), but the real behavior is different in a way that np.ma.MaskedArray is inherited from np.ndarray and np.gradient() doesn't do anything special with np.ndarray or its subclasses. So, in the case of mx = ma.masked_array(x, mask=[1, 0, 0, 0, 0, 0, 0, 1]), the gradient will be computed on array:

[--, 2, 3, 5, 5, 5, 10, --]

First it will try to calculate the gradient for the first element with by default first order difference and step h=1(np.gradient by default treats steps as unary in each dimension of the input. ):

gradient[0] = (mx[1] - mx[0]) / h

Because, it depend on mx[0] which is not allowed to use by mask, the value of gradient[0] will be masked with 'True`.

When it will try to compute gradient at index 1 for the element that you percept as the new left boundary of an array; that element is not in fact a boundary element of ndarray. When it computes the gradient for the element which you think is the new left boundary, it will actually use central differences formula with homogeneous steps h, gradient[1] = (mx[2] - mx[0]) / 2h, but because mx[0] is masked as the one that cannot use, the value for the gradient[1] cannot be retrieved either, so it is masked with True. Same happens on the other end of an array.

Now with regards to masking something in the middle, assume:

x = np.array([1, 2, 3, 4, 5])
mask = [0, 0, 1, 0, 0]
mx = ma.masked_array(x, mask=mask)

Gradient function applied to mx again will use central difference formula for homogeneous steps, and when we calculate:

masked_gradient[1] = (mx[2] - mx[0]) / 2h
masked_gradient[3] = (mx[4] - mx[2]) / 2h

Neither of these values can be calculated, because mx[2] is masked with True. And in the same time:

masked_gradient[2] = (mx[3] - mx[1]) / 2h

Can be evaluated, because all of values it depends on are masked as False, so the result will have a mask:

[0, 1, 0, 1, 0]

And values:

[1.0, --, 1.0, --, 1.0]

edited Apr 29, 2018 at 19:34

answered Apr 29, 2018 at 19:21

York's

1741 silver badge11 bronze badges

Sign up to request clarification or add additional context in comments.

5 Comments

Sebastian Beyer Over a year ago

That explains it very nicely! Thank you! Do you also have an idea how I could achieve my expected behavior? Just do a subset (slice) before I use np.gradient? Like np.gradient(x[1:-1]) and then put that back into the original array? But that will not work for the interior points...

York's Over a year ago

@SebastianBeyer for 1D array you can use MaskedArray.compressed() method which produces flattened 1D array of all non masked values of the MaskedArray and then you can compute np.gradient() on it. This should work well for the edge and interior points of 1D array: @SebastianBeyer for 1D array you can use MaskedArray.compressed()` method which produces flattened 1D array of all non masked values of the MaskedArray and then you can compute np.gradient() on it. This should work well for the edge and interior points of 1D array: x = np.array([100, 2, 3, 5, 5, 5, 10, 100])

York's Over a year ago

Continue: mx = ma.masked_array(x, mask=[1, 0, 0, 0, 1, 0, 0, 1]) np.gradient(mx.compressed()) and the result is: array([ 1. , 1.5, 1. , 2.5, 5. ]) which is computed on 5 non masked points. This want work for arrays of rank higher than 1, because it flattens everything to 1D. I think the reason it flattens everything to 1D array is because for any point at each dimension we can either mask value or not, and thus we have a question how to present it and work with it?

York's Over a year ago

So, to sum up: if the rank of the array is 1 then I would suggest using above solution; rank > 1, then you should define what you want more precisely.

Sebastian Beyer Over a year ago

Thank you, but I need to preserve the shape of the array. I have updated my question to make it clear what I want to do :-)

Davis Herring · Accepted Answer · 2018-04-29 19:06:22Z

0

It’s computing a central difference, which cannot be constructed adjacent to a boundary (or missing value). numpy.gradient switches to a one-sided difference at the ends of the array, but the masked array’s holes don’t count as ends.

Your second example shows the hole in the stencil: you can compute it at an isolated hole in the data.

answered Apr 29, 2018 at 19:06

Davis Herring

41.9k4 gold badges58 silver badges91 bronze badges

Collectives™ on Stack Overflow

Numpy gradient in masked array

2 Answers 2

5 Comments

Comments

Your Answer

Hot Network Questions

Collectives™ on Stack Overflow

2 Answers 2

5 Comments

Comments

Your Answer

Sign up or log in

Post as a guest

Related