Proper treatment of NaN values (ignoring them during convolution and
replacing NaN pixels with interpolated values)

A single function for 1-D, 2-D, and 3-D convolution

Improved options for the treatment of edges

Both direct and Fast Fourier Transform (FFT) versions

Built-in kernels that are commonly used in Astronomy

The following thumbnails show the difference between Scipy’s and
Astropy’s convolve functions on an Astronomical image that contains NaN
values. Scipy’s function essentially returns NaN for all pixels that are
within a kernel of any NaN value, which is often not the desired result.

importnumpyasnpimportmatplotlib.pyplotaspltfromastropy.ioimportfitsfromastropy.utils.dataimportget_pkg_data_filenamefromastropy.convolutionimportGaussian2DKernelfromscipy.signalimportconvolveasscipy_convolvefromastropy.convolutionimportconvolve# Load the data from data.astropy.orgfilename=get_pkg_data_filename('galactic_center/gc_msx_e.fits')hdu=fits.open(filename)[0]# Scale the file to have reasonable numbers# (this is mostly so that colorbars don't have too many digits)# Also, we crop it so you can see individual pixelsimg=hdu.data[50:90,60:100]*1e5# This example is intended to demonstrate how astropy.convolve and# scipy.convolve handle missing data, so we start by setting the# brightest pixels to NaN to simulate a "saturated" data setimg[img>2e1]=np.nan# We also create a copy of the data and set those NaNs to zero. We'll# use this for the scipy convolutionimg_zerod=img.copy()img_zerod[np.isnan(img)]=0# We smooth with a Gaussian kernel with x_stddev=1 (and y_stddev=1)# It is a 9x9 arraykernel=Gaussian2DKernel(x_stddev=1)# Convolution: scipy's direct convolution mode spreads out NaNs (see# panel 2 below)scipy_conv=scipy_convolve(img,kernel,mode='same',method='direct')# scipy's direct convolution mode run on the 'zero'd' image will not# have NaNs, but will have some very low value zones where the NaNs were# (see panel 3 below)scipy_conv_zerod=scipy_convolve(img_zerod,kernel,mode='same',method='direct')# astropy's convolution replaces the NaN pixels with a kernel-weighted# interpolation from their neighborsastropy_conv=convolve(img,kernel)# Now we do a bunch of plots. In the first two plots, the originally masked# values are marked with red X'splt.figure(1,figsize=(12,12)).clf()ax1=plt.subplot(2,2,1)im=ax1.imshow(img,vmin=-2.,vmax=2.e1,origin='lower',interpolation='nearest',cmap='viridis')y,x=np.where(np.isnan(img))ax1.set_autoscale_on(False)ax1.plot(x,y,'rx',markersize=4)ax1.set_title("Original")ax1.set_xticklabels([])ax1.set_yticklabels([])ax2=plt.subplot(2,2,2)im=ax2.imshow(scipy_conv,vmin=-2.,vmax=2.e1,origin='lower',interpolation='nearest',cmap='viridis')ax2.set_autoscale_on(False)ax2.plot(x,y,'rx',markersize=4)ax2.set_title("Scipy")ax2.set_xticklabels([])ax2.set_yticklabels([])ax3=plt.subplot(2,2,3)im=ax3.imshow(scipy_conv_zerod,vmin=-2.,vmax=2.e1,origin='lower',interpolation='nearest',cmap='viridis')ax3.set_title("Scipy nan->zero")ax3.set_xticklabels([])ax3.set_yticklabels([])ax4=plt.subplot(2,2,4)im=ax4.imshow(astropy_conv,vmin=-2.,vmax=2.e1,origin='lower',interpolation='nearest',cmap='viridis')ax4.set_title("Default astropy")ax4.set_xticklabels([])ax4.set_yticklabels([])# we make a second plot of the amplitudes vs offset position to more# clearly illustrate the value differencesplt.figure(2).clf()plt.plot(img[:,25],label='input',drawstyle='steps-mid',linewidth=2,alpha=0.5)plt.plot(scipy_conv[:,25],label='scipy',drawstyle='steps-mid',linewidth=2,alpha=0.5,marker='s')plt.plot(scipy_conv_zerod[:,25],label='scipy nan->zero',drawstyle='steps-mid',linewidth=2,alpha=0.5,marker='s')plt.plot(astropy_conv[:,25],label='astropy',drawstyle='steps-mid',linewidth=2,alpha=0.5)plt.ylabel("Amplitude")plt.ylabel("Position Offset")plt.legend(loc='best')plt.show()

convolve() is implemented as a
direct convolution algorithm, while
convolve_fft() uses a fast Fourier
transform (FFT). Thus, the former is better for small kernels, while the latter
is much more efficient for larger kernels.

For example, to convolve a 1-d dataset with a user-specified kernel, you can do:

Notice that the end points are set to zero - by default, points that are too
close to the boundary to have a convolved value calculated are set to zero.
However, the convolve() function allows for a
boundary argument that can be used to specify alternate behaviors. For
example, setting boundary='extend' causes values near the edges to be
computed, assuming the original data is simply extended using a constant
extrapolation beyond the boundary:

The values at the end are computed assuming that any value below the first
point is 1, and any value above the last point is 8. For a more
detailed discussion of boundary treatment, see Using the convolution functions.

This module also includes built-in kernels that can be imported as e.g.:

>>> fromastropy.convolutionimportGaussian1DKernel

To use a kernel, first create a specific instance of the kernel:

>>> gauss=Gaussian1DKernel(stddev=2)

