Problems with running astrosource on STELLA telescope data

[turchis]

Hello,

I'm trying to run inside a directory with already reduced and SEP embedded data (acquired from STELLA telescope) the following command:
'astrosource --ra 330.210 --dec 7.060 --indir /export/work/marcotu/COBIPULSE-North/STELLA/DATA/reduced/J2212/gfilter/ --format fits --full --lowcounts 1000 --hicounts 50000 --thresholdcounts 1000000 --closerejectd 5.0'
but I get the error message:
'CRITICAL | No files of type ".fits" found in /export/work/marcotu/COBIPULSE-North/STELLA/DATA/reduced/J2212/gfilter/'

Since this dataset was already reduced, before running astrosource I didn't run the entire BANZAI pipeline, but just the commands for performing the source extraction and append the catalogue of detected sources in all the images. These commands are found in the script photometry.py of the BANZAI package https://github.com/LCOGT/banzai, and I just readjusted these in the following script:
'import logging
from urllib.parse import urljoin

import numpy as np
from astropy.table import Table
from astropy.io import fits
from astropy.utils.data import get_pkg_data_filename
import sep
from requests import HTTPError

from banzai.utils import stats, array_utils
from banzai.utils.photometry_utils import get_reference_sources, match_catalogs, to_magnitude, fit_photometry
from banzai.data import DataTable

from skimage import measure

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)

def photometry(image,image_file,threshold,min_area):

        # Increase the internal buffer size in sep. This is most necessary for crowded fields.
        ny, nx = image[0].data.shape
        sep.set_extract_pixstack(int(nx * ny - 1))

        data = image[0].data

        # Fits can be backwards byte order, so fix that if need be and subtract
        # the background
        try:
            bkg = sep.Background(data, bw=64, bh=64, fw=3, fh=3)
        except ValueError:
            data = data.byteswap(True).newbyteorder()
            bkg = sep.Background(data, bw=64, bh=64, fw=3, fh=3)
        bkg.subfrom(data)
        bkg_rms = bkg.rms()

        # Do an initial source detection
        sources = sep.extract(data, threshold, minarea=min_area, deblend_cont=0.005, err=bkg_rms)


        # 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]

        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=bkg_rms)
        sources['flux'] = flux
        sources['fluxerr'] = fluxerr
        sources['flag'] |= flag

        # Do circular aperture photometry for diameters of 1" to 6"
        pixel_scale = 0.32
        for diameter in [1, 2, 3, 4, 5, 6]:
            flux, fluxerr, flag = sep.sum_circle(data, sources['x'], sources['y'],
                                                 diameter / 2.0 / pixel_scale, gain=1.0, err=bkg_rms)
            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]

        # Calculate the FWHMs of the stars:
        sources['fwhm'] = np.nan
        sources['fwtm'] = np.nan
        # Here we estimate contours
        for source in sources:
            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],
                                                     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]
                        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)'
        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()

        # Save some background statistics in the header
        mean_background = stats.sigma_clipped_mean(bkg.back(), 5.0)
        image[0].header['L1MEAN'] = (mean_background,
                                '[counts] Sigma clipped mean of frame background')

        median_background = np.median(bkg.back())
        image[0].header['L1MEDIAN'] = (median_background,
                                  '[counts] Median of frame background')

        std_background = stats.robust_standard_deviation(bkg.back())
        image[0].header['L1sigma'] = (std_background,
                                 '[counts] Robust std dev of frame background')

        # Save some image statistics to the header
        good_objects = catalog['flag'] == 0
        for quantity in ['fwhm', 'ellipticity', 'theta']:
            good_objects = np.logical_and(good_objects, np.logical_not(np.isnan(catalog[quantity])))
        if good_objects.sum() == 0:
            image[0].header['L1FWHM'] = ('NaN', '[arcsec] Frame FWHM in arcsec')
            image[0].header['L1FWTM'] = ('NaN', 'Ratio of FWHM to Full-Width Tenth Max')

            image[0].header['L1ELLIP'] = ('NaN', 'Mean image ellipticity (1-B/A)')
            image[0].header['L1ELLIPA'] = ('NaN', '[deg] PA of mean image ellipticity')
        else:
            seeing = np.nanmedian(catalog['fwhm'][good_objects]) * pixel_scale
            image[0].header['L1FWHM'] = (seeing, '[arcsec] Frame FWHM in arcsec')
            image[0].header['L1FWTM'] = (np.nanmedian(catalog['fwtm'][good_objects] / catalog['fwhm'][good_objects]),
                                    'Ratio of FWHM to Full-Width Tenth Max')

            mean_ellipticity = stats.sigma_clipped_mean(catalog['ellipticity'][good_objects], 3.0)
            image[0].header['L1ELLIP'] = (mean_ellipticity, 'Mean image ellipticity (1-B/A)')

            mean_position_angle = stats.sigma_clipped_mean(catalog['theta'][good_objects], 3.0)
            image[0].header['L1ELLIPA'] = (mean_position_angle,'[deg] PA of mean image ellipticity')

        

        fits.append(image_file,np.array(catalog))
        with fits.open(image_file, mode='update') as filehandle:
            filehandle[0].header['EXTNAME'] = 'SCI'
            filehandle[1].header['EXTNAME'] = 'CAT'
            
        return image

