diff --git a/.github/workflows/build-docker.yml b/.github/workflows/build-docker.yml new file mode 100644 index 00000000..b9285031 --- /dev/null +++ b/.github/workflows/build-docker.yml @@ -0,0 +1,58 @@ +name: Docker Image + +concurrency: ${{ github.workflow }}-${{ github.ref }} + +on: + push: + branches: + - "*" + tags: + - "*" + pull_request: + branches: + - master + +jobs: + build-push: + + runs-on: ubuntu-latest + + permissions: + contents: write + packages: write + + env: + REGISTRY: ghcr.io + IMAGE_NAME: ${{ github.repository }} + + steps: + - name: "Checkout" + uses: actions/checkout@v3 + + - name: Set up Docker Buildx + id: buildx + uses: docker/setup-buildx-action@v2 + with: + install: true + + - name: Login to Container Registry + uses: docker/login-action@v2 + with: + registry: ${{ env.REGISTRY }} + username: ${{ github.actor }} + password: ${{ secrets.GITHUB_TOKEN }} + + - name: Extract Docker metadata + id: meta + uses: docker/metadata-action@v4 + with: + images: ${{ env.REGISTRY }}/${{ env.IMAGE_NAME }} + + - name: Build and also push Dockerimage + id: build-and-push + uses: docker/build-push-action@v3 + with: + context: . + push: true + tags: ${{ steps.meta.outputs.tags }} + labels: ${{ steps.meta.outputs.labels }} diff --git a/.github/workflows/e2e.yaml b/.github/workflows/e2e.yaml index e71aa7a2..d43efe84 100644 --- a/.github/workflows/e2e.yaml +++ b/.github/workflows/e2e.yaml @@ -45,7 +45,7 @@ jobs: - name: Test Master Bias Creation run: | - kubectl exec banzai-e2e-test -c banzai-listener -- pytest -o log_cli=true -s --pyargs banzai.tests --durations=0 --junitxml=/archive/engineering/pytest-master-bias.xml -m master_bias + kubectl exec banzai-e2e-test -c banzai-listener -- pytest -o log_cli=true -s --pyargs banzai --durations=0 --junitxml=/archive/engineering/pytest-master-bias.xml -m master_bias - name: Cleanup run: | diff --git a/CHANGES.md b/CHANGES.md index c84f005d..7be773cb 100644 --- a/CHANGES.md +++ b/CHANGES.md @@ -1,3 +1,7 @@ +1.15.0 (2024-02-14) +------------------- +- Migrated photometry extraction to be done by astropy's photutils instead of SEP. + 1.14.1 (2024-01-29) ------------------- - Fix for null instrument/instrument size to allow processing to continue @@ -7,7 +11,7 @@ - Split worker queues for small and large images - Update e2e tests to run on github actions -1.13.1 (2023-10-21) +1.13.1 (2023-10-31) ------------------- - bugfix to the docker build diff --git a/MANIFEST.in b/MANIFEST.in index 390ebf70..37d47f2f 100644 --- a/MANIFEST.in +++ b/MANIFEST.in @@ -13,3 +13,4 @@ prune docs/_build prune docs/api include pyproject.toml +include pytest.ini diff --git a/banzai/bpm.py b/banzai/bpm.py index 473f3354..59d2ac14 100644 --- a/banzai/bpm.py +++ b/banzai/bpm.py @@ -39,6 +39,6 @@ def calibration_type(self): class SaturatedPixelFlagger(Stage): def do_stage(self, image): for image_extension in image.ccd_hdus: - image_extension.mask[image_extension.data > image_extension.saturate] |= 2 + image_extension.mask[image_extension.data >= image_extension.saturate] |= 2 return image diff --git a/banzai/photometry.py b/banzai/photometry.py index afb4eaec..e48fce32 100755 --- a/banzai/photometry.py +++ b/banzai/photometry.py @@ -2,7 +2,6 @@ import numpy as np from astropy.table import Table -import sep from requests import HTTPError from banzai.utils import stats, array_utils @@ -11,25 +10,65 @@ from banzai.data import DataTable from banzai import logs +from photutils.background import Background2D from skimage import measure +from photutils.segmentation import make_2dgaussian_kernel, detect_sources, deblend_sources, SourceCatalog +from astropy.convolution import convolve +from astropy.convolution.kernels import CustomKernel + logger = logs.get_logger() -sep.set_sub_object_limit(int(1e6)) def radius_of_contour(contour, source): x = contour[:, 1] y = contour[:, 0] - x_center = (source['xmax'] - source['xmin'] + 1) / 2.0 - 0.5 - y_center = (source['ymax'] - source['ymin'] + 1) / 2.0 - 0.5 - - return np.percentile(np.sqrt((x - x_center)**2.0 + (y - y_center)** 2.0), 90) + x_center = (source.bbox_xmax - source.bbox_xmin + 1) / 2.0 - 0.5 + y_center = (source.bbox_ymax - source.bbox_ymin + 1) / 2.0 - 0.5 + + return np.percentile(np.sqrt((x - x_center)**2.0 + (y - y_center) ** 2.0), 90) + + +def flag_sources(sources, source_labels, segmentation_map, mask, flag, mask_value): + affected_sources = np.unique(segmentation_map.data[mask == mask_value]) + sources['flag'][np.in1d(source_labels, affected_sources)] |= flag + + +def flag_deblended(sources, catalog, segmentation_map, deblended_seg_map, flag_value=2): + # By default deblending appends labels instead of reassigning them so we can just use the + # extras in the deblended map + deblended_sources = np.unique(deblended_seg_map.data[deblended_seg_map > np.max(segmentation_map)]) + # Get the sources that were originally blended + original_blends = np.unique(segmentation_map.data[deblended_seg_map > np.max(segmentation_map)]) + deblended_sources = np.hstack([deblended_sources, original_blends]) + sources['flag'][np.in1d(catalog.labels, deblended_sources)] |= flag_value + + +def flag_edge_sources(image, sources, flag_value=8): + ny, nx = image.shape + # Check 4 points on the kron aperture, one on each side of the major and minor axis + minor_xmin = sources['x'] - sources['b'] * sources['kronrad'] * np.sin(np.deg2rad(sources['theta'])) + minor_xmax = sources['x'] + sources['b'] * sources['kronrad'] * np.sin(np.deg2rad(sources['theta'])) + minor_ymin = sources['y'] - sources['b'] * sources['kronrad'] * np.cos(np.deg2rad(sources['theta'])) + minor_ymax = sources['y'] + sources['b'] * sources['kronrad'] * np.cos(np.deg2rad(sources['theta'])) + major_ymin = sources['y'] - sources['a'] * sources['kronrad'] * np.sin(np.deg2rad(sources['theta'])) + major_ymax = sources['y'] + sources['a'] * sources['kronrad'] * np.sin(np.deg2rad(sources['theta'])) + major_xmin = sources['x'] - sources['a'] * sources['kronrad'] * np.cos(np.deg2rad(sources['theta'])) + major_xmax = sources['x'] + sources['a'] * sources['kronrad'] * np.cos(np.deg2rad(sources['theta'])) + + # Note we are already 1 indexed here + sources_off = np.logical_or(minor_xmin < 1, major_xmin < 1) + sources_off = np.logical_or(sources_off, minor_ymin < 1) + sources_off = np.logical_or(sources_off, major_ymin < 1) + sources_off = np.logical_or(sources_off, minor_xmax > nx) + sources_off = np.logical_or(sources_off, major_xmax > nx) + sources_off = np.logical_or(sources_off, minor_ymax > ny) + sources_off = np.logical_or(sources_off, major_ymax > ny) + sources[sources_off]['flag'] |= flag_value class SourceDetector(Stage): - # Note that threshold is number of sigma, not an absolute number because we provide the error - # array to SEP. - threshold = 10.0 + threshold = 2.5 min_area = 9 def __init__(self, runtime_context): @@ -37,81 +76,77 @@ def __init__(self, runtime_context): def do_stage(self, image): try: - # Increase the internal buffer size in sep. This is most necessary for crowded fields. - ny, nx = image.shape - sep.set_extract_pixstack(int(nx * ny - 1)) - data = image.data.copy() error = image.uncertainty - mask = image.mask > 0 - - # Fits can be backwards byte order, so fix that if need be and subtract - # the background - try: - bkg = sep.Background(data, mask=mask, bw=32, bh=32, fw=3, fh=3) - except ValueError: - data = data.byteswap(True).newbyteorder() - bkg = sep.Background(data, mask=mask, bw=32, bh=32, fw=3, fh=3) - bkg.subfrom(data) + # From what I can piece together, the background estimator makes a low resolution mesh set by box size + # (32, 32) here and then applies a filter to the low resolution image. The filter size is 3x3 here. + # The defaults we use here are a mesh creator is from source extractor which is a mode estimator. + # The default filter that works on the mesh image is a median filter. + bkg = Background2D(data, (32, 32), filter_size=(3, 3)) + data -= bkg.background + + # Convolve the image with a 2D Guassian, but with the normalization SEP uses as + # that is correct. + # The default kernel used by Source Extractor is [[1,2,1], [2,4,2], [1,2,1]] + # The kernel corresponds to fwhm = 1.9 which we adopt here + kernel = make_2dgaussian_kernel(1.9, size=3) + convolved_data = convolve(data / (error * error), kernel) + + # We include the correct match filter normalization here that is not included in + # vanilla source extractor + kernel_squared = CustomKernel(kernel.array * kernel.array) + normalization = np.sqrt(convolve(1 / (error * error), kernel_squared)) + convolved_data /= normalization + logger.info('Running image segmentation', image=image) # Do an initial source detection - try: - sources = sep.extract(data, self.threshold, mask=mask, minarea=self.min_area, - err=error, deblend_cont=0.005) - except Exception: - logger.error(logs.format_exception(), image=image) - return image - - # Convert the detections into a table - sources = Table(sources) - - # We remove anything with a detection flag >= 8 - # This includes memory overflows and objects that are too close the edge - sources = sources[sources['flag'] < 8] + segmentation_map = detect_sources(convolved_data, self.threshold, npixels=self.min_area) + + logger.info('Deblending sources', image=image) + # Note that nlevels here is DEBLEND_NTHRESH in source extractor which is 32 by default + deblended_seg_map = deblend_sources(convolved_data, segmentation_map, + npixels=self.min_area, nlevels=32, + contrast=0.005, progress_bar=False, + nproc=1, mode='sinh') + logger.info('Finished deblending. Saving catalog', image=image) + # Convert the segmentation map to a source catalog + catalog = SourceCatalog(data, deblended_seg_map, convolved_data=convolved_data, error=error, + background=bkg.background) + + sources = Table({'x': catalog.xcentroid + 1.0, 'y': catalog.ycentroid + 1.0, + 'xwin': catalog.xcentroid_win + 1.0, 'ywin': catalog.ycentroid_win + 1.0, + 'xpeak': catalog.maxval_xindex + 1, 'ypeak': catalog.maxval_yindex + 1, + 'peak': catalog.max_value, + 'a': catalog.semimajor_sigma.value, 'b': catalog.semiminor_sigma.value, + 'theta': catalog.orientation.to('deg').value, 'ellipticity': catalog.ellipticity.value, + 'kronrad': catalog.kron_radius.value, + 'flux': catalog.kron_flux, 'fluxerr': catalog.kron_fluxerr, + 'x2': catalog.covar_sigx2.value, 'y2': catalog.covar_sigy2.value, + 'xy': catalog.covar_sigxy.value, + 'background': catalog.background_mean}) + + for r in range(1, 7): + radius_arcsec = r / image.pixel_scale + sources[f'fluxaper{r}'], sources[f'fluxerr{r}'] = catalog.circular_photometry(radius_arcsec) + + for r in [0.25, 0.5, 0.75]: + sources['fluxrad' + f'{r:.2f}'.lstrip("0.")] = catalog.fluxfrac_radius(r) + + sources['flag'] = 0 + + # Flag = 1 for sources with bad pixels + flag_sources(sources, catalog.labels, deblended_seg_map, image.mask, flag=1, mask_value=1) + # Flag = 2 for sources that are deblended + flag_deblended(sources, catalog, segmentation_map, deblended_seg_map, flag_value=2) + # Flag = 4 for sources that have saturated pixels + flag_sources(sources, catalog.labels, deblended_seg_map, image.mask, flag=4, mask_value=2) + # Flag = 8 if kron aperture falls off the image + flag_edge_sources(image, sources, flag_value=8) + # Flag = 16 if source has cosmic ray pixels + flag_sources(sources, catalog.labels, deblended_seg_map, image.mask, flag=16, mask_value=8) sources = array_utils.prune_nans_from_table(sources) - # Calculate the ellipticity - sources['ellipticity'] = 1.0 - (sources['b'] / sources['a']) - - # Fix any value of theta that are invalid due to floating point rounding - # -pi / 2 < theta < pi / 2 - sources['theta'][sources['theta'] > (np.pi / 2.0)] -= np.pi - sources['theta'][sources['theta'] < (-np.pi / 2.0)] += np.pi - - # Calculate the kron radius - kronrad, krflag = sep.kron_radius(data, sources['x'], sources['y'], - sources['a'], sources['b'], - sources['theta'], 6.0) - sources['flag'] |= krflag - sources['kronrad'] = kronrad - - # Calcuate the equivilent of flux_auto - flux, fluxerr, flag = sep.sum_ellipse(data, sources['x'], sources['y'], - sources['a'], sources['b'], - np.pi / 2.0, 2.5 * kronrad, - subpix=1, err=error) - sources['flux'] = flux - sources['fluxerr'] = fluxerr - sources['flag'] |= flag - - # Do circular aperture photometry for diameters of 1" to 6" - for diameter in [1, 2, 3, 4, 5, 6]: - flux, fluxerr, flag = sep.sum_circle(data, sources['x'], sources['y'], - diameter / 2.0 / image.pixel_scale, gain=1.0, err=error) - sources['fluxaper{0}'.format(diameter)] = flux - sources['fluxerr{0}'.format(diameter)] = fluxerr - sources['flag'] |= flag - - # Measure the flux profile - flux_radii, flag = sep.flux_radius(data, sources['x'], sources['y'], - 6.0 * sources['a'], [0.25, 0.5, 0.75], - normflux=sources['flux'], subpix=5) - sources['flag'] |= flag - sources['fluxrad25'] = flux_radii[:, 0] - sources['fluxrad50'] = flux_radii[:, 1] - sources['fluxrad75'] = flux_radii[:, 2] - # Cut individual bright pixels. Often cosmic rays sources = sources[sources['fluxrad50'] > 0.5] @@ -119,130 +154,91 @@ def do_stage(self, image): sources['fwhm'] = np.nan sources['fwtm'] = np.nan # Here we estimate contours - for source in sources: + for source, row in zip(sources, catalog): if source['flag'] == 0: for ratio, keyword in zip([0.5, 0.1], ['fwhm', 'fwtm']): - contours = measure.find_contours(data[source['ymin']: source['ymax'] + 1, - source['xmin']: source['xmax'] + 1], + contours = measure.find_contours(data[row.bbox_ymin: row.bbox_ymax + 1, + row.bbox_xmin: row.bbox_xmax + 1], ratio * source['peak']) if contours: # If there are multiple contours like a donut might have take the outer - contour_radii = [radius_of_contour(contour, source) for contour in contours] + contour_radii = [radius_of_contour(contour, row) for contour in contours] source[keyword] = 2.0 * np.nanmax(contour_radii) - # Calculate the windowed positions - sig = 2.0 / 2.35 * sources['fwhm'] - xwin, ywin, flag = sep.winpos(data, sources['x'], sources['y'], sig) - sources['flag'] |= flag - sources['xwin'] = xwin - sources['ywin'] = ywin - - # Calculate the average background at each source - bkgflux, fluxerr, flag = sep.sum_ellipse(bkg.back(), sources['x'], sources['y'], - sources['a'], sources['b'], np.pi / 2.0, - 2.5 * sources['kronrad'], subpix=1) - # masksum, fluxerr, flag = sep.sum_ellipse(mask, sources['x'], sources['y'], - # sources['a'], sources['b'], np.pi / 2.0, - # 2.5 * kronrad, subpix=1) - - background_area = (2.5 * sources['kronrad']) ** 2.0 * sources['a'] * sources['b'] * np.pi # - masksum - sources['background'] = bkgflux - sources['background'][background_area > 0] /= background_area[background_area > 0] - # Update the catalog to match fits convention instead of python array convention - sources['x'] += 1.0 - sources['y'] += 1.0 - - sources['xpeak'] += 1 - sources['ypeak'] += 1 - - sources['xwin'] += 1.0 - sources['ywin'] += 1.0 - - sources['theta'] = np.degrees(sources['theta']) - - catalog = sources['x', 'y', 'xwin', 'ywin', 'xpeak', 'ypeak', - 'flux', 'fluxerr', 'peak', 'fluxaper1', 'fluxerr1', - 'fluxaper2', 'fluxerr2', 'fluxaper3', 'fluxerr3', - 'fluxaper4', 'fluxerr4', 'fluxaper5', 'fluxerr5', - 'fluxaper6', 'fluxerr6', 'background', 'fwhm', 'fwtm', - 'a', 'b', 'theta', 'kronrad', 'ellipticity', - 'fluxrad25', 'fluxrad50', 'fluxrad75', - 'x2', 'y2', 'xy', 'flag'] - # Add the units and description to the catalogs - catalog['x'].unit = 'pixel' - catalog['x'].description = 'X coordinate of the object' - catalog['y'].unit = 'pixel' - catalog['y'].description = 'Y coordinate of the object' - catalog['xwin'].unit = 'pixel' - catalog['xwin'].description = 'Windowed X coordinate of the object' - catalog['ywin'].unit = 'pixel' - catalog['ywin'].description = 'Windowed Y coordinate of the object' - catalog['xpeak'].unit = 'pixel' - catalog['xpeak'].description = 'X coordinate of the peak' - catalog['ypeak'].unit = 'pixel' - catalog['ypeak'].description = 'Windowed Y coordinate of the peak' - catalog['flux'].unit = 'count' - catalog['flux'].description = 'Flux within a Kron-like elliptical aperture' - catalog['fluxerr'].unit = 'count' - catalog['fluxerr'].description = 'Error on the flux within Kron aperture' - catalog['peak'].unit = 'count' - catalog['peak'].description = 'Peak flux (flux at xpeak, ypeak)' + sources['x'].unit = 'pixel' + sources['x'].description = 'X coordinate of the object' + sources['y'].unit = 'pixel' + sources['y'].description = 'Y coordinate of the object' + sources['xwin'].unit = 'pixel' + sources['xwin'].description = 'Windowed X coordinate of the object' + sources['ywin'].unit = 'pixel' + sources['ywin'].description = 'Windowed Y coordinate of the object' + sources['xpeak'].unit = 'pixel' + sources['xpeak'].description = 'X coordinate of the peak' + sources['ypeak'].unit = 'pixel' + sources['ypeak'].description = 'Windowed Y coordinate of the peak' + sources['flux'].unit = 'count' + sources['flux'].description = 'Flux within a Kron-like elliptical aperture' + sources['fluxerr'].unit = 'count' + sources['fluxerr'].description = 'Error on the flux within Kron aperture' + sources['peak'].unit = 'count' + sources['peak'].description = 'Peak flux (flux at xpeak, ypeak)' for diameter in [1, 2, 3, 4, 5, 6]: - catalog['fluxaper{0}'.format(diameter)].unit = 'count' - catalog['fluxaper{0}'.format(diameter)].description = 'Flux from fixed circular aperture: {0}" diameter'.format(diameter) - catalog['fluxerr{0}'.format(diameter)].unit = 'count' - catalog['fluxerr{0}'.format(diameter)].description = 'Error on Flux from circular aperture: {0}"'.format(diameter) - - catalog['background'].unit = 'count' - catalog['background'].description = 'Average background value in the aperture' - catalog['fwhm'].unit = 'pixel' - catalog['fwhm'].description = 'FWHM of the object' - catalog['fwtm'].unit = 'pixel' - catalog['fwtm'].description = 'Full-Width Tenth Maximum' - catalog['a'].unit = 'pixel' - catalog['a'].description = 'Semi-major axis of the object' - catalog['b'].unit = 'pixel' - catalog['b'].description = 'Semi-minor axis of the object' - catalog['theta'].unit = 'degree' - catalog['theta'].description = 'Position angle of the object' - catalog['kronrad'].unit = 'pixel' - catalog['kronrad'].description = 'Kron radius used for extraction' - catalog['ellipticity'].description = 'Ellipticity' - catalog['fluxrad25'].unit = 'pixel' - catalog['fluxrad25'].description = 'Radius containing 25% of the flux' - catalog['fluxrad50'].unit = 'pixel' - catalog['fluxrad50'].description = 'Radius containing 50% of the flux' - catalog['fluxrad75'].unit = 'pixel' - catalog['fluxrad75'].description = 'Radius containing 75% of the flux' - catalog['x2'].unit = 'pixel^2' - catalog['x2'].description = 'Variance on X coordinate of the object' - catalog['y2'].unit = 'pixel^2' - catalog['y2'].description = 'Variance on Y coordinate of the object' - catalog['xy'].unit = 'pixel^2' - catalog['xy'].description = 'XY covariance of the object' - catalog['flag'].description = 'Bit mask of extraction/photometry flags' - - catalog.sort('flux') - catalog.reverse() + sources['fluxaper{0}'.format(diameter)].unit = 'count' + sources['fluxaper{0}'.format(diameter)].description = 'Flux from fixed circular aperture: {0}" diameter'.format(diameter) + sources['fluxerr{0}'.format(diameter)].unit = 'count' + sources['fluxerr{0}'.format(diameter)].description = 'Error on Flux from circular aperture: {0}"'.format(diameter) + + sources['background'].unit = 'count' + sources['background'].description = 'Average background value in the aperture' + sources['fwhm'].unit = 'pixel' + sources['fwhm'].description = 'FWHM of the object' + sources['fwtm'].unit = 'pixel' + sources['fwtm'].description = 'Full-Width Tenth Maximum' + sources['a'].unit = 'pixel' + sources['a'].description = 'Semi-major axis of the object' + sources['b'].unit = 'pixel' + sources['b'].description = 'Semi-minor axis of the object' + sources['theta'].unit = 'degree' + sources['theta'].description = 'Position angle of the object' + sources['kronrad'].unit = 'pixel' + sources['kronrad'].description = 'Kron radius used for extraction' + sources['ellipticity'].description = 'Ellipticity' + sources['fluxrad25'].unit = 'pixel' + sources['fluxrad25'].description = 'Radius containing 25% of the flux' + sources['fluxrad50'].unit = 'pixel' + sources['fluxrad50'].description = 'Radius containing 50% of the flux' + sources['fluxrad75'].unit = 'pixel' + sources['fluxrad75'].description = 'Radius containing 75% of the flux' + sources['x2'].unit = 'pixel^2' + sources['x2'].description = 'Variance on X coordinate of the object' + sources['y2'].unit = 'pixel^2' + sources['y2'].description = 'Variance on Y coordinate of the object' + sources['xy'].unit = 'pixel^2' + sources['xy'].description = 'XY covariance of the object' + sources['flag'].description = 'Bit mask of extraction/photometry flags' + + sources.sort('flux') + sources.reverse() # Save some background statistics in the header - mean_background = stats.sigma_clipped_mean(bkg.back(), 5.0) + mean_background = stats.sigma_clipped_mean(bkg.background, 5.0) image.meta['L1MEAN'] = (mean_background, '[counts] Sigma clipped mean of frame background') - median_background = np.median(bkg.back()) + median_background = np.median(bkg.background) image.meta['L1MEDIAN'] = (median_background, '[counts] Median of frame background') - std_background = stats.robust_standard_deviation(bkg.back()) + std_background = stats.robust_standard_deviation(bkg.background) image.meta['L1SIGMA'] = (std_background, '[counts] Robust std dev of frame background') # Save some image statistics to the header - good_objects = catalog['flag'] == 0 + good_objects = sources['flag'] == 0 for quantity in ['fwhm', 'ellipticity', 'theta']: - good_objects = np.logical_and(good_objects, np.logical_not(np.isnan(catalog[quantity]))) + good_objects = np.logical_and(good_objects, np.logical_not(np.isnan(sources[quantity]))) if good_objects.sum() == 0: image.meta['L1FWHM'] = ('NaN', '[arcsec] Frame FWHM in arcsec') image.meta['L1FWTM'] = ('NaN', 'Ratio of FWHM to Full-Width Tenth Max') @@ -250,23 +246,23 @@ def do_stage(self, image): image.meta['L1ELLIP'] = ('NaN', 'Mean image ellipticity (1-B/A)') image.meta['L1ELLIPA'] = ('NaN', '[deg] PA of mean image ellipticity') else: - seeing = np.nanmedian(catalog['fwhm'][good_objects]) * image.pixel_scale + seeing = np.nanmedian(sources['fwhm'][good_objects]) * image.pixel_scale image.meta['L1FWHM'] = (seeing, '[arcsec] Frame FWHM in arcsec') - image.meta['L1FWTM'] = (np.nanmedian(catalog['fwtm'][good_objects] / catalog['fwhm'][good_objects]), + image.meta['L1FWTM'] = (np.nanmedian(sources['fwtm'][good_objects] / sources['fwhm'][good_objects]), 'Ratio of FWHM to Full-Width Tenth Max') - mean_ellipticity = stats.sigma_clipped_mean(catalog['ellipticity'][good_objects], 3.0) + mean_ellipticity = stats.sigma_clipped_mean(sources['ellipticity'][good_objects], 3.0) image.meta['L1ELLIP'] = (mean_ellipticity, 'Mean image ellipticity (1-B/A)') - mean_position_angle = stats.sigma_clipped_mean(catalog['theta'][good_objects], 3.0) - image.meta['L1ELLIPA'] = (mean_position_angle,'[deg] PA of mean image ellipticity') + mean_position_angle = stats.sigma_clipped_mean(sources['theta'][good_objects], 3.0) + image.meta['L1ELLIPA'] = (mean_position_angle, '[deg] PA of mean image ellipticity') logging_tags = {key: float(image.meta[key]) for key in ['L1MEAN', 'L1MEDIAN', 'L1SIGMA', 'L1FWHM', 'L1ELLIP', 'L1ELLIPA']} logger.info('Extracted sources', image=image, extra_tags=logging_tags) # adding catalog (a data table) to the appropriate images attribute. - image.add_or_update(DataTable(catalog, name='CAT')) + image.add_or_update(DataTable(sources, name='CAT')) except Exception: logger.error(logs.format_exception(), image=image) return image diff --git a/banzai/settings.py b/banzai/settings.py index 63953e30..fdfb98b5 100644 --- a/banzai/settings.py +++ b/banzai/settings.py @@ -80,7 +80,7 @@ 'elp': {'minute': 0, 'hour': 23}, 'ogg': {'minute': 0, 'hour': 3}} -ASTROMETRY_SERVICE_URL = os.getenv('ASTROMETRY_SERVICE_URL', 'http://astrometry.lco.gtn/catalog/') +ASTROMETRY_SERVICE_URL = os.getenv('ASTROMETRY_SERVICE_URL', ' ') CALIBRATION_FILENAME_FUNCTIONS = {'BIAS': ('banzai.utils.file_utils.config_to_filename', 'banzai.utils.file_utils.ccdsum_to_filename'), diff --git a/banzai/tests/test_end_to_end.py b/banzai/tests/test_end_to_end.py index f748aeb0..b879ca69 100644 --- a/banzai/tests/test_end_to_end.py +++ b/banzai/tests/test_end_to_end.py @@ -124,7 +124,7 @@ def get_expected_number_of_calibrations(raw_filename_pattern, calibration_type): if site in frame['filename'] and instrument in frame['filename'] and dayobs in frame['filename'] and raw_filename_pattern in frame['filename'] ] - if 'calibration_type.lower()' == 'skyflat': + if calibration_type.lower() == 'skyflat': # Group by filter observed_filters = [] for frame in raw_frames_for_this_dayobs: @@ -244,7 +244,7 @@ def test_if_stacked_flat_frame_was_created(self): @pytest.mark.e2e @pytest.mark.science_files class TestScienceFileCreation: - @pytest.fixture(autouse=True) + @pytest.fixture(autouse=True, scope='class') @mock.patch('banzai.utils.observation_utils.requests.get', side_effect=observation_portal_side_effect) def reduce_science_frames(self, mock_observation_portal): run_reduce_individual_frames('e00.fits') @@ -254,7 +254,7 @@ def test_if_science_frames_were_created(self): created_files = [] for day_obs in DAYS_OBS: expected_files += [filename.replace('e00', 'e91') - for filename in TEST_FRAMES['filename']] + for filename in TEST_FRAMES['filename'] if 'e00.fits' in filename] created_files += [os.path.basename(filename) for filename in glob(os.path.join(DATA_ROOT, day_obs, 'processed', '*e91*'))] assert len(expected_files) > 0 @@ -266,6 +266,6 @@ def test_that_photometric_calibration_succeeded(self): for day_obs in DAYS_OBS: science_files += [filepath for filepath in glob(os.path.join(DATA_ROOT, day_obs, 'processed', '*e91*'))] - zeropoints = [fits.open(file)['SCI'].header.get('L1ZP') for file in science_files] + zeropoints = [fits.open(filename)['SCI'].header.get('L1ZP') for filename in science_files] # check that at least one of our images contains a zeropoint assert zeropoints.count(None) != len(zeropoints) diff --git a/docs/index.rst b/docs/index.rst index ac03b777..79ed3331 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -107,20 +107,18 @@ our models are 94% complete for our training data. Pixels flagged as cosmic-ray Source Detection ================ -Source detection uses the "Source Extraction in Python" (SEP; https://github.com/kbarbary/sep). -This is similar to Source Extractor, but is written purely in Python and Cython. This allows more -customization. +Source detection uses the "astropy.photutils" image segmentation. +This is similar to Source Extractor, but is written purely in Python. This allows more customization than the original SExtractor. We estimate the background by taking a 3x3 median filter of the image and the doing a 32x32 block average of the image. -We use the default match filter for source detection that is provided by SEP. +We use the default match filter for source detection that is provided in Source Extractor. We include the proper match filter normalization though so the extraction will be slightly different. We do aperture photometry using an elliptical aperture that is set by 2.5 times the Kron radius. This produces approximately the same results as ``FLUX_AUTO`` from SExtractor. -We set the source detection limit at 3 times the global rms of the image. ``MINAREA`` is set to 5, -the default. This should minimize false detections, but may miss the faintest sources. +We set the source detection limit at 2.5 times the uncertainty image. ``MINAREA`` is set to 9. This should minimize false detections, but may miss the faintest sources. To assess the image quality, we estimate the full-width half maximum (FWHM) of the stars in the image. We reject any sources that have a FWHM of less than a pixel to ensure that they do not bias our results. The PSFs for LCO are @@ -132,6 +130,13 @@ estimate of the FWHM. This ensures that we do not underestimate the FWHM when de we take the robust standard deviation (see below) to estimate the overall FWHM of the image. This value is recorded in the header under the L1FWHM keyword. +Flags are as follows: +- 1: Source has bad pixels in the image segmentation +- 2: Object is deblended +- 4: Source has saturated pixels in the image segmentation +- 8: Source kron aperture falls off the image +- 16: Source has cosmic ray pixels in the image segmentation + Astrometry ========== diff --git a/helm-chart/banzai/templates/listener.yaml b/helm-chart/banzai/templates/listener.yaml index d4afa66b..9d3429c9 100644 --- a/helm-chart/banzai/templates/listener.yaml +++ b/helm-chart/banzai/templates/listener.yaml @@ -98,7 +98,6 @@ spec: - "--db-address=$(DB_ADDRESS)" - "--broker-url=$(FITS_BROKER)" - "--queue-name=$(QUEUE_NAME)" - - "--no-file-cache" env: {{- include "banzai.Env" . | nindent 12 }} resources: @@ -128,7 +127,6 @@ spec: - "--db-address=$(DB_ADDRESS)" - "--broker-url=$(FITS_BROKER)" - "--log-level=info" - - "--no-file-cache" env: {{- include "banzai.Env" . | nindent 12 }} resources: