Why does np.ma.median work for only certain sizes of nan matrices and give the MaskError('Cannot alter the masked element.') for other sizes

Question

I am trying to build a code that can be given an matrix of nans and not break. I wrote a simplified version here that recreates the problem. I am using numpy version 1.23.5

Basically the code takes a grid of nans and creates a annulus mask to apply to the grid and then take the median of what ever is in the annulus. This code does not produce the error as is.

from math import floor
from scipy.signal import convolve
import numpy as np

in_radius = 10
out_radius = 31
box_size = 5

box=np.ones((box_size,box_size))
padding=floor(box_size/2)

y,x=np.ogrid[-1*out_radius:out_radius+1,-1*out_radius:out_radius+1]
grid_mask=np.logical_or(x**2+y**2<in_radius**2,x**2+y**2>out_radius**2)

grid = np.empty((67,67),)
grid[:] = np.nan

smoothed=convolve(grid,box,method='direct',mode='valid') / 2

noisepix=np.ma.masked_array(smoothed,grid_mask)
median_signal = np.ma.median(noisepix)

The last line will produce a nan, which is fine.

However if I change out_radius to anything smaller than 28, numpy is unable to perform the calculation of the median. Here is the code that produces the error:

in_radius = 10
out_radius = 28
box_size = 5

box=np.ones((box_size,box_size))
padding=floor(box_size/2)

y,x=np.ogrid[-1*out_radius:out_radius+1,-1*out_radius:out_radius+1]
grid_mask=np.logical_or(x**2+y**2<in_radius**2,x**2+y**2>out_radius**2)

grid = np.empty((61,61),)
grid[:] = np.nan

smoothed=convolve(grid,box,method='direct',mode='valid') / 2

noisepix=np.ma.masked_array(smoothed,grid_mask)
median_signal = np.ma.median(noisepix)

The full error message is:

MaskError                                 Traceback (most recent call last)
/var/folders/89/zgv9nv5563n89pzs55bv6mzm0000gq/T/ipykernel_17229/2791888929.py in <module>
     21 
     22 noisepix=np.ma.masked_array(smoothed,grid_mask)
---> 23 median_signal = np.ma.median(noisepix)

~/opt/anaconda3/lib/python3.9/site-packages/numpy/ma/extras.py in median(a, axis, out, overwrite_input, keepdims)
    733             return m
    734 
--> 735     r, k = _ureduce(a, func=_median, axis=axis, out=out,
    736                     overwrite_input=overwrite_input)
    737     if keepdims:

~/opt/anaconda3/lib/python3.9/site-packages/numpy/lib/function_base.py in _ureduce(a, func, **kwargs)
   3723         keepdim = (1,) * a.ndim
   3724 
-> 3725     r = func(a, **kwargs)
   3726     return r, keepdim
   3727 

~/opt/anaconda3/lib/python3.9/site-packages/numpy/ma/extras.py in _median(a, axis, out, overwrite_input)
    779             if not odd:
    780                 s = np.true_divide(s, 2., casting='safe', out=out)
--> 781             s = np.lib.utils._median_nancheck(asorted, s, axis)
    782         else:
    783             s = mid.mean(out=out)

~/opt/anaconda3/lib/python3.9/site-packages/numpy/lib/utils.py in _median_nancheck(data, result, axis)
   1052             return data.dtype.type(np.nan)
   1053 
-> 1054         result[n] = np.nan
   1055     return result
   1056 

~/opt/anaconda3/lib/python3.9/site-packages/numpy/ma/core.py in __setitem__(self, indx, value)
   3344         """
   3345         if self is masked:
-> 3346             raise MaskError('Cannot alter the masked element.')
   3347         _data = self._data
   3348         _mask = self._mask

MaskError: Cannot alter the masked element.

I am trying to understand why this is happening for one size of outer radius and not the other.

I know that I can work around this by replacing the array of nans with zeros or something and just compensating for it. I am more after why this is happening in some case and not in others