Satellite trail detection in ACS/WFC data using acstools.findsat_mrt#
This notebook provides examples of how to find and create masks for satellite trails in ACS/WFC imaging data using acstools.findsat_mrt, which is based on the method described in ACS ISR 2022-08. Many of the tools presented here should be applicable to any imaging data.
Table of Contents:#
Introduction
Imports, Setup, and Data
Example 1: Step-by-step guide to find trails in an FLC image
Example 2: Quick run on an FLC image
Example 3: Find trails in an FLC image using the WFC wrapper
Example 4: Step-by-step guide to find trails in a DRC image
Example 5: Find trails in a DRC image using the WFC wrapper
Example 6: Create a new kernel for detection
About this Notebook#
Author: David V. Stark, ACS Instrument Team, Space Telescope Science Institute
First Published On: 5/13/2023
Updated On: 5/15/2023
Introduction#
Despite being in orbit, HST imaging data still suffers from contamination by artificial satellites that can compromise science data unless they are identified and masked. This notebook presents examples of how to identify satellite trails in ACS/WFC data. The routine is also effective at identifying other linear features duch as diffraction spikes and glint (see Section 4.5 of the ACS DHB for further discussion on these artifacts).
A full description of the algorithm is provided in ACS ISR 2022-08. To briefly summarize, the Median Radon Transform (MRT) is calculated for an input image and used to identify linear signals in the data. The MRT is similar to the standard Radon Transform except that it calculates the median, rather than the sum, of data along all possible paths through an image. This modification makes the algorithm more robust against false signals from localized sources (e.g., stars, galaxies) but still very sensitive to persistent linear features, even well-below the background noise level.
Additional post-processing is done to filter out spurious detections, primarily eliminating them based on trail S/N, width, and persistence across the image. These parameters, especially the maximum allowed trail width, are tuned for ACS/WFC data binned 2x2 and may be different for images from other instruments. Once the final set of trails is identified and characterized, a mask can be created. The routine provides numerous ways of visualizing the results, as will be demonstrated below.
The following examples illustrate how to use acstools.findsat_mrt to identify satellite trails and then create masks for them. Examples 1 and 4 go through the analysis step by step, including how to preprocess data and run individual routines inside findsat_mrt. Examples 2, 3, and 5 demonstrate how to automate many of these steps. Our demonstrations stop at the creation of the masks. We leave it to the user to decide the best way to apply the masks to their own analysis.
Imports, setup, and data#
It is recommended that you use the latest stenv python environment when using this notebook. In particular, you must use acstools v3.6.0 or greater in order to run this notebook. You can check you version with
conda list acstools
and update if necessary with
conda update acstools
Set your working directory and import the needed packages with the following
# Import modules and setup
import matplotlib.pyplot as plt
import numpy as np
from astroquery.mast import Observations
from astropy.io import fits
from astropy.nddata import bitmask, block_reduce, block_replicate
from acstools.findsat_mrt import TrailFinder, WfcWrapper
import os
from acstools.utils_findsat_mrt import create_mrt_line_kernel
import shutil
# Check your own working directory
print('Current working directory is {}'.format(os.getcwd()))
# Define working directory if needed
# os.chdir('Insert your working directory here')
# These are optional configurations
%matplotlib inline
plt.rcParams["figure.figsize"] = (8, 6)
plt.rcParams['font.serif'] = "Georgia"
plt.rcParams['font.family'] = "serif"
Download the example data needed and place it in the working directory that you defined above. Examples 1-3 use jc8m32j5q_flc.fits, while examples 4-5 use hst_13498_32_acs_wfc_f606w_jc8m32j5_drc.fits.
# Download data files
obs_table = Observations.query_criteria(proposal_id=13498, obs_id='JC8M32010')
dl_table = Observations.download_products(obs_table['obsid'],
dataURI=['mast:HST/product/hst_13498_32_acs_wfc_f606w_jc8m32j5_drc.fits',
'mast:HST/product/jc8m32j5q_flc.fits'])
for row in dl_table:
oldfname = row['Local Path']
newfname = os.path.basename(oldfname)
os.rename(oldfname, newfname)
shutil.rmtree('mastDownload')
Example 1: Finding trails in an FLC image#
FLC images are individual exposures processed by the CALACS pipeline. The data contain two chips, but we only analyze one here.
We start by reading in an image and doing some pre-processing to remove bad pixels, subtract a median background, and make the image a bit smaller (to speed up the calculation of the MRT).
# Read in the image files and header information
image_file = 'jc8m32j5q_flc.fits'
ext = 4 # ACS image data are in extensions 1 or 4, we'll just use 4 for now (chip 1)
with fits.open(image_file) as h:
image = h[ext].data # image data
dq = h[ext+2].data # data quality bitmasks
header = h[0].header # primary header
image_header = h[1].header # image header
Below, we make a mask for bad pixels. We’re ignoring cosmic rays here because routines to make them often partially (but not fully) mask trails. By default, any masked pixels are set to NaN.
mask = bitmask.bitfield_to_boolean_mask(dq, ignore_flags=[4096, 8192, 16384])
image[mask] = np.nan
Below we subtract Subtract the background from the image. Here we just do a simple median.
image = image - np.nanmedian(image)
The MRT is computationally demanding and WFC images are big. To help things a bit, let’s rebin the images.
binsize = 2 # adjust this as needed
image_rebin = block_reduce(image, binsize, func=np.nansum)
We now set up TrailFinder. Many of the parameters in the call below are optional (and set to their current values by default) but we show them to illustrate the setup. Of note is that I’m explicitly defining the image header keys to save. These can be useful later when analyzing trail population properties. The keywords being saved here were chosen to ensure we know the original exposure ippsoot and which chip was analyzed. Additional keywords are saved that store information about the orientation of the telescope when the image was taken. In principle, the user can save any header keywords they like. We have also set plot=False in this example, so we can demonstrate how to manually create plots. Setting plot=True will automatically generate plots after specific processes are finished. Be aware that not all possible keyword parameters are defined below. See the documentation for complete information.
# Now we can set up TrailFinder
s = TrailFinder(image=image_rebin,
header=header,
image_header=image_header,
save_image_header_keys=['ROOTNAME', 'CCDCHIP', 'CRPIX1', 'CRPIX2', 'CRVAL1', 'CRVAL2',
'ORIENTAT', 'RA_APER', 'DEC_APER', 'PA_APER'],
processes=8,
plot=False,
threshold=5,
max_width=75,
check_persistence=True,
min_persistence=0.5,
output_root='example1')
Before we actually run anything, let’s plot the image we are analyzing. You should see two satellite trails in this example.
s.plot_image()
If you’re having trouble seeing the trails, you can adjust the scale keyword (the min and max values to show given as multiples of the image standard deviation)
s.plot_image(scale=[-1, 1])
Next we run the Median Radon Transform. This step can take some time depending on the image size and number of processes being used. This tutorial assumes you can run 8 processes at the same time, but adjust as needed. If you’re not sure how many processes you can run, you can see how many CPU cores are available and adjust based on that.
os.cpu_count()
s.processes = 8 # adjust this if necessary
s.run_mrt()
Now we will plot the MRT. You may be able to spot the signals from the satellite trails as two somewhat wide point-like sources.
s.plot_mrt()
Note that the x axis in in pixels, not degrees or radians. The theta array ranges from 0 to 180 with a spacing of 0.5 degrees, hence 360 pixels.
We next run the source finder on the MRT. You can create your own detection kernels, or use the defaults provided (see Example 6 for how to create detection kernels). Depending on the settings, this can pick up a lot more than the actual trails we’re interested in. There are additional steps we’ll take later to filter these false detections out. The ones we found and their location on the MRT are shown below.
The threshold in this case refers to the signal-to-noise ratio of a feature found in the MRT. The default is 5.
s.threshold = 5 # detection threshold
s.find_mrt_sources() # finds the sources
s.plot_mrt(show_sources=True) # overplots the sources on top of the MRT
We filter the sources further based on a reassessment of their S/N, width, and persistence. The default parameters (namely width) have been chosen for ACS data binned by 2 pixels in each direction. It’s possible different defaults will be better for different imaging data.
# Parameters that affect how the filtering works
s.threshold = 5
s.max_width = 75
s.check_persistence = True
s.min_persistence = 0.5
# now filter
s.filter_sources()
# note: some extra columns have been added to the source list
s.source_list
Several columns have been added to the source list that characterize the observed streak. Also, the status array has values of 0, 1, and 2 now (it just had 0 before). Those with status=2 are sources that passed all filtering stages (checks for SNR and width, then persistence). Those with status=1 are sources that passed the first filtering stage (checks for SNR and width), but not the second (persistence check). And status=0 are sources that did not pass the filtering steps.
The plot_mrt command will overplot the different statuses
s.plot_mrt(show_sources=True)
Now we can make the mask itself. By default it only uses sources in the MRT with status=2. We make two types of masks, one a simple boolean mask, and one a segementation mask where pixels corresponding to each streak are assigned the ID number. We create these below.
# make the mask
s.mask_include_status = [2]
s.make_mask()
# plot the mask and segmentation map
s.plot_mask()
s.plot_segment()
We can also overlay the mask on top of the image to make sure it makes sense.
s.plot_image(overlay_mask=True)
We can save the results now. You have the optional of saving the catalog, mask, MRT, and a diagnostic image that shows the results. In this example we’ll just save everything.
# define what to save
s.save_mask = True
s.save_mrt = True
s.save_catalog = True
s.save_diagnostic = True
s.save_output()
Keep in mind that the mask we have created is applicable to the rebinned image. To convert it into a mask that can be applied to the original unbinned image, we need to resample it using the block_replicate function. The rescaled mask is plotted below. Note the difference in image size, but the mask pattern remains the same.
full_mask = block_replicate(s.mask, binsize, conserve_sum=False)
fig, ax = plt.subplots()
ax.imshow(full_mask, origin='lower')
#
Example 2: Quick run of TrailFinder on an flc image#
Example 1 thoroughly demonstrated the steps to read in an FLC file, pre-process it, and identify trails. This example demonstrates how one can run many of the steps simultaneously once a file is read in an all parameters set.
First, we read in and preprocess the data file exactly as before.
# Read in the image files and header information
image_file = 'jc8m32j5q_flc.fits'
ext = 4 # ACS image data are in extensions 1 or 4, we'll just use 1 for now
with fits.open(image_file) as h:
image = h[ext].data # image data
dq = h[ext+2].data # data quality bitmasks
header = h[0].header # primary header
image_header = h[1].header # image header
# make a mask for bad pixels.
mask = bitmask.bitfield_to_boolean_mask(dq, ignore_flags=[4096, 8192, 16384])
image[mask] = np.nan
# Subtract the background from the image.
image = image - np.nanmedian(image)
# Rebin the image to speed up calculation
image_rebin = block_reduce(image, 2, func=np.nansum)
print(image)
And initialize trail finder as before
s2 = TrailFinder(image=image_rebin,
header=header,
image_header=image_header,
save_image_header_keys=['ROOTNAME', 'CCDCHIP', 'CRPIX1', 'CRPIX2', 'CRVAL1', 'CRVAL2',
'ORIENTAT', 'RA_APER', 'DEC_APER', 'PA_APER'],
processes=8,
plot=False,
threshold=5,
max_width=75,
check_persistence=True,
min_persistence=0.5,
output_root='example2')
If you’re feeling ok about the setup, run all the subsequent steps together with the run_all command (this calculates the MRT, finds MRT sources, filters the sources, and saves the output)
s2.run_all()
If we plot the mask, it should look identical to the one in the previous example.
s2.plot_mask()
#
Example 3: find trails in an FLC image using the WFC wrapper#
The approaches shown in examples 1 and 2 can be useful for imaging data from any telescope, not just ACS/WFC data. However, for ACS/WFC data, we provide a convenience wrapper that performs even more of the steps all together, including reading the image and pre-processing it.
The WfcWrapper class has the same properties as the TrailFinder class, but with a few additional keywords. It also contains the additional routines that read the image, rebin, mask, and subtract the background. By default, these will be run automatically when WfcWrapper is initialized, although this can be turned off. In most cases, you probably will only need to adjust the binsize keyword. The specific value of binsize is up to the user. Larger values speed up the MRT calculation, but keep in mind that the parameters to filter out spurious trails (e.g., max_width) are tuned to WFC data binned 2x2. A user may want to start with a larger value for binsize and reduce it once they get a sense for the computation time.
w = WfcWrapper('jc8m32j5q_flc.fits', binsize=2, extension=4, processes=8, output_root='example3')
We can plot the image to see that it looks like the one from the last example after preprocessing.
w.plot_image()
From here, everything is the same as the last example:
w.run_mrt()
w.find_mrt_sources()
w.filter_sources()
Below is the resulting MRT and sources
w.plot_mrt(show_sources=True)
Lastly, we generate the mask
w.make_mask()
w.plot_mask()
If you’re really feeling very confident, you can run everything in a single line by setting execute=True.
w = WfcWrapper('jc8m32j5q_flc.fits', binsize=2, extension=4, output_root='example3', processes=8,
execute=True)
We’ll plot the image and mask together to check that everything looks ok
w.plot_image(overlay_mask=True)
Example 4: Finding trails in a DRC image#
Applying TrailFinder to a DRC image (that shows both chips together) can boost sensitivity by increasing the number of pixels over which we search for trails. The DRC files also remove the distortion in the original FLC files (though this does not appear to create signficant curvature to most trails).
Here, we demonstrate the steps that go into preparing a DRC image to be analyzed. The subsequent example will illustrate how to do all of this in a single line.
There are no DQ arrays for the DRC files, so we ignore the pre-processing steps that incorporated those.
# Read in the image files and header information
image_file = 'hst_13498_32_acs_wfc_f606w_jc8m32j5_drc.fits'
ext = 1
with fits.open(image_file) as h:
image = h[ext].data # image data
wht = h[ext+1].data
image = image*wht # wht is effective exposure time, so this turns it into counts
header = h[0].header # primary header
image_header = h[1].header # image header
# Flag anything with wht == 0 as bad
image[wht == 0] = np.nan
# Subtract the background from the image.
median = np.nanmedian(image)
image = image - np.nanmedian(image)
# Let's rebin the images
binsize = 2
image_rebin = block_reduce(image, binsize, func=np.nansum)
Setting up TrailFinder is essentially the same as earlier examples at this point. We’ll use the default settings. In fact, about all the steps from here on out are the same.
s4 = TrailFinder(image=image_rebin, processes=8, output_root='example4')
We can do a quick plot of our image to make sure things look ok
s4.plot_image()
Now run the MRT calculation and plot the results
s4.run_mrt()
s4.plot_mrt(scale=[-1, 5]) # adjusted scale manually due to varying background in image
This example has a clear gradient in the background due to the cluster. This causes some large scale variation in the RT, but you can see the “point source” signals from the satellite trails around x,y = (90,700) and x,y = (300,700). This is a case where we may have wanted to explore some different background subtraction methods, but we’ll proceed with the simpler approach here. Now we’ll try to pull the sources out.
s4.find_mrt_sources()
And below we plot the MRT with the sources overlaid
s4.plot_mrt(show_sources=True)
It’s clearly shredding those large-scale features quite a bit, but we’ll try to filter these out.
s4.filter_sources()
Let’s re-plot the MRT with sources to see what made it through
s4.plot_mrt(show_sources=True)
That seems to have worked! Let’s make the map to confirm
s4.make_mask()
s4.plot_mask()
s4.plot_segment()
Let’s make a version plotting the mask on top of the original image
s4.plot_image(overlay_mask=True)
Example 5: Finding trails in a DRC image using the WFC Wrapper#
All of the setup from the last example can be streamlined using the WfcWrapper class.
w2 = WfcWrapper('hst_13498_32_acs_wfc_f606w_jc8m32j5_drc.fits', binsize=2, extension=1, processes=8,
output_root='example5')
Run full pipeline now
w2.run_all()
Let’s plot the final mask to ensure it looks the same as the earlier examples.
w2.plot_mask()
And there you go!
Example 6: Create a new kernel for trail detection#
We include a function called create_mrt_line_kernel that can be used to generate kernels for detecting trails of s specified size in the MRT. Note that kernels with widths of 1, 3, 7, and 15 pixels (convolved with a simple Gaussian HST/ACS psf model) are included already, but perhaps you want to generate a kernel with a new width, or convolved with a different PSF.
Let’s generate a kernel for a trail that has an inherent width of 5 pixels and is convolved with a Gaussian PSF with sigma=3.
out = create_mrt_line_kernel(5, 3, processes=8, plot=True)
The first plot show the model streak. The second plot shows the resulting MRT. The kernsl is created by taking a cutout around the signal in the MRT. The third double-plot shows 1D slices of the signal in the MRT, with orange lines showing the location of the maximum values. These serve as first guesses of the center, after which the center is redetermined using a Guassian fit and the cutout extracted with the kernel perfectly centered. The 4th plot above shows the final kernel.
The kernel can be saved by defining the outfile keyword in create_mrt_line_kernel. By adding this file path into the kernels keyword in TrailFinder or WfcWrapper, it will be used for source detection when running find_mrt_sources.
For more help:#
More details may be found on the ACS website and in the ACS Instrument and Data Handbooks.
Please visit the HST Help Desk. Through the help desk portal, you can explore the HST Knowledge Base and request additional help from experts.
Top of Page
