Removing Cosmic Ray Hits

This example illustrates how to remove cosmic ray hits from a LASCO C2 FITS file. using astroscrappy.detect_cosmics. Astroscrappy is a separate Python package and can be installed separately using pip or conda.

import astroscrappy
import matplotlib.pyplot as plt
from astropy.io import fits
from sunpy.map import Map
from sunpy.net import Fido
from sunpy.net import attrs as a

For more details on how to download and plot LASCO FITS file see SunPy’s example Downloading and plotting LASCO C3 data. To make this example work you need to have SunPy with all the “net” dependencies installed.

In order to download the required FITS file, we use Fido, SunPy’s downloader client. We need to define two search variables: a time range and the instrument.

time_range = a.Time("2000/11/09 00:06", "2000/11/09 00:07")
instrument = a.Instrument("LASCO")
detector = a.Detector("C2")
result = Fido.search(time_range, instrument)

downloaded_files = Fido.fetch(result[0])
data, header = fits.open(downloaded_files[0])[0].data, fits.open(downloaded_files[0])[0].header

# Add the missing meta information to the header
header["CUNIT1"] = "arcsec"
header["CUNIT2"] = "arcsec"

Files Downloaded:   0%|          | 0/1 [00:00<?, ?file/s]

22080376.fts:   0%|          | 0.00/2.11M [00:00<?, ?B/s]

22080376.fts:  18%|█▊        | 381k/2.11M [00:00<00:00, 3.75MB/s]

22080376.fts:  79%|███████▉  | 1.67M/2.11M [00:00<00:00, 6.25MB/s]


Files Downloaded: 100%|██████████| 1/1 [00:00<00:00,  1.23file/s]
Files Downloaded: 100%|██████████| 1/1 [00:00<00:00,  1.22file/s]

With this fix we can load it into a map and plot the results.

lasco_map = Map(data, header)
fig1 = plt.figure()
lasco_map.plot()

INFO: Missing metadata for solar radius: assuming the standard radius of the photosphere. [sunpy.map.mapbase]
/home/docs/checkouts/readthedocs.org/user_builds/sunkit-image/conda/latest/lib/python3.12/site-packages/sunpy/map/mapbase.py:632: SunpyMetadataWarning: Missing metadata for observer: assuming Earth-based observer.
For frame 'heliographic_stonyhurst' the following metadata is missing: hgln_obs,hglt_obs,dsun_obs
For frame 'heliographic_carrington' the following metadata is missing: crlt_obs,crln_obs,dsun_obs

  obs_coord = self.observer_coordinate

<matplotlib.image.AxesImage object at 0x7fdf23008590>

Now we will call the astroscrappy.detect_cosmics to remove the cosmic ray hits. This algorithm can perform well with both high and low noise levels in the original data. The function takes a ndarray as input so we only pass the map data. This particular image has lots of high intensity cosmic ray hits which cannot be effectively removed by using the default set of parameters. So we reduce sigclip, the Laplacian to noise ratio from 4.5 to 2 to mark more hits. We also reduce objlim, the contrast between the Laplacian image and the fine structured image to clean the high intensity bright cosmic ray hits. We also modify the readnoise parameter to obtain better results.

mask, clean_data = astroscrappy.detect_cosmics(lasco_map.data, sigclip=2, objlim=2, readnoise=4, verbose=True)

Starting 4 L.A.Cosmic iterations
Iteration 1:
169952 cosmic pixels this iteration
Iteration 2:
42674 cosmic pixels this iteration
Iteration 3:
13946 cosmic pixels this iteration
Iteration 4:
6735 cosmic pixels this iteration

This returns two variables - mask is a boolean array depicting whether there is a cosmic ray hit at that pixel, clean_data is the cleaned image after removing those hits.

Now we can now plot the cleaned image.

clean_map1 = Map(clean_data, lasco_map.meta)
fig2 = plt.figure()
clean_map1.plot()

plt.show()

INFO: Missing metadata for solar radius: assuming the standard radius of the photosphere. [sunpy.map.mapbase]
/home/docs/checkouts/readthedocs.org/user_builds/sunkit-image/conda/latest/lib/python3.12/site-packages/sunpy/map/mapbase.py:632: SunpyMetadataWarning: Missing metadata for observer: assuming Earth-based observer.
For frame 'heliographic_stonyhurst' the following metadata is missing: hgln_obs,hglt_obs,dsun_obs
For frame 'heliographic_carrington' the following metadata is missing: crlt_obs,crln_obs,dsun_obs

  obs_coord = self.observer_coordinate