python - Median filter produces unexpected result on FITS file -

this based on couple of other questions haven't quite been answered, i've started new post. i'm working on finding median of masked array in 50-pixel patches. image , mask both 901x877 telescope images.

import numpy np import matplotlib.pyplot plt astropy.io import fits  #   use fits files input image , mask hdulist = fits.open('xbulge-w1.fits') w1data = hdulist[0].data hdulist3 = fits.open('xbulge-mask.fits') mask = 1 - hdulist3[0].data w1masked = np.ma.array(w1data, mask = mask)  #   use general arrays input image , mask #w1data = np.arange(790177).reshape(901,877) #w1masked = np.ma.masked_inside(w1data, 30000, 60000)   side = 50                                                                                                                                                                                 w, h = w1data.shape                                                                                                                                                                                             width_index   = np.array(range(w//side)) * side                                                                                                                                                              height_index  = np.array(range(h//side)) * side                                                                                                                                                               def assign_patch(patch, median, side):                                                                                                                                                                          """break loop out prevent 4 nested 'for' loops"""                                                                                                                                                   j in range(side):                                                                                                                                                                                           in range(side):                                                                                                                                                                                           patch[i,j] = median                                                                                                                                                                                 return patch                                                                                                                                                                                             width in width_index:                                                                                                                                                                                       height in height_index:                                                                                                                                                                                     patch  = w1masked[width:width+side, height:height+side]                                                                                                                                                           median = np.median(patch)                                                                                                                                                                                   assign_patch(patch, median, side)                                                                             plt.imshow(w1masked) plt.show() 

the problem is, when use general arrays input image , mask (the commented out section), works fine, when use fits files, produces 'side'-sized patches on output image. can't figure out what's going on this.

i don't know how fits files there several things standing out:

• np.median doesn't take mask account. in fact in recent numpy releases (correctly) prints warning if attempted. should using np.ma.median instead. if update numpy you'll see this:

userwarning: warning: 'partition' ignore 'mask' of maskedarray.

• the assign_patch function unnecessary when know can use slice assignment:

  w1masked[width:width+side, height:height+side] = median # instead of "assign_patch(patch, median, side)" 

that's faster doing double loop replace each value.

i assume issue in fact because use np.median instead of np.ma.median. there lots of values masked pixel have including nan, 0, inf, ... if these taken account (when should ignored) produce kind of problems, if median starts returning nans or similar.

---

more if wanted median filter can't calculate median of patch , replace values in patch median. should using median filter takes mask account. unfortunately i've never seen such filter implemented in wide-spread python package. if have numba checkout (very experimental!) package of mine numbamisc contains median_filter takes masks account.