diff --git a/stellarphot/core.py b/stellarphot/core.py index 5265800c..525b0cb5 100644 --- a/stellarphot/core.py +++ b/stellarphot/core.py @@ -28,6 +28,9 @@ class Camera: The dark current of the camera in units such that, when multiplied by exposure time, the unit matches the unit of the `read_noise`. + pixel_scale : `astropy.quantity.Quantity` + The pixel scale of the camera in units of arcseconds per pixel. + Attributes ---------- @@ -42,6 +45,9 @@ class Camera: The dark current of the camera in units such that, when multiplied by exposure time, the unit matches the unit of the `read_noise`. + pixel_scale : `astropy.quantity.Quantity` + The pixel scale of the camera in units of arcseconds per pixel. + Notes ----- The gain, read noise, and dark current are all assumed to be constant @@ -53,18 +59,22 @@ class Camera: >>> from stellarphot import Camera >>> camera = Camera(gain=1.0 * u.electron / u.adu, ... read_noise=1.0 * u.electron, - ... dark_current=0.01 * u.electron / u.second) + ... dark_current=0.01 * u.electron / u.second, + ... pixel_scale=0.563 * u.arcsec / u.pixel) >>> camera.gain >>> camera.read_noise >>> camera.dark_current + >>> camera.pixel_scale + """ def __init__(self, gain=1.0 * u.electron / u.adu, read_noise=1.0 * u.electron, - dark_current=0.01 * u.electron / u.second): + dark_current=0.01 * u.electron / u.second, + pixel_scale=1.0 * u.arcsec / u.pixel): super().__init__() # Check that input has units and if not crap out. @@ -82,11 +92,16 @@ def __init__(self, gain=1.0 * u.electron / u.adu, self.dark_current = dark_current else: raise TypeError("dark_current must have units.") + if isinstance(pixel_scale, u.Quantity): + self.pixel_scale = pixel_scale + else: + raise TypeError("pixel_scale must have units.") def copy(self): return Camera(gain=self.gain, read_noise=self.read_noise, - dark_current=self.dark_current) + dark_current=self.dark_current, + pixel_scale=self.pixel_scale) def __copy__(self): return self.copy() @@ -104,13 +119,6 @@ class BaseEnhancedTable(QTable): Parameters ---------- - table_description: dict or dict-like (Default: None) - This is a dictionary where each key is a required table column name - and the value corresponding to each key is the required dtype - (can be None). This is used to check the format of the input data - table. Columns will be output in the order of the keys in the dictionary - with any additional columns tacked on the end. - input_data: `astropy.table.QTable` (Default: None) A table containing astronomical data of interest. This table will be checked to make sure all columns listed in table_description @@ -120,6 +128,13 @@ class BaseEnhancedTable(QTable): validation will not only affect the data attribute of the instance, the original input data is left unchanged. + table_description: dict or dict-like (Default: None) + This is a dictionary where each key is a required table column name + and the value corresponding to each key is the required dtype + (can be None). This is used to check the format of the input data + table. Columns will be output in the order of the keys in the dictionary + with any additional columns tacked on the end. + colname_map: dict, optional (Default: None) A dictionary containing old column names as keys and new column names as values. This is used to automatically update the column @@ -134,12 +149,12 @@ class BaseEnhancedTable(QTable): validation of the inputs is performed and the resulting table is returned. """ - def __init__(self, table_description=None, input_data=None, colname_map=None, + def __init__(self, *args, input_data=None, table_description=None, colname_map=None, **kwargs): if (table_description is None) and (input_data is None): # Assume user is trying to create an empty table and let QTable # handle it - super().__init__(**kwargs) + super().__init__(*args, **kwargs) else: # Confirm a proper table description is passed (that is dict-like with keys # and values) @@ -186,7 +201,6 @@ def __init__(self, table_description=None, input_data=None, colname_map=None, # Call QTable initializer to finish up super().__init__(data=data, **kwargs) - def _validate_columns(self, data): # Check the format of the data table matches the table_description by # checking each column listed in table_description exists and is the @@ -211,7 +225,6 @@ def _validate_columns(self, data): raise ValueError(f"data['{this_col}'] is missing from input " "data.") - def _update_colnames(self, colname_map, data): # Change column names as desired, done before validating the columns, # which is why we work on _orig_data @@ -241,16 +254,16 @@ class PhotometryData(BaseEnhancedTable): Parameters ---------- + input_data: `astropy.table.QTable`, optional (Default: None) + A table containing all the instrumental aperture photometry results + to be validated. + observatory: `astropy.coordinates.EarthLocation`, optional (Default: None) The location of the observatory. camera: `stellarphot.Camera`, optional (Default: None) A description of the CCD used to perform the photometry. - input_data: `astropy.table.QTable`, optional (Default: None) - A table containing all the instrumental aperture photometry results - to be validated. - colname_map: dict, optional (Default: None) A dictionary containing old column names as keys and new column names as values. This is used to automatically update the column @@ -268,25 +281,11 @@ class PhotometryData(BaseEnhancedTable): Attributes ---------- - gain : float - The gain of the camera in units such that the unit of the product `gain` - times the image data matches the unit of the `read_noise`. - - read_noise : float - The read noise of the camera with units. - - dark_current : float - The dark current of the camera in units such that, when multiplied by - exposure time, the unit matches the unit of the `read_noise`. - - lat: float - The location of the observatory in degrees North. - - lon: float - The location of the observatory in degrees East. + camera: `stellarphot.Camera` + A description of the CCD used to perform the photometry. - height: float - Altitute of the observatory in meters. + observatory: `astropy.coordinates.EarthLocation` + The location of the observatory. Notes ----- @@ -377,21 +376,20 @@ class PhotometryData(BaseEnhancedTable): observatory = None camera = None - def __init__(self, observatory=None, camera=None, input_data=None, + def __init__(self, *args, input_data=None, observatory=None, camera=None, colname_map=None, passband_map=None, retain_user_computed=False, **kwargs): if (observatory is None) and (camera is None) and (input_data is None): - super().__init__(**kwargs) + super().__init__(*args, **kwargs) else: - # Save passband filter name mapping - self._passband_map = passband_map.copy() - # Perform input validation if not isinstance(observatory, EarthLocation): raise TypeError("observatory must be an " - "astropy.coordinates.EarthLocation object.") + "astropy.coordinates.EarthLocation object instead " + f"of type {type(observatory)}.") if not isinstance(camera, Camera): - raise TypeError("camera must be a stellarphot.Camera object.") + raise TypeError("camera must be a stellarphot.Camera object instead " + f"of type {type(camera)}.") # Check the time column is correct format and scale try: @@ -405,9 +403,8 @@ def __init__(self, observatory=None, camera=None, input_data=None, "astropy.time.Time entries.") # Convert input data to QTable (while also checking for required columns) - super().__init__(table_description=self.phot_descript, - input_data=input_data, colname_map=colname_map, - **kwargs) + super().__init__(input_data=input_data, table_description=self.phot_descript, + colname_map=colname_map, **kwargs) # Add the TableAttributes directly to meta (and adding attribute # functions below) since using TableAttributes results in a @@ -419,6 +416,7 @@ def __init__(self, observatory=None, camera=None, input_data=None, self.meta['gain'] = camera.gain self.meta['read_noise'] = camera.read_noise self.meta['dark_current'] = camera.dark_current + self.meta['pixel_scale'] = camera.pixel_scale # Check for consistency of counts-related columns counts_columns = ['aperture_sum', 'annulus_sum', 'sky_per_pix_avg', @@ -467,7 +465,7 @@ def __init__(self, observatory=None, camera=None, input_data=None, - LocalTime.ymdhms.second*u.s shift = Column(data = delta, name='shift') # Compute MJD at local noon before the evening of this - # observation + # observation. self['night'] = Column(data= np.array((Time(self['date-obs']) + shift).to_value('mjd'), @@ -480,8 +478,10 @@ def __init__(self, observatory=None, camera=None, input_data=None, # Apply the filter/passband name update if passband_map is not None: + self._passband_map = passband_map.copy() self._update_passbands() + def add_bjd_col(self, observatory): """ Returns a astropy column of barycentric Julian date times corresponding to @@ -502,28 +502,16 @@ def add_bjd_col(self, observatory): return Time(time_barycenter + self['exposure'] / 2, scale='tdb') @property - def lat(self): - return self.meta['lat'] - - @property - def lon(self): - return self.meta['lon'] - - @property - def height(self): - return self.meta['height'] - - @property - def gain(self): - return self.meta['gain'] - - @property - def read_noise(self): - return self.meta['read_noise'] + def camera(self): + return Camera(gain=self.meta['gain'], + read_noise=self.meta['read_noise'], + dark_current=self.meta['dark_current'], + pixel_scale=self.meta['pixel_scale']) @property - def dark_current(self): - return self.meta['dark_current'] + def observatory(self): + return EarthLocation(lat=self.meta['lat'], lon=self.meta['lon'], + height=self.meta['height']) class CatalogData(BaseEnhancedTable): """ @@ -559,6 +547,15 @@ class CatalogData(BaseEnhancedTable): AAVSO passband names as values. This is used to automatically update the passband column to AAVSO standard names if desired. + Attributes + ---------- + catalog_name: str + User readable name for the catalog. + + catalog_source: str + User readable designation for the source of the catalog (could be a + URL or a journal reference). + Notes ----- For validation of inputs, you must provide input_data, catalog_name, and @@ -596,11 +593,10 @@ class CatalogData(BaseEnhancedTable): catalog_name = None catalog_source = None - def __init__(self, input_data=None, catalog_name=None, catalog_source=None, - colname_map=None, passband_map=None, - **kwargs): + def __init__(self, *args, input_data=None, catalog_name=None, catalog_source=None, + colname_map=None, passband_map=None, **kwargs): if (input_data is None) and (catalog_name is None) and (catalog_source is None): - super().__init__(**kwargs) + super().__init__(*args, **kwargs) else: self._passband_map = passband_map @@ -617,13 +613,14 @@ def __init__(self, input_data=None, catalog_name=None, catalog_source=None, self.meta['catalog_name'] = str(catalog_name) self.meta['catalog_source'] = str(catalog_source) + # Apply the filter/passband name update + if passband_map is not None: + self._passband_map = passband_map.copy() + self._update_passbands() + else: raise ValueError("You must provide input_data to CatalogData.") - # Apply the filter/passband name update - if passband_map is not None: - self._update_passbands() - @property def catalog_name(self): return self.meta['catalog_name'] @@ -700,9 +697,9 @@ class SourceListData(BaseEnhancedTable): 'ycenter' : u.pix } - def __init__(self, input_data=None, colname_map=None, **kwargs): + def __init__(self, *args, input_data=None, colname_map=None, **kwargs): if (input_data is None): - super().__init__(**kwargs) + super().__init__(*args, **kwargs) else: # Check data before copying to avoid recusive loop and non-QTable # data input. diff --git a/stellarphot/notebooks/photometry/03-photometry-template.ipynb b/stellarphot/notebooks/photometry/03-photometry-template.ipynb index 78976dc1..8ddcba4c 100755 --- a/stellarphot/notebooks/photometry/03-photometry-template.ipynb +++ b/stellarphot/notebooks/photometry/03-photometry-template.ipynb @@ -81,7 +81,10 @@ "inner_annulus = aperture + 15\n", "outer_annulus = inner_annulus + 15\n", "\n", - "feder_cg_16m = Camera(gain=1.5, read_noise=10.0, dark_current=0.01)\n", + "feder_cg_16m = Camera(gain = 1.5 * u.electron / u.adu,\n", + " read_noise = 10.0 * u.electron,\n", + " dark_current=0.01 * u.electron / u.second,\n", + " pixel_scale = 0.563 * u.arcsec / u.pix)\n", "\n", "max_adu = 50000" ] diff --git a/stellarphot/photometry/source_detection.py b/stellarphot/photometry/source_detection.py index a5c55710..dc84386f 100644 --- a/stellarphot/photometry/source_detection.py +++ b/stellarphot/photometry/source_detection.py @@ -286,7 +286,7 @@ def source_detection(ccd, fwhm=8, sigma=3.0, iters=5, 'xcentroid' : 'xcenter', 'ycentroid' : 'ycenter'} - sl_data = SourceListData(sources, colname_map=colnamemap) + sl_data = SourceListData(input_data=sources, colname_map=colnamemap) return sl_data diff --git a/stellarphot/tests/test_core.py b/stellarphot/tests/test_core.py index 0a758e62..143784e8 100644 --- a/stellarphot/tests/test_core.py +++ b/stellarphot/tests/test_core.py @@ -15,18 +15,43 @@ def test_camera_attributes(): gain = 2.0 * u.electron / u.adu read_noise = 10 * u.electron dark_current = 0.01 * u.electron / u.second - c = Camera(gain=gain, read_noise=read_noise, dark_current=dark_current) + pixel_scale = 0.563 * u.arcsec / u.pix + c = Camera(gain=gain, + read_noise=read_noise, + dark_current=dark_current, + pixel_scale=pixel_scale) assert c.gain == gain assert c.dark_current == dark_current assert c.read_noise == read_noise + assert c.pixel_scale == pixel_scale def test_camera_unitscheck(): - gain = 2.0 - read_noise = 10 - dark_current = 0.01 + gain = 2.0 * u.electron / u.adu + read_noise = 10 * u.electron + dark_current = 0.01 * u.electron / u.second + pixel_scale = 0.563 * u.arcsec / u.pix + with pytest.raises(TypeError): - c = Camera(gain=gain, read_noise=read_noise, dark_current=dark_current) + c = Camera(gain=gain.value, + read_noise=read_noise, + dark_current=dark_current, + pixel_scale=pixel_scale) + with pytest.raises(TypeError): + c = Camera(gain=gain, + read_noise=read_noise.value, + dark_current=dark_current, + pixel_scale=pixel_scale) + with pytest.raises(TypeError): + c = Camera(gain=gain, + read_noise=read_noise, + dark_current=dark_current.value, + pixel_scale=pixel_scale) + with pytest.raises(TypeError): + c = Camera(gain=gain, + read_noise=read_noise, + dark_current=dark_current, + pixel_scale=pixel_scale.value) # Create several test descriptions for use in base_enhanced_table tests. @@ -54,10 +79,57 @@ def test_camera_unitscheck(): # Define some configuration information assuming Feder telescope feder_cg_16m = Camera(gain = 1.5 * u.electron / u.adu, read_noise = 10.0 * u.electron, - dark_current=0.01 * u.electron / u.second) + dark_current=0.01 * u.electron / u.second, + pixel_scale = 0.563 * u.arcsec / u.pix) feder_passbands = {'up':'SU', 'gp':'SG', 'rp':'SR', 'zp':'SZ', 'ip':'SI'} feder_obs = EarthLocation(lat = 46.86678,lon=-96.45328, height=311) + +def test_base_enhanced_table_blank(): + # This should just return a blank BaseEnhancedTable + test_base = BaseEnhancedTable() + assert isinstance(test_base, BaseEnhancedTable) + assert len(test_base) == 0 + + +def test_base_enhanced_table_from_existing_table(): + # Should create a populated dataset properly and display the astropy data + test_base2 = BaseEnhancedTable(table_description=test_descript, input_data=testdata) + assert len(test_base2['ra']) == 1 + assert len(test_base2['dec']) == 1 + + +def test_base_enhanced_table_missing_column(): + # Should raise exception because the RA data is missing from input data + testdata_nora = testdata.copy() + testdata_nora.remove_column('ra') + with pytest.raises(ValueError): + test_base = BaseEnhancedTable(table_description=test_descript, + input_data=testdata_nora) + + +def test_base_enhanced_table_missing_badunits(): + # This will fail due to RA being in units of hours + bad_ra_descript = test_descript.copy() + bad_ra_descript[1,2] = u.hr + + with pytest.raises(ValueError): + test_base = BaseEnhancedTable(table_description=bad_ra_descript, + input_data=testdata) + + +def test_base_enhanced_table_recursive(): + # Should create a populated dataset properly and display the astropy data + test_base2 = BaseEnhancedTable(table_description=test_descript, input_data=testdata) + assert len(test_base2['ra']) == 1 + assert len(test_base2['dec']) == 1 + + # Attempt recursive call + with pytest.raises(TypeError): + test_base3 = BaseEnhancedTable(table_description=test_descript, + input_data=test_base2) + + # Define a realistic table of photometry data (a bit corrupted) photdata = np.array([[1, 2049.145245206124, 2054.0849947477964, 109070.60831212997, 154443.9371254444, 78.17278712191924, 22.505771480719375, @@ -137,65 +209,11 @@ def test_camera_unitscheck(): for this_col in computed_columns: del testphot_clean[this_col] -# Load test catalog -test_cat = ascii.read(get_pkg_data_filename('data/test_vsx_table.ecsv'), format='ecsv', - fast_reader=False) - -# Load test apertures -test_sl_data = ascii.read(get_pkg_data_filename('data/test_sourcelist.ecsv'), - format='ecsv', - fast_reader=False) - - -def test_base_enhanced_table_blank(): - # This should just return a blank BaseEnhancedTable - test_base = BaseEnhancedTable() - assert isinstance(test_base, BaseEnhancedTable) - assert len(test_base) == 0 - - -def test_base_enhanced_table_from_existing_table(): - # Should create a populated dataset properly and display the astropy data - test_base2 = BaseEnhancedTable(table_description=test_descript, input_data=testdata) - assert len(test_base2['ra']) == 1 - assert len(test_base2['dec']) == 1 - - -def test_base_enhanced_table_missing_column(): - # Should raise exception because the RA data is missing from input data - testdata_nora = testdata.copy() - testdata_nora.remove_column('ra') - with pytest.raises(ValueError): - test_base = BaseEnhancedTable(table_description=test_descript, - input_data=testdata_nora) - - -def test_base_enhanced_table_missing_badunits(): - # This will fail due to RA being in units of hours - bad_ra_descript = test_descript.copy() - bad_ra_descript[1,2] = u.hr - - with pytest.raises(ValueError): - test_base = BaseEnhancedTable(table_description=bad_ra_descript, - input_data=testdata) - - -def test_base_enhanced_table_recursive(): - # Should create a populated dataset properly and display the astropy data - test_base2 = BaseEnhancedTable(table_description=test_descript, input_data=testdata) - assert len(test_base2['ra']) == 1 - assert len(test_base2['dec']) == 1 - - # Attempt recursive call - with pytest.raises(TypeError): - test_base3 = BaseEnhancedTable(table_description=test_descript, - input_data=test_base2) - def test_photometry_blank(): # This should just return a blank PhotometryData test_base = PhotometryData() - assert type(test_base) == PhotometryData + assert isinstance(test_base, PhotometryData) assert len(test_base) == 0 @@ -205,15 +223,16 @@ def test_photometry_data(): passband_map=feder_passbands, input_data=testphot_clean) # Check some aspects of that data are sound - assert phot_data.gain == 1.5 * u.electron / u.adu - assert phot_data.read_noise == 10.0 * u.electron - assert phot_data.dark_current == 0.01 * u.electron / u.second - assert phot_data.lat.value == 46.86678 - assert phot_data.lat.unit == u.deg - assert phot_data.lon.value == -96.45328 - assert phot_data.lon.unit == u.deg - assert round(phot_data.height.value) == 311 - assert phot_data.height.unit == u.m + assert phot_data.camera.gain == 1.5 * u.electron / u.adu + assert phot_data.camera.read_noise == 10.0 * u.electron + assert phot_data.camera.dark_current == 0.01 * u.electron / u.second + assert phot_data.camera.pixel_scale == 0.563 * u.arcsec / u.pix + assert phot_data.observatory.lat.value == 46.86678 + assert phot_data.observatory.lat.unit == u.deg + assert phot_data.observatory.lon.value == -96.45328 + assert phot_data.observatory.lon.unit == u.deg + assert round(phot_data.observatory.height.value) == 311 + assert phot_data.observatory.height.unit == u.m assert phot_data['night'][0] == 59909 # Checking the BJD computation against Ohio State online calculator for @@ -228,6 +247,25 @@ def test_photometry_data(): assert (phot_data['bjd'][0].value - 2459910.775405664)*86400 < 0.05 +def test_photometry_slicing(): + # Create photometry data instance + phot_data = PhotometryData(observatory=feder_obs, camera=feder_cg_16m, + passband_map=feder_passbands, input_data=testphot_clean) + + # Test slicing works as expected, leaving attributes intact + just_cols = phot_data[['ra','dec']] + assert just_cols.camera.gain == 1.5 * u.electron / u.adu + assert just_cols.camera.read_noise == 10.0 * u.electron + assert just_cols.camera.dark_current == 0.01 * u.electron / u.second + assert just_cols.camera.pixel_scale == 0.563 * u.arcsec / u.pix + assert just_cols.observatory.lat.value == 46.86678 + assert just_cols.observatory.lat.unit == u.deg + assert just_cols.observatory.lon.value == -96.45328 + assert just_cols.observatory.lon.unit == u.deg + assert round(just_cols.observatory.height.value) == 311 + assert just_cols.observatory.height.unit == u.m + + def test_photometry_recursive(): # Create photometry data instance phot_data = PhotometryData(observatory=feder_obs, camera=feder_cg_16m, @@ -277,6 +315,11 @@ def test_photometry_inconsistent_computed_col_exists(): assert np.abs(phot_data['snr'][0].value - 46.795229859903905) < 1e-6 +# Load test catalog +test_cat = ascii.read(get_pkg_data_filename('data/test_vsx_table.ecsv'), format='ecsv', + fast_reader=False) + + def test_catalog_missing_col(): # Fails with ValueError due to not having 'ra' column with pytest.raises(ValueError): @@ -326,6 +369,12 @@ def test_catalog_recursive(): catalog_source="Vizier", colname_map=vsx_colname_map) +# Load test apertures +test_sl_data = ascii.read(get_pkg_data_filename('data/test_sourcelist.ecsv'), + format='ecsv', + fast_reader=False) + + def test_sourcelist(): sl_test = SourceListData(input_data=test_sl_data, colname_map=None) assert sl_test['star_id'][0] == 0