From 3ca9efef9272e554149df571efb4b8f33d580742 Mon Sep 17 00:00:00 2001 From: Emir Date: Tue, 23 Nov 2021 16:18:02 +0100 Subject: [PATCH 01/10] Add OB type checking Checks whether it is FixedOffset or Autojitter. --- flows/load_image.py | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/flows/load_image.py b/flows/load_image.py index 3693582..681272f 100644 --- a/flows/load_image.py +++ b/flows/load_image.py @@ -149,6 +149,14 @@ def load_image(FILENAME, target_coord=None): s = [hdul[k].header['NAXIS2'], hdul[k].header['NAXIS1']] pix = w.all_world2pix(target_radec, 0).flatten() if pix[0] >= -0.5 and pix[1] >= -0.5 and pix[0] <= s[1]-0.5 and pix[1] <= s[0]-0.5: + ob_type = hdul[k].header["HIERARCH ESO OCS DET1 IMGNAME"].split('_')[-1] + if "Auto" in ob_type: + image.ob_type = 'Autojitter' + elif "Fixed" in ob_type: + image.ob_type = 'FixedOffset' + # Should we use a default instead of raising? + else: + raise RuntimeError("Image OB Type not AutoJitter or FixedOffset") image.image = np.asarray(hdul[k].data, dtype='float64') image.shape = image.image.shape image.wcs = w From 3084d3bdfe4ce08d7649c5fd7371daa8014be77a Mon Sep 17 00:00:00 2001 From: Emir Date: Wed, 24 Nov 2021 19:44:44 +0100 Subject: [PATCH 02/10] add Instrument and Image classes to load_image --- flows/load_image.py | 738 +++++++++++++++++++++++++++++++++++++++++--- 1 file changed, 693 insertions(+), 45 deletions(-) diff --git a/flows/load_image.py b/flows/load_image.py index 681272f..c5434c3 100644 --- a/flows/load_image.py +++ b/flows/load_image.py @@ -3,6 +3,7 @@ """ Flows photometry code. +.. codeauthor:: Emir Karamehmetoglu .. codeauthor:: Rasmus Handberg """ @@ -14,9 +15,503 @@ from astropy.io import fits from astropy.time import Time from astropy.wcs import WCS, FITSFixedWarning -from . import api +from flows import api +from dataclasses import dataclass, field +import typing +from abc import ABC, abstractmethod + + +# -------------------------------------------------------------------------------------------------- +@dataclass +class FlowsImage: + image: np.ndarray + header: typing.Dict + mask: np.ndarray = None + peakmax: float = None + + def __post_init__(self): + self.shape = self.image.shape + self.wcs = self.create_wcs() + # Make empty mask + if not self.mask: + self.mask = np.zeros_like(self.image, dtype='bool') + self.check_finite() + + def check_finite(self): + self.mask |= ~np.isfinite(self.image) + + def update_mask(self, mask): + self.mask = mask + self.check_finite() + + # def get_observatory_params(self, instrument=None): + # if not instrument: + # instrument = self.instrument + # site = instrument.get_site() + # obstime = instrument.get_obstime() + # photfilter = instrument.get_photfilter() + + def create_wcs(self) -> WCS: + with warnings.catch_warnings(): + warnings.simplefilter('ignore', category=FITSFixedWarning) + return WCS(header=self.header, relax=True) + + def create_masked_image(self): + """Warning: this is destructive and will overwrite image data setting masked values to NaN""" + self.image[self.mask] = np.NaN + self.clean = np.ma.masked_array(data=self.image, mask=self.mask, copy=False) + + +class AbstractInstrument(ABC): + peakmax: int = None + + @abstractmethod + def __init__(self): + pass + + @abstractmethod + def get_site(self): + pass + + @abstractmethod + def get_exptime(self): + pass + + @abstractmethod + def get_obstime(self): + pass + + @abstractmethod + def get_photfilter(self): + pass + + +class Instrument(AbstractInstrument): + peakmax: int = None + siteid: int = None + + def __init__(self, image: FlowsImage): + self.image = image + self.image.peakmax = self.peakmax + self.image.site = self.get_site() + self.image.exptime = self.get_exptime() + self.image.obstime = self.get_obstime() + self.image.photfilter = self.get_photfilter() + self.image.create_masked_image() + + def get_site(self): + if self.siteid is not None: + return api.get_site(self.siteid) + + def get_exptime(self): + exptime = self.image.header.get('EXPTIME', None) + if exptime is None: + raise ValueError("Image exposure time could not be extracted") + return exptime + + def get_obstime(self): + '''Default for JD, jd, utc.''' + return Time(self.image.header['JD'], format='jd', scale='utc', + location=self.image.site['EarthLocation']) + + def get_photfilter(self): + return self.image.header['FILTER'] + + +class LCOGT(Instrument): + peakmax: int = 60000 + + def get_site(self): + sites = api.sites.get_all_sites() + site_keywords = {s['site_keyword']: s for s in sites} + site = site_keywords.get(self.image.header['SITE'], None) + return site + + def get_obstime(self): + observatory = coords.EarthLocation.from_geodetic(lat=self.image.header['LATITUDE'], + lon=self.image.header['LONGITUD'], + height=self.image.header['HEIGHT']) + obstime = Time(self.image.header['DATE-OBS'], format='isot', scale='utc', location=observatory) + obstime += 0.5 * self.image.exptime * u.second # Make time centre of exposure + return obstime + + def get_photfilter(self): + photfilter = {'zs': 'zp'}.get(self.image.header['FILTER'], self.image.header['FILTER']) + return photfilter + + +class HAWKI(Instrument): + siteid = 2 # Hard-coded the siteid for ESO Paranal, VLT, UT4 + + def __init__(self, image: FlowsImage): + super().__init__(image) + self.get_obtype() + + def get_obstime(self): + obstime = Time(self.image.header['DATE-OBS'], format='isot', scale='utc', + location=self.image.site['EarthLocation']) + obstime += 0.5 * self.image.exptime * u.second # Make time centre of exposure + return obstime + + def get_obtype(self): + ob_type = self.image.header["HIERARCH ESO OCS DET1 IMGNAME"].split('_')[-1] + if "Auto" in ob_type: + self.image.ob_type = 'Autojitter' + elif "Fixed" in ob_type: + self.image.ob_type = 'FixedOffset' + else: + raise RuntimeError("Image OB Type not AutoJitter or FixedOffset") + + +def get_image_extension(hdul: fits.HDUList, target_coord: coords.SkyCoord = None, fallback_extension: int = None): + # For HAWKI multi-extension images we search the extensions for which one contains + # the target, Create Image from that extension. + target_radec = [[target_coord.icrs.ra.deg, target_coord.icrs.dec.deg]] + + for k in range(1, 5): + w = WCS(header=hdul[k].header, relax=True) + s = [hdul[k].header['NAXIS2'], hdul[k].header['NAXIS1']] + pix = w.all_world2pix(target_radec, 0).flatten() + if -0.5 <= pix[0] <= s[1] - 0.5 and -0.5 <= pix[1] <= s[0] - 0.5: + return k + if not fallback_extension is None: + return fallback_extension + else: + raise RuntimeError(f"Could not find image extension that target is on!") + + +class ALFOSC(Instrument): + # Obtained from http://www.not.iac.es/instruments/detectors/CCD14/LED-linearity/20181026-200-1x1.pdf + peakmax = 80000 # For ALFOSC D, 1x1, 200; the standard for SNe. + siteid = 5 + + def get_obstime(self): + return Time(self.image.header['DATE-AVG'], format='isot', scale='utc', + location=self.image.site['EarthLocation']) + + def get_photfilter(self): + # Sometimes data from NOT does not have the FILTER keyword, + # in which case we have to try to figure out which filter + # was used based on some of the other headers: + if 'FILTER' in self.image.header: + photfilter = { + 'B Bes': 'B', + 'V Bes': 'V', + 'R Bes': 'R', + 'g SDSS': 'gp', + 'r SDSS': 'rp', + 'i SDSS': 'ip', + 'i int': 'ip', # Interference filter + 'u SDSS': 'up', + 'z SDSS': 'zp' + }.get(self.image.header['FILTER'].replace('_', ' '), self.image.header['FILTER']) + else: + filters_used = [] + for check_headers in ('ALFLTNM', 'FAFLTNM', 'FBFLTNM'): + isopen = self.image.header.get(check_headers).strip().lower() != 'open' + if self.image.header.get(check_headers) and isopen: + filters_used.append(self.image.header.get(check_headers).strip()) + if len(filters_used) == 1: + photfilter = { + 'V_Bes 530_80': 'V', + 'R_Bes 650_130': 'R', + "g'_SDSS 480_145": 'gp', + "r'_SDSS 618_148": 'rp', + "i'_SDSS 771_171": 'ip', + 'i_int 797_157': 'ip', # Interference filter + "z'_SDSS 832_LP": 'zp' + }.get(filters_used[0].replace(' ', ' '), filters_used[0]) + else: + raise RuntimeError("Could not determine filter used.") + + return photfilter + + +class NOTCAM(Instrument): + siteid = 5 + + def get_obstime(self): + return Time(self.image.header['DATE-AVG'], format='isot', scale='utc', + location=self.image.site['EarthLocation']) + + def get_photfilter(self): + # Does NOTCAM data sometimes contain a FILTER header? + # if not we have to try to figure out which filter + # was used based on some of the other headers: + if 'FILTER' in self.image.header: + raise RuntimeError("NOTCAM: Filter keyword defined") + filters_used = [] + for check_headers in ('NCFLTNM1', 'NCFLTNM2'): + isopen = self.image.header.get(check_headers).strip().lower() != 'open' + if self.image.header.get(check_headers) and isopen: + filters_used.append(self.image.header.get(check_headers).strip()) + if len(filters_used) == 1: + photfilter = {}.get(filters_used[0], filters_used[0]) + else: + raise RuntimeError("Could not determine filter used.") + return photfilter + + +class PS1(Instrument): + siteid = 6 + + def get_obstime(self): + return Time(self.image.header['MJD-OBS'], format='mjd', scale='utc', + location=self.image.site['EarthLocation']) + + def get_photfilter(self): + photfilter = { + 'g.00000': 'gp', + 'r.00000': 'rp', + 'i.00000': 'ip', + 'z.00000': 'zp' + }.get(self.image.header['FPA.FILTER'], self.image.header['FPA.FILTER']) + return photfilter + + +class Liverpool(Instrument): + siteid = 8 + + def get_obstime(self): + obstime = Time(self.image.header['DATE-OBS'], format='isot', scale='utc', + location=self.image.site['EarthLocation']) + obstime += 0.5 * self.image.exptime * u.second # Make time centre of exposure + return obstime + + def get_photfilter(self): + photfilter = { + 'Bessel-B': 'B', + 'Bessell-B': 'B', + 'Bessel-V': 'V', + 'Bessell-V': 'V', + 'SDSS-U': 'up', + 'SDSS-G': 'gp', + 'SDSS-R': 'rp', + 'SDSS-I': 'ip', + 'SDSS-Z': 'zp' + }.get(self.image.header['FILTER1'], self.image.header['FILTER1']) + return photfilter + + +class Omega2000(Instrument): + siteid = 9 + + def get_obstime(self): + obstime = Time(self.image.header['MJD-OBS'], format='mjd', scale='utc', + location=self.image.site['EarthLocation']) + obstime += 0.5 * self.image.exptime * u.second + return obstime + + +class Swope(Instrument): + siteid = 10 + + def get_photfilter(self): + photfilter = { + 'u': 'up', + 'g': 'gp', + 'r': 'rp', + 'i': 'ip', + }.get(self.image.header['FILTER'], self.image.header['FILTER']) + return photfilter + + +class Dupont(Instrument): + siteid = 14 + + def get_photfilter(self): + photfilter = { + 'u': 'up', + 'g': 'gp', + 'r': 'rp', + 'i': 'ip', + }.get(self.image.header['FILTER'], self.image.header['FILTER']) + return photfilter + + +class RetroCam(Instrument): + siteid = 16 + + def get_photfilter(self): + photfilter = { + 'Yc': 'Y', + 'Hc': 'H', + 'Jo': 'J', + }.get(self.image.header['FILTER'], self.image.header['FILTER']) + return photfilter + + +class Baade(Instrument): + siteid = 11 + + def get_exptime(self): + exptime = super().get_exptime() + exptime *= int(self.image.header['NCOMBINE']) # EXPTIME is only for a single exposure + return exptime + + def get_photfilter(self): + photfilter = { + 'Ks': 'K', + 'J1': 'Y', + }.get(self.image.header['FILTER'], self.image.header['FILTER']) + return photfilter + + +class Sofi(Instrument): + siteid = 12 + + def get_obstime(self): + if 'TMID' in self.image.header: + obstime = Time(self.image.header['TMID'], format='mjd', scale='utc', + location=self.image.site['EarthLocation']) + else: + obstime = Time(self.image.header['MJD-OBS'], format='mjd', scale='utc', + location=self.image.site['EarthLocation']) + obstime += 0.5 * self.image.exptime * u.second # Make time centre of exposure + return obstime + + def get_photfilter(self): + hdr = self.image.header + photfilter_translate = { + 'Ks': 'K' + } + if 'FILTER' in hdr: + photfilter = photfilter_translate.get(hdr['FILTER'], hdr['FILTER']) + else: + filters_used = [] + for check_headers in ('ESO INS FILT1 ID', 'ESO INS FILT2 ID'): + if hdr.get(check_headers) and hdr.get(check_headers).strip().lower() != 'open': + filters_used.append(hdr.get(check_headers).strip()) + if len(filters_used) == 1: + photfilter = photfilter_translate.get(filters_used[0], filters_used[0]) + else: + raise RuntimeError("Could not determine filter used.") + return photfilter + + +class EFOSC(Instrument): + siteid = 15 + + def get_obstime(self): + obstime = Time(self.image.header['DATE-OBS'], format='isot', scale='utc', + location=self.image.site['EarthLocation']) + obstime += 0.5 * self.image.exptime * u.second # Make time centre of exposure + return obstime + + def get_photfilter(self): + hdr = self.image.header + photfilter = { + 'g782': 'gp', + 'r784': 'rp', + 'i705': 'ip', + 'B639': 'B', + 'V641': 'V' + }.get(hdr['FILTER'], hdr['FILTER']) + return photfilter + + +class AstroNIRCam(Instrument): + siteid = 13 + + def get_exptime(self): + return self.image.header.get('FULL_EXP', super().get_exptime()) + + def get_obstime(self): + hdr = self.image.header + if 'MIDPOINT' in hdr: + obstime = Time(hdr['MIDPOINT'], format='isot', scale='utc', location=self.image.site['EarthLocation']) + else: + obstime = Time(hdr['MJD-AVG'], format='mjd', scale='utc', location=self.image.site['EarthLocation']) + return obstime + + def get_photfilter(self): + hdr = self.image.header + photfilter = { + 'H_Open': 'H', + 'K_Open': 'K', + }.get(hdr['FILTER'], hdr['FILTER']) + return photfilter + + +class OmegaCam(Instrument): + siteid = 18 # Hard-coded the siteid for ESO VLT Survey telescope + + def get_obstime(self): + obstime = Time(self.image.header['MJD-OBS'], format='mjd', scale='utc', + location=self.image.site['EarthLocation']) + obstime += 0.5 * self.image.exptime * u.second # Make time centre of exposure + return obstime + + def get_photfilter(self): + hdr = self.image.header + photfilter = { + 'i_SDSS': 'ip' + }.get(hdr['ESO INS FILT1 NAME'], hdr['ESO INS FILT1 NAME']) + return photfilter + + +class AndiCam(Instrument): + siteid = 20 # Hard-coded the siteid for ANDICAM at Cerro Tololo Interamerican Observatory (CTIO) + + def get_obstime(self): + obstime = super().get_obstime() + obstime += 0.5 * self.image.exptime * u.second + return obstime + + def get_photfilter(self): + return self.image.header['CCDFLTID'] + + +class PairTel(Instrument): + siteid = 21 + + def get_obstime(self): + hdr = self.image.header + time_start = Time(hdr['STRT_CPU'], format='iso', scale='utc', location=self.image.site['EarthLocation']) + time_stop = Time(hdr['STOP_CPU'], format='iso', scale='utc', location=self.image.site['EarthLocation']) + obstime = time_start + 0.5 * (time_stop - time_start) + return obstime + + def get_photfilter(self): + hdr = self.image.header + photfilter = { + 'j': 'J', + 'h': 'H', + 'k': 'K', + }.get(hdr['FILTER'], hdr['FILTER']) + return photfilter + +class TJO(Instrument): + siteid = 22 + + def get_obstime(self): + obstime = super().get_obstime() + obstime += 0.5 * self.image.exptime * u.second + return obstime + + +instruments = { + 'LCOGT': LCOGT, + 'HAWKI': HAWKI, + 'ALFOSC': ALFOSC, + 'NOTCAM': NOTCAM, + 'PS1': PS1, + 'Liverpool': Liverpool, + 'Omega2000': Omega2000, + 'Swope': Swope, + 'Dupont': Dupont, + 'Retrocam': RetroCam, + 'Baade': Baade, + 'Sofi': Sofi, + 'EFOSC': EFOSC, + 'AstroNIRCam': AstroNIRCam, + 'OmegaCam': OmegaCam, + 'AndiCam': AndiCam, + 'PairTel': PairTel, + 'TJO': TJO +} + -#-------------------------------------------------------------------------------------------------- def edge_mask(img, value=0): """ Create boolean mask of given value near edge of image. @@ -58,7 +553,142 @@ def edge_mask(img, value=0): return mask -#-------------------------------------------------------------------------------------------------- + +def new_load_image(filename: str, target_coord: coords.SkyCoord = None): + """ + Load FITS image using FlowsImage class and Instrument Classes. + + Parameters: + filename (str): Path to FITS file to be loaded. + target_coord (:class:`astropy.coordinates.SkyCoord`): Coordinates of target. + Only used for HAWKI images to determine which image extension to load, + for all other images it is ignored. + + Returns: + FlowsImage: instance of FlowsImage with valuues populated based on instrument. + + .. codeauthor:: Emir Karamehmetoglu + .. codeauthor:: Rasmus Handberg + """ + logger = logging.getLogger(__name__) + ext = 0 # Default extension in HDUList, individual instruments may override this. + mask = None # Instrument can override, default is to only mask all non-finite values, override is additive. + + # Read fits image, Structural Pattern Match to specific instrument. + with fits.open(filename, mode='readonly') as hdul: + hdr = hdul[ext].header + origin = hdr.get('ORIGIN', '') + telescope = hdr.get('TELESCOP', '') + instrument = hdr.get('INSTRUME', '') + + instrument_name = None + # Pattern matching begins here, ideally we use 3.10 pattern matching, or a dictionary lookup. + while instrument_name is None: + # LCOGT + if origin == "LCOGT": + instrument_name = 'LCOGT' + if 'BPM' in hdul: + mask = np.asarray(hdul['BPM'].data, dtype='bool') + else: + logger.warning('LCOGT image does not contain bad pixel map. Not applying mask.') + break + # HAWKI + elif origin == 'ESO-PARANAL' and telescope == 'ESO-VLT-U4' and instrument == 'HAWKI' and hdr.get( + 'PRODCATG') == 'SCIENCE.MEFIMAGE': + instrument_name = 'HAWKI' + + if target_coord is None: + raise ValueError("TARGET_COORD is needed for HAWKI images to find the correct extension") + if not isinstance(target_coord, coords.SkyCoord): + if len(target_coord) == 2: + target_coord = coords.SkyCoord(ra=target_coord[0] * u.deg, dec=target_coord[1] * u.deg, + frame='icrs') + else: + raise ValueError("TARGET_COORD is needed for HAWKI images to find the correct extension") + ext = get_image_extension(hdul, target_coord) # Find the one with the SN in it. + hdr = hdul[ext].header + hdul[0].header + break + + # NOT - ALFOSC + elif telescope == 'NOT' and instrument in ('ALFOSC FASU','ALFOSC_FASU') and hdr.get('OBS_MODE','').lower() == 'imaging': + instrument_name = 'ALFOSC' + break + + elif telescope == 'NOT' and instrument == 'NOTCAM' and hdr.get('OBS_MODE', '').lower() == 'imaging': + instrument_name = 'NOTCAM' + break + + elif hdr.get('FPA.TELESCOPE') == 'PS1' and hdr.get('FPA.INSTRUMENT') == 'GPC1': + instrument_name = 'PS1' + break + + elif telescope == 'Liverpool Telescope': + instrument_name = 'Liverpool' + break + + elif telescope == 'CA 3.5m' and instrument == 'Omega2000': + instrument_name = 'Omega2000' + break + + elif telescope == 'SWO' and hdr.get('SITENAME') == 'LCO': + instrument_name = 'Swope' + break + + elif telescope == 'DUP' and hdr.get('SITENAME') == 'LCO' and instrument == 'Direct/SITe2K-1': + instrument_name = 'Dupont' + break + + elif telescope == 'DUP' and instrument == 'RetroCam': + instrument_name = 'RetroCam' + break + + elif telescope == 'Baade' and hdr.get('SITENAME') == 'LCO' and instrument == 'FourStar': + instrument_name = 'Baade' + break + + elif instrument == 'SOFI' and telescope in ('ESO-NTT', 'other') and ( + origin == 'ESO' or origin.startswith('NOAO-IRAF')): + instrument_name = 'Sofi' + break + + elif telescope == 'ESO-NTT' and instrument == 'EFOSC' and (origin == 'ESO' or origin.startswith('NOAO-IRAF')): + instrument_name = 'EFOSC' + break + + elif telescope == 'SAI-2.5' and instrument == 'ASTRONIRCAM': + instrument_name = 'AstroNIRCam' + break + + elif instrument == 'OMEGACAM' and (origin == 'ESO' or origin.startswith('NOAO-IRAF')): + instrument_name = 'OmegaCam' + break + + elif instrument == 'ANDICAM-CCD' and hdr.get('OBSERVAT') == 'CTIO': + instrument_name = 'AndiCam' + break + + elif telescope == '1.3m PAIRITEL' and instrument == '2MASS Survey cam': + instrument_name = 'PairTel' + break + + elif (origin == 'OAdM' or origin.startswith('NOAO-IRAF')) and telescope == 'TJO' and instrument in ( + 'MEIA3', 'MEIA2'): + instrument_name = 'TJO' + break + + else: + raise RuntimeError("Could not determine origin of image") + + image = FlowsImage(image=np.asarray(hdul[ext].data, dtype='float64'), header=hdr, mask=mask) + ins = instruments[instrument_name](image) + + return ins.image + + +# get image and wcs solution +# with fits.open(FILENAME, mode='readonly') as hdul: + +# -------------------------------------------------------------------------------------------------- def load_image(FILENAME, target_coord=None): """ Load FITS image. @@ -78,7 +708,7 @@ def load_image(FILENAME, target_coord=None): logger = logging.getLogger(__name__) # Get image and WCS, find stars, remove galaxies - image = type('image', (object,), dict()) # image container + image = type('image', (object,), dict()) # image container # get image and wcs solution with fits.open(FILENAME, mode='readonly') as hdul: @@ -111,8 +741,8 @@ def load_image(FILENAME, target_coord=None): image.wcs = WCS(header=hdr, relax=True) # Values which will be filled out below, depending on the instrument: - image.exptime = hdr.get('EXPTIME', None) # Exposure time * u.second - image.peakmax = None # Maximum value above which data is not to be trusted + image.exptime = hdr.get('EXPTIME', None) # Exposure time * u.second + image.peakmax = None # Maximum value above which data is not to be trusted # Timestamp: if origin == 'LCOGT': @@ -120,9 +750,10 @@ def load_image(FILENAME, target_coord=None): site_keywords = {s['site_keyword']: s for s in sites} image.site = site_keywords.get(hdr['SITE'], None) - observatory = coords.EarthLocation.from_geodetic(lat=hdr['LATITUDE'], lon=hdr['LONGITUD'], height=hdr['HEIGHT']) + observatory = coords.EarthLocation.from_geodetic(lat=hdr['LATITUDE'], lon=hdr['LONGITUD'], + height=hdr['HEIGHT']) image.obstime = Time(hdr['DATE-OBS'], format='isot', scale='utc', location=observatory) - image.obstime += 0.5*image.exptime * u.second # Make time centre of exposure + image.obstime += 0.5 * image.exptime * u.second # Make time centre of exposure image.photfilter = { 'zs': 'zp' @@ -130,13 +761,14 @@ def load_image(FILENAME, target_coord=None): # Get non-linear limit # TODO: Use actual or some fraction of the non-linearity limit - #image.peakmax = hdr.get('MAXLIN') # Presumed non-linearity limit from header - image.peakmax = 60000 # From experience, this one is better. + # image.peakmax = hdr.get('MAXLIN') # Presumed non-linearity limit from header + image.peakmax = 60000 # From experience, this one is better. - elif origin == 'ESO-PARANAL' and telescope == 'ESO-VLT-U4' and instrument == 'HAWKI' and hdr.get('PRODCATG') == 'SCIENCE.MEFIMAGE': - image.site = api.get_site(2) # Hard-coded the siteid for ESO Paranal, VLT, UT4 + elif origin == 'ESO-PARANAL' and telescope == 'ESO-VLT-U4' and instrument == 'HAWKI' and hdr.get( + 'PRODCATG') == 'SCIENCE.MEFIMAGE': + image.site = api.get_site(2) # Hard-coded the siteid for ESO Paranal, VLT, UT4 image.obstime = Time(hdr['DATE-OBS'], format='isot', scale='utc', location=image.site['EarthLocation']) - image.obstime += 0.5*image.exptime * u.second # Make time centre of exposure + image.obstime += 0.5 * image.exptime * u.second # Make time centre of exposure image.photfilter = hdr['FILTER'] # For HAWKI multi-extension images we search the extensions for which one contains @@ -148,7 +780,7 @@ def load_image(FILENAME, target_coord=None): w = WCS(header=hdul[k].header, relax=True) s = [hdul[k].header['NAXIS2'], hdul[k].header['NAXIS1']] pix = w.all_world2pix(target_radec, 0).flatten() - if pix[0] >= -0.5 and pix[1] >= -0.5 and pix[0] <= s[1]-0.5 and pix[1] <= s[0]-0.5: + if pix[0] >= -0.5 and pix[1] >= -0.5 and pix[0] <= s[1] - 0.5 and pix[1] <= s[0] - 0.5: ob_type = hdul[k].header["HIERARCH ESO OCS DET1 IMGNAME"].split('_')[-1] if "Auto" in ob_type: image.ob_type = 'Autojitter' @@ -165,8 +797,9 @@ def load_image(FILENAME, target_coord=None): else: raise RuntimeError("Could not find image extension that target is on") - elif telescope == 'NOT' and instrument in ('ALFOSC FASU', 'ALFOSC_FASU') and hdr.get('OBS_MODE', '').lower() == 'imaging': - image.site = api.get_site(5) # Hard-coded the siteid for NOT + elif telescope == 'NOT' and instrument in ('ALFOSC FASU', 'ALFOSC_FASU') and hdr.get('OBS_MODE', + '').lower() == 'imaging': + image.site = api.get_site(5) # Hard-coded the siteid for NOT image.obstime = Time(hdr['DATE-AVG'], format='isot', scale='utc', location=image.site['EarthLocation']) # Sometimes data from NOT does not have the FILTER keyword, @@ -180,7 +813,7 @@ def load_image(FILENAME, target_coord=None): 'g SDSS': 'gp', 'r SDSS': 'rp', 'i SDSS': 'ip', - 'i int': 'ip', # Interference filter + 'i int': 'ip', # Interference filter 'u SDSS': 'up', 'z SDSS': 'zp' }.get(hdr['FILTER'].replace('_', ' '), hdr['FILTER']) @@ -197,7 +830,7 @@ def load_image(FILENAME, target_coord=None): "g'_SDSS 480_145": 'gp', "r'_SDSS 618_148": 'rp', "i'_SDSS 771_171": 'ip', - 'i_int 797_157': 'ip', # Interference filter + 'i_int 797_157': 'ip', # Interference filter "z'_SDSS 832_LP": 'zp' }.get(filters_used[0].replace(' ', ' '), filters_used[0]) else: @@ -206,10 +839,10 @@ def load_image(FILENAME, target_coord=None): # Get non-linear limit # Obtained from http://www.not.iac.es/instruments/detectors/CCD14/LED-linearity/20181026-200-1x1.pdf # TODO: grab these from a table for all detector setups of ALFOSC - image.peakmax = 80000 # For ALFOSC D, 1x1, 200; the standard for SNe. + image.peakmax = 80000 # For ALFOSC D, 1x1, 200; the standard for SNe. elif telescope == 'NOT' and instrument == 'NOTCAM' and hdr.get('OBS_MODE', '').lower() == 'imaging': - image.site = api.get_site(5) # Hard-coded the siteid for NOT + image.site = api.get_site(5) # Hard-coded the siteid for NOT image.obstime = Time(hdr['DATE-AVG'], format='isot', scale='utc', location=image.site['EarthLocation']) # Does NOTCAM data sometimes contain a FILTER header? @@ -232,7 +865,7 @@ def load_image(FILENAME, target_coord=None): image.mask |= edge_mask(image.image, value=0) elif hdr.get('FPA.TELESCOPE') == 'PS1' and hdr.get('FPA.INSTRUMENT') == 'GPC1': - image.site = api.get_site(6) # Hard-coded the siteid for Pan-STARRS1 + image.site = api.get_site(6) # Hard-coded the siteid for Pan-STARRS1 image.obstime = Time(hdr['MJD-OBS'], format='mjd', scale='utc', location=image.site['EarthLocation']) image.photfilter = { @@ -244,9 +877,9 @@ def load_image(FILENAME, target_coord=None): elif telescope == 'Liverpool Telescope': # Liverpool telescope - image.site = api.get_site(8) # Hard-coded the siteid for Liverpool Telescope + image.site = api.get_site(8) # Hard-coded the siteid for Liverpool Telescope image.obstime = Time(hdr['DATE-OBS'], format='isot', scale='utc', location=image.site['EarthLocation']) - image.obstime += 0.5*image.exptime * u.second # Make time centre of exposure + image.obstime += 0.5 * image.exptime * u.second # Make time centre of exposure image.photfilter = { 'Bessel-B': 'B', 'Bessell-B': 'B', @@ -261,13 +894,13 @@ def load_image(FILENAME, target_coord=None): elif telescope == 'CA 3.5m' and instrument == 'Omega2000': # Calar Alto 3.5m (Omege2000) - image.site = api.get_site(9) # Hard-coded the siteid for Calar Alto 3.5m + image.site = api.get_site(9) # Hard-coded the siteid for Calar Alto 3.5m image.obstime = Time(hdr['MJD-OBS'], format='mjd', scale='utc', location=image.site['EarthLocation']) - image.obstime += 0.5*image.exptime * u.second # Make time centre of exposure + image.obstime += 0.5 * image.exptime * u.second # Make time centre of exposure image.photfilter = hdr['FILTER'] elif telescope == 'SWO' and hdr.get('SITENAME') == 'LCO': - image.site = api.get_site(10) # Hard-coded the siteid for Swope, Las Campanas Observatory + image.site = api.get_site(10) # Hard-coded the siteid for Swope, Las Campanas Observatory image.obstime = Time(hdr['JD'], format='jd', scale='utc', location=image.site['EarthLocation']) image.photfilter = { 'u': 'up', @@ -277,7 +910,7 @@ def load_image(FILENAME, target_coord=None): }.get(hdr['FILTER'], hdr['FILTER']) elif telescope == 'DUP' and hdr.get('SITENAME') == 'LCO' and instrument == 'Direct/SITe2K-1': - image.site = api.get_site(14) # Hard-coded the siteid for Du Pont, Las Campanas Observatory + image.site = api.get_site(14) # Hard-coded the siteid for Du Pont, Las Campanas Observatory image.obstime = Time(hdr['JD'], format='jd', scale='utc', location=image.site['EarthLocation']) image.photfilter = { 'u': 'up', @@ -287,7 +920,7 @@ def load_image(FILENAME, target_coord=None): }.get(hdr['FILTER'], hdr['FILTER']) elif telescope == 'DUP' and instrument == 'RetroCam': - image.site = api.get_site(16) # Hard-coded the siteid for Du Pont, Las Campanas Observatory + image.site = api.get_site(16) # Hard-coded the siteid for Du Pont, Las Campanas Observatory image.obstime = Time(hdr['JD'], format='jd', scale='utc', location=image.site['EarthLocation']) image.photfilter = { 'Yc': 'Y', @@ -296,21 +929,22 @@ def load_image(FILENAME, target_coord=None): }.get(hdr['FILTER'], hdr['FILTER']) elif telescope == 'Baade' and hdr.get('SITENAME') == 'LCO' and instrument == 'FourStar': - image.site = api.get_site(11) # Hard-coded the siteid for Swope, Las Campanas Observatory + image.site = api.get_site(11) # Hard-coded the siteid for Swope, Las Campanas Observatory image.obstime = Time(hdr['JD'], format='jd', scale='utc', location=image.site['EarthLocation']) image.photfilter = { 'Ks': 'K', 'J1': 'Y', }.get(hdr['FILTER'], hdr['FILTER']) - image.exptime *= int(hdr['NCOMBINE']) # EXPTIME is only for a single exposure + image.exptime *= int(hdr['NCOMBINE']) # EXPTIME is only for a single exposure - elif instrument == 'SOFI' and telescope in ('ESO-NTT', 'other') and (origin == 'ESO' or origin.startswith('NOAO-IRAF')): - image.site = api.get_site(12) # Hard-coded the siteid for SOFT, ESO NTT + elif instrument == 'SOFI' and telescope in ('ESO-NTT', 'other') and ( + origin == 'ESO' or origin.startswith('NOAO-IRAF')): + image.site = api.get_site(12) # Hard-coded the siteid for SOFT, ESO NTT if 'TMID' in hdr: image.obstime = Time(hdr['TMID'], format='mjd', scale='utc', location=image.site['EarthLocation']) else: image.obstime = Time(hdr['MJD-OBS'], format='mjd', scale='utc', location=image.site['EarthLocation']) - image.obstime += 0.5*image.exptime * u.second # Make time centre of exposure + image.obstime += 0.5 * image.exptime * u.second # Make time centre of exposure # Photometric filter: photfilter_translate = { @@ -332,9 +966,9 @@ def load_image(FILENAME, target_coord=None): image.mask |= edge_mask(image.image, value=0) elif telescope == 'ESO-NTT' and instrument == 'EFOSC' and (origin == 'ESO' or origin.startswith('NOAO-IRAF')): - image.site = api.get_site(15) # Hard-coded the siteid for EFOSC, ESO NTT + image.site = api.get_site(15) # Hard-coded the siteid for EFOSC, ESO NTT image.obstime = Time(hdr['DATE-OBS'], format='isot', scale='utc', location=image.site['EarthLocation']) - image.obstime += 0.5*image.exptime * u.second # Make time centre of exposure + image.obstime += 0.5 * image.exptime * u.second # Make time centre of exposure image.photfilter = { 'g782': 'gp', 'r784': 'rp', @@ -344,7 +978,7 @@ def load_image(FILENAME, target_coord=None): }.get(hdr['FILTER'], hdr['FILTER']) elif telescope == 'SAI-2.5' and instrument == 'ASTRONIRCAM': - image.site = api.get_site(13) # Hard-coded the siteid for Caucasus Mountain Observatory + image.site = api.get_site(13) # Hard-coded the siteid for Caucasus Mountain Observatory if 'MIDPOINT' in hdr: image.obstime = Time(hdr['MIDPOINT'], format='isot', scale='utc', location=image.site['EarthLocation']) else: @@ -356,24 +990,25 @@ def load_image(FILENAME, target_coord=None): image.exptime = hdr.get('FULL_EXP', image.exptime) elif instrument == 'OMEGACAM' and (origin == 'ESO' or origin.startswith('NOAO-IRAF')): - image.site = api.get_site(18) # Hard-coded the siteid for ESO VLT Survey telescope + image.site = api.get_site(18) # Hard-coded the siteid for ESO VLT Survey telescope image.obstime = Time(hdr['MJD-OBS'], format='mjd', scale='utc', location=image.site['EarthLocation']) - image.obstime += 0.5*image.exptime * u.second # Make time centre of exposure + image.obstime += 0.5 * image.exptime * u.second # Make time centre of exposure image.photfilter = { 'i_SDSS': 'ip' }.get(hdr['ESO INS FILT1 NAME'], hdr['ESO INS FILT1 NAME']) elif instrument == 'ANDICAM-CCD' and hdr.get('OBSERVAT') == 'CTIO': - image.site = api.get_site(20) # Hard-coded the siteid for ANDICAM at Cerro Tololo Interamerican Observatory (CTIO) + image.site = api.get_site( + 20) # Hard-coded the siteid for ANDICAM at Cerro Tololo Interamerican Observatory (CTIO) image.obstime = Time(hdr['JD'], format='jd', scale='utc', location=image.site['EarthLocation']) - image.obstime += 0.5*image.exptime * u.second # Make time centre of exposure + image.obstime += 0.5 * image.exptime * u.second # Make time centre of exposure image.photfilter = hdr['CCDFLTID'] elif telescope == '1.3m PAIRITEL' and instrument == '2MASS Survey cam': - image.site = api.get_site(21) # Hard-coded the siteid for Peters Automated InfraRed Imaging TELescope + image.site = api.get_site(21) # Hard-coded the siteid for Peters Automated InfraRed Imaging TELescope time_start = Time(hdr['STRT_CPU'], format='iso', scale='utc', location=image.site['EarthLocation']) time_stop = Time(hdr['STOP_CPU'], format='iso', scale='utc', location=image.site['EarthLocation']) - image.obstime = time_start + 0.5*(time_stop - time_start) + image.obstime = time_start + 0.5 * (time_stop - time_start) image.photfilter = { 'j': 'J', 'h': 'H', @@ -383,10 +1018,12 @@ def load_image(FILENAME, target_coord=None): # Mask out "halo" of pixels with zero value along edge of image: image.mask |= edge_mask(image.image, value=0) - elif (origin == 'OAdM' or origin.startswith('NOAO-IRAF')) and telescope == 'TJO' and instrument in ('MEIA3', 'MEIA2'): - image.site = api.get_site(22) # Hard-coded the siteid for Telescopi Joan Oró (TJO) at Observatori Astronòmic del Montsec + elif (origin == 'OAdM' or origin.startswith('NOAO-IRAF')) and telescope == 'TJO' and instrument in ( + 'MEIA3', 'MEIA2'): + image.site = api.get_site( + 22) # Hard-coded the siteid for Telescopi Joan Oró (TJO) at Observatori Astronòmic del Montsec image.obstime = Time(hdr['JD'], format='jd', scale='utc', location=image.site['EarthLocation']) - image.obstime += 0.5*image.exptime * u.second # Make time centre of exposure + image.obstime += 0.5 * image.exptime * u.second # Make time centre of exposure image.photfilter = hdr['FILTER'] else: @@ -401,3 +1038,14 @@ def load_image(FILENAME, target_coord=None): image.clean = np.ma.masked_array(data=image.image, mask=image.mask, copy=False) return image + + +def main(): + filename = "../ADP.2021-11-08T12_34_46.383.fits" + img = new_load_image(filename, (325.772590825, -17.54410925)) + print(img.peakmax) + return img + + +if __name__ == '__main__': + main() From 84c32d47b2d43aad740c0945b0e2e7afc665871a Mon Sep 17 00:00:00 2001 From: Emir Date: Tue, 30 Nov 2021 11:34:54 +0100 Subject: [PATCH 03/10] Implement strategy design pattern to load_image.py This is done to help in adding in Edge correction and make the code more manageable and extendable. --- flows/load_image.py | 107 ++++++++++++++++++++++++++++++++++++++------ 1 file changed, 94 insertions(+), 13 deletions(-) diff --git a/flows/load_image.py b/flows/load_image.py index c5434c3..e605d2c 100644 --- a/flows/load_image.py +++ b/flows/load_image.py @@ -15,11 +15,15 @@ from astropy.io import fits from astropy.time import Time from astropy.wcs import WCS, FITSFixedWarning +from typing import Tuple + from flows import api from dataclasses import dataclass, field import typing from abc import ABC, abstractmethod +logger = logging.getLogger(__name__) # Singleton logger instance + # -------------------------------------------------------------------------------------------------- @dataclass @@ -44,13 +48,6 @@ def update_mask(self, mask): self.mask = mask self.check_finite() - # def get_observatory_params(self, instrument=None): - # if not instrument: - # instrument = self.instrument - # site = instrument.get_site() - # obstime = instrument.get_obstime() - # photfilter = instrument.get_photfilter() - def create_wcs(self) -> WCS: with warnings.catch_warnings(): warnings.simplefilter('ignore', category=FITSFixedWarning) @@ -61,6 +58,84 @@ def create_masked_image(self): self.image[self.mask] = np.NaN self.clean = np.ma.masked_array(data=self.image, mask=self.mask, copy=False) + def set_edge_rows_to_value(self, y: Tuple[float] = None, value: typing.Union[int, float, np.NaN] = 0): + if y is None: + pass + for row in y: + self.image[row] = value + + def set_edge_columns_to_value(self, x: Tuple[float] = None, value: typing.Union[int, float, np.NaN] = 0): + if x is None: + pass + for col in x: + self.image[:, col] = value + + @staticmethod + def get_edge_mask(img: np.ndarray, value: typing.Union[int, float, np.NaN] = 0): + """ + Create boolean mask of given value near edge of image. + + Parameters: + img (ndarray): image with values for masking. + value (float): Value to detect near edge. Default=0. + + Returns: + ndarray: Pixel mask with given values on the edge of image. + + .. codeauthor:: Rasmus Handberg + """ + + mask1 = (img == value) + mask = np.zeros_like(img, dtype='bool') + + # Mask entire rows and columns which are only the value: + mask[np.all(mask1, axis=1), :] = True + mask[:, np.all(mask1, axis=0)] = True + + # Detect "uneven" edges column-wise in image: + a = np.argmin(mask1, axis=0) + b = np.argmin(np.flipud(mask1), axis=0) + for col in range(img.shape[1]): + if mask1[0, col]: + mask[:a[col], col] = True + if mask1[-1, col]: + mask[-b[col]:, col] = True + + # Detect "uneven" edges row-wise in image: + a = np.argmin(mask1, axis=1) + b = np.argmin(np.fliplr(mask1), axis=1) + for row in range(img.shape[0]): + if mask1[row, 0]: + mask[row, :a[row]] = True + if mask1[row, -1]: + mask[row, -b[row]:] = True + + return mask + + def apply_edge_mask(self, y: Tuple[int] = None, x: Tuple[int] = None, apply_existing_mask_first: bool = False): + """ + Masks given rows and columns of image but will replace the current mask! Set apply_existing_mask_first to True + if the current mask should be kept. + :param y: Tuple[int] of rows to mask + :param x: Tuple[int] of columns to mask + :param apply_existing_mask_first: Whether to apply the existing mask to image first, before overwriting mask. + :return: None + """ + if y is None and x is None: + logger.debug("(y,x) was None when applying edge mask. Edge was not actually masked.") + + if apply_existing_mask_first: + self.create_masked_image() + + if y is not None: + self.set_edge_rows_to_value(y=y) + + if x is not None: + self.set_edge_columns_to_value(x=x) + + self.mask = self.get_edge_mask(self.image) + self.create_masked_image() + class AbstractInstrument(ABC): peakmax: int = None @@ -149,7 +224,7 @@ def __init__(self, image: FlowsImage): def get_obstime(self): obstime = Time(self.image.header['DATE-OBS'], format='isot', scale='utc', - location=self.image.site['EarthLocation']) + location=self.image.site['EarthLocation']) obstime += 0.5 * self.image.exptime * u.second # Make time centre of exposure return obstime @@ -186,8 +261,12 @@ class ALFOSC(Instrument): siteid = 5 def get_obstime(self): - return Time(self.image.header['DATE-AVG'], format='isot', scale='utc', - location=self.image.site['EarthLocation']) + return Time( + self.image.header['DATE-AVG'], + format='isot', + scale='utc', + location=self.image.site['EarthLocation'] + ) def get_photfilter(self): # Sometimes data from NOT does not have the FILTER keyword, @@ -481,6 +560,7 @@ def get_photfilter(self): }.get(hdr['FILTER'], hdr['FILTER']) return photfilter + class TJO(Instrument): siteid = 22 @@ -570,7 +650,6 @@ def new_load_image(filename: str, target_coord: coords.SkyCoord = None): .. codeauthor:: Emir Karamehmetoglu .. codeauthor:: Rasmus Handberg """ - logger = logging.getLogger(__name__) ext = 0 # Default extension in HDUList, individual instruments may override this. mask = None # Instrument can override, default is to only mask all non-finite values, override is additive. @@ -610,7 +689,8 @@ def new_load_image(filename: str, target_coord: coords.SkyCoord = None): break # NOT - ALFOSC - elif telescope == 'NOT' and instrument in ('ALFOSC FASU','ALFOSC_FASU') and hdr.get('OBS_MODE','').lower() == 'imaging': + elif telescope == 'NOT' and instrument in ('ALFOSC FASU', 'ALFOSC_FASU') and hdr.get('OBS_MODE', + '').lower() == 'imaging': instrument_name = 'ALFOSC' break @@ -651,7 +731,8 @@ def new_load_image(filename: str, target_coord: coords.SkyCoord = None): instrument_name = 'Sofi' break - elif telescope == 'ESO-NTT' and instrument == 'EFOSC' and (origin == 'ESO' or origin.startswith('NOAO-IRAF')): + elif telescope == 'ESO-NTT' and instrument == 'EFOSC' and ( + origin == 'ESO' or origin.startswith('NOAO-IRAF')): instrument_name = 'EFOSC' break From d8f36ec3f98d4cd19e29929da149445bdf7ddbd2 Mon Sep 17 00:00:00 2001 From: Emir Date: Tue, 21 Dec 2021 09:48:46 +0100 Subject: [PATCH 04/10] change pattern matching implementation --- flows/load_image.py | 73 ++++++++++++-------- tests/input/ADP.2021-11-08T12_34_46.383.fits | 3 + tests/input/ADP.2021-11-08T12_46_22.832.fits | 3 + 3 files changed, 50 insertions(+), 29 deletions(-) create mode 100644 tests/input/ADP.2021-11-08T12_34_46.383.fits create mode 100644 tests/input/ADP.2021-11-08T12_46_22.832.fits diff --git a/flows/load_image.py b/flows/load_image.py index e605d2c..316ec90 100644 --- a/flows/load_image.py +++ b/flows/load_image.py @@ -160,19 +160,16 @@ def get_obstime(self): def get_photfilter(self): pass + @abstractmethod + def process_image(self): + class Instrument(AbstractInstrument): peakmax: int = None siteid: int = None - def __init__(self, image: FlowsImage): + def __init__(self, image: FlowsImage=None): self.image = image - self.image.peakmax = self.peakmax - self.image.site = self.get_site() - self.image.exptime = self.get_exptime() - self.image.obstime = self.get_obstime() - self.image.photfilter = self.get_photfilter() - self.image.create_masked_image() def get_site(self): if self.siteid is not None: @@ -192,6 +189,23 @@ def get_obstime(self): def get_photfilter(self): return self.image.header['FILTER'] + def _get_clean_image(self): + self.image.peakmax = self.peakmax + self.image.site = self.get_site() + self.image.exptime = self.get_exptime() + self.image.obstime = self.get_obstime() + self.image.photfilter = self.get_photfilter() + self.image.create_masked_image() + + def process_image(self, image: FlowsImage=None): + """Process existing or new image.""" + if image is not None: + self.image = image + if self.image is None: + raise AttributeError('No FlowsImage to be processed. Self.image was None') + + self._get_clean_image() + class LCOGT(Instrument): peakmax: int = 60000 @@ -218,9 +232,10 @@ def get_photfilter(self): class HAWKI(Instrument): siteid = 2 # Hard-coded the siteid for ESO Paranal, VLT, UT4 - def __init__(self, image: FlowsImage): + def __init__(self, image: FlowsImage=None): super().__init__(image) - self.get_obtype() + if self.image is not None: + self.get_obtype() def get_obstime(self): obstime = Time(self.image.header['DATE-OBS'], format='isot', scale='utc', @@ -665,7 +680,7 @@ def new_load_image(filename: str, target_coord: coords.SkyCoord = None): while instrument_name is None: # LCOGT if origin == "LCOGT": - instrument_name = 'LCOGT' + instrument_name = LCOGT() if 'BPM' in hdul: mask = np.asarray(hdul['BPM'].data, dtype='bool') else: @@ -674,7 +689,7 @@ def new_load_image(filename: str, target_coord: coords.SkyCoord = None): # HAWKI elif origin == 'ESO-PARANAL' and telescope == 'ESO-VLT-U4' and instrument == 'HAWKI' and hdr.get( 'PRODCATG') == 'SCIENCE.MEFIMAGE': - instrument_name = 'HAWKI' + instrument_name = HAWKI() if target_coord is None: raise ValueError("TARGET_COORD is needed for HAWKI images to find the correct extension") @@ -691,79 +706,79 @@ def new_load_image(filename: str, target_coord: coords.SkyCoord = None): # NOT - ALFOSC elif telescope == 'NOT' and instrument in ('ALFOSC FASU', 'ALFOSC_FASU') and hdr.get('OBS_MODE', '').lower() == 'imaging': - instrument_name = 'ALFOSC' + instrument_name = ALFOSC() break elif telescope == 'NOT' and instrument == 'NOTCAM' and hdr.get('OBS_MODE', '').lower() == 'imaging': - instrument_name = 'NOTCAM' + instrument_name = NOTCAM() break elif hdr.get('FPA.TELESCOPE') == 'PS1' and hdr.get('FPA.INSTRUMENT') == 'GPC1': - instrument_name = 'PS1' + instrument_name = PS1() break elif telescope == 'Liverpool Telescope': - instrument_name = 'Liverpool' + instrument_name = Liverpool() break elif telescope == 'CA 3.5m' and instrument == 'Omega2000': - instrument_name = 'Omega2000' + instrument_name = Omega2000() break elif telescope == 'SWO' and hdr.get('SITENAME') == 'LCO': - instrument_name = 'Swope' + instrument_name = Swope() break elif telescope == 'DUP' and hdr.get('SITENAME') == 'LCO' and instrument == 'Direct/SITe2K-1': - instrument_name = 'Dupont' + instrument_name = Dupont() break elif telescope == 'DUP' and instrument == 'RetroCam': - instrument_name = 'RetroCam' + instrument_name = RetroCam() break elif telescope == 'Baade' and hdr.get('SITENAME') == 'LCO' and instrument == 'FourStar': - instrument_name = 'Baade' + instrument_name = Baade() break elif instrument == 'SOFI' and telescope in ('ESO-NTT', 'other') and ( origin == 'ESO' or origin.startswith('NOAO-IRAF')): - instrument_name = 'Sofi' + instrument_name = Sofi() break elif telescope == 'ESO-NTT' and instrument == 'EFOSC' and ( origin == 'ESO' or origin.startswith('NOAO-IRAF')): - instrument_name = 'EFOSC' + instrument_name = EFOSC() break elif telescope == 'SAI-2.5' and instrument == 'ASTRONIRCAM': - instrument_name = 'AstroNIRCam' + instrument_name = AstroNIRCam() break elif instrument == 'OMEGACAM' and (origin == 'ESO' or origin.startswith('NOAO-IRAF')): - instrument_name = 'OmegaCam' + instrument_name = OmegaCam() break elif instrument == 'ANDICAM-CCD' and hdr.get('OBSERVAT') == 'CTIO': - instrument_name = 'AndiCam' + instrument_name = AndiCam() break elif telescope == '1.3m PAIRITEL' and instrument == '2MASS Survey cam': - instrument_name = 'PairTel' + instrument_name = PairTel() break elif (origin == 'OAdM' or origin.startswith('NOAO-IRAF')) and telescope == 'TJO' and instrument in ( 'MEIA3', 'MEIA2'): - instrument_name = 'TJO' + instrument_name = TJO() break else: raise RuntimeError("Could not determine origin of image") image = FlowsImage(image=np.asarray(hdul[ext].data, dtype='float64'), header=hdr, mask=mask) - ins = instruments[instrument_name](image) + instrument_name.process_image(image) - return ins.image + return instrument_name.image # get image and wcs solution diff --git a/tests/input/ADP.2021-11-08T12_34_46.383.fits b/tests/input/ADP.2021-11-08T12_34_46.383.fits new file mode 100644 index 0000000..2e642f0 --- /dev/null +++ b/tests/input/ADP.2021-11-08T12_34_46.383.fits @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:180dc91ea06017e38c9af76cd8cf83a1aefc0e018932ec56296566bc2a638f29 +size 71887680 diff --git a/tests/input/ADP.2021-11-08T12_46_22.832.fits b/tests/input/ADP.2021-11-08T12_46_22.832.fits new file mode 100644 index 0000000..2e3a783 --- /dev/null +++ b/tests/input/ADP.2021-11-08T12_46_22.832.fits @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:e223925768ebde91ab5bc76cace959c1462b98a683707bc94f9f036b523b2559 +size 73353600 From c61852a286f557c675b8ae13fa30c0b287c5b4b2 Mon Sep 17 00:00:00 2001 From: Emir Date: Mon, 14 Feb 2022 17:40:30 +0100 Subject: [PATCH 05/10] finish v1 OO rewrite of load_img with tests passing --- flows/load_image.py | 1876 ++++++++++++++++++-------------------- tests/test_load_image.py | 23 +- 2 files changed, 881 insertions(+), 1018 deletions(-) diff --git a/flows/load_image.py b/flows/load_image.py index 316ec90..c08f250 100644 --- a/flows/load_image.py +++ b/flows/load_image.py @@ -18,7 +18,7 @@ from typing import Tuple from flows import api -from dataclasses import dataclass, field +from dataclasses import dataclass # , field import typing from abc import ABC, abstractmethod @@ -28,1120 +28,982 @@ # -------------------------------------------------------------------------------------------------- @dataclass class FlowsImage: - image: np.ndarray - header: typing.Dict - mask: np.ndarray = None - peakmax: float = None - - def __post_init__(self): - self.shape = self.image.shape - self.wcs = self.create_wcs() - # Make empty mask - if not self.mask: - self.mask = np.zeros_like(self.image, dtype='bool') - self.check_finite() - - def check_finite(self): - self.mask |= ~np.isfinite(self.image) - - def update_mask(self, mask): - self.mask = mask - self.check_finite() - - def create_wcs(self) -> WCS: - with warnings.catch_warnings(): - warnings.simplefilter('ignore', category=FITSFixedWarning) - return WCS(header=self.header, relax=True) - - def create_masked_image(self): - """Warning: this is destructive and will overwrite image data setting masked values to NaN""" - self.image[self.mask] = np.NaN - self.clean = np.ma.masked_array(data=self.image, mask=self.mask, copy=False) - - def set_edge_rows_to_value(self, y: Tuple[float] = None, value: typing.Union[int, float, np.NaN] = 0): - if y is None: - pass - for row in y: - self.image[row] = value - - def set_edge_columns_to_value(self, x: Tuple[float] = None, value: typing.Union[int, float, np.NaN] = 0): - if x is None: - pass - for col in x: - self.image[:, col] = value - - @staticmethod - def get_edge_mask(img: np.ndarray, value: typing.Union[int, float, np.NaN] = 0): - """ - Create boolean mask of given value near edge of image. - - Parameters: - img (ndarray): image with values for masking. - value (float): Value to detect near edge. Default=0. - - Returns: - ndarray: Pixel mask with given values on the edge of image. - - .. codeauthor:: Rasmus Handberg - """ - - mask1 = (img == value) - mask = np.zeros_like(img, dtype='bool') - - # Mask entire rows and columns which are only the value: - mask[np.all(mask1, axis=1), :] = True - mask[:, np.all(mask1, axis=0)] = True - - # Detect "uneven" edges column-wise in image: - a = np.argmin(mask1, axis=0) - b = np.argmin(np.flipud(mask1), axis=0) - for col in range(img.shape[1]): - if mask1[0, col]: - mask[:a[col], col] = True - if mask1[-1, col]: - mask[-b[col]:, col] = True - - # Detect "uneven" edges row-wise in image: - a = np.argmin(mask1, axis=1) - b = np.argmin(np.fliplr(mask1), axis=1) - for row in range(img.shape[0]): - if mask1[row, 0]: - mask[row, :a[row]] = True - if mask1[row, -1]: - mask[row, -b[row]:] = True - - return mask - - def apply_edge_mask(self, y: Tuple[int] = None, x: Tuple[int] = None, apply_existing_mask_first: bool = False): - """ - Masks given rows and columns of image but will replace the current mask! Set apply_existing_mask_first to True - if the current mask should be kept. - :param y: Tuple[int] of rows to mask - :param x: Tuple[int] of columns to mask - :param apply_existing_mask_first: Whether to apply the existing mask to image first, before overwriting mask. - :return: None - """ - if y is None and x is None: - logger.debug("(y,x) was None when applying edge mask. Edge was not actually masked.") - - if apply_existing_mask_first: - self.create_masked_image() - - if y is not None: - self.set_edge_rows_to_value(y=y) - - if x is not None: - self.set_edge_columns_to_value(x=x) - - self.mask = self.get_edge_mask(self.image) - self.create_masked_image() + image: np.ndarray + header: typing.Dict + mask: np.ndarray = None + peakmax: float = None + exptime: float = None + clean: np.ma.MaskedArray = None + + def __post_init__(self): + self.shape = self.image.shape + self.wcs = self.create_wcs() + # Make empty mask + if not self.mask: + self.mask = np.zeros_like(self.image, dtype='bool') + self.check_finite() + + def check_finite(self): + self.mask |= ~np.isfinite(self.image) + + def update_mask(self, mask): + self.mask = mask + self.check_finite() + + def create_wcs(self) -> WCS: + with warnings.catch_warnings(): + warnings.simplefilter('ignore', category=FITSFixedWarning) + return WCS(header=self.header, relax=True) + + def create_masked_image(self): + """Warning: this is destructive and will overwrite image data setting masked values to NaN""" + self.image[self.mask] = np.NaN + self.clean = np.ma.masked_array(data=self.image, mask=self.mask, copy=False) + + def set_edge_rows_to_value(self, y: Tuple[float] = None, value: typing.Union[int, float, np.float64] = 0): + if y is None: + pass + for row in y: + self.image[row] = value + + def set_edge_columns_to_value(self, x: Tuple[float] = None, value: typing.Union[int, float, np.float64] = 0): + if x is None: + pass + for col in x: + self.image[:, col] = value + + @staticmethod + def get_edge_mask(img: np.ndarray, value: typing.Union[int, float, np.float64] = 0): + """ + Create boolean mask of given value near edge of image. + + Parameters: + img (ndarray): image with values for masking. + value (float): Value to detect near edge. Default=0. + + Returns: + ndarray: Pixel mask with given values on the edge of image. + + .. codeauthor:: Rasmus Handberg + """ + + mask1 = (img == value) + mask = np.zeros_like(img, dtype='bool') + + # Mask entire rows and columns which are only the value: + mask[np.all(mask1, axis=1), :] = True + mask[:, np.all(mask1, axis=0)] = True + + # Detect "uneven" edges column-wise in image: + a = np.argmin(mask1, axis=0) + b = np.argmin(np.flipud(mask1), axis=0) + for col in range(img.shape[1]): + if mask1[0, col]: + mask[:a[col], col] = True + if mask1[-1, col]: + mask[-b[col]:, col] = True + + # Detect "uneven" edges row-wise in image: + a = np.argmin(mask1, axis=1) + b = np.argmin(np.fliplr(mask1), axis=1) + for row in range(img.shape[0]): + if mask1[row, 0]: + mask[row, :a[row]] = True + if mask1[row, -1]: + mask[row, -b[row]:] = True + + return mask + + def apply_edge_mask(self, y: Tuple[int] = None, x: Tuple[int] = None, apply_existing_mask_first: bool = False): + """ + Masks given rows and columns of image but will replace the current mask! Set apply_existing_mask_first to True + if the current mask should be kept. + :param y: Tuple[int] of rows to mask + :param x: Tuple[int] of columns to mask + :param apply_existing_mask_first: Whether to apply the existing mask to image first, before overwriting mask. + :return: None + """ + if y is None and x is None: + logger.debug("(y,x) was None when applying edge mask. Edge was not actually masked.") + + if apply_existing_mask_first: + self.create_masked_image() + + if y is not None: + self.set_edge_rows_to_value(y=y) + + if x is not None: + self.set_edge_columns_to_value(x=x) + + self.mask = self.get_edge_mask(self.image) + self.create_masked_image() class AbstractInstrument(ABC): - peakmax: int = None + peakmax: int = None - @abstractmethod - def __init__(self): - pass + @abstractmethod + def __init__(self): + pass - @abstractmethod - def get_site(self): - pass + @abstractmethod + def get_site(self): + pass - @abstractmethod - def get_exptime(self): - pass + @abstractmethod + def get_exptime(self): + pass - @abstractmethod - def get_obstime(self): - pass + @abstractmethod + def get_obstime(self): + pass - @abstractmethod - def get_photfilter(self): - pass + @abstractmethod + def get_photfilter(self): + pass - @abstractmethod - def process_image(self): + @abstractmethod + def process_image(self): + pass class Instrument(AbstractInstrument): - peakmax: int = None - siteid: int = None + peakmax: int = None + siteid: int = None - def __init__(self, image: FlowsImage=None): - self.image = image + def __init__(self, image: FlowsImage = None): + self.image = image - def get_site(self): - if self.siteid is not None: - return api.get_site(self.siteid) + def get_site(self): + if self.siteid is not None: + return api.get_site(self.siteid) - def get_exptime(self): - exptime = self.image.header.get('EXPTIME', None) - if exptime is None: - raise ValueError("Image exposure time could not be extracted") - return exptime + def get_exptime(self): + exptime = self.image.header.get('EXPTIME', None) + if exptime is None: + raise ValueError("Image exposure time could not be extracted") + return exptime - def get_obstime(self): - '''Default for JD, jd, utc.''' - return Time(self.image.header['JD'], format='jd', scale='utc', - location=self.image.site['EarthLocation']) + def get_obstime(self): + """Default for JD, jd, utc.""" + return Time(self.image.header['JD'], format='jd', scale='utc', location=self.image.site['EarthLocation']) - def get_photfilter(self): - return self.image.header['FILTER'] + def get_photfilter(self): + return self.image.header['FILTER'] - def _get_clean_image(self): - self.image.peakmax = self.peakmax - self.image.site = self.get_site() - self.image.exptime = self.get_exptime() - self.image.obstime = self.get_obstime() - self.image.photfilter = self.get_photfilter() - self.image.create_masked_image() + def _get_clean_image(self): + self.image.peakmax = self.peakmax + self.image.site = self.get_site() + self.image.exptime = self.get_exptime() + self.image.obstime = self.get_obstime() + self.image.photfilter = self.get_photfilter() + self.image.create_masked_image() - def process_image(self, image: FlowsImage=None): - """Process existing or new image.""" - if image is not None: - self.image = image - if self.image is None: - raise AttributeError('No FlowsImage to be processed. Self.image was None') + def process_image(self, image: FlowsImage = None): + """Process existing or new image.""" + if image is not None: + self.image = image + if self.image is None: + raise AttributeError('No FlowsImage to be processed. Self.image was None') - self._get_clean_image() + self._get_clean_image() class LCOGT(Instrument): - peakmax: int = 60000 + peakmax: int = 60000 - def get_site(self): - sites = api.sites.get_all_sites() - site_keywords = {s['site_keyword']: s for s in sites} - site = site_keywords.get(self.image.header['SITE'], None) - return site + def get_site(self): + sites = api.sites.get_all_sites() + site_keywords = {s['site_keyword']: s for s in sites} + site = site_keywords.get(self.image.header['SITE'], None) + return site - def get_obstime(self): - observatory = coords.EarthLocation.from_geodetic(lat=self.image.header['LATITUDE'], - lon=self.image.header['LONGITUD'], - height=self.image.header['HEIGHT']) - obstime = Time(self.image.header['DATE-OBS'], format='isot', scale='utc', location=observatory) - obstime += 0.5 * self.image.exptime * u.second # Make time centre of exposure - return obstime + def get_obstime(self): + observatory = coords.EarthLocation.from_geodetic(lat=self.image.header['LATITUDE'], + lon=self.image.header['LONGITUD'], + height=self.image.header['HEIGHT']) + obstime = Time(self.image.header['DATE-OBS'], format='isot', scale='utc', location=observatory) + obstime += 0.5 * self.image.exptime * u.second # Make time centre of exposure + return obstime - def get_photfilter(self): - photfilter = {'zs': 'zp'}.get(self.image.header['FILTER'], self.image.header['FILTER']) - return photfilter + def get_photfilter(self): + photfilter = {'zs': 'zp'}.get(self.image.header['FILTER'], self.image.header['FILTER']) + return photfilter class HAWKI(Instrument): - siteid = 2 # Hard-coded the siteid for ESO Paranal, VLT, UT4 + siteid = 2 # Hard-coded the siteid for ESO Paranal, VLT, UT4 - def __init__(self, image: FlowsImage=None): - super().__init__(image) - if self.image is not None: - self.get_obtype() + def __init__(self, image: FlowsImage = None): + super().__init__(image) + if self.image is not None: + self.get_obtype() - def get_obstime(self): - obstime = Time(self.image.header['DATE-OBS'], format='isot', scale='utc', - location=self.image.site['EarthLocation']) - obstime += 0.5 * self.image.exptime * u.second # Make time centre of exposure - return obstime + def get_obstime(self): + obstime = Time(self.image.header['DATE-OBS'], format='isot', scale='utc', + location=self.image.site['EarthLocation']) + obstime += 0.5 * self.image.exptime * u.second # Make time centre of exposure + return obstime - def get_obtype(self): - ob_type = self.image.header["HIERARCH ESO OCS DET1 IMGNAME"].split('_')[-1] - if "Auto" in ob_type: - self.image.ob_type = 'Autojitter' - elif "Fixed" in ob_type: - self.image.ob_type = 'FixedOffset' - else: - raise RuntimeError("Image OB Type not AutoJitter or FixedOffset") + def get_obtype(self): + ob_type = self.image.header["HIERARCH ESO OCS DET1 IMGNAME"].split('_')[-1] + if "Auto" in ob_type: + self.image.ob_type = 'Autojitter' + elif "Fixed" in ob_type: + self.image.ob_type = 'FixedOffset' + else: + raise RuntimeError("Image OB Type not AutoJitter or FixedOffset") def get_image_extension(hdul: fits.HDUList, target_coord: coords.SkyCoord = None, fallback_extension: int = None): - # For HAWKI multi-extension images we search the extensions for which one contains - # the target, Create Image from that extension. - target_radec = [[target_coord.icrs.ra.deg, target_coord.icrs.dec.deg]] - - for k in range(1, 5): - w = WCS(header=hdul[k].header, relax=True) - s = [hdul[k].header['NAXIS2'], hdul[k].header['NAXIS1']] - pix = w.all_world2pix(target_radec, 0).flatten() - if -0.5 <= pix[0] <= s[1] - 0.5 and -0.5 <= pix[1] <= s[0] - 0.5: - return k - if not fallback_extension is None: - return fallback_extension - else: - raise RuntimeError(f"Could not find image extension that target is on!") + # For HAWKI multi-extension images we search the extensions for which one contains + # the target, Create Image from that extension. + target_radec = [[target_coord.icrs.ra.deg, target_coord.icrs.dec.deg]] + + for k in range(1, 5): + w = WCS(header=hdul[k].header, relax=True) + s = [hdul[k].header['NAXIS2'], hdul[k].header['NAXIS1']] + pix = w.all_world2pix(target_radec, 0).flatten() + if -0.5 <= pix[0] <= s[1] - 0.5 and -0.5 <= pix[1] <= s[0] - 0.5: + return k + if fallback_extension is not None: + return fallback_extension + else: + raise RuntimeError(f"Could not find image extension that target is on!") class ALFOSC(Instrument): - # Obtained from http://www.not.iac.es/instruments/detectors/CCD14/LED-linearity/20181026-200-1x1.pdf - peakmax = 80000 # For ALFOSC D, 1x1, 200; the standard for SNe. - siteid = 5 - - def get_obstime(self): - return Time( - self.image.header['DATE-AVG'], - format='isot', - scale='utc', - location=self.image.site['EarthLocation'] - ) - - def get_photfilter(self): - # Sometimes data from NOT does not have the FILTER keyword, - # in which case we have to try to figure out which filter - # was used based on some of the other headers: - if 'FILTER' in self.image.header: - photfilter = { - 'B Bes': 'B', - 'V Bes': 'V', - 'R Bes': 'R', - 'g SDSS': 'gp', - 'r SDSS': 'rp', - 'i SDSS': 'ip', - 'i int': 'ip', # Interference filter - 'u SDSS': 'up', - 'z SDSS': 'zp' - }.get(self.image.header['FILTER'].replace('_', ' '), self.image.header['FILTER']) - else: - filters_used = [] - for check_headers in ('ALFLTNM', 'FAFLTNM', 'FBFLTNM'): - isopen = self.image.header.get(check_headers).strip().lower() != 'open' - if self.image.header.get(check_headers) and isopen: - filters_used.append(self.image.header.get(check_headers).strip()) - if len(filters_used) == 1: - photfilter = { - 'V_Bes 530_80': 'V', - 'R_Bes 650_130': 'R', - "g'_SDSS 480_145": 'gp', - "r'_SDSS 618_148": 'rp', - "i'_SDSS 771_171": 'ip', - 'i_int 797_157': 'ip', # Interference filter - "z'_SDSS 832_LP": 'zp' - }.get(filters_used[0].replace(' ', ' '), filters_used[0]) - else: - raise RuntimeError("Could not determine filter used.") - - return photfilter + # Obtained from http://www.not.iac.es/instruments/detectors/CCD14/LED-linearity/20181026-200-1x1.pdf + peakmax = 80000 # For ALFOSC D, 1x1, 200; the standard for SNe. + siteid = 5 + + def get_obstime(self): + return Time(self.image.header['DATE-AVG'], format='isot', scale='utc', + location=self.image.site['EarthLocation']) + + def get_photfilter(self): + # Sometimes data from NOT does not have the FILTER keyword, + # in which case we have to try to figure out which filter + # was used based on some of the other headers: + if 'FILTER' in self.image.header: + photfilter = {'B Bes': 'B', 'V Bes': 'V', 'R Bes': 'R', 'g SDSS': 'gp', 'r SDSS': 'rp', 'i SDSS': 'ip', + 'i int': 'ip', # Interference filter + 'u SDSS': 'up', 'z SDSS': 'zp'}.get(self.image.header['FILTER'].replace('_', ' '), + self.image.header['FILTER']) + else: + filters_used = [] + for check_headers in ('ALFLTNM', 'FAFLTNM', 'FBFLTNM'): + isopen = self.image.header.get(check_headers).strip().lower() != 'open' + if self.image.header.get(check_headers) and isopen: + filters_used.append(self.image.header.get(check_headers).strip()) + if len(filters_used) == 1: + photfilter = {'B_Bes 440_100': 'B', 'V_Bes 530_80': 'V', 'R_Bes 650_130': 'R', "g'_SDSS 480_145": 'gp', + "r'_SDSS 618_148": 'rp', "i'_SDSS 771_171": 'ip', 'i_int 797_157': 'ip', + # Interference filter + "z'_SDSS 832_LP": 'zp'}.get(filters_used[0].replace(' ', ' '), filters_used[0]) + else: + raise RuntimeError("Could not determine filter used.") + + return photfilter class NOTCAM(Instrument): - siteid = 5 - - def get_obstime(self): - return Time(self.image.header['DATE-AVG'], format='isot', scale='utc', - location=self.image.site['EarthLocation']) - - def get_photfilter(self): - # Does NOTCAM data sometimes contain a FILTER header? - # if not we have to try to figure out which filter - # was used based on some of the other headers: - if 'FILTER' in self.image.header: - raise RuntimeError("NOTCAM: Filter keyword defined") - filters_used = [] - for check_headers in ('NCFLTNM1', 'NCFLTNM2'): - isopen = self.image.header.get(check_headers).strip().lower() != 'open' - if self.image.header.get(check_headers) and isopen: - filters_used.append(self.image.header.get(check_headers).strip()) - if len(filters_used) == 1: - photfilter = {}.get(filters_used[0], filters_used[0]) - else: - raise RuntimeError("Could not determine filter used.") - return photfilter + siteid = 5 + + def get_obstime(self): + return Time(self.image.header['DATE-AVG'], format='isot', scale='utc', + location=self.image.site['EarthLocation']) + + def get_photfilter(self): + # Does NOTCAM data sometimes contain a FILTER header? + # if not we have to try to figure out which filter + # was used based on some of the other headers: + if 'FILTER' in self.image.header: + raise RuntimeError("NOTCAM: Filter keyword defined") + filters_used = [] + for check_headers in ('NCFLTNM1', 'NCFLTNM2'): + isopen = self.image.header.get(check_headers).strip().lower() != 'open' + if self.image.header.get(check_headers) and isopen: + filters_used.append(self.image.header.get(check_headers).strip()) + if len(filters_used) == 1: + photfilter = {'Ks': 'K'}.get(filters_used[0], filters_used[0]) + else: + raise RuntimeError("Could not determine filter used.") + return photfilter class PS1(Instrument): - siteid = 6 + siteid = 6 - def get_obstime(self): - return Time(self.image.header['MJD-OBS'], format='mjd', scale='utc', - location=self.image.site['EarthLocation']) + def get_obstime(self): + return Time(self.image.header['MJD-OBS'], format='mjd', scale='utc', location=self.image.site['EarthLocation']) - def get_photfilter(self): - photfilter = { - 'g.00000': 'gp', - 'r.00000': 'rp', - 'i.00000': 'ip', - 'z.00000': 'zp' - }.get(self.image.header['FPA.FILTER'], self.image.header['FPA.FILTER']) - return photfilter + def get_photfilter(self): + photfilter = {'g.00000': 'gp', 'r.00000': 'rp', 'i.00000': 'ip', 'z.00000': 'zp'}.get( + self.image.header['FPA.FILTER'], self.image.header['FPA.FILTER']) + return photfilter class Liverpool(Instrument): - siteid = 8 - - def get_obstime(self): - obstime = Time(self.image.header['DATE-OBS'], format='isot', scale='utc', - location=self.image.site['EarthLocation']) - obstime += 0.5 * self.image.exptime * u.second # Make time centre of exposure - return obstime - - def get_photfilter(self): - photfilter = { - 'Bessel-B': 'B', - 'Bessell-B': 'B', - 'Bessel-V': 'V', - 'Bessell-V': 'V', - 'SDSS-U': 'up', - 'SDSS-G': 'gp', - 'SDSS-R': 'rp', - 'SDSS-I': 'ip', - 'SDSS-Z': 'zp' - }.get(self.image.header['FILTER1'], self.image.header['FILTER1']) - return photfilter + siteid = 8 + + def get_obstime(self): + obstime = Time(self.image.header['DATE-OBS'], format='isot', scale='utc', + location=self.image.site['EarthLocation']) + obstime += 0.5 * self.image.exptime * u.second # Make time centre of exposure + return obstime + + def get_photfilter(self): + photfilter = {'Bessel-B': 'B', 'Bessell-B': 'B', 'Bessel-V': 'V', 'Bessell-V': 'V', 'SDSS-U': 'up', + 'SDSS-G': 'gp', 'SDSS-R': 'rp', 'SDSS-I': 'ip', 'SDSS-Z': 'zp'}.get(self.image.header['FILTER1'], + self.image.header['FILTER1']) + return photfilter class Omega2000(Instrument): - siteid = 9 + siteid = 9 - def get_obstime(self): - obstime = Time(self.image.header['MJD-OBS'], format='mjd', scale='utc', - location=self.image.site['EarthLocation']) - obstime += 0.5 * self.image.exptime * u.second - return obstime + def get_obstime(self): + obstime = Time(self.image.header['MJD-OBS'], format='mjd', scale='utc', + location=self.image.site['EarthLocation']) + obstime += 0.5 * self.image.exptime * u.second + return obstime class Swope(Instrument): - siteid = 10 + siteid = 10 - def get_photfilter(self): - photfilter = { - 'u': 'up', - 'g': 'gp', - 'r': 'rp', - 'i': 'ip', - }.get(self.image.header['FILTER'], self.image.header['FILTER']) - return photfilter + def get_photfilter(self): + photfilter = {'u': 'up', 'g': 'gp', 'r': 'rp', 'i': 'ip', }.get(self.image.header['FILTER'], + self.image.header['FILTER']) + return photfilter class Dupont(Instrument): - siteid = 14 + siteid = 14 - def get_photfilter(self): - photfilter = { - 'u': 'up', - 'g': 'gp', - 'r': 'rp', - 'i': 'ip', - }.get(self.image.header['FILTER'], self.image.header['FILTER']) - return photfilter + def get_photfilter(self): + photfilter = {'u': 'up', 'g': 'gp', 'r': 'rp', 'i': 'ip', }.get(self.image.header['FILTER'], + self.image.header['FILTER']) + return photfilter class RetroCam(Instrument): - siteid = 16 + siteid = 16 - def get_photfilter(self): - photfilter = { - 'Yc': 'Y', - 'Hc': 'H', - 'Jo': 'J', - }.get(self.image.header['FILTER'], self.image.header['FILTER']) - return photfilter + def get_photfilter(self): + photfilter = {'Yc': 'Y', 'Hc': 'H', 'Jo': 'J', }.get(self.image.header['FILTER'], self.image.header['FILTER']) + return photfilter class Baade(Instrument): - siteid = 11 + siteid = 11 - def get_exptime(self): - exptime = super().get_exptime() - exptime *= int(self.image.header['NCOMBINE']) # EXPTIME is only for a single exposure - return exptime + def get_exptime(self): + exptime = super().get_exptime() + exptime *= int(self.image.header['NCOMBINE']) # EXPTIME is only for a single exposure + return exptime - def get_photfilter(self): - photfilter = { - 'Ks': 'K', - 'J1': 'Y', - }.get(self.image.header['FILTER'], self.image.header['FILTER']) - return photfilter + def get_photfilter(self): + photfilter = {'Ks': 'K', 'J1': 'Y', }.get(self.image.header['FILTER'], self.image.header['FILTER']) + return photfilter class Sofi(Instrument): - siteid = 12 - - def get_obstime(self): - if 'TMID' in self.image.header: - obstime = Time(self.image.header['TMID'], format='mjd', scale='utc', - location=self.image.site['EarthLocation']) - else: - obstime = Time(self.image.header['MJD-OBS'], format='mjd', scale='utc', - location=self.image.site['EarthLocation']) - obstime += 0.5 * self.image.exptime * u.second # Make time centre of exposure - return obstime - - def get_photfilter(self): - hdr = self.image.header - photfilter_translate = { - 'Ks': 'K' - } - if 'FILTER' in hdr: - photfilter = photfilter_translate.get(hdr['FILTER'], hdr['FILTER']) - else: - filters_used = [] - for check_headers in ('ESO INS FILT1 ID', 'ESO INS FILT2 ID'): - if hdr.get(check_headers) and hdr.get(check_headers).strip().lower() != 'open': - filters_used.append(hdr.get(check_headers).strip()) - if len(filters_used) == 1: - photfilter = photfilter_translate.get(filters_used[0], filters_used[0]) - else: - raise RuntimeError("Could not determine filter used.") - return photfilter + siteid = 12 + + def get_obstime(self): + if 'TMID' in self.image.header: + obstime = Time(self.image.header['TMID'], format='mjd', scale='utc', + location=self.image.site['EarthLocation']) + else: + obstime = Time(self.image.header['MJD-OBS'], format='mjd', scale='utc', + location=self.image.site['EarthLocation']) + obstime += 0.5 * self.image.exptime * u.second # Make time centre of exposure + return obstime + + def get_photfilter(self): + hdr = self.image.header + photfilter_translate = {'Ks': 'K'} + if 'FILTER' in hdr: + photfilter = photfilter_translate.get(hdr['FILTER'], hdr['FILTER']) + else: + filters_used = [] + for check_headers in ('ESO INS FILT1 ID', 'ESO INS FILT2 ID'): + if hdr.get(check_headers) and hdr.get(check_headers).strip().lower() != 'open': + filters_used.append(hdr.get(check_headers).strip()) + if len(filters_used) == 1: + photfilter = photfilter_translate.get(filters_used[0], filters_used[0]) + else: + raise RuntimeError("Could not determine filter used.") + return photfilter class EFOSC(Instrument): - siteid = 15 - - def get_obstime(self): - obstime = Time(self.image.header['DATE-OBS'], format='isot', scale='utc', - location=self.image.site['EarthLocation']) - obstime += 0.5 * self.image.exptime * u.second # Make time centre of exposure - return obstime - - def get_photfilter(self): - hdr = self.image.header - photfilter = { - 'g782': 'gp', - 'r784': 'rp', - 'i705': 'ip', - 'B639': 'B', - 'V641': 'V' - }.get(hdr['FILTER'], hdr['FILTER']) - return photfilter + siteid = 15 + + def get_obstime(self): + obstime = Time(self.image.header['DATE-OBS'], format='isot', scale='utc', + location=self.image.site['EarthLocation']) + obstime += 0.5 * self.image.exptime * u.second # Make time centre of exposure + return obstime + + def get_photfilter(self): + hdr = self.image.header + photfilter = {'g782': 'gp', 'r784': 'rp', 'i705': 'ip', 'B639': 'B', 'V641': 'V'}.get(hdr['FILTER'], + hdr['FILTER']) + return photfilter class AstroNIRCam(Instrument): - siteid = 13 + siteid = 13 - def get_exptime(self): - return self.image.header.get('FULL_EXP', super().get_exptime()) + def get_exptime(self): + return self.image.header.get('FULL_EXP', super().get_exptime()) - def get_obstime(self): - hdr = self.image.header - if 'MIDPOINT' in hdr: - obstime = Time(hdr['MIDPOINT'], format='isot', scale='utc', location=self.image.site['EarthLocation']) - else: - obstime = Time(hdr['MJD-AVG'], format='mjd', scale='utc', location=self.image.site['EarthLocation']) - return obstime + def get_obstime(self): + hdr = self.image.header + if 'MIDPOINT' in hdr: + obstime = Time(hdr['MIDPOINT'], format='isot', scale='utc', location=self.image.site['EarthLocation']) + else: + obstime = Time(hdr['MJD-AVG'], format='mjd', scale='utc', location=self.image.site['EarthLocation']) + return obstime - def get_photfilter(self): - hdr = self.image.header - photfilter = { - 'H_Open': 'H', - 'K_Open': 'K', - }.get(hdr['FILTER'], hdr['FILTER']) - return photfilter + def get_photfilter(self): + hdr = self.image.header + photfilter = {'H_Open': 'H', 'K_Open': 'K', }.get(hdr['FILTER'], hdr['FILTER']) + return photfilter class OmegaCam(Instrument): - siteid = 18 # Hard-coded the siteid for ESO VLT Survey telescope + siteid = 18 # Hard-coded the siteid for ESO VLT Survey telescope - def get_obstime(self): - obstime = Time(self.image.header['MJD-OBS'], format='mjd', scale='utc', - location=self.image.site['EarthLocation']) - obstime += 0.5 * self.image.exptime * u.second # Make time centre of exposure - return obstime + def get_obstime(self): + obstime = Time(self.image.header['MJD-OBS'], format='mjd', scale='utc', + location=self.image.site['EarthLocation']) + obstime += 0.5 * self.image.exptime * u.second # Make time centre of exposure + return obstime - def get_photfilter(self): - hdr = self.image.header - photfilter = { - 'i_SDSS': 'ip' - }.get(hdr['ESO INS FILT1 NAME'], hdr['ESO INS FILT1 NAME']) - return photfilter + def get_photfilter(self): + hdr = self.image.header + photfilter = {'i_SDSS': 'ip'}.get(hdr['ESO INS FILT1 NAME'], hdr['ESO INS FILT1 NAME']) + return photfilter class AndiCam(Instrument): - siteid = 20 # Hard-coded the siteid for ANDICAM at Cerro Tololo Interamerican Observatory (CTIO) + siteid = 20 # Hard-coded the siteid for ANDICAM at Cerro Tololo Interamerican Observatory (CTIO) - def get_obstime(self): - obstime = super().get_obstime() - obstime += 0.5 * self.image.exptime * u.second - return obstime + def get_obstime(self): + obstime = super().get_obstime() + obstime += 0.5 * self.image.exptime * u.second + return obstime - def get_photfilter(self): - return self.image.header['CCDFLTID'] + def get_photfilter(self): + return self.image.header['CCDFLTID'] class PairTel(Instrument): - siteid = 21 + siteid = 21 - def get_obstime(self): - hdr = self.image.header - time_start = Time(hdr['STRT_CPU'], format='iso', scale='utc', location=self.image.site['EarthLocation']) - time_stop = Time(hdr['STOP_CPU'], format='iso', scale='utc', location=self.image.site['EarthLocation']) - obstime = time_start + 0.5 * (time_stop - time_start) - return obstime + def get_obstime(self): + hdr = self.image.header + time_start = Time(hdr['STRT_CPU'], format='iso', scale='utc', location=self.image.site['EarthLocation']) + time_stop = Time(hdr['STOP_CPU'], format='iso', scale='utc', location=self.image.site['EarthLocation']) + obstime = time_start + 0.5 * (time_stop - time_start) + return obstime - def get_photfilter(self): - hdr = self.image.header - photfilter = { - 'j': 'J', - 'h': 'H', - 'k': 'K', - }.get(hdr['FILTER'], hdr['FILTER']) - return photfilter + def get_photfilter(self): + hdr = self.image.header + photfilter = {'j': 'J', 'h': 'H', 'k': 'K', }.get(hdr['FILTER'], hdr['FILTER']) + return photfilter class TJO(Instrument): - siteid = 22 - - def get_obstime(self): - obstime = super().get_obstime() - obstime += 0.5 * self.image.exptime * u.second - return obstime - - -instruments = { - 'LCOGT': LCOGT, - 'HAWKI': HAWKI, - 'ALFOSC': ALFOSC, - 'NOTCAM': NOTCAM, - 'PS1': PS1, - 'Liverpool': Liverpool, - 'Omega2000': Omega2000, - 'Swope': Swope, - 'Dupont': Dupont, - 'Retrocam': RetroCam, - 'Baade': Baade, - 'Sofi': Sofi, - 'EFOSC': EFOSC, - 'AstroNIRCam': AstroNIRCam, - 'OmegaCam': OmegaCam, - 'AndiCam': AndiCam, - 'PairTel': PairTel, - 'TJO': TJO -} + siteid = 22 + def get_obstime(self): + obstime = super().get_obstime() + obstime += 0.5 * self.image.exptime * u.second + return obstime -def edge_mask(img, value=0): - """ - Create boolean mask of given value near edge of image. - - Parameters: - img (ndarray): Image of - value (float): Value to detect near edge. Default=0. - - Returns: - ndarray: Pixel mask with given values on the edge of image. - - .. codeauthor:: Rasmus Handberg - """ - - mask1 = (img == value) - mask = np.zeros_like(img, dtype='bool') - - # Mask entire rows and columns which are only the value: - mask[np.all(mask1, axis=1), :] = True - mask[:, np.all(mask1, axis=0)] = True - - # Detect "uneven" edges column-wise in image: - a = np.argmin(mask1, axis=0) - b = np.argmin(np.flipud(mask1), axis=0) - for col in range(img.shape[1]): - if mask1[0, col]: - mask[:a[col], col] = True - if mask1[-1, col]: - mask[-b[col]:, col] = True - - # Detect "uneven" edges row-wise in image: - a = np.argmin(mask1, axis=1) - b = np.argmin(np.fliplr(mask1), axis=1) - for row in range(img.shape[0]): - if mask1[row, 0]: - mask[row, :a[row]] = True - if mask1[row, -1]: - mask[row, -b[row]:] = True - - return mask - - -def new_load_image(filename: str, target_coord: coords.SkyCoord = None): - """ - Load FITS image using FlowsImage class and Instrument Classes. - - Parameters: - filename (str): Path to FITS file to be loaded. - target_coord (:class:`astropy.coordinates.SkyCoord`): Coordinates of target. - Only used for HAWKI images to determine which image extension to load, - for all other images it is ignored. - - Returns: - FlowsImage: instance of FlowsImage with valuues populated based on instrument. - - .. codeauthor:: Emir Karamehmetoglu - .. codeauthor:: Rasmus Handberg - """ - ext = 0 # Default extension in HDUList, individual instruments may override this. - mask = None # Instrument can override, default is to only mask all non-finite values, override is additive. - - # Read fits image, Structural Pattern Match to specific instrument. - with fits.open(filename, mode='readonly') as hdul: - hdr = hdul[ext].header - origin = hdr.get('ORIGIN', '') - telescope = hdr.get('TELESCOP', '') - instrument = hdr.get('INSTRUME', '') - - instrument_name = None - # Pattern matching begins here, ideally we use 3.10 pattern matching, or a dictionary lookup. - while instrument_name is None: - # LCOGT - if origin == "LCOGT": - instrument_name = LCOGT() - if 'BPM' in hdul: - mask = np.asarray(hdul['BPM'].data, dtype='bool') - else: - logger.warning('LCOGT image does not contain bad pixel map. Not applying mask.') - break - # HAWKI - elif origin == 'ESO-PARANAL' and telescope == 'ESO-VLT-U4' and instrument == 'HAWKI' and hdr.get( - 'PRODCATG') == 'SCIENCE.MEFIMAGE': - instrument_name = HAWKI() - - if target_coord is None: - raise ValueError("TARGET_COORD is needed for HAWKI images to find the correct extension") - if not isinstance(target_coord, coords.SkyCoord): - if len(target_coord) == 2: - target_coord = coords.SkyCoord(ra=target_coord[0] * u.deg, dec=target_coord[1] * u.deg, - frame='icrs') - else: - raise ValueError("TARGET_COORD is needed for HAWKI images to find the correct extension") - ext = get_image_extension(hdul, target_coord) # Find the one with the SN in it. - hdr = hdul[ext].header + hdul[0].header - break - - # NOT - ALFOSC - elif telescope == 'NOT' and instrument in ('ALFOSC FASU', 'ALFOSC_FASU') and hdr.get('OBS_MODE', - '').lower() == 'imaging': - instrument_name = ALFOSC() - break - - elif telescope == 'NOT' and instrument == 'NOTCAM' and hdr.get('OBS_MODE', '').lower() == 'imaging': - instrument_name = NOTCAM() - break - - elif hdr.get('FPA.TELESCOPE') == 'PS1' and hdr.get('FPA.INSTRUMENT') == 'GPC1': - instrument_name = PS1() - break - - elif telescope == 'Liverpool Telescope': - instrument_name = Liverpool() - break - - elif telescope == 'CA 3.5m' and instrument == 'Omega2000': - instrument_name = Omega2000() - break - - elif telescope == 'SWO' and hdr.get('SITENAME') == 'LCO': - instrument_name = Swope() - break - - elif telescope == 'DUP' and hdr.get('SITENAME') == 'LCO' and instrument == 'Direct/SITe2K-1': - instrument_name = Dupont() - break - - elif telescope == 'DUP' and instrument == 'RetroCam': - instrument_name = RetroCam() - break - - elif telescope == 'Baade' and hdr.get('SITENAME') == 'LCO' and instrument == 'FourStar': - instrument_name = Baade() - break - - elif instrument == 'SOFI' and telescope in ('ESO-NTT', 'other') and ( - origin == 'ESO' or origin.startswith('NOAO-IRAF')): - instrument_name = Sofi() - break - - elif telescope == 'ESO-NTT' and instrument == 'EFOSC' and ( - origin == 'ESO' or origin.startswith('NOAO-IRAF')): - instrument_name = EFOSC() - break - - elif telescope == 'SAI-2.5' and instrument == 'ASTRONIRCAM': - instrument_name = AstroNIRCam() - break - - elif instrument == 'OMEGACAM' and (origin == 'ESO' or origin.startswith('NOAO-IRAF')): - instrument_name = OmegaCam() - break - - elif instrument == 'ANDICAM-CCD' and hdr.get('OBSERVAT') == 'CTIO': - instrument_name = AndiCam() - break - - elif telescope == '1.3m PAIRITEL' and instrument == '2MASS Survey cam': - instrument_name = PairTel() - break - - elif (origin == 'OAdM' or origin.startswith('NOAO-IRAF')) and telescope == 'TJO' and instrument in ( - 'MEIA3', 'MEIA2'): - instrument_name = TJO() - break - - else: - raise RuntimeError("Could not determine origin of image") - - image = FlowsImage(image=np.asarray(hdul[ext].data, dtype='float64'), header=hdr, mask=mask) - instrument_name.process_image(image) - - return instrument_name.image +instruments = {'LCOGT': LCOGT, 'HAWKI': HAWKI, 'ALFOSC': ALFOSC, 'NOTCAM': NOTCAM, 'PS1': PS1, 'Liverpool': Liverpool, + 'Omega2000': Omega2000, 'Swope': Swope, 'Dupont': Dupont, 'Retrocam': RetroCam, 'Baade': Baade, + 'Sofi': Sofi, 'EFOSC': EFOSC, 'AstroNIRCam': AstroNIRCam, 'OmegaCam': OmegaCam, 'AndiCam': AndiCam, + 'PairTel': PairTel, 'TJO': TJO} -# get image and wcs solution -# with fits.open(FILENAME, mode='readonly') as hdul: -# -------------------------------------------------------------------------------------------------- -def load_image(FILENAME, target_coord=None): - """ - Load FITS image. - - Parameters: - FILENAME (str): Path to FITS file to be loaded. - target_coord (:class:`astropy.coordinates.SkyCoord`): Coordinates of target. - Only used for HAWKI images to determine which image extension to load, - for all other images it is ignored. - - Returns: - object: Image constainer. - - .. codeauthor:: Rasmus Handberg - """ - - logger = logging.getLogger(__name__) - - # Get image and WCS, find stars, remove galaxies - image = type('image', (object,), dict()) # image container - - # get image and wcs solution - with fits.open(FILENAME, mode='readonly') as hdul: - - hdr = hdul[0].header - image.header = hdr - origin = hdr.get('ORIGIN', '') - telescope = hdr.get('TELESCOP', '') - instrument = hdr.get('INSTRUME', '') - - # Load image data: - image.image = np.asarray(hdul[0].data, dtype='float64') - image.shape = image.image.shape - - # Load image mask: - if origin == 'LCOGT': - if 'BPM' in hdul: - image.mask = np.asarray(hdul['BPM'].data, dtype='bool') - else: - logger.warning('LCOGT image does not contain bad pixel map. Not applying mask.') - image.mask = np.zeros_like(image.image, dtype='bool') - else: - image.mask = np.zeros_like(image.image, dtype='bool') - - image.mask |= ~np.isfinite(image.image) - - # World Coordinate System: - with warnings.catch_warnings(): - warnings.simplefilter('ignore', category=FITSFixedWarning) - image.wcs = WCS(header=hdr, relax=True) - - # Values which will be filled out below, depending on the instrument: - image.exptime = hdr.get('EXPTIME', None) # Exposure time * u.second - image.peakmax = None # Maximum value above which data is not to be trusted - - # Timestamp: - if origin == 'LCOGT': - sites = api.sites.get_all_sites() - site_keywords = {s['site_keyword']: s for s in sites} - image.site = site_keywords.get(hdr['SITE'], None) - - observatory = coords.EarthLocation.from_geodetic(lat=hdr['LATITUDE'], lon=hdr['LONGITUD'], - height=hdr['HEIGHT']) - image.obstime = Time(hdr['DATE-OBS'], format='isot', scale='utc', location=observatory) - image.obstime += 0.5 * image.exptime * u.second # Make time centre of exposure - - image.photfilter = { - 'zs': 'zp' - }.get(hdr['FILTER'], hdr['FILTER']) - - # Get non-linear limit - # TODO: Use actual or some fraction of the non-linearity limit - # image.peakmax = hdr.get('MAXLIN') # Presumed non-linearity limit from header - image.peakmax = 60000 # From experience, this one is better. - - elif origin == 'ESO-PARANAL' and telescope == 'ESO-VLT-U4' and instrument == 'HAWKI' and hdr.get( - 'PRODCATG') == 'SCIENCE.MEFIMAGE': - image.site = api.get_site(2) # Hard-coded the siteid for ESO Paranal, VLT, UT4 - image.obstime = Time(hdr['DATE-OBS'], format='isot', scale='utc', location=image.site['EarthLocation']) - image.obstime += 0.5 * image.exptime * u.second # Make time centre of exposure - image.photfilter = hdr['FILTER'] - - # For HAWKI multi-extension images we search the extensions for which one contains - # the target, and overwrites the image data with that: - if target_coord is None: - raise ValueError("TARGET_COORD is needed for HAWKI images to find the correct extension") - target_radec = [[target_coord.icrs.ra.deg, target_coord.icrs.dec.deg]] - for k in range(1, 5): - w = WCS(header=hdul[k].header, relax=True) - s = [hdul[k].header['NAXIS2'], hdul[k].header['NAXIS1']] - pix = w.all_world2pix(target_radec, 0).flatten() - if pix[0] >= -0.5 and pix[1] >= -0.5 and pix[0] <= s[1] - 0.5 and pix[1] <= s[0] - 0.5: - ob_type = hdul[k].header["HIERARCH ESO OCS DET1 IMGNAME"].split('_')[-1] - if "Auto" in ob_type: - image.ob_type = 'Autojitter' - elif "Fixed" in ob_type: - image.ob_type = 'FixedOffset' - # Should we use a default instead of raising? - else: - raise RuntimeError("Image OB Type not AutoJitter or FixedOffset") - image.image = np.asarray(hdul[k].data, dtype='float64') - image.shape = image.image.shape - image.wcs = w - image.mask = ~np.isfinite(image.image) - break - else: - raise RuntimeError("Could not find image extension that target is on") - - elif telescope == 'NOT' and instrument in ('ALFOSC FASU', 'ALFOSC_FASU') and hdr.get('OBS_MODE', - '').lower() == 'imaging': - image.site = api.get_site(5) # Hard-coded the siteid for NOT - image.obstime = Time(hdr['DATE-AVG'], format='isot', scale='utc', location=image.site['EarthLocation']) - - # Sometimes data from NOT does not have the FILTER keyword, - # in which case we have to try to figure out which filter - # was used based on some of the other headers: - if 'FILTER' in hdr: - image.photfilter = { - 'B Bes': 'B', - 'V Bes': 'V', - 'R Bes': 'R', - 'g SDSS': 'gp', - 'r SDSS': 'rp', - 'i SDSS': 'ip', - 'i int': 'ip', # Interference filter - 'u SDSS': 'up', - 'z SDSS': 'zp' - }.get(hdr['FILTER'].replace('_', ' '), hdr['FILTER']) - else: - filters_used = [] - for check_headers in ('ALFLTNM', 'FAFLTNM', 'FBFLTNM'): - if hdr.get(check_headers) and hdr.get(check_headers).strip().lower() != 'open': - filters_used.append(hdr.get(check_headers).strip()) - if len(filters_used) == 1: - image.photfilter = { - 'B_Bes 440_100': 'B', - 'V_Bes 530_80': 'V', - 'R_Bes 650_130': 'R', - "g'_SDSS 480_145": 'gp', - "r'_SDSS 618_148": 'rp', - "i'_SDSS 771_171": 'ip', - 'i_int 797_157': 'ip', # Interference filter - "z'_SDSS 832_LP": 'zp' - }.get(filters_used[0].replace(' ', ' '), filters_used[0]) - else: - raise RuntimeError("Could not determine filter used.") - - # Get non-linear limit - # Obtained from http://www.not.iac.es/instruments/detectors/CCD14/LED-linearity/20181026-200-1x1.pdf - # TODO: grab these from a table for all detector setups of ALFOSC - image.peakmax = 80000 # For ALFOSC D, 1x1, 200; the standard for SNe. - - elif telescope == 'NOT' and instrument == 'NOTCAM' and hdr.get('OBS_MODE', '').lower() == 'imaging': - image.site = api.get_site(5) # Hard-coded the siteid for NOT - image.obstime = Time(hdr['DATE-AVG'], format='isot', scale='utc', location=image.site['EarthLocation']) - - # Does NOTCAM data sometimes contain a FILTER header? - # if not we have to try to figure out which filter - # was used based on some of the other headers: - if 'FILTER' in hdr: - raise RuntimeError("NOTCAM: Filter keyword defined") - filters_used = [] - for check_headers in ('NCFLTNM1', 'NCFLTNM2'): - if hdr.get(check_headers) and hdr.get(check_headers).strip().lower() != 'open': - filters_used.append(hdr.get(check_headers).strip()) - if len(filters_used) == 1: - image.photfilter = { - 'Ks': 'K' - }.get(filters_used[0], filters_used[0]) - else: - raise RuntimeError("Could not determine filter used.") - - # Mask out "halo" of pixels with zero value along edge of image: - image.mask |= edge_mask(image.image, value=0) - - elif hdr.get('FPA.TELESCOPE') == 'PS1' and hdr.get('FPA.INSTRUMENT') == 'GPC1': - image.site = api.get_site(6) # Hard-coded the siteid for Pan-STARRS1 - image.obstime = Time(hdr['MJD-OBS'], format='mjd', scale='utc', location=image.site['EarthLocation']) - - image.photfilter = { - 'g.00000': 'gp', - 'r.00000': 'rp', - 'i.00000': 'ip', - 'z.00000': 'zp' - }.get(hdr['FPA.FILTER'], hdr['FPA.FILTER']) - - elif telescope == 'Liverpool Telescope': - # Liverpool telescope - image.site = api.get_site(8) # Hard-coded the siteid for Liverpool Telescope - image.obstime = Time(hdr['DATE-OBS'], format='isot', scale='utc', location=image.site['EarthLocation']) - image.obstime += 0.5 * image.exptime * u.second # Make time centre of exposure - image.photfilter = { - 'Bessel-B': 'B', - 'Bessell-B': 'B', - 'Bessel-V': 'V', - 'Bessell-V': 'V', - 'SDSS-U': 'up', - 'SDSS-G': 'gp', - 'SDSS-R': 'rp', - 'SDSS-I': 'ip', - 'SDSS-Z': 'zp' - }.get(hdr['FILTER1'], hdr['FILTER1']) - - elif telescope == 'CA 3.5m' and instrument == 'Omega2000': - # Calar Alto 3.5m (Omege2000) - image.site = api.get_site(9) # Hard-coded the siteid for Calar Alto 3.5m - image.obstime = Time(hdr['MJD-OBS'], format='mjd', scale='utc', location=image.site['EarthLocation']) - image.obstime += 0.5 * image.exptime * u.second # Make time centre of exposure - image.photfilter = hdr['FILTER'] - - elif telescope == 'SWO' and hdr.get('SITENAME') == 'LCO': - image.site = api.get_site(10) # Hard-coded the siteid for Swope, Las Campanas Observatory - image.obstime = Time(hdr['JD'], format='jd', scale='utc', location=image.site['EarthLocation']) - image.photfilter = { - 'u': 'up', - 'g': 'gp', - 'r': 'rp', - 'i': 'ip', - }.get(hdr['FILTER'], hdr['FILTER']) - - elif telescope == 'DUP' and hdr.get('SITENAME') == 'LCO' and instrument == 'Direct/SITe2K-1': - image.site = api.get_site(14) # Hard-coded the siteid for Du Pont, Las Campanas Observatory - image.obstime = Time(hdr['JD'], format='jd', scale='utc', location=image.site['EarthLocation']) - image.photfilter = { - 'u': 'up', - 'g': 'gp', - 'r': 'rp', - 'i': 'ip', - }.get(hdr['FILTER'], hdr['FILTER']) - - elif telescope == 'DUP' and instrument == 'RetroCam': - image.site = api.get_site(16) # Hard-coded the siteid for Du Pont, Las Campanas Observatory - image.obstime = Time(hdr['JD'], format='jd', scale='utc', location=image.site['EarthLocation']) - image.photfilter = { - 'Yc': 'Y', - 'Hc': 'H', - 'Jo': 'J', - }.get(hdr['FILTER'], hdr['FILTER']) - - elif telescope == 'Baade' and hdr.get('SITENAME') == 'LCO' and instrument == 'FourStar': - image.site = api.get_site(11) # Hard-coded the siteid for Swope, Las Campanas Observatory - image.obstime = Time(hdr['JD'], format='jd', scale='utc', location=image.site['EarthLocation']) - image.photfilter = { - 'Ks': 'K', - 'J1': 'Y', - }.get(hdr['FILTER'], hdr['FILTER']) - image.exptime *= int(hdr['NCOMBINE']) # EXPTIME is only for a single exposure - - elif instrument == 'SOFI' and telescope in ('ESO-NTT', 'other') and ( - origin == 'ESO' or origin.startswith('NOAO-IRAF')): - image.site = api.get_site(12) # Hard-coded the siteid for SOFT, ESO NTT - if 'TMID' in hdr: - image.obstime = Time(hdr['TMID'], format='mjd', scale='utc', location=image.site['EarthLocation']) - else: - image.obstime = Time(hdr['MJD-OBS'], format='mjd', scale='utc', location=image.site['EarthLocation']) - image.obstime += 0.5 * image.exptime * u.second # Make time centre of exposure - - # Photometric filter: - photfilter_translate = { - 'Ks': 'K' - } - if 'FILTER' in hdr: - image.photfilter = photfilter_translate.get(hdr['FILTER'], hdr['FILTER']) - else: - filters_used = [] - for check_headers in ('ESO INS FILT1 ID', 'ESO INS FILT2 ID'): - if hdr.get(check_headers) and hdr.get(check_headers).strip().lower() != 'open': - filters_used.append(hdr.get(check_headers).strip()) - if len(filters_used) == 1: - image.photfilter = photfilter_translate.get(filters_used[0], filters_used[0]) - else: - raise RuntimeError("Could not determine filter used.") - - # Mask out "halo" of pixels with zero value along edge of image: - image.mask |= edge_mask(image.image, value=0) - - elif telescope == 'ESO-NTT' and instrument == 'EFOSC' and (origin == 'ESO' or origin.startswith('NOAO-IRAF')): - image.site = api.get_site(15) # Hard-coded the siteid for EFOSC, ESO NTT - image.obstime = Time(hdr['DATE-OBS'], format='isot', scale='utc', location=image.site['EarthLocation']) - image.obstime += 0.5 * image.exptime * u.second # Make time centre of exposure - image.photfilter = { - 'g782': 'gp', - 'r784': 'rp', - 'i705': 'ip', - 'B639': 'B', - 'V641': 'V' - }.get(hdr['FILTER'], hdr['FILTER']) - - elif telescope == 'SAI-2.5' and instrument == 'ASTRONIRCAM': - image.site = api.get_site(13) # Hard-coded the siteid for Caucasus Mountain Observatory - if 'MIDPOINT' in hdr: - image.obstime = Time(hdr['MIDPOINT'], format='isot', scale='utc', location=image.site['EarthLocation']) - else: - image.obstime = Time(hdr['MJD-AVG'], format='mjd', scale='utc', location=image.site['EarthLocation']) - image.photfilter = { - 'H_Open': 'H', - 'K_Open': 'K', - }.get(hdr['FILTER'], hdr['FILTER']) - image.exptime = hdr.get('FULL_EXP', image.exptime) - - elif instrument == 'OMEGACAM' and (origin == 'ESO' or origin.startswith('NOAO-IRAF')): - image.site = api.get_site(18) # Hard-coded the siteid for ESO VLT Survey telescope - image.obstime = Time(hdr['MJD-OBS'], format='mjd', scale='utc', location=image.site['EarthLocation']) - image.obstime += 0.5 * image.exptime * u.second # Make time centre of exposure - image.photfilter = { - 'i_SDSS': 'ip' - }.get(hdr['ESO INS FILT1 NAME'], hdr['ESO INS FILT1 NAME']) - - elif instrument == 'ANDICAM-CCD' and hdr.get('OBSERVAT') == 'CTIO': - image.site = api.get_site( - 20) # Hard-coded the siteid for ANDICAM at Cerro Tololo Interamerican Observatory (CTIO) - image.obstime = Time(hdr['JD'], format='jd', scale='utc', location=image.site['EarthLocation']) - image.obstime += 0.5 * image.exptime * u.second # Make time centre of exposure - image.photfilter = hdr['CCDFLTID'] - - elif telescope == '1.3m PAIRITEL' and instrument == '2MASS Survey cam': - image.site = api.get_site(21) # Hard-coded the siteid for Peters Automated InfraRed Imaging TELescope - time_start = Time(hdr['STRT_CPU'], format='iso', scale='utc', location=image.site['EarthLocation']) - time_stop = Time(hdr['STOP_CPU'], format='iso', scale='utc', location=image.site['EarthLocation']) - image.obstime = time_start + 0.5 * (time_stop - time_start) - image.photfilter = { - 'j': 'J', - 'h': 'H', - 'k': 'K', - }.get(hdr['FILTER'], hdr['FILTER']) - - # Mask out "halo" of pixels with zero value along edge of image: - image.mask |= edge_mask(image.image, value=0) - - elif (origin == 'OAdM' or origin.startswith('NOAO-IRAF')) and telescope == 'TJO' and instrument in ( - 'MEIA3', 'MEIA2'): - image.site = api.get_site( - 22) # Hard-coded the siteid for Telescopi Joan Oró (TJO) at Observatori Astronòmic del Montsec - image.obstime = Time(hdr['JD'], format='jd', scale='utc', location=image.site['EarthLocation']) - image.obstime += 0.5 * image.exptime * u.second # Make time centre of exposure - image.photfilter = hdr['FILTER'] - - else: - raise RuntimeError("Could not determine origin of image") - - # Sanity checks: - if image.exptime is None: - raise ValueError("Image exposure time could not be extracted") - - # Create masked version of image: - image.image[image.mask] = np.NaN - image.clean = np.ma.masked_array(data=image.image, mask=image.mask, copy=False) - - return image +def edge_mask(img, value=0): + """ + Create boolean mask of given value near edge of image. + + Parameters: + img (ndarray): Image of + value (float): Value to detect near edge. Default=0. + + Returns: + ndarray: Pixel mask with given values on the edge of image. + + .. codeauthor:: Rasmus Handberg + """ + + mask1 = (img == value) + mask = np.zeros_like(img, dtype='bool') + + # Mask entire rows and columns which are only the value: + mask[np.all(mask1, axis=1), :] = True + mask[:, np.all(mask1, axis=0)] = True + + # Detect "uneven" edges column-wise in image: + a = np.argmin(mask1, axis=0) + b = np.argmin(np.flipud(mask1), axis=0) + for col in range(img.shape[1]): + if mask1[0, col]: + mask[:a[col], col] = True + if mask1[-1, col]: + mask[-b[col]:, col] = True + + # Detect "uneven" edges row-wise in image: + a = np.argmin(mask1, axis=1) + b = np.argmin(np.fliplr(mask1), axis=1) + for row in range(img.shape[0]): + if mask1[row, 0]: + mask[row, :a[row]] = True + if mask1[row, -1]: + mask[row, -b[row]:] = True + + return mask + + +def load_image(filename: str, target_coord: typing.Union[coords.SkyCoord, typing.Tuple[float, float]] = None): + """ + Load FITS image using FlowsImage class and Instrument Classes. + + Parameters: + filename (str): Path to FITS file to be loaded. + target_coord (:class:`astropy.coordinates.SkyCoord`): Coordinates of target. + Only used for HAWKI images to determine which image extension to load, + for all other images it is ignored. + + Returns: + FlowsImage: instance of FlowsImage with valuues populated based on instrument. + + .. codeauthor:: Emir Karamehmetoglu + .. codeauthor:: Rasmus Handberg + """ + ext = 0 # Default extension in HDUList, individual instruments may override this. + mask = None # Instrument can override, default is to only mask all non-finite values, override is additive. + + # Read fits image, Structural Pattern Match to specific instrument. + with fits.open(filename, mode='readonly') as hdul: + hdr = hdul[ext].header + origin = hdr.get('ORIGIN', '') + telescope = hdr.get('TELESCOP', '') + instrument = hdr.get('INSTRUME', '') + + instrument_name = None + # Pattern matching begins here, ideally we use 3.10 pattern matching, or a dictionary lookup. + while instrument_name is None: + # LCOGT + if origin == "LCOGT": + instrument_name = LCOGT() + if 'BPM' in hdul: + mask = np.asarray(hdul['BPM'].data, dtype='bool') + else: + logger.warning('LCOGT image does not contain bad pixel map. Not applying mask.') + break + # HAWKI + elif origin == 'ESO-PARANAL' and telescope == 'ESO-VLT-U4' and instrument == 'HAWKI' and hdr.get( + 'PRODCATG') == 'SCIENCE.MEFIMAGE': + instrument_name = HAWKI() + + if target_coord is None: + raise ValueError("TARGET_COORD is needed for HAWKI images to find the correct extension") + if not isinstance(target_coord, coords.SkyCoord): + if len(target_coord) == 2: + target_coord = coords.SkyCoord(ra=target_coord[0] * u.deg, dec=target_coord[1] * u.deg, + frame='icrs') + else: + raise ValueError("TARGET_COORD is needed for HAWKI images to find the correct extension") + ext = get_image_extension(hdul, target_coord) # Find the one with the SN in it. + hdr = hdul[ext].header + hdul[0].header + break + + # NOT - ALFOSC + elif telescope == 'NOT' and instrument in ('ALFOSC FASU', 'ALFOSC_FASU') \ + and hdr.get('OBS_MODE', '').lower() == 'imaging': + instrument_name = ALFOSC() + break + + elif telescope == 'NOT' and instrument == 'NOTCAM' and hdr.get('OBS_MODE', '').lower() == 'imaging': + instrument_name = NOTCAM() + break + + elif hdr.get('FPA.TELESCOPE') == 'PS1' and hdr.get('FPA.INSTRUMENT') == 'GPC1': + instrument_name = PS1() + break + + elif telescope == 'Liverpool Telescope': + instrument_name = Liverpool() + break + + elif telescope == 'CA 3.5m' and instrument == 'Omega2000': + instrument_name = Omega2000() + break + + elif telescope == 'SWO' and hdr.get('SITENAME') == 'LCO': + instrument_name = Swope() + break + + elif telescope == 'DUP' and hdr.get('SITENAME') == 'LCO' and instrument == 'Direct/SITe2K-1': + instrument_name = Dupont() + break + + elif telescope == 'DUP' and instrument == 'RetroCam': + instrument_name = RetroCam() + break + + elif telescope == 'Baade' and hdr.get('SITENAME') == 'LCO' and instrument == 'FourStar': + instrument_name = Baade() + break + + elif instrument == 'SOFI' and telescope in ('ESO-NTT', 'other') and ( + origin == 'ESO' or origin.startswith('NOAO-IRAF')): + instrument_name = Sofi() + break + + elif telescope == 'ESO-NTT' and instrument == 'EFOSC' and ( + origin == 'ESO' or origin.startswith('NOAO-IRAF')): + instrument_name = EFOSC() + break + + elif telescope == 'SAI-2.5' and instrument == 'ASTRONIRCAM': + instrument_name = AstroNIRCam() + break + + elif instrument == 'OMEGACAM' and (origin == 'ESO' or origin.startswith('NOAO-IRAF')): + instrument_name = OmegaCam() + break + + elif instrument == 'ANDICAM-CCD' and hdr.get('OBSERVAT') == 'CTIO': + instrument_name = AndiCam() + break + + elif telescope == '1.3m PAIRITEL' and instrument == '2MASS Survey cam': + instrument_name = PairTel() + break + + elif (origin == 'OAdM' or origin.startswith('NOAO-IRAF')) and telescope == 'TJO' and instrument in ( + 'MEIA3', 'MEIA2'): + instrument_name = TJO() + break + + else: + raise RuntimeError("Could not determine origin of image") + + image = FlowsImage(image=np.asarray(hdul[ext].data, dtype='float64'), header=hdr, mask=mask) + instrument_name.process_image(image) + + return instrument_name.image def main(): - filename = "../ADP.2021-11-08T12_34_46.383.fits" - img = new_load_image(filename, (325.772590825, -17.54410925)) - print(img.peakmax) - return img + filename = "../tests/input/ADP.2021-11-08T12_34_46.383.fits" + img = load_image(filename, (325.772590825, -17.54410925)) + print(img.peakmax) + return img if __name__ == '__main__': - main() + img = main() + +# get image and wcs solution +# with fits.open(FILENAME, mode='readonly') as hdul: + +# -------------------------------------------------------------------------------------------------- +# def load_image(FILENAME, target_coord=None): +# """ +# Load FITS image. +# +# Parameters: +# FILENAME (str): Path to FITS file to be loaded. +# target_coord (:class:`astropy.coordinates.SkyCoord`): Coordinates of target. +# Only used for HAWKI images to determine which image extension to load, +# for all other images it is ignored. +# +# Returns: +# object: Image constainer. +# +# .. codeauthor:: Rasmus Handberg +# """ +# +# logger = logging.getLogger(__name__) +# +# # Get image and WCS, find stars, remove galaxies +# image = type('image', (object,), dict()) # image container +# +# # get image and wcs solution +# with fits.open(FILENAME, mode='readonly') as hdul: +# +# hdr = hdul[0].header +# image.header = hdr +# origin = hdr.get('ORIGIN', '') +# telescope = hdr.get('TELESCOP', '') +# instrument = hdr.get('INSTRUME', '') +# +# # Load image data: +# image.image = np.asarray(hdul[0].data, dtype='float64') +# image.shape = image.image.shape +# +# # Load image mask: +# if origin == 'LCOGT': +# if 'BPM' in hdul: +# image.mask = np.asarray(hdul['BPM'].data, dtype='bool') +# else: +# logger.warning('LCOGT image does not contain bad pixel map. Not applying mask.') +# image.mask = np.zeros_like(image.image, dtype='bool') +# else: +# image.mask = np.zeros_like(image.image, dtype='bool') +# +# image.mask |= ~np.isfinite(image.image) +# +# # World Coordinate System: +# with warnings.catch_warnings(): +# warnings.simplefilter('ignore', category=FITSFixedWarning) +# image.wcs = WCS(header=hdr, relax=True) +# +# # Values which will be filled out below, depending on the instrument: +# image.exptime = hdr.get('EXPTIME', None) # Exposure time * u.second +# image.peakmax = None # Maximum value above which data is not to be trusted +# +# # Timestamp: +# if origin == 'LCOGT': +# sites = api.sites.get_all_sites() +# site_keywords = {s['site_keyword']: s for s in sites} +# image.site = site_keywords.get(hdr['SITE'], None) +# +# observatory = coords.EarthLocation.from_geodetic(lat=hdr['LATITUDE'], lon=hdr['LONGITUD'], +# height=hdr['HEIGHT']) +# image.obstime = Time(hdr['DATE-OBS'], format='isot', scale='utc', location=observatory) +# image.obstime += 0.5 * image.exptime * u.second # Make time centre of exposure +# +# image.photfilter = {'zs': 'zp'}.get(hdr['FILTER'], hdr['FILTER']) +# +# # Get non-linear limit +# # TODO: Use actual or some fraction of the non-linearity limit +# # image.peakmax = hdr.get('MAXLIN') # Presumed non-linearity limit from header +# image.peakmax = 60000 # From experience, this one is better. +# +# elif origin == 'ESO-PARANAL' and telescope == 'ESO-VLT-U4' and instrument == 'HAWKI' and hdr.get( +# 'PRODCATG') == 'SCIENCE.MEFIMAGE': +# image.site = api.get_site(2) # Hard-coded the siteid for ESO Paranal, VLT, UT4 +# image.obstime = Time(hdr['DATE-OBS'], format='isot', scale='utc', location=image.site['EarthLocation']) +# image.obstime += 0.5 * image.exptime * u.second # Make time centre of exposure +# image.photfilter = hdr['FILTER'] +# +# # For HAWKI multi-extension images we search the extensions for which one contains +# # the target, and overwrites the image data with that: +# if target_coord is None: +# raise ValueError("TARGET_COORD is needed for HAWKI images to find the correct extension") +# target_radec = [[target_coord.icrs.ra.deg, target_coord.icrs.dec.deg]] +# for k in range(1, 5): +# w = WCS(header=hdul[k].header, relax=True) +# s = [hdul[k].header['NAXIS2'], hdul[k].header['NAXIS1']] +# pix = w.all_world2pix(target_radec, 0).flatten() +# if pix[0] >= -0.5 and pix[1] >= -0.5 and pix[0] <= s[1] - 0.5 and pix[1] <= s[0] - 0.5: +# ob_type = hdul[k].header["HIERARCH ESO OCS DET1 IMGNAME"].split('_')[-1] +# if "Auto" in ob_type: +# image.ob_type = 'Autojitter' +# elif "Fixed" in ob_type: +# image.ob_type = 'FixedOffset' +# # Should we use a default instead of raising? +# else: +# raise RuntimeError("Image OB Type not AutoJitter or FixedOffset") +# image.image = np.asarray(hdul[k].data, dtype='float64') +# image.shape = image.image.shape +# image.wcs = w +# image.mask = ~np.isfinite(image.image) +# break +# else: +# raise RuntimeError("Could not find image extension that target is on") +# +# elif telescope == 'NOT' and instrument in ('ALFOSC FASU', 'ALFOSC_FASU') and hdr.get('OBS_MODE', +# '').lower() == 'imaging': +# image.site = api.get_site(5) # Hard-coded the siteid for NOT +# image.obstime = Time(hdr['DATE-AVG'], format='isot', scale='utc', location=image.site['EarthLocation']) +# +# # Sometimes data from NOT does not have the FILTER keyword, +# # in which case we have to try to figure out which filter +# # was used based on some of the other headers: +# if 'FILTER' in hdr: +# image.photfilter = {'B Bes': 'B', 'V Bes': 'V', 'R Bes': 'R', 'g SDSS': 'gp', 'r SDSS': 'rp', +# 'i SDSS': 'ip', 'i int': 'ip', # Interference filter +# 'u SDSS': 'up', 'z SDSS': 'zp'}.get(hdr['FILTER'].replace('_', ' '), hdr['FILTER']) +# else: +# filters_used = [] +# for check_headers in ('ALFLTNM', 'FAFLTNM', 'FBFLTNM'): +# if hdr.get(check_headers) and hdr.get(check_headers).strip().lower() != 'open': +# filters_used.append(hdr.get(check_headers).strip()) +# if len(filters_used) == 1: +# image.photfilter = {'V_Bes 530_80': 'V', 'R_Bes 650_130': 'R', "g'_SDSS 480_145": 'gp', +# "r'_SDSS 618_148": 'rp', "i'_SDSS 771_171": 'ip', 'i_int 797_157': 'ip', +# # Interference filter +# "z'_SDSS 832_LP": 'zp'}.get(filters_used[0].replace(' ', ' '), filters_used[0]) +# else: +# raise RuntimeError("Could not determine filter used.") +# +# # Get non-linear limit +# # Obtained from http://www.not.iac.es/instruments/detectors/CCD14/LED-linearity/20181026-200-1x1.pdf +# # TODO: grab these from a table for all detector setups of ALFOSC +# image.peakmax = 80000 # For ALFOSC D, 1x1, 200; the standard for SNe. +# +# elif telescope == 'NOT' and instrument == 'NOTCAM' and hdr.get('OBS_MODE', '').lower() == 'imaging': +# image.site = api.get_site(5) # Hard-coded the siteid for NOT +# image.obstime = Time(hdr['DATE-AVG'], format='isot', scale='utc', location=image.site['EarthLocation']) +# +# # Does NOTCAM data sometimes contain a FILTER header? +# # if not we have to try to figure out which filter +# # was used based on some of the other headers: +# if 'FILTER' in hdr: +# raise RuntimeError("NOTCAM: Filter keyword defined") +# filters_used = [] +# for check_headers in ('NCFLTNM1', 'NCFLTNM2'): +# if hdr.get(check_headers) and hdr.get(check_headers).strip().lower() != 'open': +# filters_used.append(hdr.get(check_headers).strip()) +# if len(filters_used) == 1: +# image.photfilter = {}.get(filters_used[0], filters_used[0]) +# else: +# raise RuntimeError("Could not determine filter used.") +# +# # Mask out "halo" of pixels with zero value along edge of image: +# image.mask |= edge_mask(image.image, value=0) +# +# elif hdr.get('FPA.TELESCOPE') == 'PS1' and hdr.get('FPA.INSTRUMENT') == 'GPC1': +# image.site = api.get_site(6) # Hard-coded the siteid for Pan-STARRS1 +# image.obstime = Time(hdr['MJD-OBS'], format='mjd', scale='utc', location=image.site['EarthLocation']) +# +# image.photfilter = {'g.00000': 'gp', 'r.00000': 'rp', 'i.00000': 'ip', 'z.00000': 'zp'}.get( +# hdr['FPA.FILTER'], hdr['FPA.FILTER']) +# +# elif telescope == 'Liverpool Telescope': +# # Liverpool telescope +# image.site = api.get_site(8) # Hard-coded the siteid for Liverpool Telescope +# image.obstime = Time(hdr['DATE-OBS'], format='isot', scale='utc', location=image.site['EarthLocation']) +# image.obstime += 0.5 * image.exptime * u.second # Make time centre of exposure +# image.photfilter = {'Bessel-B': 'B', 'Bessell-B': 'B', 'Bessel-V': 'V', 'Bessell-V': 'V', 'SDSS-U': 'up', +# 'SDSS-G': 'gp', 'SDSS-R': 'rp', 'SDSS-I': 'ip', 'SDSS-Z': 'zp'}.get(hdr['FILTER1'], +# hdr['FILTER1']) +# +# elif telescope == 'CA 3.5m' and instrument == 'Omega2000': +# # Calar Alto 3.5m (Omege2000) +# image.site = api.get_site(9) # Hard-coded the siteid for Calar Alto 3.5m +# image.obstime = Time(hdr['MJD-OBS'], format='mjd', scale='utc', location=image.site['EarthLocation']) +# image.obstime += 0.5 * image.exptime * u.second # Make time centre of exposure +# image.photfilter = hdr['FILTER'] +# +# elif telescope == 'SWO' and hdr.get('SITENAME') == 'LCO': +# image.site = api.get_site(10) # Hard-coded the siteid for Swope, Las Campanas Observatory +# image.obstime = Time(hdr['JD'], format='jd', scale='utc', location=image.site['EarthLocation']) +# image.photfilter = {'u': 'up', 'g': 'gp', 'r': 'rp', 'i': 'ip', }.get(hdr['FILTER'], hdr['FILTER']) +# +# elif telescope == 'DUP' and hdr.get('SITENAME') == 'LCO' and instrument == 'Direct/SITe2K-1': +# image.site = api.get_site(14) # Hard-coded the siteid for Du Pont, Las Campanas Observatory +# image.obstime = Time(hdr['JD'], format='jd', scale='utc', location=image.site['EarthLocation']) +# image.photfilter = {'u': 'up', 'g': 'gp', 'r': 'rp', 'i': 'ip', }.get(hdr['FILTER'], hdr['FILTER']) +# +# elif telescope == 'DUP' and instrument == 'RetroCam': +# image.site = api.get_site(16) # Hard-coded the siteid for Du Pont, Las Campanas Observatory +# image.obstime = Time(hdr['JD'], format='jd', scale='utc', location=image.site['EarthLocation']) +# image.photfilter = {'Yc': 'Y', 'Hc': 'H', 'Jo': 'J', }.get(hdr['FILTER'], hdr['FILTER']) +# +# elif telescope == 'Baade' and hdr.get('SITENAME') == 'LCO' and instrument == 'FourStar': +# image.site = api.get_site(11) # Hard-coded the siteid for Swope, Las Campanas Observatory +# image.obstime = Time(hdr['JD'], format='jd', scale='utc', location=image.site['EarthLocation']) +# image.photfilter = {'Ks': 'K', 'J1': 'Y', }.get(hdr['FILTER'], hdr['FILTER']) +# image.exptime *= int(hdr['NCOMBINE']) # EXPTIME is only for a single exposure +# +# elif instrument == 'SOFI' and telescope in ('ESO-NTT', 'other') and ( +# origin == 'ESO' or origin.startswith('NOAO-IRAF')): +# image.site = api.get_site(12) # Hard-coded the siteid for SOFT, ESO NTT +# if 'TMID' in hdr: +# image.obstime = Time(hdr['TMID'], format='mjd', scale='utc', location=image.site['EarthLocation']) +# else: +# image.obstime = Time(hdr['MJD-OBS'], format='mjd', scale='utc', location=image.site['EarthLocation']) +# image.obstime += 0.5 * image.exptime * u.second # Make time centre of exposure +# +# # Photometric filter: +# photfilter_translate = {'Ks': 'K'} +# if 'FILTER' in hdr: +# image.photfilter = photfilter_translate.get(hdr['FILTER'], hdr['FILTER']) +# else: +# filters_used = [] +# for check_headers in ('ESO INS FILT1 ID', 'ESO INS FILT2 ID'): +# if hdr.get(check_headers) and hdr.get(check_headers).strip().lower() != 'open': +# filters_used.append(hdr.get(check_headers).strip()) +# if len(filters_used) == 1: +# image.photfilter = photfilter_translate.get(filters_used[0], filters_used[0]) +# else: +# raise RuntimeError("Could not determine filter used.") +# +# # Mask out "halo" of pixels with zero value along edge of image: +# image.mask |= edge_mask(image.image, value=0) +# +# elif telescope == 'ESO-NTT' and instrument == 'EFOSC' and (origin == 'ESO' or origin.startswith('NOAO-IRAF')): +# image.site = api.get_site(15) # Hard-coded the siteid for EFOSC, ESO NTT +# image.obstime = Time(hdr['DATE-OBS'], format='isot', scale='utc', location=image.site['EarthLocation']) +# image.obstime += 0.5 * image.exptime * u.second # Make time centre of exposure +# image.photfilter = {'g782': 'gp', 'r784': 'rp', 'i705': 'ip', 'B639': 'B', 'V641': 'V'}.get(hdr['FILTER'], +# hdr['FILTER']) +# +# elif telescope == 'SAI-2.5' and instrument == 'ASTRONIRCAM': +# image.site = api.get_site(13) # Hard-coded the siteid for Caucasus Mountain Observatory +# if 'MIDPOINT' in hdr: +# image.obstime = Time(hdr['MIDPOINT'], format='isot', scale='utc', location=image.site['EarthLocation']) +# else: +# image.obstime = Time(hdr['MJD-AVG'], format='mjd', scale='utc', location=image.site['EarthLocation']) +# image.photfilter = {'H_Open': 'H', 'K_Open': 'K', }.get(hdr['FILTER'], hdr['FILTER']) +# image.exptime = hdr.get('FULL_EXP', image.exptime) +# +# elif instrument == 'OMEGACAM' and (origin == 'ESO' or origin.startswith('NOAO-IRAF')): +# image.site = api.get_site(18) # Hard-coded the siteid for ESO VLT Survey telescope +# image.obstime = Time(hdr['MJD-OBS'], format='mjd', scale='utc', location=image.site['EarthLocation']) +# image.obstime += 0.5 * image.exptime * u.second # Make time centre of exposure +# image.photfilter = {'i_SDSS': 'ip'}.get(hdr['ESO INS FILT1 NAME'], hdr['ESO INS FILT1 NAME']) +# +# elif instrument == 'ANDICAM-CCD' and hdr.get('OBSERVAT') == 'CTIO': +# image.site = api.get_site( +# 20) # Hard-coded the siteid for ANDICAM at Cerro Tololo Interamerican Observatory (CTIO) +# image.obstime = Time(hdr['JD'], format='jd', scale='utc', location=image.site['EarthLocation']) +# image.obstime += 0.5 * image.exptime * u.second # Make time centre of exposure +# image.photfilter = hdr['CCDFLTID'] +# +# elif telescope == '1.3m PAIRITEL' and instrument == '2MASS Survey cam': +# image.site = api.get_site(21) # Hard-coded the siteid for Peters Automated InfraRed Imaging TELescope +# time_start = Time(hdr['STRT_CPU'], format='iso', scale='utc', location=image.site['EarthLocation']) +# time_stop = Time(hdr['STOP_CPU'], format='iso', scale='utc', location=image.site['EarthLocation']) +# image.obstime = time_start + 0.5 * (time_stop - time_start) +# image.photfilter = {'j': 'J', 'h': 'H', 'k': 'K', }.get(hdr['FILTER'], hdr['FILTER']) +# +# # Mask out "halo" of pixels with zero value along edge of image: +# image.mask |= edge_mask(image.image, value=0) +# +# elif (origin == 'OAdM' or origin.startswith('NOAO-IRAF')) and telescope == 'TJO' and instrument in ( +# 'MEIA3', 'MEIA2'): +# image.site = api.get_site( +# 22) # Hard-coded the siteid for Telescopi Joan Oró (TJO) at Observatori Astronòmic del Montsec +# image.obstime = Time(hdr['JD'], format='jd', scale='utc', location=image.site['EarthLocation']) +# image.obstime += 0.5 * image.exptime * u.second # Make time centre of exposure +# image.photfilter = hdr['FILTER'] +# +# else: +# raise RuntimeError("Could not determine origin of image") +# +# # Sanity checks: +# if image.exptime is None: +# raise ValueError("Image exposure time could not be extracted") +# +# # Create masked version of image: +# image.image[image.mask] = np.NaN +# image.clean = np.ma.masked_array(data=image.image, mask=image.mask, copy=False) +# +# return image diff --git a/tests/test_load_image.py b/tests/test_load_image.py index 3dd3927..2f202a6 100644 --- a/tests/test_load_image.py +++ b/tests/test_load_image.py @@ -15,6 +15,7 @@ import conftest # noqa: F401 from flows.api import get_filters from flows.load_image import load_image +from flows.api import get_filters #-------------------------------------------------------------------------------------------------- @pytest.mark.parametrize('fpath,siteid', [ @@ -32,18 +33,18 @@ def test_load_image(fpath, siteid): # Get list of all available filters: all_filters = set(get_filters().keys()) - # The test input directory containing the test-images: - INPUT_DIR = os.path.join(os.path.abspath(os.path.dirname(__file__)), 'input') + # The test input directory containing the test-images: + INPUT_DIR = os.path.join(os.path.abspath(os.path.dirname(__file__)), 'input') - # Target coordinates, only used for HAWKI image: - target_coord = SkyCoord( - ra=347.6230189, - dec=7.5888196, - unit='deg', - frame='icrs') + # Target coordinates, only used for HAWKI image: + target_coord = SkyCoord( + ra=347.6230189, + dec=7.5888196, + unit='deg', + frame='icrs') - # Load the image from the test-set: - img = load_image(os.path.join(INPUT_DIR, fpath), target_coord=target_coord) + # Load the image from the test-set: + img = load_image(os.path.join(INPUT_DIR, fpath), target_coord=target_coord) # Check the attributes of the image object: assert isinstance(img.image, np.ndarray) @@ -63,4 +64,4 @@ def test_load_image(fpath, siteid): #-------------------------------------------------------------------------------------------------- if __name__ == '__main__': - pytest.main([__file__]) + pytest.main([__file__]) From f28fb62d889cc2d696c655673bd457091c301a3a Mon Sep 17 00:00:00 2001 From: Emir Date: Mon, 14 Feb 2022 18:33:58 +0100 Subject: [PATCH 06/10] fix double import --- tests/test_load_image.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/test_load_image.py b/tests/test_load_image.py index 2f202a6..62b1219 100644 --- a/tests/test_load_image.py +++ b/tests/test_load_image.py @@ -15,7 +15,7 @@ import conftest # noqa: F401 from flows.api import get_filters from flows.load_image import load_image -from flows.api import get_filters + #-------------------------------------------------------------------------------------------------- @pytest.mark.parametrize('fpath,siteid', [ From ad6a9cf29db073a1edc384e174bf00f74f032e78 Mon Sep 17 00:00:00 2001 From: Emir Karamehmetoglu Date: Mon, 7 Mar 2022 11:36:56 +0100 Subject: [PATCH 07/10] add type hint backwards compat --- flows/load_image.py | 9 ++------- 1 file changed, 2 insertions(+), 7 deletions(-) diff --git a/flows/load_image.py b/flows/load_image.py index b9e3a61..573a4e2 100644 --- a/flows/load_image.py +++ b/flows/load_image.py @@ -1,12 +1,7 @@ -#!/usr/bin/env python3 -# -*- coding: utf-8 -*- """ -Flows photometry code. - -.. codeauthor:: Emir Karamehmetoglu -.. codeauthor:: Rasmus Handberg +Load image code. """ - +from __future__ import annotations import numpy as np import warnings import logging From b0b6839af26f4f4573d612319b6e702aad2cb537 Mon Sep 17 00:00:00 2001 From: Emir Karamehmetoglu Date: Wed, 9 Mar 2022 17:49:13 +0100 Subject: [PATCH 08/10] Remove old code --- flows/load_image.py | 302 -------------------------------------------- 1 file changed, 302 deletions(-) diff --git a/flows/load_image.py b/flows/load_image.py index 573a4e2..76282fa 100644 --- a/flows/load_image.py +++ b/flows/load_image.py @@ -700,305 +700,3 @@ def load_image(filename: str, target_coord: typing.Union[coords.SkyCoord, typing instrument_name.process_image(image) return instrument_name.image - - -def main(): - filename = "../tests/input/ADP.2021-11-08T12_34_46.383.fits" - img = load_image(filename, (325.772590825, -17.54410925)) - print(img.peakmax) - return img - - -if __name__ == '__main__': - img = main() - -# get image and wcs solution -# with fits.open(FILENAME, mode='readonly') as hdul: - -# -------------------------------------------------------------------------------------------------- -# def load_image(FILENAME, target_coord=None): -# """ -# Load FITS image. -# -# Parameters: -# FILENAME (str): Path to FITS file to be loaded. -# target_coord (:class:`astropy.coordinates.SkyCoord`): Coordinates of target. -# Only used for HAWKI images to determine which image extension to load, -# for all other images it is ignored. -# -# Returns: -# object: Image constainer. -# -# .. codeauthor:: Rasmus Handberg -# """ -# -# logger = logging.getLogger(__name__) -# -# # Get image and WCS, find stars, remove galaxies -# image = type('image', (object,), dict()) # image container -# -# # get image and wcs solution -# with fits.open(FILENAME, mode='readonly') as hdul: -# -# hdr = hdul[0].header -# image.header = hdr -# origin = hdr.get('ORIGIN', '') -# telescope = hdr.get('TELESCOP', '') -# instrument = hdr.get('INSTRUME', '') -# -# # Load image data: -# image.image = np.asarray(hdul[0].data, dtype='float64') -# image.shape = image.image.shape -# -# # Load image mask: -# if origin == 'LCOGT': -# if 'BPM' in hdul: -# image.mask = np.asarray(hdul['BPM'].data, dtype='bool') -# else: -# logger.warning('LCOGT image does not contain bad pixel map. Not applying mask.') -# image.mask = np.zeros_like(image.image, dtype='bool') -# else: -# image.mask = np.zeros_like(image.image, dtype='bool') -# -# image.mask |= ~np.isfinite(image.image) -# -# # World Coordinate System: -# with warnings.catch_warnings(): -# warnings.simplefilter('ignore', category=FITSFixedWarning) -# image.wcs = WCS(header=hdr, relax=True) -# -# # Values which will be filled out below, depending on the instrument: -# image.exptime = hdr.get('EXPTIME', None) # Exposure time * u.second -# image.peakmax = None # Maximum value above which data is not to be trusted -# -# # Timestamp: -# if origin == 'LCOGT': -# sites = api.sites.get_all_sites() -# site_keywords = {s['site_keyword']: s for s in sites} -# image.site = site_keywords.get(hdr['SITE'], None) -# -# observatory = coords.EarthLocation.from_geodetic(lat=hdr['LATITUDE'], lon=hdr['LONGITUD'], -# height=hdr['HEIGHT']) -# image.obstime = Time(hdr['DATE-OBS'], format='isot', scale='utc', location=observatory) -# image.obstime += 0.5 * image.exptime * u.second # Make time centre of exposure -# -# image.photfilter = {'zs': 'zp'}.get(hdr['FILTER'], hdr['FILTER']) -# -# # Get non-linear limit -# # TODO: Use actual or some fraction of the non-linearity limit -# # image.peakmax = hdr.get('MAXLIN') # Presumed non-linearity limit from header -# image.peakmax = 60000 # From experience, this one is better. -# -# elif origin == 'ESO-PARANAL' and telescope == 'ESO-VLT-U4' and instrument == 'HAWKI' and hdr.get( -# 'PRODCATG') == 'SCIENCE.MEFIMAGE': -# image.site = api.get_site(2) # Hard-coded the siteid for ESO Paranal, VLT, UT4 -# image.obstime = Time(hdr['DATE-OBS'], format='isot', scale='utc', location=image.site['EarthLocation']) -# image.obstime += 0.5 * image.exptime * u.second # Make time centre of exposure -# image.photfilter = hdr['FILTER'] -# -# # For HAWKI multi-extension images we search the extensions for which one contains -# # the target, and overwrites the image data with that: -# if target_coord is None: -# raise ValueError("TARGET_COORD is needed for HAWKI images to find the correct extension") -# target_radec = [[target_coord.icrs.ra.deg, target_coord.icrs.dec.deg]] -# for k in range(1, 5): -# w = WCS(header=hdul[k].header, relax=True) -# s = [hdul[k].header['NAXIS2'], hdul[k].header['NAXIS1']] -# pix = w.all_world2pix(target_radec, 0).flatten() -# if pix[0] >= -0.5 and pix[1] >= -0.5 and pix[0] <= s[1] - 0.5 and pix[1] <= s[0] - 0.5: -# ob_type = hdul[k].header["HIERARCH ESO OCS DET1 IMGNAME"].split('_')[-1] -# if "Auto" in ob_type: -# image.ob_type = 'Autojitter' -# elif "Fixed" in ob_type: -# image.ob_type = 'FixedOffset' -# # Should we use a default instead of raising? -# else: -# raise RuntimeError("Image OB Type not AutoJitter or FixedOffset") -# image.image = np.asarray(hdul[k].data, dtype='float64') -# image.shape = image.image.shape -# image.wcs = w -# image.mask = ~np.isfinite(image.image) -# break -# else: -# raise RuntimeError("Could not find image extension that target is on") -# -# elif telescope == 'NOT' and instrument in ('ALFOSC FASU', 'ALFOSC_FASU') and hdr.get('OBS_MODE', -# '').lower() == 'imaging': -# image.site = api.get_site(5) # Hard-coded the siteid for NOT -# image.obstime = Time(hdr['DATE-AVG'], format='isot', scale='utc', location=image.site['EarthLocation']) -# -# # Sometimes data from NOT does not have the FILTER keyword, -# # in which case we have to try to figure out which filter -# # was used based on some of the other headers: -# if 'FILTER' in hdr: -# image.photfilter = {'B Bes': 'B', 'V Bes': 'V', 'R Bes': 'R', 'g SDSS': 'gp', 'r SDSS': 'rp', -# 'i SDSS': 'ip', 'i int': 'ip', # Interference filter -# 'u SDSS': 'up', 'z SDSS': 'zp'}.get(hdr['FILTER'].replace('_', ' '), hdr['FILTER']) -# else: -# filters_used = [] -# for check_headers in ('ALFLTNM', 'FAFLTNM', 'FBFLTNM'): -# if hdr.get(check_headers) and hdr.get(check_headers).strip().lower() != 'open': -# filters_used.append(hdr.get(check_headers).strip()) -# if len(filters_used) == 1: -# image.photfilter = {'V_Bes 530_80': 'V', 'R_Bes 650_130': 'R', "g'_SDSS 480_145": 'gp', -# "r'_SDSS 618_148": 'rp', "i'_SDSS 771_171": 'ip', 'i_int 797_157': 'ip', -# # Interference filter -# "z'_SDSS 832_LP": 'zp'}.get(filters_used[0].replace(' ', ' '), filters_used[0]) -# else: -# raise RuntimeError("Could not determine filter used.") -# -# # Get non-linear limit -# # Obtained from http://www.not.iac.es/instruments/detectors/CCD14/LED-linearity/20181026-200-1x1.pdf -# # TODO: grab these from a table for all detector setups of ALFOSC -# image.peakmax = 80000 # For ALFOSC D, 1x1, 200; the standard for SNe. -# -# elif telescope == 'NOT' and instrument == 'NOTCAM' and hdr.get('OBS_MODE', '').lower() == 'imaging': -# image.site = api.get_site(5) # Hard-coded the siteid for NOT -# image.obstime = Time(hdr['DATE-AVG'], format='isot', scale='utc', location=image.site['EarthLocation']) -# -# # Does NOTCAM data sometimes contain a FILTER header? -# # if not we have to try to figure out which filter -# # was used based on some of the other headers: -# if 'FILTER' in hdr: -# raise RuntimeError("NOTCAM: Filter keyword defined") -# filters_used = [] -# for check_headers in ('NCFLTNM1', 'NCFLTNM2'): -# if hdr.get(check_headers) and hdr.get(check_headers).strip().lower() != 'open': -# filters_used.append(hdr.get(check_headers).strip()) -# if len(filters_used) == 1: -# image.photfilter = {}.get(filters_used[0], filters_used[0]) -# else: -# raise RuntimeError("Could not determine filter used.") -# -# # Mask out "halo" of pixels with zero value along edge of image: -# image.mask |= edge_mask(image.image, value=0) -# -# elif hdr.get('FPA.TELESCOPE') == 'PS1' and hdr.get('FPA.INSTRUMENT') == 'GPC1': -# image.site = api.get_site(6) # Hard-coded the siteid for Pan-STARRS1 -# image.obstime = Time(hdr['MJD-OBS'], format='mjd', scale='utc', location=image.site['EarthLocation']) -# -# image.photfilter = {'g.00000': 'gp', 'r.00000': 'rp', 'i.00000': 'ip', 'z.00000': 'zp'}.get( -# hdr['FPA.FILTER'], hdr['FPA.FILTER']) -# -# elif telescope == 'Liverpool Telescope': -# # Liverpool telescope -# image.site = api.get_site(8) # Hard-coded the siteid for Liverpool Telescope -# image.obstime = Time(hdr['DATE-OBS'], format='isot', scale='utc', location=image.site['EarthLocation']) -# image.obstime += 0.5 * image.exptime * u.second # Make time centre of exposure -# image.photfilter = {'Bessel-B': 'B', 'Bessell-B': 'B', 'Bessel-V': 'V', 'Bessell-V': 'V', 'SDSS-U': 'up', -# 'SDSS-G': 'gp', 'SDSS-R': 'rp', 'SDSS-I': 'ip', 'SDSS-Z': 'zp'}.get(hdr['FILTER1'], -# hdr['FILTER1']) -# -# elif telescope == 'CA 3.5m' and instrument == 'Omega2000': -# # Calar Alto 3.5m (Omege2000) -# image.site = api.get_site(9) # Hard-coded the siteid for Calar Alto 3.5m -# image.obstime = Time(hdr['MJD-OBS'], format='mjd', scale='utc', location=image.site['EarthLocation']) -# image.obstime += 0.5 * image.exptime * u.second # Make time centre of exposure -# image.photfilter = hdr['FILTER'] -# -# elif telescope == 'SWO' and hdr.get('SITENAME') == 'LCO': -# image.site = api.get_site(10) # Hard-coded the siteid for Swope, Las Campanas Observatory -# image.obstime = Time(hdr['JD'], format='jd', scale='utc', location=image.site['EarthLocation']) -# image.photfilter = {'u': 'up', 'g': 'gp', 'r': 'rp', 'i': 'ip', }.get(hdr['FILTER'], hdr['FILTER']) -# -# elif telescope == 'DUP' and hdr.get('SITENAME') == 'LCO' and instrument == 'Direct/SITe2K-1': -# image.site = api.get_site(14) # Hard-coded the siteid for Du Pont, Las Campanas Observatory -# image.obstime = Time(hdr['JD'], format='jd', scale='utc', location=image.site['EarthLocation']) -# image.photfilter = {'u': 'up', 'g': 'gp', 'r': 'rp', 'i': 'ip', }.get(hdr['FILTER'], hdr['FILTER']) -# -# elif telescope == 'DUP' and instrument == 'RetroCam': -# image.site = api.get_site(16) # Hard-coded the siteid for Du Pont, Las Campanas Observatory -# image.obstime = Time(hdr['JD'], format='jd', scale='utc', location=image.site['EarthLocation']) -# image.photfilter = {'Yc': 'Y', 'Hc': 'H', 'Jo': 'J', }.get(hdr['FILTER'], hdr['FILTER']) -# -# elif telescope == 'Baade' and hdr.get('SITENAME') == 'LCO' and instrument == 'FourStar': -# image.site = api.get_site(11) # Hard-coded the siteid for Swope, Las Campanas Observatory -# image.obstime = Time(hdr['JD'], format='jd', scale='utc', location=image.site['EarthLocation']) -# image.photfilter = {'Ks': 'K', 'J1': 'Y', }.get(hdr['FILTER'], hdr['FILTER']) -# image.exptime *= int(hdr['NCOMBINE']) # EXPTIME is only for a single exposure -# -# elif instrument == 'SOFI' and telescope in ('ESO-NTT', 'other') and ( -# origin == 'ESO' or origin.startswith('NOAO-IRAF')): -# image.site = api.get_site(12) # Hard-coded the siteid for SOFT, ESO NTT -# if 'TMID' in hdr: -# image.obstime = Time(hdr['TMID'], format='mjd', scale='utc', location=image.site['EarthLocation']) -# else: -# image.obstime = Time(hdr['MJD-OBS'], format='mjd', scale='utc', location=image.site['EarthLocation']) -# image.obstime += 0.5 * image.exptime * u.second # Make time centre of exposure -# -# # Photometric filter: -# photfilter_translate = {'Ks': 'K'} -# if 'FILTER' in hdr: -# image.photfilter = photfilter_translate.get(hdr['FILTER'], hdr['FILTER']) -# else: -# filters_used = [] -# for check_headers in ('ESO INS FILT1 ID', 'ESO INS FILT2 ID'): -# if hdr.get(check_headers) and hdr.get(check_headers).strip().lower() != 'open': -# filters_used.append(hdr.get(check_headers).strip()) -# if len(filters_used) == 1: -# image.photfilter = photfilter_translate.get(filters_used[0], filters_used[0]) -# else: -# raise RuntimeError("Could not determine filter used.") -# -# # Mask out "halo" of pixels with zero value along edge of image: -# image.mask |= edge_mask(image.image, value=0) -# -# elif telescope == 'ESO-NTT' and instrument == 'EFOSC' and (origin == 'ESO' or origin.startswith('NOAO-IRAF')): -# image.site = api.get_site(15) # Hard-coded the siteid for EFOSC, ESO NTT -# image.obstime = Time(hdr['DATE-OBS'], format='isot', scale='utc', location=image.site['EarthLocation']) -# image.obstime += 0.5 * image.exptime * u.second # Make time centre of exposure -# image.photfilter = {'g782': 'gp', 'r784': 'rp', 'i705': 'ip', 'B639': 'B', 'V641': 'V'}.get(hdr['FILTER'], -# hdr['FILTER']) -# -# elif telescope == 'SAI-2.5' and instrument == 'ASTRONIRCAM': -# image.site = api.get_site(13) # Hard-coded the siteid for Caucasus Mountain Observatory -# if 'MIDPOINT' in hdr: -# image.obstime = Time(hdr['MIDPOINT'], format='isot', scale='utc', location=image.site['EarthLocation']) -# else: -# image.obstime = Time(hdr['MJD-AVG'], format='mjd', scale='utc', location=image.site['EarthLocation']) -# image.photfilter = {'H_Open': 'H', 'K_Open': 'K', }.get(hdr['FILTER'], hdr['FILTER']) -# image.exptime = hdr.get('FULL_EXP', image.exptime) -# -# elif instrument == 'OMEGACAM' and (origin == 'ESO' or origin.startswith('NOAO-IRAF')): -# image.site = api.get_site(18) # Hard-coded the siteid for ESO VLT Survey telescope -# image.obstime = Time(hdr['MJD-OBS'], format='mjd', scale='utc', location=image.site['EarthLocation']) -# image.obstime += 0.5 * image.exptime * u.second # Make time centre of exposure -# image.photfilter = {'i_SDSS': 'ip'}.get(hdr['ESO INS FILT1 NAME'], hdr['ESO INS FILT1 NAME']) -# -# elif instrument == 'ANDICAM-CCD' and hdr.get('OBSERVAT') == 'CTIO': -# image.site = api.get_site( -# 20) # Hard-coded the siteid for ANDICAM at Cerro Tololo Interamerican Observatory (CTIO) -# image.obstime = Time(hdr['JD'], format='jd', scale='utc', location=image.site['EarthLocation']) -# image.obstime += 0.5 * image.exptime * u.second # Make time centre of exposure -# image.photfilter = hdr['CCDFLTID'] -# -# elif telescope == '1.3m PAIRITEL' and instrument == '2MASS Survey cam': -# image.site = api.get_site(21) # Hard-coded the siteid for Peters Automated InfraRed Imaging TELescope -# time_start = Time(hdr['STRT_CPU'], format='iso', scale='utc', location=image.site['EarthLocation']) -# time_stop = Time(hdr['STOP_CPU'], format='iso', scale='utc', location=image.site['EarthLocation']) -# image.obstime = time_start + 0.5 * (time_stop - time_start) -# image.photfilter = {'j': 'J', 'h': 'H', 'k': 'K', }.get(hdr['FILTER'], hdr['FILTER']) -# -# # Mask out "halo" of pixels with zero value along edge of image: -# image.mask |= edge_mask(image.image, value=0) -# -# elif (origin == 'OAdM' or origin.startswith('NOAO-IRAF')) and telescope == 'TJO' and instrument in ( -# 'MEIA3', 'MEIA2'): -# image.site = api.get_site( -# 22) # Hard-coded the siteid for Telescopi Joan Oró (TJO) at Observatori Astronòmic del Montsec -# image.obstime = Time(hdr['JD'], format='jd', scale='utc', location=image.site['EarthLocation']) -# image.obstime += 0.5 * image.exptime * u.second # Make time centre of exposure -# image.photfilter = hdr['FILTER'] -# -# else: -# raise RuntimeError("Could not determine origin of image") -# -# # Sanity checks: -# if image.exptime is None: -# raise ValueError("Image exposure time could not be extracted") -# -# # Create masked version of image: -# image.image[image.mask] = np.NaN -# image.clean = np.ma.masked_array(data=image.image, mask=image.mask, copy=False) -# -# return image From cce81483d2c5bdfc0fe538e32bae703b539d4d29 Mon Sep 17 00:00:00 2001 From: Emir Karamehmetoglu Date: Wed, 9 Mar 2022 19:56:46 +0100 Subject: [PATCH 09/10] change tests to spaces --- tests/test_load_image.py | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/tests/test_load_image.py b/tests/test_load_image.py index 0dc3ffa..f480de4 100644 --- a/tests/test_load_image.py +++ b/tests/test_load_image.py @@ -15,15 +15,15 @@ #-------------------------------------------------------------------------------------------------- @pytest.mark.parametrize('fpath,siteid', [ - ['SN2020aatc_K_20201213_495s.fits.gz', 13], - ['ADP.2021-10-15T11_40_06.553.fits.gz', 2], - #['TJO2459406.56826_V_imc.fits.gz', 22], - #['lsc1m009-fa04-20210704-0044-e91_v1.fits.gz', 4], - #['SN2021rcp_59409.931159242_B.fits.gz', 22], - #['SN2021rhu_59465.86130221_B.fits.gz', 22], - #['20200613_SN2020lao_u_stacked_meandiff.fits.gz', 1], - #['2021aess_20220104_K.fits.gz', 5], - #['2021aess_B01_20220207v1.fits.gz', 5], + ['SN2020aatc_K_20201213_495s.fits.gz', 13], + ['ADP.2021-10-15T11_40_06.553.fits.gz', 2], + #['TJO2459406.56826_V_imc.fits.gz', 22], + #['lsc1m009-fa04-20210704-0044-e91_v1.fits.gz', 4], + #['SN2021rcp_59409.931159242_B.fits.gz', 22], + #['SN2021rhu_59465.86130221_B.fits.gz', 22], + #['20200613_SN2020lao_u_stacked_meandiff.fits.gz', 1], + #['2021aess_20220104_K.fits.gz', 5], + #['2021aess_B01_20220207v1.fits.gz', 5], ]) def test_load_image(fpath, siteid): # Get list of all available filters: From 7d264eceff87e8db0980f8f6db4742f2d7edf6f0 Mon Sep 17 00:00:00 2001 From: Emir Karamehmetoglu Date: Thu, 10 Mar 2022 11:05:44 +0100 Subject: [PATCH 10/10] Update tests, remove repeat code --- flows/load_image.py | 46 +--------------------------------------- tests/test_load_image.py | 40 +++++++++++++++++----------------- 2 files changed, 21 insertions(+), 65 deletions(-) diff --git a/flows/load_image.py b/flows/load_image.py index 76282fa..bc72a25 100644 --- a/flows/load_image.py +++ b/flows/load_image.py @@ -528,48 +528,6 @@ def get_obstime(self): 'PairTel': PairTel, 'TJO': TJO} -def edge_mask(img, value=0): - """ - Create boolean mask of given value near edge of image. - - Parameters: - img (ndarray): Image of - value (float): Value to detect near edge. Default=0. - - Returns: - ndarray: Pixel mask with given values on the edge of image. - - .. codeauthor:: Rasmus Handberg - """ - - mask1 = (img == value) - mask = np.zeros_like(img, dtype='bool') - - # Mask entire rows and columns which are only the value: - mask[np.all(mask1, axis=1), :] = True - mask[:, np.all(mask1, axis=0)] = True - - # Detect "uneven" edges column-wise in image: - a = np.argmin(mask1, axis=0) - b = np.argmin(np.flipud(mask1), axis=0) - for col in range(img.shape[1]): - if mask1[0, col]: - mask[:a[col], col] = True - if mask1[-1, col]: - mask[-b[col]:, col] = True - - # Detect "uneven" edges row-wise in image: - a = np.argmin(mask1, axis=1) - b = np.argmin(np.fliplr(mask1), axis=1) - for row in range(img.shape[0]): - if mask1[row, 0]: - mask[row, :a[row]] = True - if mask1[row, -1]: - mask[row, -b[row]:] = True - - return mask - - def load_image(filename: str, target_coord: typing.Union[coords.SkyCoord, typing.Tuple[float, float]] = None): """ Load FITS image using FlowsImage class and Instrument Classes. @@ -581,10 +539,8 @@ def load_image(filename: str, target_coord: typing.Union[coords.SkyCoord, typing for all other images it is ignored. Returns: - FlowsImage: instance of FlowsImage with valuues populated based on instrument. + FlowsImage: instance of FlowsImage with values populated based on instrument. - .. codeauthor:: Emir Karamehmetoglu - .. codeauthor:: Rasmus Handberg """ ext = 0 # Default extension in HDUList, individual instruments may override this. mask = None # Instrument can override, default is to only mask all non-finite values, override is additive. diff --git a/tests/test_load_image.py b/tests/test_load_image.py index f480de4..4616687 100644 --- a/tests/test_load_image.py +++ b/tests/test_load_image.py @@ -10,24 +10,15 @@ import os.path import conftest # noqa: F401 from tendrils import api -from flows.load_image import load_image +from flows.load_image import load_image, FlowsImage, instruments +# Get list of all available filters: +ALL_FILTERS = set(api.get_filters().keys()) -#-------------------------------------------------------------------------------------------------- @pytest.mark.parametrize('fpath,siteid', [ ['SN2020aatc_K_20201213_495s.fits.gz', 13], - ['ADP.2021-10-15T11_40_06.553.fits.gz', 2], - #['TJO2459406.56826_V_imc.fits.gz', 22], - #['lsc1m009-fa04-20210704-0044-e91_v1.fits.gz', 4], - #['SN2021rcp_59409.931159242_B.fits.gz', 22], - #['SN2021rhu_59465.86130221_B.fits.gz', 22], - #['20200613_SN2020lao_u_stacked_meandiff.fits.gz', 1], - #['2021aess_20220104_K.fits.gz', 5], - #['2021aess_B01_20220207v1.fits.gz', 5], -]) + ['ADP.2021-10-15T11_40_06.553.fits.gz', 2],]) def test_load_image(fpath, siteid): - # Get list of all available filters: - all_filters = set(api.get_filters().keys()) # The test input directory containing the test-images: INPUT_DIR = os.path.join(os.path.abspath(os.path.dirname(__file__)), 'input') @@ -45,18 +36,27 @@ def test_load_image(fpath, siteid): # Check the attributes of the image object: assert isinstance(img.image, np.ndarray) assert img.image.dtype in ('float32', 'float64') - assert isinstance(img.mask, np.ndarray) - assert img.mask.dtype == 'bool' + if img.mask is not None: + assert isinstance(img.mask, np.ndarray) + assert img.mask.dtype == 'bool' assert isinstance(img.clean, np.ma.MaskedArray) assert img.clean.dtype == img.image.dtype - assert isinstance(img.obstime, Time) + assert isinstance(img.exptime, float) - assert img.exptime > 0 + assert img.exptime > 0. assert isinstance(img.wcs, WCS) - assert isinstance(img.site, dict) - assert img.site['siteid'] == siteid + assert isinstance(img.photfilter, str) - assert img.photfilter in all_filters + assert img.photfilter in ALL_FILTERS + + +def test_instruments(): + for instrument_name, instrument_class in instruments.items(): + instrument = instrument_class() + # get site: + site = api.get_site(instrument.siteid) + + assert site['siteid'] == instrument.siteid # --------------------------------------------------------------------------------------------------