diff --git a/README.md b/README.md index 3c8de144..f04b8846 100644 --- a/README.md +++ b/README.md @@ -44,6 +44,7 @@ Required: * scipy * astropy >= 4.0 * h5py +* pyuvdata Optional: diff --git a/ci/min_deps.yaml b/ci/min_deps.yaml index 3e2c08df..82a84144 100644 --- a/ci/min_deps.yaml +++ b/ci/min_deps.yaml @@ -9,4 +9,5 @@ dependencies: - scipy - coverage - pytest-cov + - pyuvdata - setuptools_scm diff --git a/environment.yaml b/environment.yaml index 23a62ea7..13af4a95 100644 --- a/environment.yaml +++ b/environment.yaml @@ -11,6 +11,7 @@ dependencies: - h5py - coverage - pytest-cov + - pyuvdata - setuptools_scm - flake8 - pip diff --git a/pyradiosky/data/mock_hera_text_2458098.27471.txt b/pyradiosky/data/mock_hera_text_2458098.27471.txt new file mode 100644 index 00000000..f7382302 --- /dev/null +++ b/pyradiosky/data/mock_hera_text_2458098.27471.txt @@ -0,0 +1,44 @@ +SOURCE_ID RA_J2000 [deg] Dec_J2000 [deg] Flux [Jy] Frequency [Hz] +src0 28.103507 -32.541308 1.00 150000000.00 +src1 25.734946 -32.672235 1.00 150000000.00 +src2 23.361369 -32.760516 1.00 150000000.00 +src3 22.172001 -32.787771 1.00 150000000.00 +src4 20.983198 -32.804964 1.00 150000000.00 +src5 18.604394 -32.806366 1.00 150000000.00 +src6 16.226085 -32.764718 1.00 150000000.00 +src7 13.852232 -32.679226 1.00 150000000.00 +src8 11.483261 -32.551071 1.00 150000000.00 +src9 28.013547 -31.546839 1.00 150000000.00 +src10 25.670235 -31.675179 1.00 150000000.00 +src11 23.321217 -31.761294 1.00 150000000.00 +src12 18.616585 -31.806457 1.00 150000000.00 +src13 17.440155 -31.791342 1.00 150000000.00 +src14 13.915314 -31.682104 1.00 150000000.00 +src15 11.571614 -31.556511 1.00 150000000.00 +src16 27.927165 -30.552168 1.00 150000000.00 +src17 26.768313 -30.620199 1.00 150000000.00 +src18 25.607956 -30.677942 1.00 150000000.00 +src19 23.283697 -30.762426 1.00 150000000.00 +src20 22.120290 -30.789115 1.00 150000000.00 +src21 18.627961 -30.806786 1.00 150000000.00 +src22 17.463988 -30.791865 1.00 150000000.00 +src23 16.300503 -30.766550 1.00 150000000.00 +src24 13.975997 -30.684803 1.00 150000000.00 +src25 12.815469 -30.628424 1.00 150000000.00 +src26 11.656418 -30.561752 1.00 150000000.00 +src27 27.844011 -29.557348 1.00 150000000.00 +src28 25.548066 -29.680626 1.00 150000000.00 +src29 23.247039 -29.763546 1.00 150000000.00 +src30 18.639216 -29.807115 1.00 150000000.00 +src31 16.335582 -29.767633 1.00 150000000.00 +src32 14.034319 -29.687427 1.00 150000000.00 +src33 11.738024 -29.566848 1.00 150000000.00 +src34 27.764133 -28.562410 1.00 150000000.00 +src35 25.490401 -28.683338 1.00 150000000.00 +src36 23.212815 -28.764195 1.00 150000000.00 +src37 22.071845 -28.790420 1.00 150000000.00 +src38 20.931473 -28.805845 1.00 150000000.00 +src39 18.649710 -28.807197 1.00 150000000.00 +src40 17.509296 -28.793122 1.00 150000000.00 +src41 16.368255 -28.768247 1.00 150000000.00 +src42 12.952166 -28.635639 1.00 150000000.00 diff --git a/pyradiosky/data/single_source.vot b/pyradiosky/data/single_source.vot new file mode 100644 index 00000000..7cb03006 --- /dev/null +++ b/pyradiosky/data/single_source.vot @@ -0,0 +1,635 @@ + + + + + + + GLEAM first year observing parameters + + + + Observation date + + + + + + RA range + + + + + RA range + + + + + Declination + + + + + Number of flagged tiles out of the 128 available + + + + + Calibrator (1) + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
2013-08-0919.5 3.5-55103C444
2013-08-1019.5 3.5-2743C444
2013-08-1819.5 3.5-727Pictor A
2013-08-2219.5 3.5-1333C444
2013-08-2519.5 3.5-4053C444
2013-11-05 0.0 8.0-134Hydra A
2013-11-06 0.0 8.0-404Hydra A
2013-11-07 0.0 8.024Hydra A
2013-11-08 0.0 8.0-554Hydra A
2013-11-11 0.0 8.0194Hydra A
2013-11-12 0.0 8.0-724Hydra A
2013-11-25 0.0 8.0-274Hydra A
2014-03-03 6.016.0-274Hydra A
2014-03-04 6.016.0-134Hydra A
2014-03-06 6.016.024Hydra A
2014-03-08 6.016.0194Hydra A
2014-03-09 6.016.0-724Virgo A
2014-03-16 6.016.0-404Virgo A
2014-03-17 6.016.0-554Hydra A
2014-06-0912.022.0-2783C444
2014-06-1012.022.0-4083C444
2014-06-1112.022.028Hercules A
2014-06-1212.018.5-5583C444
2014-06-1312.019.0-138Centaurus A
2014-06-1412.022.0-729Hercules A
2014-06-1512.022.01813Virgo A
2014-06-1618.522.0-1383C444
2014-06-1818.522.0-5583C444
+ + + + Right ascension (FK5, Equinox=J2000.0) (computed by VizieR, not + part of the original data. The format may include more digits + than the original data because of internal accuracy requirements + in VizieR and across other CDS services) + + + + + Declination (FK5, Equinox=J2000.0) (computed by VizieR, not part + of the original data. The format may include more digits than the + original data because of internal accuracy requirements in VizieR + and across other CDS services) + + + + + GLEAM name (JHHMMSS+DDMMSS) (Name) + + + + + Right Ascension J2000 (RAJ2000) + + + + + Declination J2000 (DEJ2000) + + + + + Peak flux in wide (170-231MHz) image (peak_flux_wide) + + + + + Integrated flux in wide (170-231MHz) image (int_flux_wide) + + + + + Percent error in absolute flux scale - all frequencies + (err_abs_flux_pct) + + + + + Percent error on internal flux scale - all frequencies + (err_fit_flux_pct) + + + + + ? Peak flux in 072-080MHz image (peak_flux_076) + + + + + + ? Integrated flux in 072-080MHz image (int_flux_076) + + + + + + ? Peak flux in 080-088MHz image (peak_flux_084) + + + + + + ? Integrated flux in 080-088MHz image (int_flux_084) + + + + + + ? Peak flux in 088-095MHz image (peak_flux_092) + + + + + + ? Integrated flux in 088-095MHz image (int_flux_092) + + + + + + ? Peak flux in 095-103MHz image (peak_flux_099) + + + + + + ? Integrated flux in 095-103MHz image (int_flux_099) + + + + + + ? Peak flux in 103-111MHz image (peak_flux_107) + + + + + + ? Integrated flux in 103-111MHz image (int_flux_107) + + + + + + ? Peak flux in 111-118MHz image (peak_flux_115) + + + + + + ? Integrated flux in 111-118MHz image (int_flux_115) + + + + + + ? Peak flux in 118-126MHz image (peak_flux_122) + + + + + + ? Integrated flux in 118-126MHz image (int_flux_122) + + + + + + ? Peak flux in 126-134MHz image (peak_flux_130) + + + + + + ? Integrated flux in 126-134MHz image (int_flux_130) + + + + + + ? Peak flux in 139-147MHz image (peak_flux_143) + + + + + + ? Integrated flux in 139-147MHz image (int_flux_143) + + + + + + ? Peak flux in 147-154MHz image (peak_flux_151) + + + + + + ? Integrated flux in 147-154MHz image (int_flux_151) + + + + + + ? Peak flux in 154-162MHz image (peak_flux_158) + + + + + + ? Integrated flux in 154-162MHz image (int_flux_158) + + + + + + ? Peak flux in 162-170MHz image (peak_flux_166) + + + + + + ? Integrated flux in 162-170MHz image (int_flux_166) + + + + + + ? Peak flux in 170-177MHz image (peak_flux_174) + + + + + + ? Integrated flux in 170-177MHz image (int_flux_174) + + + + + + ? Peak flux in 177-185MHz image (peak_flux_181) + + + + + + ? Integrated flux in 177-185MHz image (int_flux_181) + + + + + + ? Peak flux in 185-193MHz image (peak_flux_189) + + + + + + ? Integrated flux in 185-193MHz image (int_flux_189) + + + + + + ? Peak flux in 193-200MHz image (peak_flux_197) + + + + + + ? Integrated flux in 193-200MHz image (int_flux_197) + + + + + + ? Peak flux in 200-208MHz image (peak_flux_204) + + + + + + ? Integrated flux in 200-208MHz image (int_flux_204) + + + + + + ? Peak flux in 208-216MHz image (peak_flux_212) + + + + + + ? Integrated flux in 208-216MHz image (int_flux_212) + + + + + + ? Peak flux in 216-223MHz image (peak_flux_220) + + + + + + ? Integrated flux in 216-223MHz image (int_flux_220) + + + + + + ? Peak flux in 223-231MHz image (peak_flux_227) + + + + + + ? Integrated flux in 223-231MHz image (int_flux_227) + + + + + + ? Fitted spectral index (alpha) + + + + + + ? Fitted 200MHz integrated flux density (int_flux_fit_200) + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
0.000000000.0000000TEST0.000000-30.0000000.000000 1.00000016959630.0000000.0000000.0000000.0000000.0000000.0000000.0000000.0000000.0000000.0000000.0000000.0000000.0000000.0000000.0000000.0000000.0000000.0000000.0000000.0000000.0000000.0000000.0000000.0000000.0000000.0000000.0000000.0000000.0000000.0000000.0000000.0000000.0000000.0000000.0000000.0000000.0000000.0000000.0000000.0000000.0000000.000000
+ + diff --git a/pyradiosky/skymodel.py b/pyradiosky/skymodel.py index ad366932..b1a4edfe 100644 --- a/pyradiosky/skymodel.py +++ b/pyradiosky/skymodel.py @@ -8,12 +8,20 @@ import numpy as np from numpy.lib import recfunctions from scipy.linalg import orthogonal_procrustes as ortho_procr -from astropy.coordinates import Angle, EarthLocation, AltAz +import scipy.io +from astropy.coordinates import Angle, EarthLocation, AltAz, Latitude, Longitude from astropy.time import Time import astropy.units as units from astropy.units import Quantity from astropy.io import votable -import scipy.io +from pyuvdata.uvbase import UVBase +from pyuvdata.parameter import UVParameter + +try: + from pyuvdata.uvbeam.cst_beam import CSTBeam +except ImportError: # pragma: no cover + # backwards compatility for older pyuvdata versions + from pyuvdata.cst_beam import CSTBeam from . import utils as skyutils from . import spherical_coords_transforms as sct @@ -27,6 +35,7 @@ "skymodel_to_array", "array_to_skymodel", "source_cuts", + "read_gleam_catalog", "read_votable_catalog", "read_text_catalog", "read_idl_catalog", @@ -61,7 +70,7 @@ class LunarTopo: # Convert stokes and coherency to Astropy quantities. -class SkyModel(object): +class SkyModel(UVBase): """ Object to hold point source and diffuse models. @@ -80,30 +89,26 @@ class SkyModel(object): ---------- name : array_like of str Unique identifier for each source component, shape (Ncomponents,). - ra : astropy Angle object + ra : astropy Longitude object source RA in J2000 (or ICRS) coordinates, shape (Ncomponents,). - dec : astropy Angle object + dec : astropy Latitude object source Dec in J2000 (or ICRS) coordinates, shape (Ncomponents,). stokes : array_like of float 4 element vector giving the source [I, Q, U, V], shape (4, Nfreqs, Ncomponents). - freq : astropy Quantity - Reference frequencies of flux values, shape (Ncomponents,). spectral_type : str Indicates how fluxes should be calculated at each frequency. - spectral_index : array_like of float - Spectral index of each source, shape (Nfreqs, Ncomponents). - None if spectral_type is not 'spectral_index'. Options: - 'flat' : Flat spectrum. - 'full' : Flux is defined by a saved value at each frequency. - 'subband' : Flux is given at a set of band centers. (TODO) - 'spectral_index' : Flux is given at a reference frequency. (TODO) - rise_lst : array_like of float - Approximate lst (radians) when the source rises, shape (Ncomponents,). Set by - coarse horizon cut in simsetup. Default is nan, meaning the source never rises. - set_lst : array_like of float - Approximate lst (radians) when the source sets, shape (Ncomponents,). - Default is None, meaning the source never sets. + freq_array : astropy Quantity + Array of frequencies that fluxes are provided for, shape (Nfreqs,). + reference_frequency : astropy Quantity + Reference frequencies of flux values, shape (Ncomponents,). + spectral_index : array_like of float + Spectral index of each source, shape (Ncomponents). + None if spectral_type is not 'spectral_index'. pos_tol : float position tolerance in degrees, defaults to minimum float in numpy position tolerance in degrees @@ -114,27 +119,8 @@ class SkyModel(object): Beam amplitude at the source position, shape (4, Nfreqs, Ncomponents). 4 element vector corresponds to [XX, YY, XY, YX] instrumental polarizations. - """ - Ncomponents = None # Number of point source components represented here. - - _Ncomp_attrs = [ - "ra", - "dec", - "coherency_radec", - "coherency_local", - "stokes", - "alt_az", - "rise_lst", - "set_lst", - "freq", - "pos_lmn", - "name", - "horizon_mask", - ] - _scalar_attrs = ["Ncomponents", "time", "pos_tol"] - - _member_funcs = ["coherency_calc", "update_positions"] + """ def __init__( self, @@ -142,85 +128,341 @@ def __init__( ra, dec, stokes, - freq_array, spectral_type, + freq_array=None, + reference_frequency=None, spectral_index=None, - rise_lst=None, - set_lst=None, pos_tol=np.finfo(float).eps, extended_model_group=None, beam_amp=None, ): + # standard angle tolerance: 1 mas in radians. + angle_tol = Angle(1, units.arcsec) + self.future_angle_tol = Angle(1e-3, units.arcsec) - if not isinstance(ra, Angle): - raise ValueError( - "ra must be an astropy Angle object. " "value was: {ra}".format(ra=ra) - ) + self._Ncomponents = UVParameter( + "Ncomponents", description="Number of components", expected_type=int + ) - if not isinstance(dec, Angle): - raise ValueError( - "dec must be an astropy Angle object. " - "value was: {dec}".format(dec=dec) - ) + desc = ( + "Number of frequencies if spectral_type is 'full' or 'subband', " + "1 otherwise." + ) + self._Nfreqs = UVParameter("Nfreqs", description=desc, expected_type=int) + + desc = "Right ascension of components in ICRS coordinates." + self._ra = UVParameter( + "ra", + description=desc, + form=("Ncomponents",), + expected_type=Longitude, + tols=angle_tol, + ) + + desc = "Declination of components in ICRS coordinates." + self._dec = UVParameter( + "dec", + description=desc, + form=("Ncomponents",), + expected_type=Latitude, + tols=angle_tol, + ) + + desc = "Component name." + self._name = UVParameter( + "name", description=desc, form=("Ncomponents",), expected_type=str + ) - self.Ncomponents = ra.size + desc = "Frequency array in Hz, only required if spectral_type is 'full' or 'subband'." + self._freq_array = UVParameter( + "freq_array", + description=desc, + form=("Nfreqs",), + expected_type=Quantity, + required=False, + ) + + desc = "Reference frequency in Hz, only required if spectral_type is 'spectral_index'." + self._reference_frequency = UVParameter( + "reference_frequency", + description=desc, + form=("Ncomponents",), + expected_type=Quantity, + required=False, + ) + + desc = "Component flux per frequency and Stokes parameter" + self._stokes = UVParameter( + "stokes", + description=desc, + form=(4, "Nfreqs", "Ncomponents"), + expected_type=float, + ) + + # The coherency is a 2x2 matrix giving electric field correlation in Jy + self._coherency_radec = UVParameter( + "coherency_radec", + description="Ra/Dec coherency per component", + form=(2, 2, "Nfreqs", "Ncomponents"), + expected_type=np.complex, + ) - self.name = np.atleast_1d(np.asarray(name)) - self.freq_array = np.atleast_1d(freq_array) - self.Nfreqs = self.freq_array.size - self.stokes = np.asarray(stokes) + desc = ( + "Type of spectral flux specification, options are: " + "'full','flat', 'subband', 'spectral_index'." + ) + self._spectral_type = UVParameter( + "spectral_type", + description=desc, + expected_type=str, + acceptable_vals=["full", "flat", "subband", "spectral_index"], + ) + + self._spectral_index = UVParameter( + "spectral_index", + description="Spectral indexm only required if spectral_type is 'spectral_index'.", + form=("Ncomponents",), + expected_type=float, + required=False, + ) + + self._beam_amp = UVParameter( + "beam_amp", + description=( + "Beam amplitude at the source position as a function " + "of instrument polarization and frequency." + ), + form=(4, "Nfreqs", "Ncomponents"), + expected_type=float, + required=False, + ) + + self._extended_model_group = UVParameter( + "extended_model_group", + description=( + "Identifier that groups components of an extended " + "source model. Set to -1 for point sources." + ), + form=("Ncomponents",), + expected_type=int, + required=False, + ) + + desc = "Time for local position calculations." + self._time = UVParameter( + "time", description=desc, expected_type=Time, required=False + ) + + desc = "Telescope Location for local position calculations." + self._telescope_location = UVParameter( + "telescope_location", + description=desc, + expected_type=EarthLocation, + required=False, + ) + if hasmoon: + self._telescope_location.expected_type = (EarthLocation, MoonLocation) + + desc = "Altitude and Azimuth of components in local coordinates." + self._alt_az = UVParameter( + "alt_az", + description=desc, + form=(2, "Ncomponents"), + expected_type=float, + tols=np.finfo(float).eps, + required=False, + ) + + desc = "Position cosines of components in local coordinates." + self._pos_lmn = UVParameter( + "pos_lmn", + description=desc, + form=(3, "Ncomponents"), + expected_type=float, + tols=np.finfo(float).eps, + required=False, + ) + + desc = ( + "Boolean indicator of whether this source is below the horizon." + "True indicates the source is below the horizon." + ) + self._horizon_mask = UVParameter( + "horizon_mask", + description=desc, + form=("Ncomponents",), + expected_type=bool, + required=False, + ) + + # initialize the underlying UVBase properties + super(SkyModel, self).__init__() + + self.name = np.atleast_1d(name) + self.Ncomponents = self.name.size self.ra = np.atleast_1d(ra) self.dec = np.atleast_1d(dec) - self.pos_tol = pos_tol - self.spectral_type = spectral_type - self.beam_amp = beam_amp - self.extended_model_group = extended_model_group - if self.spectral_type == "spectral_index": - self.spectral_index = spectral_index - - self.has_rise_set_lsts = False - if (rise_lst is not None) and (set_lst is not None): - self.rise_lst = np.asarray(rise_lst) - self.set_lst = np.asarray(set_lst) - self.has_rise_set_lsts = True + # handle old parameter order + # (use to be: name, ra, dec, stokes, freq_array spectral_type) + if isinstance(spectral_type, (np.ndarray, list, float, Quantity)): + warnings.warn( + "The input parameters to SkyModel.__init__ have changed. Please " + "update the call.", + category=DeprecationWarning, + ) + freqs_use = spectral_type + spectral_type = freq_array - self.alt_az = np.zeros((2, self.Ncomponents), dtype=float) - self.pos_lmn = np.zeros((3, self.Ncomponents), dtype=float) + if spectral_type == "flat" and np.unique(freqs_use).size == 1: + reference_frequency = np.zeros((self.Ncomponents), dtype=np.float) + reference_frequency.fill(freqs_use[0]) + freq_array = None + else: + freq_array = freqs_use + reference_frequency = None + + self.set_spectral_type_params(spectral_type) + + if freq_array is not None: + if not isinstance(freq_array, (Quantity,)): + warnings.warn( + "In the future, the freq_array will be required to be an " + "astropy Quantity with units that are convertable to Hz. " + "Currently, floats are assumed to be in Hz.", + category=DeprecationWarning, + ) + freq_array = freq_array * units.Hz + self.freq_array = np.atleast_1d(freq_array) + self.Nfreqs = self.freq_array.size + else: + self.Nfreqs = 1 + + if reference_frequency is not None: + if not isinstance(reference_frequency, (Quantity,)): + warnings.warn( + "In the future, the reference_frequency will be required to be an " + "astropy Quantity with units that are convertable to Hz. " + "Currently, floats are assumed to be in Hz.", + category=DeprecationWarning, + ) + reference_frequency = reference_frequency * units.Hz + self.reference_frequency = np.atleast_1d(reference_frequency) - self.horizon_mask = np.zeros(self.Ncomponents).astype( - bool - ) # If true, source component is below horizon. + if spectral_index is not None: + self.spectral_index = np.atleast_1d(spectral_index) + self.stokes = np.asarray(stokes, dtype=np.float) if self.Ncomponents == 1: self.stokes = self.stokes.reshape(4, self.Nfreqs, 1) - # The coherency is a 2x2 matrix giving electric field correlation in Jy - # Multiply by .5 to ensure that Trace sums to I not 2*I - # Shape = (2,2,Ncomponents) - self.coherency_radec = skyutils.stokes_to_coherency(self.stokes) - self.time = None + self.check() - assert np.all( - [ - self.Ncomponents == l - for l in [self.ra.size, self.dec.size, self.stokes.shape[2]] - ] - ), "Inconsistent quantity dimensions." + def set_spectral_type_params(self, spectral_type): + """Set parameters depending on spectral_type.""" + self.spectral_type = spectral_type - def _calc_average_rotation_matrix(self, telescope_location): + if spectral_type == "spectral_index": + self._spectral_index.required = True + self._reference_frequency.required = True + self._Nfreqs.acceptable_vals = [1] + self._freq_array.required = False + elif spectral_type in ["full", "subband"]: + self._freq_array.required = True + self._spectral_index.required = False + self._reference_frequency.required = False + self._Nfreqs.acceptable_vals = None + elif spectral_type == "flat": + self._freq_array.required = False + self._spectral_index.required = False + self._reference_frequency.required = False + self._Nfreqs.acceptable_vals = [1] + + def check(self, check_extra=True, run_check_acceptability=True): """ - Calculate the "average" rotation matrix from RA/Dec to AltAz. + Check that all required parameters are set reasonably. - This gets us close to the right value, then need to calculate a correction - for each source separately. + Check that required parameters exist and have appropriate shapes. + Optionally check if the values are acceptable. Parameters ---------- - telescope_location : astropy EarthLocation or lunarsky MoonLocation object - Location of the telescope. + check_extra : bool + Option to check optional parameters as well as required ones. + run_check_acceptability : bool + Option to check if values in required parameters are acceptable. + + """ + # first make sure the required parameters and forms are set properly + # for the spectral_type + self.set_spectral_type_params(self.spectral_type) + + # make sure only one of freq_array and reference_frequency is defined + if self.freq_array is not None and self.reference_frequency is not None: + raise ValueError( + "Only one of freq_array and reference_frequency can be specified, not both." + ) + + # Run the basic check from UVBase + super(SkyModel, self).check( + check_extra=check_extra, run_check_acceptability=run_check_acceptability + ) + + # make sure freq_array or reference_frequency if present is compatible with Hz + if self.freq_array is not None: + try: + self.freq_array.to("Hz") + except (units.UnitConversionError) as e: + raise ValueError( + "freq_array must have a unit that can be converted to Hz." + ) from e + + if self.reference_frequency is not None: + try: + self.reference_frequency.to("Hz") + except (units.UnitConversionError) as e: + raise ValueError( + "reference_frequency must have a unit that can be converted to Hz." + ) from e + + def __eq__(self, other, check_extra=True): + """Check for equality, check for future equality.""" + # Run the basic __eq__ from UVBase + equal = super(SkyModel, self).__eq__(other, check_extra=check_extra) + + # Issue deprecation warning if ra/decs aren't close to future_angle_tol levels + if not np.allclose(self.ra, other.ra, rtol=0, atol=self.future_angle_tol): + warnings.warn( + "The _ra parameters are not within the future tolerance. " + f"Left is {self.ra}, right is {other.ra}", + category=DeprecationWarning, + ) + + if not np.allclose(self.dec, other.dec, rtol=0, atol=self.future_angle_tol): + warnings.warn( + "The _dec parameters are not within the future tolerance. " + f"Left is {self.dec}, right is {other.dec}", + category=DeprecationWarning, + ) + + if not equal: + warnings.warn( + "Future equality does not pass, probably because the " + "frequencies were not checked in the deprecated equality checking.", + category=DeprecationWarning, + ) + equal = super(SkyModel, self).__eq__(other, check_extra=False) + + return equal + + def _calc_average_rotation_matrix(self): + """ + Calculate the "average" rotation matrix from RA/Dec to AltAz. + + This gets us close to the right value, then need to calculate a correction + for each source separately. Returns ------- @@ -238,7 +480,7 @@ def _calc_average_rotation_matrix(self, telescope_location): y=y_c, z=z_c, obstime=self.time, - location=telescope_location, + location=self.telescope_location, frame="icrs", representation_type="cartesian", ) @@ -257,15 +499,10 @@ def _calc_average_rotation_matrix(self, telescope_location): return R_really_orthogonal - def _calc_rotation_matrix(self, telescope_location): + def _calc_rotation_matrix(self): """ Calculate the true rotation matrix from RA/Dec to AltAz for each component. - Parameters - ---------- - telescope_location : astropy EarthLocation or lunarsky MoonLocation object - Location of the telescope. - Returns ------- array of floats @@ -284,7 +521,7 @@ def _calc_rotation_matrix(self, telescope_location): altaz_vec = sct.r_hat(theta_altaz, phi_altaz) assert altaz_vec.shape == (3, self.Ncomponents) - R_avg = self._calc_average_rotation_matrix(telescope_location) + R_avg = self._calc_average_rotation_matrix() R_exact = np.zeros((3, 3, self.Ncomponents), dtype=np.float) @@ -297,22 +534,17 @@ def _calc_rotation_matrix(self, telescope_location): return R_exact - def _calc_coherency_rotation(self, telescope_location): + def _calc_coherency_rotation(self): """ Calculate the rotation matrix to apply to the RA/Dec coherency to get it into alt/az. - Parameters - ---------- - telescope_location : astropy EarthLocation or lunarsky MoonLocation object - Location of the telescope. - Returns ------- array of floats Rotation matrix that takes the coherency from (RA,Dec) --> (Alt,Az), shape (2, 2, Ncomponents). """ - basis_rotation_matrix = self._calc_rotation_matrix(telescope_location) + basis_rotation_matrix = self._calc_rotation_matrix() # Find mathematical points and vectors for RA/Dec theta_radec = np.pi / 2.0 - self.dec.rad @@ -336,7 +568,7 @@ def _calc_coherency_rotation(self, telescope_location): return coherency_rot_matrix - def coherency_calc(self, telescope_location): + def coherency_calc(self, deprecated_location=None): """ Calculate the local coherency in alt/az basis for this source at a time & location. @@ -346,22 +578,31 @@ def coherency_calc(self, telescope_location): Parameters ---------- - telescope_location : astropy EarthLocation or lunarsky MoonLocation object - location of the telescope. + deprecated_location : EarthLocation object + This keyword is deprecated. It is preserved to maintain backwards + compatibility and sets the EarthLocation on this SkyModel object. Returns ------- array of float - local coherency in alt/az basis, shape (2, 2, Ncomponents) + local coherency in alt/az basis, shape (2, 2, Nfreqs, Ncomponents) """ - if not isinstance(telescope_location, (EarthLocation, MoonLocation)): + if deprecated_location is not None: + warnings.warn( + "Passing telescope_location to SkyModel.coherency_calc is " + "deprecated. Set the telescope_location via SkyModel.update_positions.", + category=DeprecationWarning, + ) + self.update_positions(self.time, deprecated_location) + + if not isinstance(self.telescope_location, (EarthLocation, MoonLocation)): errm = "telescope_location must be an astropy EarthLocation object" if hasmoon: errm += " or a lunarsky MoonLocation object " errm += ". " raise ValueError( - errm + "value was: {al}".format(al=str(telescope_location)) + errm + "value was: {al}".format(al=str(self.telescope_location)) ) Ionly_mask = np.sum(self.stokes[1:, :, :], axis=0) == 0.0 @@ -372,7 +613,7 @@ def coherency_calc(self, telescope_location): if NstokesI < self.Ncomponents: # If there are any polarized sources, do rotation. - rotation_matrix = self._calc_coherency_rotation(telescope_location) + rotation_matrix = self._calc_coherency_rotation() polarized_sources = np.where(~Ionly_mask)[0] @@ -412,7 +653,7 @@ def update_positions(self, time, telescope_location): """ if not isinstance(time, Time): raise ValueError( - "time must be an astropy Time object. " "value was: {t}".format(t=time) + "time must be an astropy Time object. value was: {t}".format(t=time) ) if not isinstance(telescope_location, (EarthLocation, MoonLocation)): @@ -425,22 +666,23 @@ def update_positions(self, time, telescope_location): errm + "value was: {al}".format(al=str(telescope_location)) ) - if self.time == time: # Don't repeat calculations + # Don't repeat calculations + if self.time == time and self.telescope_location == telescope_location: return + self.time = time + self.telescope_location = telescope_location + skycoord_use = SkyCoord(self.ra, self.dec, frame="icrs") - if isinstance(telescope_location, MoonLocation): + if isinstance(self.telescope_location, MoonLocation): source_altaz = skycoord_use.transform_to( - LunarTopo(obstime=time, location=telescope_location) + LunarTopo(obstime=self.time, location=self.telescope_location) ) else: source_altaz = skycoord_use.transform_to( - AltAz(obstime=time, location=telescope_location) + AltAz(obstime=self.time, location=self.telescope_location) ) - time.location = telescope_location - - self.time = time alt_az = np.array([source_altaz.alt.rad, source_altaz.az.rad]) self.alt_az = alt_az @@ -449,6 +691,8 @@ def update_positions(self, time, telescope_location): pos_m = np.cos(alt_az[1, :]) * np.cos(alt_az[0, :]) pos_n = np.sin(alt_az[0, :]) + if self.pos_lmn is None: + self.pos_lmn = np.zeros((3, self.Ncomponents), dtype=float) self.pos_lmn[0, :] = pos_l self.pos_lmn[1, :] = pos_m self.pos_lmn[2, :] = pos_n @@ -456,18 +700,6 @@ def update_positions(self, time, telescope_location): # Horizon mask: self.horizon_mask = self.alt_az[0, :] < 0.0 - def __eq__(self, other): - """Test equality with another sky model.""" - time_check = self.time is None and other.time is None - if not time_check: - time_check = np.isclose(self.time, other.time) - return ( - np.allclose(self.ra.deg, other.ra.deg, atol=self.pos_tol) - and np.allclose(self.stokes, other.stokes) - and np.all(self.name == other.name) - and time_check - ) - def read_healpix_hdf5(hdf5_filename): """ @@ -616,13 +848,13 @@ def healpix_to_sky(hpmap, indices, freqs): stokes[0] = (hpmap.T / skyutils.jy_to_ksr(freq)).T stokes[0] = stokes[0] * astropy_healpix.nside_to_pixel_area(nside) - sky = SkyModel(indices.astype("str"), ra, dec, stokes, freq, "full") + sky = SkyModel(indices.astype("str"), ra, dec, stokes, "full", freq_array=freq) return sky def skymodel_to_array(sky): """ - Make a recarrayof source components from a SkyModel object. + Make a recarray of source components from a SkyModel object. Parameters ---------- @@ -635,20 +867,74 @@ def skymodel_to_array(sky): recarray to turn into a SkyModel object. """ - fieldtypes = ["U10", "f8", "f8", "f8", "f8"] - fieldnames = ["source_id", "ra_j2000", "dec_j2000", "flux_density_I", "frequency"] - fieldshapes = [()] * 3 + [(sky.Nfreqs,)] * 2 + sky.check() + max_name_len = np.max([len(name) for name in sky.name]) + fieldtypes = ["U" + str(max_name_len), "f8", "f8", "f8"] + fieldnames = ["source_id", "ra_j2000", "dec_j2000", "flux_density_I"] + fieldshapes = [()] * 3 + if sky.freq_array is not None: + if sky.spectral_type == "subband": + fieldnames.append("subband_frequency") + else: + fieldnames.append("frequency") + fieldtypes.append("f8") + fieldshapes.extend([(sky.Nfreqs,)] * 2) + elif sky.reference_frequency is not None: + # add frequency field (a copy of reference_frequency) for backwards + # compatibility. + warnings.warn( + "The frequency field is included in the recarray as a copy of the " + "reference_frequency field for backwards compatibility. In future " + "only the reference_frequency will be included.", + category=DeprecationWarning, + ) + fieldnames.extend(["frequency", "reference_frequency"]) + fieldtypes.extend(["f8"] * 2) + fieldshapes.extend([()] * 3) + if sky.spectral_index is not None: + fieldnames.append("spectral_index") + fieldtypes.append("f8") + fieldshapes.append(()) + else: + # flat spectrum, no freq info + fieldshapes.append(()) + + if hasattr(sky, "_rise_lst"): + fieldnames.append("rise_lst") + fieldtypes.append("f8") + fieldshapes.append(()) + + if hasattr(sky, "_set_lst"): + fieldnames.append("set_lst") + fieldtypes.append("f8") + fieldshapes.append(()) dt = np.dtype(list(zip(fieldnames, fieldtypes, fieldshapes))) - if sky.Nfreqs == 1: - sky.freq_array = sky.freq_array[:, None] arr = np.empty(sky.Ncomponents, dtype=dt) arr["source_id"] = sky.name arr["ra_j2000"] = sky.ra.deg arr["dec_j2000"] = sky.dec.deg - arr["flux_density_I"] = sky.stokes[0, :, :].T # Swaps component and frequency axes - arr["frequency"] = sky.freq_array + if sky.freq_array is not None: + if sky.spectral_type == "subband": + arr["subband_frequency"] = sky.freq_array + else: + arr["frequency"] = sky.freq_array + arr["flux_density_I"] = sky.stokes[0, :, :].T + elif sky.reference_frequency is not None: + arr["reference_frequency"] = sky.reference_frequency + arr["frequency"] = sky.reference_frequency + arr["flux_density_I"] = np.squeeze(sky.stokes[0, :, :]) + if sky.spectral_index is not None: + arr["spectral_index"] = sky.spectral_index + else: + # flat spectral type, no freq info + arr["flux_density_I"] = np.squeeze(sky.stokes[0, :, :]) + + if hasattr(sky, "_rise_lst"): + arr["rise_lst"] = sky._rise_lst + if hasattr(sky, "_set_lst"): + arr["set_lst"] = sky._set_lst return arr @@ -667,37 +953,81 @@ def array_to_skymodel(catalog_table): :class:`pyradiosky.SkyModel` """ - ra = Angle(catalog_table["ra_j2000"], units.deg) - dec = Angle(catalog_table["dec_j2000"], units.deg) - ids = catalog_table["source_id"] + ra = Longitude(catalog_table["ra_j2000"], units.deg) + dec = Latitude(catalog_table["dec_j2000"], units.deg) + ids = np.asarray(catalog_table["source_id"]).astype(str) flux_I = np.atleast_1d(catalog_table["flux_density_I"]) - if flux_I.ndim == 1: - flux_I = flux_I[:, None] - stokes = np.pad(np.expand_dims(flux_I, 2), ((0, 0), (0, 0), (0, 3)), "constant").T rise_lst = None set_lst = None - source_freqs = np.atleast_1d(catalog_table["frequency"][0]) + + fieldnames = catalog_table.dtype.names + if "reference_frequency" in fieldnames: + reference_frequency = Quantity( + np.atleast_1d(catalog_table["reference_frequency"]), "hertz" + ) + if "spectral_index" in fieldnames: + spectral_index = np.atleast_1d(catalog_table["spectral_index"]) + spectral_type = "spectral_index" + else: + spectral_type = "flat" + spectral_index = None + freq_array = None + elif "frequency" in fieldnames or "subband_frequency" in fieldnames: + if "frequency" in fieldnames: + freq_array = Quantity(np.atleast_1d(catalog_table["frequency"]), "hertz") + else: + spectral_type = "subband" + freq_array = Quantity( + np.atleast_1d(catalog_table["subband_frequency"]), "hertz" + ) + # freq_array gets copied for every component, so its zeroth axis is + # length Ncomponents. Just take the first one. + freq_array = freq_array[0, :] + if freq_array.size > 1: + if "subband_frequency" not in fieldnames: + spectral_type = "full" + else: + spectral_type = "flat" + reference_frequency = None + spectral_index = None + else: + # flat spectrum, no freq info + spectral_type = "flat" + freq_array = None + reference_frequency = None + spectral_index = None if "rise_lst" in catalog_table.dtype.names: rise_lst = catalog_table["rise_lst"] set_lst = catalog_table["set_lst"] - spectral_type = "flat" - if source_freqs.size > 1: - spectral_type = "full" - sourcelist = SkyModel( + if freq_array is not None: + flux_I = flux_I.T + + if flux_I.ndim == 1: + flux_I = flux_I[np.newaxis, :] + + # pad non-I flux values with zeros + stokes = np.pad(flux_I[np.newaxis, :, :], ((0, 3), (0, 0), (0, 0))) + + skymodel = SkyModel( ids, ra, dec, stokes, - source_freqs, spectral_type, - rise_lst=rise_lst, - set_lst=set_lst, + freq_array=freq_array, + reference_frequency=reference_frequency, + spectral_index=spectral_index, ) - return sourcelist + if rise_lst is not None: + skymodel._rise_lst = rise_lst + if set_lst is not None: + skymodel._set_lst = set_lst + + return skymodel def source_cuts( @@ -706,6 +1036,7 @@ def source_cuts( horizon_buffer=0.04364, min_flux=None, max_flux=None, + freq_range=None, ): """ Perform flux and horizon selections on recarray of source components. @@ -714,7 +1045,9 @@ def source_cuts( ---------- catalog_table : recarray recarray of source catalog information. Must have the columns: - 'source_id', 'ra_j2000', 'dec_j2000', 'flux_density_I', 'frequency' + 'source_id', 'ra_j2000', 'dec_j2000', 'flux_density_I' + may also have the colums: + 'frequency' or 'reference_frequency' latitude_deg : float Latitude of telescope in degrees. Used to estimate rise/set lst. horizon_buffer : float @@ -730,6 +1063,9 @@ def source_cuts( Minimum stokes I flux to select [Jy] max_flux : float Maximum stokes I flux to select [Jy] + freq_range : astropy Quantity + Frequency range over which the min and max flux tests should be performed. + Must be length 2. If None, use the range over which the object is defined. Returns ------- @@ -742,14 +1078,74 @@ def source_cuts( lat_rad = np.radians(latitude_deg) buff = horizon_buffer - if min_flux: - catalog_table = catalog_table[catalog_table["flux_density_I"] > min_flux] - - if max_flux: - catalog_table = catalog_table[catalog_table["flux_density_I"] < max_flux] - - ra = Angle(catalog_table["ra_j2000"], units.deg) - dec = Angle(catalog_table["dec_j2000"], units.deg) + if freq_range is not None: + if not isinstance(freq_range, (Quantity,)): + raise ValueError("freq_range must be an astropy Quantity.") + if not np.atleast_1d(freq_range).size == 2: + raise ValueError("freq_range must have 2 elements.") + + fieldnames = catalog_table.dtype.names + if min_flux or max_flux: + if "reference_frequency" in fieldnames: + if "spectral_index" in fieldnames: + raise NotImplementedError( + "Flux cuts with spectral index type objects is not supported yet." + ) + else: + # flat spectral type + if min_flux: + catalog_table = catalog_table[ + catalog_table["flux_density_I"] > min_flux + ] + if max_flux: + catalog_table = catalog_table[ + catalog_table["flux_density_I"] < max_flux + ] + elif "frequency" in fieldnames or "subband_frequency" in fieldnames: + if "frequency" in fieldnames: + freq_array = Quantity( + np.atleast_1d(catalog_table["frequency"]), "hertz" + ) + else: + freq_array = Quantity( + np.atleast_1d(catalog_table["subband_frequency"]), "hertz" + ) + # freq_array gets copied for every component, so its zeroth axis is + # length Ncomponents. Just take the first one. + freq_array = freq_array[0, :] + if freq_range is not None: + freqs_inds_use = np.where( + (freq_array >= np.min(freq_range)) + & (freq_array <= np.max(freq_range)) + )[0] + if freqs_inds_use.size == 0: + raise ValueError("No frequencies in freq_range.") + else: + freqs_inds_use = np.arange(freq_array.size) + flux_data = np.atleast_1d(catalog_table["flux_density_I"]) + if flux_data.ndim > 1: + flux_data = flux_data[:, freqs_inds_use] + if min_flux: + row_inds = np.where(np.min(flux_data, axis=1) > min_flux) + catalog_table = catalog_table[row_inds] + flux_data = flux_data[row_inds] + if max_flux: + catalog_table = catalog_table[ + np.where(np.max(flux_data, axis=1) < max_flux) + ] + else: + # flat spectral type + if min_flux: + catalog_table = catalog_table[ + catalog_table["flux_density_I"] > min_flux + ] + if max_flux: + catalog_table = catalog_table[ + catalog_table["flux_density_I"] < max_flux + ] + + ra = Longitude(catalog_table["ra_j2000"], units.deg) + dec = Latitude(catalog_table["dec_j2000"], units.deg) if coarse_horizon_cut: tans = np.tan(lat_rad) * np.tan(dec.rad) @@ -777,22 +1173,87 @@ def source_cuts( catalog_table, ["rise_lst", "set_lst"], [rise_lst, set_lst], usemask=False ) + if len(catalog_table) == 0: + warnings.warn("All sources eliminated by cuts.") + return catalog_table -def read_votable_catalog(gleam_votable, source_select_kwds={}, return_table=False): +def _get_matching_fields( + name_to_match, name_list, exclude_start_pattern=None, brittle=True +): + match_list = [name for name in name_list if name_to_match.lower() in name.lower()] + if len(match_list) > 1: + # try requiring exact match + match_list_temp = [ + name for name in match_list if name_to_match.lower() == name.lower() + ] + if len(match_list_temp) == 1: + match_list = match_list_temp + elif exclude_start_pattern is not None: + # try excluding columns which start with exclude_start_pattern + match_list_temp = [ + name + for name in match_list + if not name.startswith(exclude_start_pattern) + ] + if len(match_list_temp) == 1: + match_list = match_list_temp + if len(match_list) > 1: + if brittle: + raise ValueError( + f"More than one match for {name_to_match} in {name_list}." + ) + else: + return match_list + elif len(match_list) == 0: + if brittle: + raise ValueError(f"No match for {name_to_match} in {name_list}.") + else: + return None + return match_list[0] + + +def read_votable_catalog( + votable_file, + table_name="GLEAM", + id_column="GLEAM", + ra_column="RAJ2000", + dec_column="DEJ2000", + flux_columns="Fintwide", + reference_frequency=200e6 * units.Hz, + freq_array=None, + spectral_index_column=None, + source_select_kwds=None, + return_table=False, +): """ - Create a list of pyradiosky source objects from a votable catalog. + Create a SkyModel object from a votable catalog. Tested on: GLEAM EGC catalog, version 2 Parameters ---------- - gleam_votable: str + votable_file : str Path to votable catalog file. - return_table: bool, optional + table_name : str + Part of expected table name. Should match only one table name in votable_file. + id_column : str + Part of expected ID column. Should match only one column in the table. + ra_column : str + Part of expected RA column. Should match only one column in the table. + dec_column : str + Part of expected Dec column. Should match only one column in the table. + flux_columns : str or list of str + Part of expected Flux column(s). Each one should match only one column in the table. + reference_frequency : astropy Quantity + Reference frequency for flux values, assumed to be the same value for all rows. + freq_array : astropy Quantity + Frequencies corresponding to flux_columns (should be same length). + Required for multiple flux columns. + return_table : bool, optional Whether to return the astropy table instead of a list of Source objects. - source_select_kwds: dict, optional + source_select_kwds : dict, optional Dictionary of keywords for source selection Valid options: * `lst_array`: For coarse RA horizon cuts, lsts used in the simulation [radians] @@ -809,55 +1270,212 @@ def read_votable_catalog(gleam_votable, source_select_kwds={}, return_table=Fals recarray or :class:`pyradiosky.SkyModel` if return_table, recarray of source parameters, otherwise :class:`pyradiosky.SkyModel` instance """ - resources = votable.parse(gleam_votable).resources + parsed_vo = votable.parse(votable_file) - class Found(Exception): - pass + tables = list(parsed_vo.iter_tables()) + table_ids = [table._ID for table in tables] - try: - for rs in resources: - for tab in rs.tables: - if "GLEAM" in tab.array.dtype.names: - raise Found - except Found: - table = tab.to_table() # Convert to astropy Table - - fieldnames = ["GLEAM", "RAJ2000", "DEJ2000", "Fintwide"] - newnames = ["source_id", "ra_j2000", "dec_j2000", "flux_density_I", "frequency"] - data = table[fieldnames] - Nsrcs = len(data) - freq = 200e6 - for t in data.colnames: - i = fieldnames.index(t) - data[t] = data[t] - data[t].name = newnames[i] - data.add_column(table.Column(np.ones(Nsrcs) * freq, name="frequency")) - data = data.as_array().data - - if len(source_select_kwds) > 0: - data = source_cuts(data, **source_select_kwds) + if None not in table_ids: + table_name_use = _get_matching_fields(table_name, table_ids) + table_match = [table for table in tables if table._ID == table_name_use][0] + else: + warnings.warn( + f"File {votable_file} contains tables with no name or ID, Support for " + "such files is deprecated.", + category=DeprecationWarning, + ) + # Find correct table using the field names + tables_match = [] + for table in tables: + id_col_use = _get_matching_fields( + id_column, table.to_table().colnames, brittle=False + ) + if id_col_use is not None: + tables_match.append(table) + if len(tables_match) > 1: + raise ValueError("More than one matching table.") + else: + table_match = tables_match[0] + + # Convert to astropy Table + astropy_table = table_match.to_table() + + # get ID column + id_col_use = _get_matching_fields(id_column, astropy_table.colnames) + + # get RA & Dec columns, if multiple matches, exclude VizieR calculated columns + # which start with an underscore + ra_col_use = _get_matching_fields( + ra_column, astropy_table.colnames, exclude_start_pattern="_" + ) + dec_col_use = _get_matching_fields( + dec_column, astropy_table.colnames, exclude_start_pattern="_" + ) + + if isinstance(flux_columns, (str)): + flux_columns = [flux_columns] + flux_cols_use = [] + for col in flux_columns: + flux_cols_use.append(_get_matching_fields(col, astropy_table.colnames)) + + if len(flux_columns) > 1 and freq_array is None: + raise ValueError("freq_array must be provided for multiple flux columns.") + + if reference_frequency is not None or len(flux_cols_use) == 1: + if reference_frequency is not None: + if not isinstance(reference_frequency, (Quantity,)): + raise ValueError("reference_frequency must be an astropy Quantity.") + reference_frequency = ( + np.array([reference_frequency.value] * len(astropy_table)) + * reference_frequency.unit + ) + if spectral_index_column is not None: + spectral_type = "spectral_index" + spec_index_col_use = _get_matching_fields( + spectral_index_column, astropy_table.colnames + ) + spectral_index = astropy_table[spec_index_col_use].data.data + else: + spectral_type = "flat" + spectral_index = None + else: + spectral_type = "subband" + spectral_index = None + + stokes = np.zeros((4, len(flux_cols_use), len(astropy_table)), dtype=np.float) + for index, col in enumerate(flux_cols_use): + stokes[0, index, :] = astropy_table[col].data.data + + skymodel_obj = SkyModel( + astropy_table[id_col_use].data.data.astype("str"), + Longitude(astropy_table[ra_col_use].data.data, astropy_table[ra_col_use].unit), + Latitude(astropy_table[dec_col_use].data.data, astropy_table[dec_col_use].unit), + stokes, + spectral_type, + freq_array=freq_array, + reference_frequency=reference_frequency, + spectral_index=spectral_index, + ) + + if source_select_kwds is not None: + sky_array = skymodel_to_array(skymodel_obj) + sky_array = source_cuts(sky_array, **source_select_kwds) + if return_table: + return sky_array + + skymodel_obj = array_to_skymodel(sky_array) if return_table: - return data + return skymodel_to_array(skymodel_obj) + + return skymodel_obj + + +def read_gleam_catalog( + votable_file, spectral_type="flat", source_select_kwds=None, return_table=False +): + """ + Create a SkyModel object from the GLEAM votable catalog. + + Tested on: GLEAM EGC catalog, version 2 + + Parameters + ---------- + votable_file : str + Path to GLEAM votable catalog file. + spectral_type : str + One of 'flat', 'subband' or 'spectral_index'. If set to 'flat', the + wide band integrated flux will be used, if set to 'spectral_index' the + fitted flux at 200 MHz will be used for the flux column. + return_table : bool, optional + Whether to return the astropy table instead of a SkyModel object. + source_select_kwds : dict, optional + Dictionary of keywords for source selection Valid options: + + * `lst_array`: For coarse RA horizon cuts, lsts used in the simulation [radians] + * `latitude_deg`: Latitude of telescope in degrees. Used for declination coarse + horizon cut. + * `horizon_buffer`: Angle (float, in radians) of buffer for coarse horizon cut. + Default is about 10 minutes of sky rotation. (See caveats in + :func:`array_to_skymodel` docstring) + * `min_flux`: Minimum stokes I flux to select [Jy] + * `max_flux`: Maximum stokes I flux to select [Jy] + + Returns + ------- + recarray or :class:`pyradiosky.SkyModel` + if return_table, recarray of source parameters, otherwise :class:`pyradiosky.SkyModel` instance + """ + spec_type_list = ["flat", "spectral_index", "subband"] + if spectral_type not in spec_type_list: + raise ValueError( + f"spectral_type {spectral_type} is not an allowed type. " + f"Allowed types are: {spec_type_list}" + ) + + if spectral_type == "flat": + flux_columns = "Fintwide" + reference_frequency = 200e6 * units.Hz + freq_array = None + spectral_index_column = None + elif spectral_type == "spectral_index": + flux_columns = "Fintfit200" + reference_frequency = 200e6 * units.Hz + spectral_index_column = "alpha" + freq_array = None + else: + # fmt: off + flux_columns = ["Fint076", "Fint084", "Fint092", "Fint099", "Fint107", + "Fint115", "Fint122", "Fint130", "Fint143", "Fint151", + "Fint158", "Fint166", "Fint174", "Fint181", "Fint189", + "Fint197", "Fint204", "Fint212", "Fint220", "Fint227"] + freq_array = [76, 84, 92, 99, 107, 115, 122, 130, 143, 151, 158, 166, + 174, 181, 189, 197, 204, 212, 220, 227] + freq_array = np.array(freq_array) * 1e6 * units.Hz + reference_frequency = None + spectral_index_column = None + # fmt: on + + gleam_skymodel = read_votable_catalog( + votable_file, + table_name="GLEAM", + id_column="GLEAM", + ra_column="RAJ2000", + dec_column="DEJ2000", + flux_columns=flux_columns, + freq_array=freq_array, + reference_frequency=reference_frequency, + spectral_index_column=spectral_index_column, + source_select_kwds=source_select_kwds, + return_table=return_table, + ) - return array_to_skymodel(data) + return gleam_skymodel -def read_text_catalog(catalog_csv, source_select_kwds={}, return_table=False): +def read_text_catalog(catalog_csv, source_select_kwds=None, return_table=False): """ Read in a text file of sources. Parameters ---------- catalog_csv: str - Path to tab separated value file with the following expected columns - (For now, all sources are flat spectrum): - + Path to tab separated value file with the following required columns: * `Source_ID`: source name as a string of maximum 10 characters * `ra_j2000`: right ascension at J2000 epoch, in decimal degrees * `dec_j2000`: declination at J2000 epoch, in decimal degrees - * `flux_density_I`: Stokes I flux density in Janskys - * `frequency`: reference frequency (for future spectral indexing) [Hz] + * `Flux [Jy]`: Stokes I flux density in Janskys + + If flux is specified at multiple frequencies (must be the same set for all + components), the frequencies must be included in each column name, + e.g. `Flux at 150 MHz [Jy]`. Recognized units are ('Hz', 'kHz', 'MHz' or 'GHz'): + + If flux is only specified at one reference frequency (can be different per + component), a frequency column should be added (note: assumed to be in Hz): + * `Frequency`: reference frequency [Hz] + + Optionally a spectral index can be specified per component with: + * `Spectral_Index`: spectral index source_select_kwds : dict, optional Dictionary of keywords for source selection. Valid options: @@ -881,25 +1499,115 @@ def read_text_catalog(catalog_csv, source_select_kwds={}, return_table=False): header = [ h.strip() for h in header.split() if not h[0] == "[" ] # Ignore units in header - dt = np.format_parser( - ["U10", "f8", "f8", "f8", "f8"], - ["source_id", "ra_j2000", "dec_j2000", "flux_density_I", "frequency"], - header, - ) + + flux_fields = [colname for colname in header if colname.lower().startswith("flux")] + flux_fields_lower = [colname.lower() for colname in flux_fields] + + header_lower = [colname.lower() for colname in header] + + expected_cols = ["source_id", "ra_j2000", "dec_j2000"] + if "frequency" in header_lower: + if len(flux_fields) != 1: + raise ValueError( + "If frequency column is present, only one flux column allowed." + ) + freq_array = None + expected_cols.extend([flux_fields_lower[0], "frequency"]) + if "spectral_index" in header_lower: + spectral_type = "spectral_index" + expected_cols.append("spectral_index") + freq_array = None + else: + spectral_type = "flat" + n_freqs = 1 + else: + frequencies = [] + for fluxname in flux_fields: + if "Hz" in fluxname: + cst_obj = CSTBeam() + freq = cst_obj.name2freq(fluxname) + frequencies.append(freq) + else: + if len(flux_fields) > 1: + raise ValueError( + "Multiple flux fields, but they do not all contain a frequency." + ) + if len(frequencies) > 0: + n_freqs = len(frequencies) + if "subband" in flux_fields[0]: + spectral_type = "subband" + else: + if len(frequencies) > 1: + spectral_type = "full" + else: + spectral_type = "flat" + # This has a freq_array + expected_cols.extend(flux_fields_lower) + freq_array = np.array(frequencies) * units.Hz + else: + # This is a flat spectrum (no freq info) + n_freqs = 1 + spectral_type = "flat" + freq_array = None + expected_cols.append("flux") + + if expected_cols != header_lower: + raise ValueError( + "Header does not match expectations. Expected columns" + f"are: {expected_cols}, header columns were: {header_lower}" + ) catalog_table = np.genfromtxt( - catalog_csv, autostrip=True, skip_header=1, dtype=dt.dtype + catalog_csv, autostrip=True, skip_header=1, dtype=None, encoding="utf-8" ) catalog_table = np.atleast_1d(catalog_table) - if len(source_select_kwds) > 0: - catalog_table = source_cuts(catalog_table, **source_select_kwds) + col_names = catalog_table.dtype.names + + names = catalog_table[col_names[0]].astype("str") + ras = Longitude(catalog_table[col_names[1]], units.deg) + decs = Latitude(catalog_table[col_names[2]], units.deg) + + stokes = np.zeros((4, n_freqs, len(catalog_table)), dtype=np.float) + for ind in np.arange(n_freqs): + stokes[0, ind, :] = catalog_table[col_names[ind + 3]] + + if "frequency" in header_lower and freq_array is None: + freq_ind = np.where(np.array(header_lower) == "frequency")[0][0] + reference_frequency = catalog_table[col_names[freq_ind]] * units.Hz + if "spectral_index" in header_lower: + si_ind = np.where(np.array(header_lower) == "spectral_index")[0][0] + spectral_index = catalog_table[col_names[si_ind]] + else: + spectral_index = None + else: + reference_frequency = None + spectral_index = None + + skymodel_obj = SkyModel( + names, + ras, + decs, + stokes, + spectral_type, + freq_array=freq_array, + reference_frequency=reference_frequency, + spectral_index=spectral_index, + ) + + if source_select_kwds is not None: + sky_array = skymodel_to_array(skymodel_obj) + sky_array = source_cuts(sky_array, **source_select_kwds) + if return_table: + return sky_array + + skymodel_obj = array_to_skymodel(sky_array) if return_table: - return catalog_table + return skymodel_to_array(skymodel_obj) - return array_to_skymodel(catalog_table) + return skymodel_obj def read_idl_catalog(filename_sav, expand_extended=True): @@ -920,7 +1628,7 @@ def read_idl_catalog(filename_sav, expand_extended=True): :class:`pyradiosky.SkyModel` """ catalog = scipy.io.readsav(filename_sav)["catalog"] - ids = catalog["id"] + ids = catalog["id"].astype(str) ra = catalog["ra"] dec = catalog["dec"] source_freqs = catalog["freq"] @@ -1007,16 +1715,16 @@ def read_idl_catalog(filename_sav, expand_extended=True): spectral_index = np.insert(spectral_index, use_index, src["alpha"]) source_inds = np.insert(source_inds, use_index, np.full(Ncomps, ext)) - ra = Angle(ra, units.deg) - dec = Angle(dec, units.deg) + ra = Longitude(ra, units.deg) + dec = Latitude(dec, units.deg) stokes = stokes[:, np.newaxis, :] # Add frequency axis sourcelist = SkyModel( ids, ra, dec, stokes, - source_freqs, spectral_type="spectral_index", + reference_frequency=Quantity(source_freqs, "hertz"), spectral_index=spectral_index, beam_amp=beam_amp, extended_model_group=extended_model_group, @@ -1024,29 +1732,69 @@ def read_idl_catalog(filename_sav, expand_extended=True): return sourcelist -def write_catalog_to_file(filename, catalog): +def write_catalog_to_file(filename, skymodel): """ Write out a catalog to a text file. - Readable with :meth:`simsetup.read_catalog_text()`. + Readable with :meth:`read_text_catalog()`. Parameters ---------- filename : str Path to output file (string) - catalog : :class:`SkyModel` + skymodel : :class:`SkyModel` The sky model to write to file. """ + header = "SOURCE_ID\tRA_J2000 [deg]\tDec_J2000 [deg]" + format_str = "{}\t{:0.8f}\t{:0.8f}" + if skymodel.reference_frequency is not None: + header += "\tFlux [Jy]\tFrequency [Hz]" + format_str += "\t{:0.8f}\t{:0.8f}" + if skymodel.spectral_index is not None: + header += "\tSpectral_Index" + format_str += "\t{:0.8f}" + elif skymodel.freq_array is not None: + for freq in skymodel.freq_array: + freq_hz_val = freq.to(units.Hz).value + if freq_hz_val > 1e9: + freq_str = "{:g}_GHz".format(freq_hz_val * 1e-9) + elif freq_hz_val > 1e6: + freq_str = "{:g}_MHz".format(freq_hz_val * 1e-6) + elif freq_hz_val > 1e3: + freq_str = "{:g}_kHz".format(freq_hz_val * 1e-3) + else: + freq_str = "{:g}_Hz".format(freq_hz_val) + + if skymodel.spectral_type == "subband": + header += f"\tFlux_subband_{freq_str} [Jy]" + else: + header += f"\tFlux_{freq_str} [Jy]" + format_str += "\t{:0.8f}" + else: + # flat spectral response, no freq info + header += "\tFlux [Jy]" + format_str += "\t{:0.8f}" + + header += "\n" + format_str += "\n" + with open(filename, "w+") as fo: - fo.write( - "SOURCE_ID\tRA_J2000 [deg]\tDec_J2000 [deg]\tFlux [Jy]\tFrequency [Hz]\n" - ) - arr = skymodel_to_array(catalog) + fo.write(header) + arr = skymodel_to_array(skymodel) for src in arr: - srcid, ra, dec, flux_i, freq = src - - fo.write( - "{}\t{:f}\t{:f}\t{:0.2f}\t{:0.2f}\n".format( - srcid, ra, dec, flux_i[0], freq[0] - ) - ) + if skymodel.reference_frequency is not None: + if skymodel.spectral_index is not None: + srcid, ra, dec, flux_i, rfreq, freq, spec_index = src + fo.write( + format_str.format(srcid, ra, dec, flux_i, freq, spec_index) + ) + else: + srcid, ra, dec, flux_i, rfreq, freq = src + fo.write(format_str.format(srcid, ra, dec, flux_i, rfreq)) + elif skymodel.freq_array is not None: + srcid, ra, dec, flux_i, freq = src + fo.write(format_str.format(srcid, ra, dec, *flux_i)) + else: + # flat spectral response, no freq info + srcid, ra, dec, flux_i = src + fo.write(format_str.format(srcid, ra, dec, flux_i)) diff --git a/pyradiosky/tests/test_skymodel.py b/pyradiosky/tests/test_skymodel.py index 0231165b..7e213b39 100644 --- a/pyradiosky/tests/test_skymodel.py +++ b/pyradiosky/tests/test_skymodel.py @@ -3,11 +3,20 @@ # Licensed under the 3-clause BSD License import os +import fileinput import pytest import numpy as np +import warnings from astropy import units -from astropy.coordinates import SkyCoord, EarthLocation, Angle, AltAz +from astropy.coordinates import ( + SkyCoord, + EarthLocation, + Angle, + AltAz, + Longitude, + Latitude, +) from astropy.time import Time, TimeDelta import scipy.io @@ -19,11 +28,111 @@ GLEAM_vot = os.path.join(SKY_DATA_PATH, "gleam_50srcs.vot") -def test_source_zenith_from_icrs(): - """Test single source position at zenith constructed using icrs.""" +@pytest.fixture +def time_location(): array_location = EarthLocation(lat="-30d43m17.5s", lon="21d25m41.9s", height=1073.0) + time = Time("2018-03-01 00:00:00", scale="utc", location=array_location) + return time, array_location + + +@pytest.fixture +def zenith_skycoord(time_location): + time, array_location = time_location + + source_coord = SkyCoord( + alt=Angle(90, unit=units.deg), + az=Angle(0, unit=units.deg), + obstime=time, + frame="altaz", + location=array_location, + ) + icrs_coord = source_coord.transform_to("icrs") + + return icrs_coord + + +@pytest.fixture +def zenith_skymodel(zenith_skycoord): + icrs_coord = zenith_skycoord + + ra = icrs_coord.ra + dec = icrs_coord.dec + + names = "zen_source" + stokes = [1.0, 0, 0, 0] + zenith_source = skymodel.SkyModel(names, ra, dec, stokes, "flat") + + return zenith_source + + +@pytest.fixture +def moonsky(): + pytest.importorskip("lunarsky") + + from lunarsky import MoonLocation, SkyCoord as SkyC + + # Tranquility base + array_location = MoonLocation(lat="00d41m15s", lon="23d26m00s", height=0.0) + + time = Time.now() + zen_coord = SkyC( + alt=Angle(90, unit=units.deg), + az=Angle(0, unit=units.deg), + obstime=time, + frame="lunartopo", + location=array_location, + ) + + icrs_coord = zen_coord.transform_to("icrs") + + ra = icrs_coord.ra + dec = icrs_coord.dec + names = "zen_source" + stokes = [1.0, 0, 0, 0] + zenith_source = skymodel.SkyModel(names, ra, dec, stokes, "flat") + + zenith_source.update_positions(time, array_location) + + yield zenith_source + + +@pytest.fixture +def healpix_data(): + pytest.importorskip("astropy_healpix") + import astropy_healpix + + nside = 32 + npix = astropy_healpix.nside_to_npix(nside) + hp_obj = astropy_healpix.HEALPix(nside=nside) + + frequencies = np.linspace(100, 110, 10) + pixel_area = astropy_healpix.nside_to_pixel_area(nside) + + # Note that the cone search includes any pixels that overlap with the search + # region. With such a low resolution, this returns some slightly different + # results from the equivalent healpy search. Subtracting(0.75 * pixres) from + # the pixel area resolves this discrepancy for the test. + + pixres = hp_obj.pixel_resolution.to("deg").value + ipix_disc = hp_obj.cone_search_lonlat( + 135 * units.deg, 0 * units.deg, radius=(10 - pixres * 0.75) * units.deg + ) + + return { + "nside": nside, + "npix": npix, + "frequencies": frequencies, + "pixel_area": pixel_area, + "ipix_disc": ipix_disc, + } + + +def test_source_zenith_from_icrs(time_location): + """Test single source position at zenith constructed using icrs.""" + time, array_location = time_location + lst = time.sidereal_time("apparent") tee_ra = lst @@ -41,52 +150,235 @@ def test_source_zenith_from_icrs(): ra = icrs_coord.ra dec = icrs_coord.dec - # Check error cases - with pytest.raises(ValueError) as cm: - skymodel.SkyModel("icrs_zen", ra.rad, dec.rad, [1, 0, 0, 0], 1e8, "flat") - assert str(cm.value).startswith( - "ra must be an astropy Angle object. " "value was: 3.14" - ) - - with pytest.raises(ValueError) as cm: - skymodel.SkyModel("icrs_zen", ra, dec.rad, [1, 0, 0, 0], 1e8, "flat") - assert str(cm.value).startswith( - "dec must be an astropy Angle object. " "value was: -0.53" - ) - zenith_source = skymodel.SkyModel("icrs_zen", ra, dec, [1, 0, 0, 0], 1e8, "flat") + zenith_source = skymodel.SkyModel("icrs_zen", ra, dec, [1.0, 0, 0, 0], "flat") zenith_source.update_positions(time, array_location) zenith_source_lmn = zenith_source.pos_lmn.squeeze() assert np.allclose(zenith_source_lmn, np.array([0, 0, 1]), atol=1e-5) -def test_source_zenith(): +def test_source_zenith(time_location, zenith_skymodel): """Test single source position at zenith.""" - time = Time("2018-03-01 00:00:00", scale="utc") + time, array_location = time_location - array_location = EarthLocation(lat="-30d43m17.5s", lon="21d25m41.9s", height=1073.0) + zenith_skymodel.update_positions(time, array_location) + zenith_source_lmn = zenith_skymodel.pos_lmn.squeeze() + assert np.allclose(zenith_source_lmn, np.array([0, 0, 1])) - source_coord = SkyCoord( - alt=Angle(90, unit=units.deg), - az=Angle(0, unit=units.deg), - obstime=time, - frame="altaz", - location=array_location, - ) - icrs_coord = source_coord.transform_to("icrs") + +def test_skymodel_init_errors(zenith_skycoord): + icrs_coord = zenith_skycoord ra = icrs_coord.ra dec = icrs_coord.dec - names = "zen_source" - stokes = [1, 0, 0, 0] - freqs = [1e8] - zenith_source = skymodel.SkyModel(names, ra, dec, stokes, freqs, "flat") + # Check error cases + with pytest.raises( + ValueError, + match=( + "UVParameter _ra is not the appropriate type. Is: float64. " + "Should be: " + ), + ): + skymodel.SkyModel("icrs_zen", ra.rad, dec, [1.0, 0, 0, 0], "flat") + + with pytest.raises( + ValueError, + match=( + "UVParameter _dec is not the appropriate type. Is: float64. " + "Should be: " + ), + ): + skymodel.SkyModel("icrs_zen", ra, dec.rad, [1.0, 0, 0, 0], "flat") + + with pytest.raises( + ValueError, + match=( + "Only one of freq_array and reference_frequency can be specified, not both." + ), + ): + skymodel.SkyModel( + "icrs_zen", + ra, + dec, + [1.0, 0, 0, 0], + "flat", + reference_frequency=[1e8] * units.Hz, + freq_array=[1e8] * units.Hz, + ) - zenith_source.update_positions(time, array_location) - zenith_source_lmn = zenith_source.pos_lmn.squeeze() - assert np.allclose(zenith_source_lmn, np.array([0, 0, 1])) + with pytest.raises( + ValueError, match=("freq_array must have a unit that can be converted to Hz.") + ): + skymodel.SkyModel( + "icrs_zen", ra, dec, [1.0, 0, 0, 0], "flat", freq_array=[1e8] * units.m + ) + + with pytest.raises( + ValueError, + match=("reference_frequency must have a unit that can be converted to Hz."), + ): + skymodel.SkyModel( + "icrs_zen", + ra, + dec, + [1.0, 0, 0, 0], + "flat", + reference_frequency=[1e8] * units.m, + ) + + +def test_skymodel_deprecated(): + """Test that old init works with deprecation.""" + source_new = skymodel.SkyModel( + "Test", + Longitude(12.0 * units.hr), + Latitude(-30.0 * units.deg), + [1.0, 0.0, 0.0, 0.0], + "flat", + reference_frequency=np.array([1e8]) * units.Hz, + ) + + with pytest.warns( + DeprecationWarning, + match="The input parameters to SkyModel.__init__ have changed", + ): + source_old = skymodel.SkyModel( + "Test", + Longitude(12.0 * units.hr), + Latitude(-30.0 * units.deg), + [1.0, 0.0, 0.0, 0.0], + np.array([1e8]) * units.Hz, + "flat", + ) + assert source_new == source_old + + with pytest.warns( + DeprecationWarning, + match="In the future, the reference_frequency will be required to be an astropy Quantity", + ): + source_old = skymodel.SkyModel( + "Test", + Longitude(12.0 * units.hr), + Latitude(-30.0 * units.deg), + [1.0, 0.0, 0.0, 0.0], + "flat", + reference_frequency=np.array([1e8]), + ) + assert source_new == source_old + + source_old = skymodel.SkyModel( + "Test", + Longitude(12.0 * units.hr), + Latitude(-30.0 * units.deg), + [1.0, 0.0, 0.0, 0.0], + "flat", + reference_frequency=np.array([1.5e8]) * units.Hz, + ) + with pytest.warns( + DeprecationWarning, + match=( + "Future equality does not pass, probably because the frequencies " + "were not checked" + ), + ): + assert source_new == source_old + + source_old = skymodel.SkyModel( + "Test", + Longitude(12.0 * units.hr), + Latitude(-30.0 * units.deg + 2e-3 * units.arcsec), + [1.0, 0.0, 0.0, 0.0], + "flat", + reference_frequency=np.array([1e8]) * units.Hz, + ) + with pytest.warns( + DeprecationWarning, + match=("The _dec parameters are not within the future tolerance"), + ): + assert source_new == source_old + + source_old = skymodel.SkyModel( + "Test", + Longitude(Longitude(12.0 * units.hr) + Longitude(2e-3 * units.arcsec)), + Latitude(-30.0 * units.deg), + [1.0, 0.0, 0.0, 0.0], + "flat", + reference_frequency=np.array([1e8]) * units.Hz, + ) + with pytest.warns( + DeprecationWarning, + match=("The _ra parameters are not within the future tolerance"), + ): + assert source_new == source_old + + stokes = np.zeros((4, 2, 1), dtype=np.float) + stokes[0, :, :] = 1.0 + source_new = skymodel.SkyModel( + "Test", + Longitude(12.0 * units.hr), + Latitude(-30.0 * units.deg), + stokes, + "subband", + freq_array=np.array([1e8, 1.5e8]) * units.Hz, + ) + with pytest.warns( + DeprecationWarning, + match="The input parameters to SkyModel.__init__ have changed", + ): + source_old = skymodel.SkyModel( + "Test", + Longitude(12.0 * units.hr), + Latitude(-30.0 * units.deg), + stokes, + np.array([1e8, 1.5e8]) * units.Hz, + "subband", + ) + assert source_new == source_old + + with pytest.warns( + DeprecationWarning, + match="In the future, the freq_array will be required to be an astropy Quantity", + ): + source_old = skymodel.SkyModel( + "Test", + Longitude(12.0 * units.hr), + Latitude(-30.0 * units.deg), + stokes, + "subband", + freq_array=np.array([1e8, 1.5e8]), + ) + assert source_new == source_old + + telescope_location = EarthLocation( + lat="-30d43m17.5s", lon="21d25m41.9s", height=1073.0 + ) + time = Time("2018-01-01 00:00") + with pytest.warns( + DeprecationWarning, + match="Passing telescope_location to SkyModel.coherency_calc is deprecated", + ): + source_new.update_positions(time, telescope_location) + source_new.coherency_calc(telescope_location) + + +def test_update_position_errors(zenith_skymodel, time_location): + time, array_location = time_location + with pytest.raises(ValueError, match=("time must be an astropy Time object.")): + zenith_skymodel.update_positions("2018-03-01 00:00:00", array_location) + + +def test_coherency_calc_errors(): + """Test that correct errors are raised when providing invalid location object.""" + coord = SkyCoord(ra=30.0 * units.deg, dec=40 * units.deg, frame="icrs") + + stokes_radec = [1, -0.2, 0.3, 0.1] + + source = skymodel.SkyModel("test", coord.ra, coord.dec, stokes_radec, "flat") + + with pytest.raises(ValueError, match="telescope_location must be an"): + source.coherency_calc().squeeze() def test_calc_basis_rotation_matrix(): @@ -102,15 +394,14 @@ def test_calc_basis_rotation_matrix(): source = skymodel.SkyModel( "Test", - Angle(12.0 * units.hr), - Angle(-30.0 * units.deg), + Longitude(12.0 * units.hr), + Latitude(-30.0 * units.deg), [1.0, 0.0, 0.0, 0.0], - 1e8, "flat", ) source.update_positions(time, telescope_location) - basis_rot_matrix = source._calc_average_rotation_matrix(telescope_location) + basis_rot_matrix = source._calc_average_rotation_matrix() assert np.allclose(np.matmul(basis_rot_matrix, basis_rot_matrix.T), np.eye(3)) assert np.allclose(np.matmul(basis_rot_matrix.T, basis_rot_matrix), np.eye(3)) @@ -129,15 +420,14 @@ def test_calc_vector_rotation(): source = skymodel.SkyModel( "Test", - Angle(12.0 * units.hr), - Angle(-30.0 * units.deg), + Longitude(12.0 * units.hr), + Latitude(-30.0 * units.deg), [1.0, 0.0, 0.0, 0.0], - 1e8, "flat", ) source.update_positions(time, telescope_location) - coherency_rotation = np.squeeze(source._calc_coherency_rotation(telescope_location)) + coherency_rotation = np.squeeze(source._calc_coherency_rotation()) assert np.isclose(np.linalg.det(coherency_rotation), 1) @@ -192,10 +482,9 @@ def test_polarized_source_visibilities(): source = skymodel.SkyModel( "icrs_zen", - zenith_icrs.ra + raoff, - zenith_icrs.dec + decoff, + Longitude(zenith_icrs.ra + raoff), + Latitude(zenith_icrs.dec + decoff), stokes_radec, - 1e8, "flat", ) @@ -210,7 +499,7 @@ def test_polarized_source_visibilities(): alts[ti] = alt azs[ti] = az - coherency_tmp = source.coherency_calc(array_location).squeeze() + coherency_tmp = source.coherency_calc().squeeze() coherency_matrix_local[:, :, ti] = coherency_tmp zas = np.pi / 2.0 - alts @@ -291,7 +580,7 @@ def test_polarized_source_smooth_visibilities(): stokes_radec = [1, -0.2, 0.3, 0.1] source = skymodel.SkyModel( - "icrs_zen", zenith_icrs.ra, zenith_icrs.dec, stokes_radec, 1e8, "flat" + "icrs_zen", zenith_icrs.ra, zenith_icrs.dec, stokes_radec, "flat" ) coherency_matrix_local = np.zeros([2, 2, ntimes], dtype="complex128") @@ -305,7 +594,7 @@ def test_polarized_source_smooth_visibilities(): alts[ti] = alt azs[ti] = az - coherency_tmp = source.coherency_calc(array_location).squeeze() + coherency_tmp = source.coherency_calc().squeeze() coherency_matrix_local[:, :, ti] = coherency_tmp zas = np.pi / 2.0 - alts @@ -338,117 +627,83 @@ def test_polarized_source_smooth_visibilities(): assert np.all(imag_stokes == 0) -def test_coherency_calc_errors(): - """Test that correct errors are raised when providing invalid location object.""" - coord = SkyCoord(ra=30.0 * units.deg, dec=40 * units.deg, frame="icrs") +def test_read_healpix_hdf5(healpix_data): + m = np.arange(healpix_data["npix"]) + m[healpix_data["ipix_disc"]] = healpix_data["npix"] - 1 - stokes_radec = [1, -0.2, 0.3, 0.1] - - source = skymodel.SkyModel("test", coord.ra, coord.dec, stokes_radec, 1e8, "flat") - time = Time.now() - array_location = None - with pytest.raises(ValueError) as err: - source.update_positions(time, telescope_location=array_location) - assert str(err.value).startswith("telescope_location must be an") - - with pytest.raises(ValueError) as err: - source.coherency_calc(array_location).squeeze() - assert str(err.value).startswith("telescope_location must be an") + indices = np.arange(healpix_data["npix"]) + hpmap, inds, freqs = skymodel.read_healpix_hdf5( + os.path.join(SKY_DATA_PATH, "healpix_disk.hdf5") + ) -class TestHealpixHdf5: - def setup(self): - # Common features of tests - pytest.importorskip("astropy_healpix") - import astropy_healpix - - self.nside = 32 - self.Npix = astropy_healpix.nside_to_npix(self.nside) - hp_obj = astropy_healpix.HEALPix(nside=self.nside) - self.frequencies = np.linspace(100, 110, 10) - self.pixel_area = astropy_healpix.nside_to_pixel_area(self.nside) - # Note that the cone search includes any pixels that overlap with the search - # region. With such a low resolution, this returns some slightly different - # results from the equivalent healpy search. Subtracting(0.75 * pixres) from - # the pixel area resolves this discrepancy for the test. - - pixres = hp_obj.pixel_resolution.to("deg").value - self.ipix_disc = hp_obj.cone_search_lonlat( - 135 * units.deg, 0 * units.deg, radius=(10 - pixres * 0.75) * units.deg - ) + assert np.allclose(hpmap[0, :], m) + assert np.allclose(inds, indices) + assert np.allclose(freqs, healpix_data["frequencies"]) - def test_read_healpix_hdf5(self): - m = np.arange(self.Npix) - m[self.ipix_disc] = self.Npix - 1 - indices = np.arange(self.Npix) +def test_healpix_to_sky(healpix_data): + hpmap, inds, freqs = skymodel.read_healpix_hdf5( + os.path.join(SKY_DATA_PATH, "healpix_disk.hdf5") + ) - hpmap, inds, freqs = skymodel.read_healpix_hdf5( - os.path.join(SKY_DATA_PATH, "healpix_disk.hdf5") - ) + hmap_orig = np.arange(healpix_data["npix"]) + hmap_orig[healpix_data["ipix_disc"]] = healpix_data["npix"] - 1 - assert np.allclose(hpmap[0, :], m) - assert np.allclose(inds, indices) - assert np.allclose(freqs, self.frequencies) + hmap_orig = np.repeat(hmap_orig[None, :], 10, axis=0) + hmap_orig = (hmap_orig.T / skyutils.jy_to_ksr(freqs)).T + hmap_orig = hmap_orig * healpix_data["pixel_area"] + sky = skymodel.healpix_to_sky(hpmap, inds, freqs) - def test_healpix_to_sky(self): - hpmap, inds, freqs = skymodel.read_healpix_hdf5( - os.path.join(SKY_DATA_PATH, "healpix_disk.hdf5") - ) + assert np.allclose(sky.stokes[0], hmap_orig.value) - hmap_orig = np.arange(self.Npix) - hmap_orig[self.ipix_disc] = self.Npix - 1 - hmap_orig = np.repeat(hmap_orig[None, :], 10, axis=0) - hmap_orig = (hmap_orig.T / skyutils.jy_to_ksr(freqs)).T - hmap_orig = hmap_orig * self.pixel_area - sky = skymodel.healpix_to_sky(hpmap, inds, freqs) +def test_units_healpix_to_sky(healpix_data): + hpmap, inds, freqs = skymodel.read_healpix_hdf5( + os.path.join(SKY_DATA_PATH, "healpix_disk.hdf5") + ) + freqs = freqs * units.Hz - assert np.allclose(sky.stokes[0], hmap_orig.value) + brightness_temperature_conv = units.brightness_temperature( + freqs, beam_area=healpix_data["pixel_area"] + ) + stokes = (hpmap.T * units.K).to(units.Jy, brightness_temperature_conv).T + sky = skymodel.healpix_to_sky(hpmap, inds, freqs) - def test_units_healpix_to_sky(self): - hpmap, inds, freqs = skymodel.read_healpix_hdf5( - os.path.join(SKY_DATA_PATH, "healpix_disk.hdf5") - ) - freqs = freqs * units.Hz + assert np.allclose(sky.stokes[0, 0], stokes.value[0]) - brightness_temperature_conv = units.brightness_temperature( - freqs, beam_area=self.pixel_area - ) - stokes = (hpmap.T * units.K).to(units.Jy, brightness_temperature_conv).T - sky = skymodel.healpix_to_sky(hpmap, inds, freqs) - assert np.allclose(sky.stokes[0, 0], stokes.value[0]) +def test_read_write_healpix(healpix_data): + hpmap, inds, freqs = skymodel.read_healpix_hdf5( + os.path.join(SKY_DATA_PATH, "healpix_disk.hdf5") + ) + freqs = freqs * units.Hz + filename = "tempfile.hdf5" + with pytest.raises(ValueError) as verr: + skymodel.write_healpix_hdf5(filename, hpmap, inds[:10], freqs.to("Hz").value) + assert str(verr.value).startswith( + "Need to provide nside if giving a subset of the map." + ) - def test_read_write_healpix(self): - hpmap, inds, freqs = skymodel.read_healpix_hdf5( - os.path.join(SKY_DATA_PATH, "healpix_disk.hdf5") + with pytest.raises(ValueError) as verr: + skymodel.write_healpix_hdf5( + filename, + hpmap, + inds[:10], + freqs.to("Hz").value, + nside=healpix_data["nside"], ) - freqs = freqs * units.Hz - filename = "tempfile.hdf5" - with pytest.raises(ValueError) as verr: - skymodel.write_healpix_hdf5( - filename, hpmap, inds[:10], freqs.to("Hz").value - ) - assert str(verr.value).startswith( - "Need to provide nside if giving a subset of the map." - ) - - with pytest.raises(ValueError) as verr: - skymodel.write_healpix_hdf5( - filename, hpmap, inds[:10], freqs.to("Hz").value, nside=self.nside - ) - assert str(verr.value).startswith("Invalid map shape") + assert str(verr.value).startswith("Invalid map shape") - skymodel.write_healpix_hdf5(filename, hpmap, inds, freqs.to("Hz").value) + skymodel.write_healpix_hdf5(filename, hpmap, inds, freqs.to("Hz").value) - hpmap_new, inds_new, freqs_new = skymodel.read_healpix_hdf5(filename) + hpmap_new, inds_new, freqs_new = skymodel.read_healpix_hdf5(filename) - os.remove(filename) + os.remove(filename) - assert np.allclose(hpmap_new, hpmap) - assert np.allclose(inds_new, inds) - assert np.allclose(freqs_new, freqs.to("Hz").value) + assert np.allclose(hpmap_new, hpmap) + assert np.allclose(inds_new, inds) + assert np.allclose(freqs_new, freqs.to("Hz").value) def test_healpix_import_err(): @@ -525,6 +780,55 @@ def test_healpix_positions(): assert np.isclose(src_lmn[2][ipix], src_n) +@pytest.mark.parametrize("spec_type", ["flat", "subband", "spectral_index", "full"]) +def test_array_to_skymodel_loop(spec_type): + if spec_type == "full": + spectral_type = "subband" + else: + spectral_type = spec_type + + sky = skymodel.read_gleam_catalog(GLEAM_vot, spectral_type=spectral_type) + if spec_type == "full": + sky.spectral_type = "full" + + # This should be removed after pyuvdata PR #790 is merged. + # GLEAM has NaNs for the spectral_index of some sources + # Currently, arrays with NaNs are never equal even if they are equal where + # they are not nan and are nan in the same locations. + # So remove those components for the round trip equality test. + if spectral_type == "spectral_index": + wh_not_nan = np.squeeze(np.argwhere(~np.isnan(sky.spectral_index))) + sky.Ncomponents = wh_not_nan.size + sky.ra = sky.ra[wh_not_nan] + sky.dec = sky.dec[wh_not_nan] + sky.name = sky.name[wh_not_nan] + sky.reference_frequency = sky.reference_frequency[wh_not_nan] + sky.spectral_index = sky.spectral_index[wh_not_nan] + sky.stokes = sky.stokes[:, :, wh_not_nan] + sky.coherency_radec = skyutils.stokes_to_coherency(sky.stokes) + + arr = skymodel.skymodel_to_array(sky) + sky2 = skymodel.array_to_skymodel(arr) + + assert sky == sky2 + + if spec_type == "flat": + # again with no reference_frequency field + reference_frequency = sky.reference_frequency + sky.reference_frequency = None + arr = skymodel.skymodel_to_array(sky) + sky2 = skymodel.array_to_skymodel(arr) + + assert sky == sky2 + + # again with flat & freq_array + sky.freq_array = np.atleast_1d(np.unique(reference_frequency)) + arr = skymodel.skymodel_to_array(sky) + sky2 = skymodel.array_to_skymodel(arr) + + assert sky == sky2 + + def test_param_flux_cuts(): # Check that min/max flux limits in test params work. @@ -537,34 +841,358 @@ def test_param_flux_cuts(): assert np.all(0.2 < sI < 1.5) -def test_point_catalog_reader(): - catfile = os.path.join(SKY_DATA_PATH, "pointsource_catalog.txt") - srcs = skymodel.read_text_catalog(catfile) +@pytest.mark.parametrize( + "spec_type, init_kwargs, cut_kwargs", + [ + ("flat", {}, {}), + ("flat", {"reference_frequency": np.ones(20) * 200e6 * units.Hz}, {}), + ("full", {"freq_array": np.array([1e8, 1.5e8]) * units.Hz}, {}), + ( + "subband", + {"freq_array": np.array([1e8, 1.5e8]) * units.Hz}, + {"freq_range": np.array([0.9e8, 2e8]) * units.Hz}, + ), + ( + "subband", + {"freq_array": np.array([1e8, 1.5e8]) * units.Hz}, + {"freq_range": np.array([1.1e8, 2e8]) * units.Hz}, + ), + ( + "flat", + {"freq_array": np.array([1e8]) * units.Hz}, + {"freq_range": np.array([0.9e8, 2e8]) * units.Hz}, + ), + ], +) +def test_flux_cuts(spec_type, init_kwargs, cut_kwargs): + Nsrcs = 20 - with open(catfile, "r") as fhandle: - header = fhandle.readline() - header = [h.strip() for h in header.split()] - dt = np.format_parser( - ["U10", "f8", "f8", "f8", "f8"], - ["source_id", "ra_j2000", "dec_j2000", "flux_density_I", "frequency"], - header, + minflux = 0.5 + maxflux = 3.0 + + ids = ["src{}".format(i) for i in range(Nsrcs)] + ras = Longitude(np.random.uniform(0, 360.0, Nsrcs), units.deg) + decs = Latitude(np.linspace(-90, 90, Nsrcs), units.deg) + stokes = np.zeros((4, 1, Nsrcs), dtype=np.float) + if spec_type == "flat": + stokes[0, :, :] = np.linspace(minflux, maxflux, Nsrcs) + else: + stokes = np.zeros((4, 2, Nsrcs), dtype=np.float) + stokes[0, 0, :] = np.linspace(minflux, maxflux / 2.0, Nsrcs) + stokes[0, 1, :] = np.linspace(minflux * 2.0, maxflux, Nsrcs) + + skymodel_obj = skymodel.SkyModel(ids, ras, decs, stokes, spec_type, **init_kwargs) + catalog_table = skymodel.skymodel_to_array(skymodel_obj) + + minI_cut = 1.0 + maxI_cut = 2.3 + + cut_sourcelist = skymodel.source_cuts( + catalog_table, + latitude_deg=30.0, + min_flux=minI_cut, + max_flux=maxI_cut, + **cut_kwargs, ) - catalog_table = np.genfromtxt( - catfile, autostrip=True, skip_header=1, dtype=dt.dtype + if "freq_range" in cut_kwargs and np.min( + cut_kwargs["freq_range"] > np.min(init_kwargs["freq_array"]) + ): + assert np.all(cut_sourcelist["flux_density_I"] < maxI_cut) + else: + assert np.all(cut_sourcelist["flux_density_I"] > minI_cut) + assert np.all(cut_sourcelist["flux_density_I"] < maxI_cut) + + +@pytest.mark.parametrize( + "spec_type, init_kwargs, cut_kwargs, error_category, error_message", + [ + ( + "spectral_index", + { + "reference_frequency": np.ones(20) * 200e6 * units.Hz, + "spectral_index": np.ones(20) * 0.8, + }, + {}, + NotImplementedError, + "Flux cuts with spectral index type objects is not supported yet.", + ), + ( + "full", + {"freq_array": np.array([1e8, 1.5e8]) * units.Hz}, + {"freq_range": [0.9e8, 2e8]}, + ValueError, + "freq_range must be an astropy Quantity.", + ), + ( + "subband", + {"freq_array": np.array([1e8, 1.5e8]) * units.Hz}, + {"freq_range": 0.9e8 * units.Hz}, + ValueError, + "freq_range must have 2 elements.", + ), + ( + "subband", + {"freq_array": np.array([1e8, 1.5e8]) * units.Hz}, + {"freq_range": np.array([1.1e8, 1.4e8]) * units.Hz}, + ValueError, + "No frequencies in freq_range.", + ), + ], +) +def test_source_cut_error( + spec_type, init_kwargs, cut_kwargs, error_category, error_message +): + Nsrcs = 20 + + minflux = 0.5 + maxflux = 3.0 + + ids = ["src{}".format(i) for i in range(Nsrcs)] + ras = Longitude(np.random.uniform(0, 360.0, Nsrcs), units.deg) + decs = Latitude(np.linspace(-90, 90, Nsrcs), units.deg) + stokes = np.zeros((4, 1, Nsrcs), dtype=np.float) + stokes = np.zeros((4, 1, Nsrcs), dtype=np.float) + if spec_type == "flat" or spec_type == "spectral_index": + stokes[0, :, :] = np.linspace(minflux, maxflux, Nsrcs) + else: + stokes = np.zeros((4, 2, Nsrcs), dtype=np.float) + stokes[0, 0, :] = np.linspace(minflux, maxflux / 2.0, Nsrcs) + stokes[0, 1, :] = np.linspace(minflux * 2.0, maxflux, Nsrcs) + + skymodel_obj = skymodel.SkyModel(ids, ras, decs, stokes, spec_type, **init_kwargs) + catalog_table = skymodel.skymodel_to_array(skymodel_obj) + + minI_cut = 1.0 + maxI_cut = 2.3 + + with pytest.raises(error_category, match=error_message): + skymodel.source_cuts( + catalog_table, + latitude_deg=30.0, + min_flux=minI_cut, + max_flux=maxI_cut, + **cut_kwargs, + ) + + +def test_circumpolar_nonrising(time_location): + # Check that the source_cut function correctly identifies sources that are circumpolar or + # won't rise. + # Working with an observatory at the HERA latitude. + + time, location = time_location + + Ntimes = 100 + Nras = 20 + Ndecs = 20 + Nsrcs = Nras * Ndecs + + times = time + TimeDelta(np.linspace(0, 1.0, Ntimes), format="jd") + lon = location.lon.deg + ra = np.linspace(lon - 90, lon + 90, Nras) + dec = np.linspace(-90, 90, Ndecs) + + ra, dec = map(np.ndarray.flatten, np.meshgrid(ra, dec)) + ra = Longitude(ra, units.deg) + dec = Latitude(dec, units.deg) + + names = ["src{}".format(i) for i in range(Nsrcs)] + stokes = np.zeros((4, 1, Nsrcs)) + stokes[0, ...] = 1.0 + + sky = skymodel.SkyModel(names, ra, dec, stokes, "flat") + + src_arr = skymodel.skymodel_to_array(sky) + src_arr = skymodel.source_cuts(src_arr, latitude_deg=location.lat.deg) + + # Boolean array identifying nonrising sources that were removed by source_cuts + nonrising = np.array( + [sky.name[ind] not in src_arr["source_id"] for ind in range(Nsrcs)] ) + is_below_horizon = np.zeros(Nsrcs).astype(bool) + is_below_horizon[nonrising] = True + + alts, azs = [], [] + for ti in range(Ntimes): + sky.update_positions(times[ti], location) + alts.append(sky.alt_az[0]) + azs.append(sky.alt_az[1]) + + # Check that sources below the horizon by coarse cut are + # indeed below the horizon. + lst = times[ti].sidereal_time("mean").rad + dt0 = lst - src_arr["rise_lst"] + dt1 = src_arr["set_lst"] - src_arr["rise_lst"] + with warnings.catch_warnings(): + warnings.filterwarnings( + "ignore", message="invalid value encountered", category=RuntimeWarning + ) - assert sorted(srcs.name) == sorted(catalog_table["source_id"]) - assert srcs.ra.deg in catalog_table["ra_j2000"] - assert srcs.dec.deg in catalog_table["dec_j2000"] - assert srcs.stokes[0] in catalog_table["flux_density_I"] + dt0[dt0 < 0] += 2 * np.pi + dt1[dt1 < 0] += 2 * np.pi + + is_below_horizon[~nonrising] = dt0 > dt1 + assert np.all(sky.alt_az[0][is_below_horizon] < 0.0) + + alts = np.degrees(alts) + + # Check that the circumpolar and nonrising sources match expectation + # from the tan(lat) * tan(dec) values. + nonrising_test = np.where(np.all(alts < 0, axis=0))[0] + circumpolar_test = np.where(np.all(alts > 0, axis=0))[0] + + tans = np.tan(location.lat.rad) * np.tan(dec.rad) + nonrising_calc = np.where(tans < -1) + circumpolar_calc = np.where(tans > 1) + assert np.all(circumpolar_calc == circumpolar_test) + assert np.all(nonrising_calc == nonrising_test) + + # Confirm that the source cuts excluded the non-rising sources. + assert np.all(np.where(nonrising)[0] == nonrising_test) + + # check that rise_lst and set_lst get added to object when converted + new_sky = skymodel.array_to_skymodel(src_arr) + assert hasattr(new_sky, "_rise_lst") + assert hasattr(new_sky, "_set_lst") + + # and that it's round tripped + src_arr2 = skymodel.skymodel_to_array(new_sky) + assert src_arr.dtype == src_arr2.dtype + assert len(src_arr) == len(src_arr2) + + for name in src_arr.dtype.names: + if isinstance(src_arr[name][0], (str,)): + assert np.array_equal(src_arr[name], src_arr2[name]) + else: + assert np.allclose(src_arr[name], src_arr2[name], equal_nan=True) + + +@pytest.mark.parametrize( + "name_to_match, name_list, kwargs, result", + [ + ("gle", ["_RAJ2000", "_DEJ2000", "RAJ2000", "DEJ2000", "GLEAM"], {}, "GLEAM"), + ( + "raj2000", + ["_RAJ2000", "_DEJ2000", "RAJ2000", "DEJ2000", "GLEAM"], + {}, + "RAJ2000", + ), + ( + "ra", + ["_RAJ2000", "_DEJ2000", "RAJ2000", "DEJ2000", "GLEAM"], + {"exclude_start_pattern": "_"}, + "RAJ2000", + ), + ( + "foo", + ["_RAJ2000", "_DEJ2000", "RAJ2000", "DEJ2000", "GLEAM"], + {"brittle": False}, + None, + ), + ], +) +def test_get_matching_fields(name_to_match, name_list, kwargs, result): + assert skymodel._get_matching_fields(name_to_match, name_list, **kwargs) == result + + +@pytest.mark.parametrize( + "name_to_match, name_list, error_message", + [ + ( + "j2000", + ["_RAJ2000", "_DEJ2000", "RAJ2000", "DEJ2000", "GLEAM"], + "More than one match for j2000 in", + ), + ( + "foo", + ["_RAJ2000", "_DEJ2000", "RAJ2000", "DEJ2000", "GLEAM"], + "No match for foo in", + ), + ], +) +def test_get_matching_fields_errors(name_to_match, name_list, error_message): + with pytest.raises(ValueError, match=error_message): + skymodel._get_matching_fields(name_to_match, name_list) + + +@pytest.mark.parametrize("spec_type", ["flat", "subband"]) +def test_read_gleam(spec_type): + skymodel_obj = skymodel.read_gleam_catalog(GLEAM_vot, spectral_type=spec_type) + + assert skymodel_obj.Ncomponents == 50 + if spec_type == "subband": + assert skymodel_obj.Nfreqs == 20 # Check cuts - source_select_kwds = {"min_flux": 1.0} - catalog = skymodel.read_text_catalog( - catfile, source_select_kwds=source_select_kwds, return_table=True + source_select_kwds = {"min_flux": 0.5} + cut_catalog = skymodel.read_gleam_catalog( + GLEAM_vot, + spectral_type=spec_type, + source_select_kwds=source_select_kwds, + return_table=True, ) - assert len(catalog) == 2 + + assert len(cut_catalog) < skymodel_obj.Ncomponents + + cut_obj = skymodel.read_gleam_catalog( + GLEAM_vot, spectral_type=spec_type, source_select_kwds=source_select_kwds + ) + + assert len(cut_catalog) == cut_obj.Ncomponents + + source_select_kwds = {"min_flux": 10.0} + with pytest.warns(UserWarning, match="All sources eliminated by cuts."): + skymodel.read_gleam_catalog( + GLEAM_vot, + spectral_type=spec_type, + source_select_kwds=source_select_kwds, + return_table=True, + ) + + +def test_read_gleam_errors(): + with pytest.raises(ValueError, match="spectral_type full is not an allowed type"): + skymodel.read_gleam_catalog(GLEAM_vot, spectral_type="full") + + +def test_read_deprecated_votable(): + votable_file = os.path.join(SKY_DATA_PATH, "single_source.vot") + + with pytest.warns( + DeprecationWarning, + match=( + f"File {votable_file} contains tables with no name or ID, " + "Support for such files is deprecated." + ), + ): + skymodel_obj = skymodel.read_votable_catalog(votable_file) + + assert skymodel_obj.Ncomponents == 1 + + with pytest.raises(ValueError, match=("More than one matching table.")): + skymodel.read_votable_catalog(votable_file, id_column="de") + + +def test_read_votable_errors(): + + # fmt: off + flux_columns = ["Fint076", "Fint084", "Fint092", "Fint099", "Fint107", + "Fint115", "Fint122", "Fint130", "Fint143", "Fint151", + "Fint158", "Fint166", "Fint174", "Fint181", "Fint189", + "Fint197", "Fint204", "Fint212", "Fint220", "Fint227"] + # fmt: on + with pytest.raises( + ValueError, match="freq_array must be provided for multiple flux columns." + ): + skymodel.read_votable_catalog( + GLEAM_vot, flux_columns=flux_columns, reference_frequency=200e6 * units.Hz + ) + + with pytest.raises( + ValueError, match="reference_frequency must be an astropy Quantity." + ): + skymodel.read_votable_catalog(GLEAM_vot, reference_frequency=200e6) def test_idl_catalog_reader(): @@ -587,90 +1215,34 @@ def test_idl_catalog_reader_extended_sources(): assert len(sourcelist.ra) == len(catalog) - len(ext_inds) + sum(ext_Ncomps) -def test_flux_cuts(): - Nsrcs = 20 +def test_point_catalog_reader(): + catfile = os.path.join(SKY_DATA_PATH, "pointsource_catalog.txt") + srcs = skymodel.read_text_catalog(catfile) + with open(catfile, "r") as fhandle: + header = fhandle.readline() + header = [h.strip() for h in header.split()] dt = np.format_parser( ["U10", "f8", "f8", "f8", "f8"], ["source_id", "ra_j2000", "dec_j2000", "flux_density_I", "frequency"], - [], - ) - - minflux = 0.5 - maxflux = 3.0 - - catalog_table = np.recarray(Nsrcs, dtype=dt.dtype) - catalog_table["source_id"] = ["src{}".format(i) for i in range(Nsrcs)] - catalog_table["ra_j2000"] = np.random.uniform(0, 360.0, Nsrcs) - catalog_table["dec_j2000"] = np.linspace(-90, 90, Nsrcs) - catalog_table["flux_density_I"] = np.linspace(minflux, maxflux, Nsrcs) - catalog_table["frequency"] = np.ones(Nsrcs) * 200e6 - - minI_cut = 1.0 - maxI_cut = 2.3 - - cut_sourcelist = skymodel.source_cuts( - catalog_table, latitude_deg=30.0, min_flux=minI_cut, max_flux=maxI_cut + header, ) - assert np.all(cut_sourcelist["flux_density_I"] > minI_cut) - assert np.all(cut_sourcelist["flux_density_I"] < maxI_cut) - -def test_circumpolar_nonrising(): - # Check that the source_cut function correctly identifies sources that are circumpolar or - # won't rise. - # Working with an observatory at the HERA latitude - - lat = -31.0 - lon = 0.0 - - Ntimes = 100 - Nsrcs = 50 - - j2000 = 2451545.0 - times = Time( - np.linspace(j2000 - 0.5, j2000 + 0.5, Ntimes), format="jd", scale="utc" + catalog_table = np.genfromtxt( + catfile, autostrip=True, skip_header=1, dtype=dt.dtype ) - ra = np.zeros(Nsrcs) - dec = np.linspace(-90, 90, Nsrcs) - - ra = Angle(ra, units.deg) - dec = Angle(dec, units.deg) - - coord = SkyCoord(ra=ra, dec=dec, frame="icrs") - alts = [] - azs = [] - - loc = EarthLocation.from_geodetic(lat=lat, lon=lon) - for i in range(Ntimes): - altaz = coord.transform_to(AltAz(obstime=times[i], location=loc)) - alts.append(altaz.alt.deg) - azs.append(altaz.az.deg) - alts = np.array(alts) - - nonrising = np.where(np.all(alts < 0, axis=0))[0] - circumpolar = np.where(np.all(alts > 0, axis=0))[0] - - tans = np.tan(np.radians(lat)) * np.tan(dec.rad) - nonrising_check = np.where(tans < -1) - circumpolar_check = np.where(tans > 1) - assert np.all(circumpolar_check == circumpolar) - assert np.all(nonrising_check == nonrising) - - -def test_read_gleam(): - sourcelist = skymodel.read_votable_catalog(GLEAM_vot) - - assert sourcelist.Ncomponents == 50 + assert sorted(srcs.name) == sorted(catalog_table["source_id"]) + assert srcs.ra.deg in catalog_table["ra_j2000"] + assert srcs.dec.deg in catalog_table["dec_j2000"] + assert srcs.stokes[0] in catalog_table["flux_density_I"] # Check cuts source_select_kwds = {"min_flux": 1.0} - catalog = skymodel.read_votable_catalog( - GLEAM_vot, source_select_kwds=source_select_kwds, return_table=True + catalog = skymodel.read_text_catalog( + catfile, source_select_kwds=source_select_kwds, return_table=True ) - - assert len(catalog) < sourcelist.Ncomponents + assert len(catalog) == 2 def test_catalog_file_writer(): @@ -690,9 +1262,8 @@ def test_catalog_file_writer(): dec = icrs_coord.dec names = "zen_source" - stokes = [1, 0, 0, 0] - freqs = [1e8] - zenith_source = skymodel.SkyModel(names, ra, dec, stokes, freqs, "flat") + stokes = [1.0, 0, 0, 0] + zenith_source = skymodel.SkyModel(names, ra, dec, stokes, "flat") fname = os.path.join(SKY_DATA_PATH, "temp_cat.txt") @@ -702,72 +1273,172 @@ def test_catalog_file_writer(): os.remove(fname) -def test_array_to_skymodel_loop(): - sky = skymodel.read_votable_catalog(GLEAM_vot) - sky.ra = Angle(sky.ra.rad, "rad") - sky.dec = Angle(sky.dec.rad, "rad") - arr = skymodel.skymodel_to_array(sky) - sky2 = skymodel.array_to_skymodel(arr) +@pytest.mark.parametrize("spec_type", ["flat", "subband", "spectral_index", "full"]) +def test_text_catalog_loop(spec_type): + if spec_type == "full": + spectral_type = "subband" + else: + spectral_type = spec_type + + sky = skymodel.read_gleam_catalog(GLEAM_vot, spectral_type=spectral_type) + if spec_type == "full": + sky.spectral_type = "full" + + # This should be removed after pyuvdata PR #790 is merged. + # GLEAM has NaNs for the spectral_index of some sources + # Currently, arrays with NaNs are never equal even if they are equal where + # they are not nan and are nan in the same locations. + # So remove those components for the round trip equality test. + if spectral_type == "spectral_index": + wh_not_nan = np.squeeze(np.argwhere(~np.isnan(sky.spectral_index))) + sky.Ncomponents = wh_not_nan.size + sky.ra = sky.ra[wh_not_nan] + sky.dec = sky.dec[wh_not_nan] + sky.name = sky.name[wh_not_nan] + sky.reference_frequency = sky.reference_frequency[wh_not_nan] + sky.spectral_index = sky.spectral_index[wh_not_nan] + sky.stokes = sky.stokes[:, :, wh_not_nan] + sky.coherency_radec = skyutils.stokes_to_coherency(sky.stokes) - assert np.allclose((sky.ra - sky2.ra).rad, 0.0) - assert np.allclose((sky.dec - sky2.dec).rad, 0.0) + fname = os.path.join(SKY_DATA_PATH, "temp_cat.txt") + skymodel.write_catalog_to_file(fname, sky) + sky2 = skymodel.read_text_catalog(fname) + sky_arr2 = skymodel.read_text_catalog(fname, return_table=True) + sky3 = skymodel.array_to_skymodel(sky_arr2) + os.remove(fname) + assert sky == sky2 + assert sky == sky3 -class TestMoon: - """ - Series of tests for Moon-based observers - """ + if spec_type == "flat": + # again with no reference_frequency field + reference_frequency = sky.reference_frequency + sky.reference_frequency = None + arr = skymodel.skymodel_to_array(sky) + sky2 = skymodel.array_to_skymodel(arr) - def setup(self): - pytest.importorskip("lunarsky") + assert sky == sky2 - from lunarsky import MoonLocation, SkyCoord as SkyC + # again with flat & freq_array + sky.freq_array = np.atleast_1d(np.unique(reference_frequency)) + arr = skymodel.skymodel_to_array(sky) + sky2 = skymodel.array_to_skymodel(arr) - # Tranquility base - self.array_location = MoonLocation(lat="00d41m15s", lon="23d26m00s", height=0.0) + assert sky == sky2 - self.time = Time.now() - self.zen_coord = SkyC( - alt=Angle(90, unit=units.deg), - az=Angle(0, unit=units.deg), - obstime=self.time, - frame="lunartopo", - location=self.array_location, - ) - icrs_coord = self.zen_coord.transform_to("icrs") +@pytest.mark.parametrize("freq_mult", [1e-6, 1e-3, 1e3]) +def test_text_catalog_loop_other_freqs(freq_mult): + sky = skymodel.read_gleam_catalog(GLEAM_vot) + sky.freq_array = np.atleast_1d(np.unique(sky.reference_frequency) * freq_mult) + sky.reference_frequency = None + + fname = os.path.join(SKY_DATA_PATH, "temp_cat.txt") + skymodel.write_catalog_to_file(fname, sky) + sky2 = skymodel.read_text_catalog(fname) + print(sky.spectral_type) + print(sky2.spectral_type) + os.remove(fname) + + assert sky == sky2 + + +@pytest.mark.parametrize("spec_type", ["flat", "subband"]) +def test_read_text_source_cuts(spec_type): + + sky = skymodel.read_gleam_catalog(GLEAM_vot, spectral_type=spec_type) + fname = os.path.join(SKY_DATA_PATH, "temp_cat.txt") + skymodel.write_catalog_to_file(fname, sky) + + source_select_kwds = {"min_flux": 0.5} + cut_catalog = skymodel.read_text_catalog( + fname, source_select_kwds=source_select_kwds, return_table=True, + ) + + assert len(cut_catalog) < sky.Ncomponents + + cut_obj = skymodel.read_text_catalog(fname, source_select_kwds=source_select_kwds) + os.remove(fname) + + assert len(cut_catalog) == cut_obj.Ncomponents + + +def test_pyuvsim_mock_catalog_read(): + mock_cat_file = os.path.join(SKY_DATA_PATH, "mock_hera_text_2458098.27471.txt") + + mock_sky = skymodel.read_text_catalog(mock_cat_file) + expected_names = ["src" + str(val) for val in np.arange(mock_sky.Ncomponents)] + assert mock_sky.name.tolist() == expected_names + - ra = icrs_coord.ra - dec = icrs_coord.dec - names = "zen_source" - stokes = [1, 0, 0, 0] - freqs = [1e8] - self.zenith_source = skymodel.SkyModel(names, ra, dec, stokes, freqs, "flat") +def test_read_text_errors(): + sky = skymodel.read_gleam_catalog(GLEAM_vot, spectral_type="subband") - self.zenith_source.update_positions(self.time, self.array_location) + fname = os.path.join(SKY_DATA_PATH, "temp_cat.txt") + skymodel.write_catalog_to_file(fname, sky) + with fileinput.input(files=fname, inplace=True) as infile: + for line in infile: + line = line.replace("Flux_subband_76_MHz [Jy]", "Frequency [Hz]") + print(line, end="") + + with pytest.raises( + ValueError, + match="If frequency column is present, only one flux column allowed.", + ): + skymodel.read_text_catalog(fname) + + skymodel.write_catalog_to_file(fname, sky) + with fileinput.input(files=fname, inplace=True) as infile: + for line in infile: + line = line.replace("Flux_subband_76_MHz [Jy]", "Flux [Jy]") + print(line, end="") + + with pytest.raises( + ValueError, + match="Multiple flux fields, but they do not all contain a frequency.", + ): + skymodel.read_text_catalog(fname) + + skymodel.write_catalog_to_file(fname, sky) + with fileinput.input(files=fname, inplace=True) as infile: + for line in infile: + line = line.replace("SOURCE_ID", "NAME") + print(line, end="") + + with pytest.raises(ValueError, match="Header does not match expectations."): + skymodel.read_text_catalog(fname) + + os.remove(fname) + + +def test_zenith_on_moon(moonsky): + """Source at zenith from the Moon.""" + + zenith_source = moonsky + zenith_source.check() + + zenith_source_lmn = zenith_source.pos_lmn.squeeze() + assert np.allclose(zenith_source_lmn, np.array([0, 0, 1])) - def test_zenith_on_moon(self): - """Source at zenith from the Moon.""" - zenith_source_lmn = self.zenith_source.pos_lmn.squeeze() - assert np.allclose(zenith_source_lmn, np.array([0, 0, 1])) +def test_source_motion(moonsky): + """ Check that period is about 28 days.""" - def test_source_motion(self): - """ Check that period is about 28 days.""" + zenith_source = moonsky - Ntimes = 500 - ets = np.linspace(0, 4 * 28 * 24 * 3600, Ntimes) - times = self.time + TimeDelta(ets, format="sec") + Ntimes = 500 + ets = np.linspace(0, 4 * 28 * 24 * 3600, Ntimes) + times = zenith_source.time + TimeDelta(ets, format="sec") - lmns = np.zeros((Ntimes, 3)) - for ti in range(Ntimes): - self.zenith_source.update_positions(times[ti], self.array_location) - lmns[ti] = self.zenith_source.pos_lmn.squeeze() - _els = np.fft.fft(lmns[:, 0]) - dt = np.diff(ets)[0] - _freqs = np.fft.fftfreq(Ntimes, d=dt) + lmns = np.zeros((Ntimes, 3)) + for ti in range(Ntimes): + zenith_source.update_positions(times[ti], zenith_source.telescope_location) + lmns[ti] = zenith_source.pos_lmn.squeeze() + _els = np.fft.fft(lmns[:, 0]) + dt = np.diff(ets)[0] + _freqs = np.fft.fftfreq(Ntimes, d=dt) - f_28d = 1 / (28 * 24 * 3600.0) + f_28d = 1 / (28 * 24 * 3600.0) - maxf = _freqs[np.argmax(np.abs(_els[_freqs > 0]) ** 2)] - assert np.isclose(maxf, f_28d, atol=2 / ets[-1]) + maxf = _freqs[np.argmax(np.abs(_els[_freqs > 0]) ** 2)] + assert np.isclose(maxf, f_28d, atol=2 / ets[-1]) diff --git a/pyradiosky/utils.py b/pyradiosky/utils.py index 17f10537..0a0c4dcb 100644 --- a/pyradiosky/utils.py +++ b/pyradiosky/utils.py @@ -82,7 +82,7 @@ def stokes_to_coherency(stokes_vector): Returns ------- coherency matrix : array of float - Array of coherencies, shape (2, 2) or (2, 2, Ncomponents) + Array of coherencies, shape (2, 2) or (2, 2, Nfreqs, Ncomponents) """ stokes_arr = np.atleast_1d(np.asarray(stokes_vector)) initial_shape = stokes_arr.shape diff --git a/setup.py b/setup.py index 7528ba3a..7c2d1293 100644 --- a/setup.py +++ b/setup.py @@ -36,7 +36,7 @@ class MockVersion: "scripts": glob.glob("scripts/*"), "version": version.version, "include_package_data": True, - "install_requires": ["numpy>=1.15", "scipy", "astropy>=4.0", "h5py"], + "install_requires": ["numpy>=1.15", "scipy", "astropy>=4.0", "h5py", "pyuvdata"], "tests_require": ["pytest"], "extras_require": { "healpix": ["astropy-healpix"],