#image_file = get_pkg_data_filename('lsc1m009-fl03-20151120-0184-e90.fits')
#image = fits.open('lsc1m009-fl03-20151120-0184-e90.fits')
#photometry(image,image_file=image_file,threshold=1.5,min_area=9)

import os
path = '.'
images = [f for f in os.listdir(path) if f.endswith('.fits')]

for i in range(0,(len(images))):
image_file = get_pkg_data_filename(images[i])
image = fits.open(images[i])
photometry(image,image_file=image_file,threshold=1.5,min_area=9)
print(i)
image.close()'

Could you please help me to understand why astrosource doesn't recognize the fits files? Thanks a lot.

[mfitzasp]

Oooooh! I think this may be a simple one to solve. The code does run with BANZAI data...... but it is expecting that data to be in .fits.fz format rather than just fits (this is the LCO default). The solution is just to incorporate this I think..... although the error you've got up the top says that it cannot find the fits files (which is more likely to be that they are not in that directory?) --- either way, perhaps could you send the data you are trying to use to [email protected] and I will take a look?

[turchis]

Hi Michael, thank you very much for your answer. It seems unlikely to me that the problem is that the files are in the fits format instead the fits.fz format. Indeed, when I try to run astrosource on some 2022 LCO fits data (extracted using the same command 'funpack' used to uncompress STELLA data) with the command:
"astrosource --ra 136.7046 --dec -21.3724 --indir /export/work/marcotu/COBIPULSE-South/LCO/2022A/J0906/gfilter/ --format fits --full --lowcounts 1000 --hicounts 50000 --thresholdcounts 1000000 --closerejectd 5.0"

I get the following error message:
"Inspecting input files
Traceback (most recent call last):
File "/home/gudrun/marcotu/python-envs/astrosource/bin/astrosource", line 8, in
sys.exit(main())
File "/home/gudrun/marcotu/python-envs/astrosource/lib64/python3.10/site-packages/click/core.py", line 1128, in call
return self.main(*args, **kwargs)
File "/home/gudrun/marcotu/python-envs/astrosource/lib64/python3.10/site-packages/click/core.py", line 1053, in main
rv = self.invoke(ctx)
File "/home/gudrun/marcotu/python-envs/astrosource/lib64/python3.10/site-packages/click/core.py", line 1395, in invoke
return ctx.invoke(self.callback, **ctx.params)
File "/home/gudrun/marcotu/python-envs/astrosource/lib64/python3.10/site-packages/click/core.py", line 754, in invoke
return __callback(*args, **kwargs)
File "/home/gudrun/marcotu/python-envs/astrosource/lib64/python3.10/site-packages/astrosource/main.py", line 75, in main
ts = TimeSeries(indir=parentPath,
File "/home/gudrun/marcotu/python-envs/astrosource/lib64/python3.10/site-packages/astrosource/astrosource.py", line 53, in init
self.files, self.filtercode = gather_files(self.paths, filelist=filelist, filetype=self.format, bjd=bjd)
File "/home/gudrun/marcotu/python-envs/astrosource/lib64/python3.10/site-packages/astrosource/identify.py", line 136, in gather_files
phot_list_temp = export_photometry_files(filelist, paths['parent'], bjd=bjd)
File "/home/gudrun/marcotu/python-envs/astrosource/lib64/python3.10/site-packages/astrosource/identify.py", line 66, in export_photometry_files
filepath = extract_photometry(fitsobj, indir, bjd=bjd)
File "/home/gudrun/marcotu/python-envs/astrosource/lib64/python3.10/site-packages/astrosource/identify.py", line 83, in extract_photometry
outfile = rename_data_file(hdulist[1].header, bjd=bjd)
File "/home/gudrun/marcotu/python-envs/astrosource/lib64/python3.10/site-packages/astrosource/identify.py", line 32, in rename_data_file
filters.append(prihdr['FILTER1'])
File "/home/gudrun/marcotu/python-envs/astrosource/lib64/python3.10/site-packages/astropy/io/fits/header.py", line 157, in getitem
card = self._cards[self._cardindex(key)]
File "/home/gudrun/marcotu/python-envs/astrosource/lib64/python3.10/site-packages/astropy/io/fits/header.py", line 1754, in _cardindex
raise KeyError(f"Keyword {keyword!r} not found.")
KeyError: "Keyword 'FILTER1' not found."

that doesn't make sense to me since in the script identify.py I find the following rows:
"def rename_data_file(prihdr, bjd=False):

prihdrkeys = prihdr.keys()

if any("OBJECT" in s for s in prihdrkeys):
    objectTemp=prihdr['OBJECT'].replace('-','d').replace('+','p').replace('.','d').replace(' ','').replace('_','').replace('=','e').replace('(','').replace(')','').replace('<','').replace('>','').replace('/','')
else:
    objectTemp="UNKNOWN"

if 'FILTER' in prihdr:
    filterOne=(prihdr['FILTER'])
else:
    filters = []
    filters.append(prihdr['FILTER1'])
    filters.append(prihdr['FILTER2'])
    filters.append(prihdr['FILTER3'])
    filter =list(set(filters))
    filter.remove('air')
    filterOne = filter[0]

"
and filter keyword for 2022 LCO data is FILTER2, that should be identified by the script.

So running astrosource on 2022 LCO data, it seems to find some fits files in the directory, in contrast to what happens with STELLA data.

Thanks again, I just sent to your email the two datasets on which I'm trying to run astrosource, respectively STELLA data and 2022 LCO data.