gauss is not an array, but a kernel object. The underlying array can be retrieved with:

Astropy’s convolution methods can be used to replace bad data with values
interpolated from their neighbors. Kernel-based interpolation is useful for
handling images with a few bad pixels or for interpolating sparsely sampled
images.

Some contexts in which you might want to use kernel-based interpolation include:

images with saturated pixels. Generally, these are the highest-intensity
regions in the imaged area, and the interpolated values are not reliable,
but this can be useful for display purposes

images with flagged pixels, e.g., with a few small regions affected by cosmic
rays or other spurious signals that require those pixels to be flagged out.
If the affected region is small enough, the resulting interpolation will have
a small effect on source statistics and may allow for robust source finding
algorithms to be run on the resulting data

sparsely sampled images such as those constructed with single-pixel
detectors. Such images will only have a few discrete points sampled across
the imaged area, but an approximation of the extended sky emission can still
be constructed.

The script below shows an example of kernel interpolation to fill in
flagged-out pixels:

importnumpyasnpimportmatplotlib.pyplotaspltfromastropy.ioimportfitsfromastropy.utils.dataimportget_pkg_data_filenamefromastropy.convolutionimportGaussian2DKernel,interpolate_replace_nans# Load the data from data.astropy.orgfilename=get_pkg_data_filename('galactic_center/gc_msx_e.fits')hdu=fits.open(filename)[0]img=hdu.data[50:90,60:100]*1e5# This example is intended to demonstrate how astropy.convolve and# scipy.convolve handle missing data, so we start by setting the brightest# pixels to NaN to simulate a "saturated" data setimg[img>2e1]=np.nan# We smooth with a Gaussian kernel with x_stddev=1 (and y_stddev=1)# It is a 9x9 arraykernel=Gaussian2DKernel(x_stddev=1)# create a "fixed" image with NaNs replaced by interpolated valuesfixed_image=interpolate_replace_nans(img,kernel)# Now we do a bunch of plots. In the first two plots, the originally masked# values are marked with red X'splt.figure(1,figsize=(12,6)).clf()plt.close(2)# close the second plot from aboveax1=plt.subplot(1,2,1)im=ax1.imshow(img,vmin=-2.,vmax=2.e1,origin='lower',interpolation='nearest',cmap='viridis')y,x=np.where(np.isnan(img))ax1.set_autoscale_on(False)ax1.plot(x,y,'rx',markersize=4)ax1.set_title("Original")ax1.set_xticklabels([])ax1.set_yticklabels([])ax2=plt.subplot(1,2,2)im=ax2.imshow(fixed_image,vmin=-2.,vmax=2.e1,origin='lower',interpolation='nearest',cmap='viridis')ax2.set_title("Fixed")ax2.set_xticklabels([])ax2.set_yticklabels([])

This script shows the power of this technique for reconstructing images from
sparse sampling. Note that the image is not perfect - the pointlike sources
are sometimes missed - but the extended structure is very well recovered (by
eye).

importnumpyasnpimportmatplotlib.pyplotaspltfromastropy.ioimportfitsfromastropy.utils.dataimportget_pkg_data_filenamefromastropy.convolutionimportGaussian2DKernel,interpolate_replace_nans# Load the data from data.astropy.orgfilename=get_pkg_data_filename('galactic_center/gc_msx_e.fits')hdu=fits.open(filename)[0]img=hdu.data[50:90,60:100]*1e5indices=np.random.randint(low=0,high=img.size,size=300)sampled_data=img.flat[indices]# Build a new, sparsely sampled version of the original imagenew_img=np.tile(np.nan,img.shape)new_img.flat[indices]=sampled_data# We smooth with a Gaussian kernel with x_stddev=1 (and y_stddev=1)# It is a 9x9 arraykernel=Gaussian2DKernel(x_stddev=1)# create a "reconstructed" image with NaNs replaced by interpolated valuesreconstructed_image=interpolate_replace_nans(new_img,kernel)# Now we do a bunch of plots. In the first two plots, the originally masked# values are marked with red X'splt.figure(1,figsize=(12,6)).clf()ax1=plt.subplot(1,3,1)im=ax1.imshow(img,vmin=-2.,vmax=2.e1,origin='lower',interpolation='nearest',cmap='viridis')y,x=np.where(np.isnan(img))ax1.set_autoscale_on(False)ax1.set_title("Original")ax1.set_xticklabels([])ax1.set_yticklabels([])ax2=plt.subplot(1,3,2)im=ax2.imshow(new_img,vmin=-2.,vmax=2.e1,origin='lower',interpolation='nearest',cmap='viridis')ax2.set_title("Sparsely Sampled")ax2.set_xticklabels([])ax2.set_yticklabels([])ax2=plt.subplot(1,3,3)im=ax2.imshow(reconstructed_image,vmin=-2.,vmax=2.e1,origin='lower',interpolation='nearest',cmap='viridis')ax2.set_title("Reconstructed")ax2.set_xticklabels([])ax2.set_yticklabels([])

The behavior of astropy’s direct convolution
(convolve()) changed in version 2.0. Generally, the
old version is undesirable. However, to recover the behavior of the old
(astropy version <2.0) direct convolution function, you can interpolate and
then convolve, e.g.:

Note that the default behavior of both convolve and
convolve_fft is to perform normalized convolution and
interpolate NaNs during that process. The example given in this note, and what
was previously done only in direct convolution in old versions of astropy, does
a two-step process: first, it replaces the NaNs with their interpolated values
while leaving all non-NaN values unchanged, then it convolves the resulting
image with the specified kernel.