diff --git a/.gitignore b/.gitignore index 82711ac6..7a67bdd4 100644 --- a/.gitignore +++ b/.gitignore @@ -114,3 +114,6 @@ karabo/examples/karabo/test/data/results/* # Other (playground, experiment, ..) .experiment/* data_download/* + +# files that are generated by OSKAR +karabo/data/meerkat.tm/element_pattern* diff --git a/karabo/simulation/coordinate_helper.py b/karabo/simulation/coordinate_helper.py index cc654747..a6cd65a7 100644 --- a/karabo/simulation/coordinate_helper.py +++ b/karabo/simulation/coordinate_helper.py @@ -19,10 +19,10 @@ def east_north_to_long_lat( """ # https://stackoverflow.com/questions/7477003/calculating-new-longitude-latitude-from-old-n-meters - r_earth = 6371000 - new_latitude = lat + (east_relative / r_earth) * (180 / np.pi) - new_longitude = long + (north_relative / r_earth) * (180 / np.pi) / np.cos( - long * np.pi / 180 + r_earth = 6378100 # from astropy (astropy.constants.R_earth.value) + new_latitude = lat + np.rad2deg(east_relative / r_earth) + new_longitude = long + np.rad2deg(north_relative / r_earth) / np.cos( + np.deg2rad(long) ) return new_longitude, new_latitude @@ -31,7 +31,7 @@ def wgs84_to_cartesian( lon: Union[float, NDArray[np.float64]], lat: Union[float, NDArray[np.float64]], alt: Union[float, NDArray[np.float64]], - radius: int = 6371000, + radius: int = 6378100, ) -> NDArray[np.float64]: """Transforms WGS84 to cartesian in meters. diff --git a/karabo/simulation/sky_model.py b/karabo/simulation/sky_model.py index ceffcccd..80f2552a 100644 --- a/karabo/simulation/sky_model.py +++ b/karabo/simulation/sky_model.py @@ -819,14 +819,14 @@ def read_from_file(cls: Type[_TSkyModel], path: str) -> _TSkyModel: if dataframe.shape[1] < 3: raise KaraboSkyModelError( - f"CSV does not have the necessary 3 basic columns (RA, DEC and " + "CSV does not have the necessary 3 basic columns (RA, DEC and " f"STOKES I), but only {dataframe.shape[1]} columns." ) if dataframe.shape[1] > cls.SOURCES_COLS: print( - f"""CSV has {dataframe.shape[1] - cls.SOURCES_COLS + 1} - too many rows. The extra rows will be cut off.""" + f"CSV has {dataframe.shape[1] - cls.SOURCES_COLS + 1} " + "rows too many. The extra rows will be cut off." ) return cls(dataframe) diff --git a/karabo/simulation/telescope.py b/karabo/simulation/telescope.py index e0269def..8307005c 100644 --- a/karabo/simulation/telescope.py +++ b/karabo/simulation/telescope.py @@ -6,7 +6,7 @@ import os import re import shutil -from itertools import product +from itertools import combinations from typing import ( Dict, List, @@ -24,21 +24,17 @@ import pandas as pd from astropy import constants as const from astropy import units as u +from astropy.coordinates import EarthLocation from numpy.typing import NDArray from oskar.telescope import Telescope as OskarTelescope -from rascil.processing_components.simulation.simulation_helpers import ( - plot_configuration, -) +from scipy.spatial.distance import pdist from ska_sdp_datamodels.configuration.config_create import create_named_configuration from ska_sdp_datamodels.configuration.config_model import Configuration from typing_extensions import assert_never import karabo.error from karabo.error import KaraboError -from karabo.simulation.coordinate_helper import ( - east_north_to_long_lat, - wgs84_to_cartesian, -) +from karabo.simulation.coordinate_helper import wgs84_to_cartesian from karabo.simulation.east_north_coordinate import EastNorthCoordinate from karabo.simulation.station import Station from karabo.simulation.telescope_versions import ( @@ -206,6 +202,7 @@ def __init__( Altitude (in meters) at the center of the telescope, default is 0. """ self.path: Optional[DirPathType] = None + self._name: Optional[str] = None self.centre_longitude = longitude self.centre_latitude = latitude self.centre_altitude = altitude @@ -317,7 +314,7 @@ def constructor( ) assert name in get_args(RASCILTelescopes) try: - configuration = create_named_configuration(name) + telescope: Telescope = cls.__convert_to_karabo_telescope(name) except ValueError as e: raise ValueError( f"""Requested telescope {name} is not supported by this backend. @@ -325,19 +322,71 @@ def constructor( https://gitlab.com/ska-telescope/sdp/ska-sdp-datamodels/-/blob/d6dcce6288a7bf6d9ce63ab16e799977723e7ae5/src/ska_sdp_datamodels/configuration/config_create.py""" # noqa ) from e - config_earth_location = configuration.location - telescope = Telescope( - longitude=config_earth_location.lon.to("deg").value, - latitude=config_earth_location.lat.to("deg").value, - altitude=config_earth_location.height.to("m").value, + # Function like _get_station_infos() and + # create_baseline_cut_telescope() need access to an OSKAR telescope + # model (.tm). This is not available for RASCIL datasets. + # Thus, we create a temporary one. + disk_cache = FileHandler().get_tmp_dir( + prefix="telescope-constructor-rascil-", + mkdir=False, ) - telescope.backend = SimulatorBackend.RASCIL - telescope.RASCIL_configuration = configuration + tm_path = os.path.join(disk_cache, f"{name}.tm") + telescope.write_to_disk(tm_path) + telescope.path = tm_path return telescope else: assert_never(backend) + @classmethod + def __convert_to_karabo_telescope(cls, instr_name: str) -> Telescope: + """Converts a site saved in RASCIl data format into a Karabo Telescope. + This function acts as an adapter to make the functionality in Telescope + class work for a RASCIL telescope. Namely the functions max_baseline() + and get_baseline_lengths(). + It derives the necessary data structures from the RASCIL_configuration + and fits them into those of the Telescope class. The resuting class is + a SimulatorBackend.RASCIL but has the stations: List[Station] + list filled as well. Nevertheless, it should only be used as a RASCIL + telescope class. + + :param instr_name: The name of the instrument to convert. + + :returns: An instance of Karabo Telescope. + :rtype: karabo.simulation.telescope.Telecope + :raises: ValueError if instr_name is not a valid RASCIL telescope + """ + config = create_named_configuration(instr_name) + + site_location_gc: EarthLocation = config.location + # this conversion returns complex type with unit + # lon,lat,alt = site_location_gc.geodetic + longitude = site_location_gc.lon.to("deg").value + latitude = site_location_gc.lat.to("deg").value + altitude = site_location_gc.height.to("m").value + + telescope = Telescope(longitude, latitude, altitude) + # This is used in some inteferometer simulations + telescope.RASCIL_configuration = config + + station_coords = config.xyz.data + for i, coord in enumerate(station_coords): + telescope.add_station( + horizontal_x=coord[0], + horizontal_y=coord[1], + horizontal_z=coord[2], + ) + + # there are only stations in the rascil files no antennas. + # we add a dummy antenna in order to avoid the creation + # of an empty file. This matches other files. See + # karabo/data/aca.all.tm/station000/layout.txt for example. + # Reason: Value not set to 0 probably to compensate + # for dish diameter. (see comment for PR #631) + telescope.add_antenna_to_station(i, 0.1, 0.1) + telescope.backend = SimulatorBackend.RASCIL + return telescope + @property def name(self) -> Optional[str]: """Gets the telescope name (if available). @@ -347,10 +396,19 @@ def name(self) -> Optional[str]: Returns: Telescope name or `None`. """ + if self._name is not None: + return self._name if self.path is None: return None return os.path.split(self.path)[-1].split(".")[0] + @name.setter + def name(self, value: str) -> None: + """Sets the name of the telescope. Usually, this is the name + of the telescope model file w/o ending + """ + self._name = value + def get_backend_specific_information(self) -> Union[DirPathType, Configuration]: if self.backend is SimulatorBackend.OSKAR: return self.path @@ -421,7 +479,7 @@ def add_antenna_to_station( :param horizontal_z_coordinate_error: altitude of antenna error :return: """ - if station_index < len(self.stations): + if station_index >= 0 and station_index < len(self.stations): station = self.stations[station_index] station.add_station_antenna( EastNorthCoordinate( @@ -433,6 +491,11 @@ def add_antenna_to_station( horizontal_z_coordinate_error, ) ) + else: + raise IndexError( + "You tried to add an antenna to a station that doesn't exist.\n" + f"station_index must be between 0 and {len(self.stations)-1}" + ) def plot_telescope(self, file: Optional[str] = None) -> None: """ @@ -442,7 +505,9 @@ def plot_telescope(self, file: Optional[str] = None) -> None: if self.backend is SimulatorBackend.OSKAR: self.plot_telescope_OSKAR(file) elif self.backend is SimulatorBackend.RASCIL: - plot_configuration(self.get_backend_specific_information(), plot_file=file) + # we can use plot_telescope_OSKAR here because we converted + # the RASCIl setup into an OSKAR setup when constructing it. + self.plot_telescope_OSKAR(file) else: logging.warning( f"""Backend {self.backend} is not valid. @@ -459,32 +524,23 @@ def plot_telescope_OSKAR( import matplotlib.pyplot as plt fig, ax = plt.subplots() - antenna_x = [] - antenna_y = [] station_x = [] station_y = [] for station in self.stations: station_x.append(station.longitude) station_y.append(station.latitude) - for antenna in station.antennas: - long, lat = east_north_to_long_lat( - antenna.x, antenna.y, station.longitude, station.latitude - ) - antenna_x.append(long) - antenna_y.append(lat) - - ax.scatter(antenna_x, antenna_y, label="Antennas") - ax.scatter(station_x, station_y, label="Stations") + # we set the colour manually in order to keep the colour scheme. + ax.scatter(station_x, station_y, label="Stations", c="tab:orange") x = np.array([self.centre_longitude]) y = np.array([self.centre_latitude]) - ax.scatter(x, y, label="Centre") + ax.scatter(x, y, label="Centre", c="tab:green") ax.ticklabel_format(useOffset=False) ax.set_xlabel("Longitude [deg]") ax.set_ylabel("Latitude [deg]") - ax.set_title("Antenna Locations") + ax.set_title(f"{self.name} Overview") ax.legend(loc="upper left", shadow=False, fontsize="medium") if file is not None: @@ -543,7 +599,7 @@ def __write_position_txt(self, position_file_path: str) -> None: position_file = open(position_file_path, "a") position_file.write( - f"{self.centre_longitude} {self.centre_latitude} {self.centre_altitude} \n" + f"{self.centre_longitude} {self.centre_latitude} {self.centre_altitude}\n" ) position_file.close() @@ -553,8 +609,8 @@ def __write_layout_txt( layout_file = open(layout_path, "a") for element in elements: layout_file.write( - f"{element.x}, {element.y}, {element.z}, {element.x_error}, " - + f"{element.y_error}, {element.z_error} \n" + f"{element.x} {element.y} {element.z} {element.x_error} " + + f"{element.y_error} {element.z_error}\n" ) layout_file.close() @@ -563,6 +619,15 @@ def get_cartesian_position(self) -> NDArray[np.float_]: @classmethod def read_OSKAR_tm_file(cls, path: DirPathType) -> Telescope: + """Reads an OSKAR telescope model from disk and + returns an object of karabo.simulation.telescope.Telescope + + :param path: Path to a valid telescope model (extemsion *.tm) + :return: A karabo.simulation.telescope.Telescope object. Importantn: + The object has the backend set to SimulatorBackend.OSKAR. + :raises: A karabo.error.KaraboError if the path does not exit, + or the data in the file cannot be read. + """ path_ = str(path) abs_station_dir_paths = [] center_position_file = None @@ -578,12 +643,14 @@ def read_OSKAR_tm_file(cls, path: DirPathType) -> Telescope: ) if center_position_file is None: - raise karabo.error.KaraboError("Missing crucial position.txt file_or_dir") + raise karabo.error.KaraboError( + f"Missing crucial position.txt in {file_or_dir}" + ) if station_layout_file is None: raise karabo.error.KaraboError( "Missing layout.txt file in station directory. " - "Only Layout.txt is support. " + "Only layout.txt is supported. " "The layout_ecef.txt and layout_wgs84.txt as " "defined in the OSKAR Telescope .tm specification are not " "supported currently." @@ -716,7 +783,7 @@ def _get_station_infos(cls, tel_path: DirPathType) -> pd.DataFrame: raise KaraboError( f"Stations found in {tel_path} are not ascending from station<0 - n>. " ) - stations = np.loadtxt(os.path.join(tel_path, "layout.txt")) + stations = np.array(cls.__read_layout_txt(os.path.join(tel_path, "layout.txt"))) if (n_stations_layout := stations.shape[0]) != (n_stations := df_tel.shape[0]): raise KaraboError( f"Number of stations mismatch of {n_stations_layout=} & {n_stations=}" @@ -733,17 +800,18 @@ def create_baseline_cut_telescope( tel: Telescope, tm_path: Optional[DirPathType] = None, ) -> Tuple[DirPathType, Dict[str, str]]: - """Cut telescope `tel` for baseline-lengths. - - Args: - lcut: Lower cut - hcut: Higher cut - tel: Telescope to cut off - tm_path: .tm file-path to save the cut-telescope. - `tm_path` will get overwritten if it already exists. - - Returns: - .tm file-path & station-name conversion (e.g. station055 -> station009) + """Returns a telescope model for telescope `tel` with baseline lengths + only between `lcut` and `hcut` metres. + + Args: + lcut: Lower cut + hcut: Higher cut + tel: Telescope to cut off + tm_path: .tm file-path to save the cut-telescope. + `tm_path` will get overwritten if it already exists. + + Returns: + .tm file-path & station-name conversion (e.g. station055 -> station009) """ if tel.path is None: raise KaraboError( @@ -754,25 +822,17 @@ def create_baseline_cut_telescope( raise KaraboError(f"{tm_path=} must end with '.tm'.") df_tel = Telescope._get_station_infos(tel_path=tel.path) n_stations = df_tel.shape[0] + baseline_idx: List[Tuple[int, int]] = list(combinations(range(n_stations), 2)) + + n_baselines = len(baseline_idx) + station_x = df_tel["x"].to_numpy() station_y = df_tel["y"].to_numpy() - baselines: List[Tuple[int, int]] = sorted( - [ # each unique combination-idx a station with another station - tuple(station_idx) # type: ignore[misc] - for station_idx in set( - map( - frozenset, product(np.arange(n_stations), np.arange(n_stations)) - ) - ) - if len(station_idx) > 1 - ] - ) - n_baselines = len(baselines) baseline_dist = np.zeros(n_baselines) - for i, (x, y) in enumerate(baselines): + for i, (x, y) in enumerate(baseline_idx): baseline_dist[i] = np.linalg.norm(station_x[x] - station_y[y]) cut_idx = np.where((baseline_dist > lcut) & (baseline_dist < hcut)) - cut_station_list = np.unique(np.array(baselines)[cut_idx]) + cut_station_list = np.unique(np.array(baseline_idx)[cut_idx]) df_tel = df_tel[df_tel["station-nr"].isin(cut_station_list)].reset_index( drop=True ) @@ -806,16 +866,16 @@ def create_baseline_cut_telescope( dst=os.path.join(tm_path, "position.txt"), ) cut_stations = df_tel[["x", "y"]].to_numpy() - np.savetxt(os.path.join(tm_path, "layout.txt"), cut_stations) + np.savetxt(os.path.join(tm_path, "layout.txt"), cut_stations, delimiter=" ") return tm_path, conversions - def get_baselines_wgs84(self) -> NDArray[np.float64]: - """Gets the interferometer baselines in WGS84. + def get_stations_wgs84(self) -> NDArray[np.float64]: + """Gets the coordinates of the interferometer stations in WGS84. This function assumes that `self.stations` provides WGS84 coordinates. Returns: - Baselines lon[deg]/lat[deg]/alt[m] (nx3). + Stations lon[deg]/lat[deg]/alt[m] (nx3). """ return np.array( [ @@ -824,14 +884,14 @@ def get_baselines_wgs84(self) -> NDArray[np.float64]: ] ) - @classmethod # cls-fun to detach instance constraint - def get_baselines_dists( + @classmethod + def get_baseline_lengths( cls, - baselines_wgs84: NDArray[np.float64], + stations_wgs84: NDArray[np.float64], ) -> NDArray[np.float64]: """Gets the interferometer baselines distances in meters. - It's euclidean distance, not geodesic. + It's euclidean distance (aka geocentric), not geodesic. Args: baselines_wgs84: nx3 wgs84 baselines. @@ -840,15 +900,14 @@ def get_baselines_dists( Interferometer baselines dists in meters. """ lon, lat, alt = ( - baselines_wgs84[:, 0], - baselines_wgs84[:, 1], - baselines_wgs84[:, 2], - ) - cart_coords = wgs84_to_cartesian(lon, lat, alt) - dists: NDArray[np.float64] = np.linalg.norm( - cart_coords[:, np.newaxis, :] - cart_coords[np.newaxis, :, :], axis=2 + stations_wgs84[:, 0], + stations_wgs84[:, 1], + stations_wgs84[:, 2], ) - return dists + + cart_coords: NDArray[np.float64] = wgs84_to_cartesian(lon, lat, alt) + + return cast(NDArray[np.float64], pdist(cart_coords)) def max_baseline(self) -> np.float64: """Gets the longest baseline in meters. @@ -856,7 +915,7 @@ def max_baseline(self) -> np.float64: Returns: Length of longest baseline. """ - dists = self.get_baselines_dists(baselines_wgs84=self.get_baselines_wgs84()) + dists = self.get_baseline_lengths(stations_wgs84=self.get_stations_wgs84()) max_distance = np.max(dists) return max_distance diff --git a/karabo/test/test_coordinate_helper.py b/karabo/test/test_coordinate_helper.py new file mode 100644 index 00000000..6fff68ce --- /dev/null +++ b/karabo/test/test_coordinate_helper.py @@ -0,0 +1,46 @@ +import numpy as np +import pytest +from numpy.typing import NDArray + +from karabo.simulation.coordinate_helper import ( + east_north_to_long_lat, + wgs84_to_cartesian, +) + + +# wgs for LOFAR +# longitude 6.86763008, latitude 52.91139459, height 50.11317741 +# expected geocentric (a.k.a. EarthLocation in astropy) (metres) +# 3826923.9, 460915.1, 5064643.2 +def test_wgs84_to_cartesian(): + cart_coord: NDArray[np.float64] = wgs84_to_cartesian( + lon=6.86763008, lat=52.91139459, alt=50.11317741 + ) + + geocentric_coords_expected = [3826923.9, 460915.1, 5064643.2] + assert np.isclose(cart_coord, geocentric_coords_expected, rtol=0.01).all() + + +# +# east,north = 1000, 1000 --> lat,lon = 52.920378,6.876678 +# east,north = -1000, -1000 --> lat,lon = 52.902411,6.858582 + +testdata = [ + (1000, 1000, 52.920378, 6.876678), # go east and north 1000m + (-1000, -1000, 52.902411, 6.858582), # go west and south 1000m +] + + +@pytest.mark.parametrize("east, north, test_lat, test_lon", testdata) +def test_east_north_to_long_lat(east, north, test_lat, test_lon): + # coords of LOFAR in NL + lon = 6.86763008 + lat = 52.91139459 + + new_coords = np.array( + east_north_to_long_lat( + east_relative=east, north_relative=north, long=lon, lat=lat + ) + ) + + assert np.isclose(new_coords, np.array([test_lon, test_lat]), atol=1e-4).all() diff --git a/karabo/test/test_rascil_telescope_setup.py b/karabo/test/test_rascil_telescope_setup.py new file mode 100644 index 00000000..73cc2736 --- /dev/null +++ b/karabo/test/test_rascil_telescope_setup.py @@ -0,0 +1,52 @@ +import pytest + +from karabo.simulation.telescope import SimulatorBackend, Telescope + +rascil_telesecopes_to_test = [ + # (site_name, num_stations) + # data files are read from + # /envs/karabo/lib/python3.9/site-packages/ska_sdp_datamodels/configuration/ \ + # example_antenna_files/ + ("LOWBD2", 512), + ("MID", 197), + ("ASKAP", 36), + ("LOFAR", 134), +] + + +def test_set_telescope_name(): + site_name = "ASKAP" + + # create dummy telescope, name will be overriden later + site: Telescope = Telescope(0.0, 0.0, 0.0) + # __init__() sets name to None + assert site.name is None + + site.name = site_name + assert site.name == site_name + + +@pytest.mark.parametrize("site_name, _", rascil_telesecopes_to_test) +def test_set_telescope_from_oskar_telescope(site_name, _): + # we must set the backend to RASCIL. Otherwise OSKAR is used by default + site: Telescope = Telescope.constructor(site_name, backend=SimulatorBackend.RASCIL) + + site.name = site_name + assert site.name == site_name + assert site.backend == SimulatorBackend.RASCIL + + +@pytest.mark.parametrize("site_name, num_stations", rascil_telesecopes_to_test) +def test_num_of_stations(site_name, num_stations): + site: Telescope = Telescope.constructor(site_name, backend=SimulatorBackend.RASCIL) + assert len(site.stations) == num_stations + + +@pytest.mark.parametrize("site_name, num_stations", rascil_telesecopes_to_test) +def test_num_of_baselines(site_name, num_stations): + site: Telescope = Telescope.constructor(site_name, backend=SimulatorBackend.RASCIL) + + # This is the expected number of baselines + num_baselines = num_stations * (num_stations - 1) // 2 + stations = site.get_stations_wgs84() + assert len(site.get_baseline_lengths(stations)) == num_baselines diff --git a/karabo/test/test_telescope.py b/karabo/test/test_telescope.py index 077b1b2b..e73a382c 100644 --- a/karabo/test/test_telescope.py +++ b/karabo/test/test_telescope.py @@ -66,7 +66,11 @@ def test_OSKAR_telescope_plot_file_created(): tel = Telescope.constructor("MeerKAT") tel.plot_telescope(temp_plot_file_name) assert os.path.exists(temp_plot_file_name) - assert os.path.getsize(temp_plot_file_name) == 28171 + # It is tedious to check a specific file size. Even + # small changes to the code creating the image will make + # this test fail. Thus, I check only if the file size + # is > 0. + assert os.path.getsize(temp_plot_file_name) > 0 def test_create_alma_telescope(): @@ -202,7 +206,11 @@ def test_RASCIL_telescope_plot_file_created(): tel = Telescope.constructor("MID", backend=SimulatorBackend.RASCIL) tel.plot_telescope(temp_plot_file_name) assert os.path.exists(temp_plot_file_name) - assert os.path.getsize(temp_plot_file_name) == 20583 + # It is tedious to check a specific file size. Even + # small changes to the code creating the image will make + # this test fail. Thus, I check only if the file size + # is > 0. + assert os.path.getsize(temp_plot_file_name) > 0 # There is an if statement in Telescope::plot_telescope for the diff --git a/karabo/test/test_telescope_baselines.py b/karabo/test/test_telescope_baselines.py index 171866da..55f45f2e 100644 --- a/karabo/test/test_telescope_baselines.py +++ b/karabo/test/test_telescope_baselines.py @@ -1,8 +1,11 @@ +import math import os import tempfile from datetime import datetime, timedelta +from typing import Dict, Tuple import numpy as np +import pytest from numpy.typing import NDArray from karabo.imaging.imager_rascil import RascilDirtyImager, RascilDirtyImagerConfig @@ -10,43 +13,78 @@ from karabo.simulation.observation import Observation from karabo.simulation.sky_model import SkyModel from karabo.simulation.telescope import Telescope +from karabo.simulator_backend import SimulatorBackend +from karabo.util.data_util import get_module_absolute_path +from karabo.util.file_handler import DirPathType -def test_baselines_based_cutoff(sky_data: NDArray[np.float64]): +@pytest.fixture +def oskar_telescope() -> Telescope: + return Telescope.constructor("MeerKAT", backend=SimulatorBackend.OSKAR) + + +@pytest.fixture +def rascil_telescope() -> Telescope: + return Telescope.constructor("LOFAR", backend=SimulatorBackend.RASCIL) + + +@pytest.fixture +def sky_model() -> SkyModel: + sky = SkyModel() + sky_data = np.array( + [ + [20.0, -30.0, 1, 0, 0, 0, 100.0e6, -0.7, 0.0, 0, 0, 0], + [20.0, -30.5, 3, 2, 2, 0, 100.0e6, -0.7, 0.0, 600, 50, 45], + [20.5, -30.5, 3, 0, 0, 2, 100.0e6, -0.7, 0.0, 700, 10, -10], + ] + ) + sky.add_point_sources(sky_data) + return sky + + +# This test only tests that the function runs without error. This includes +# - counting number of stations in telescope after cut +# - assert that output measurement set was created +# - assert that FITS image calculated with cut baseline was written +# However, it doesn't test if the result is correct. You would expect a +# different image quality. This is because a cut baseline reduces the +# resolution of the instrument. +def test_baselines_based_cutoff(oskar_telescope: Telescope, sky_data: SkyModel): + # Max. baselength of MeerKAT is 7500 m. Thus, we cut somewhere + # inbetween, lcut = 5000 hcut = 10000 # Lower cut off and higher cut-off in meters - tel = Telescope.constructor("MeerKAT") with tempfile.TemporaryDirectory() as tmpdir: tm_path = os.path.join(tmpdir, "tel-cut.tm") - telescope_path, _ = Telescope.create_baseline_cut_telescope( + baseline_cut: Tuple[ + DirPathType, Dict[str, str] + ] = Telescope.create_baseline_cut_telescope( lcut, hcut, - tel, + oskar_telescope, tm_path=tm_path, ) + + telescope_path, _ = baseline_cut telescope = Telescope.read_OSKAR_tm_file(telescope_path) + # There are 64 stations fpr MeerKAT. After baseline cut there are + # 11 left. + assert len(telescope.get_stations_wgs84()) == 11 + sky = SkyModel() sky.add_point_sources(sky_data) simulation = InterferometerSimulation( channel_bandwidth_hz=1e6, time_average_sec=1, - noise_enable=False, - noise_seed="time", - noise_freq="Range", - noise_rms="Range", - noise_start_freq=1.0e9, - noise_inc_freq=1.0e8, - noise_number_freq=24, - noise_rms_start=5000, - noise_rms_end=10000, ) + observation = Observation( phase_centre_ra_deg=20.0, + phase_centre_dec_deg=-30.5, start_date_and_time=datetime(2022, 1, 1, 11, 00, 00, 521489), length=timedelta(hours=0, minutes=0, seconds=1, milliseconds=0), - phase_centre_dec_deg=-30.5, number_of_time_steps=1, - start_frequency_hz=1.0e9, + start_frequency_hz=100e6, frequency_increment_hz=1e6, number_of_channels=1, ) @@ -67,5 +105,47 @@ def test_baselines_based_cutoff(sky_data: NDArray[np.float64]): ) ) dirty = dirty_imager.create_dirty_image(visibility) + assert os.path.isdir(visibility.path) + dirty.write_to_file(os.path.join(tmpdir, "baseline_cut.fits"), overwrite=True) - dirty.plot(title="Flux Density (Jy)") + + assert os.path.isfile(os.path.join(tmpdir, "baseline_cut.fits")) + + dirty.plot( + title="Flux Density (Jy)", + filename=os.path.join( + get_module_absolute_path(), "test/data/image_cut.png" + ), + ) + + +def test_telescope_max_baseline_length( + oskar_telescope: Telescope, rascil_telescope: Telescope +): + max_length_oskar: np.float64 = oskar_telescope.max_baseline() + # Should be the same +/- 1 m + assert math.isclose(max_length_oskar - 7500.0, 0, abs_tol=1) + + max_length_rascil: np.float64 = rascil_telescope.max_baseline() + # Should be the same +/- 1 m + assert math.isclose(max_length_rascil - 995242.0, 0, abs_tol=1) + + freq_Hz = 100e6 + angular_res: float = Telescope.ang_res(freq_Hz, max_length_oskar) + assert math.isclose(angular_res, 1.44, rel_tol=1e-2) + + +def test_telescope_stations(oskar_telescope: Telescope, rascil_telescope: Telescope): + # station has 30 stations according to *.tm file + baseline_wgs: NDArray[np.float64] = oskar_telescope.get_stations_wgs84() + assert len(baseline_wgs) == 64 + + baseline_wgs: NDArray[np.float64] = rascil_telescope.get_stations_wgs84() + assert len(baseline_wgs) == 134 + + +def test_telescope_baseline_length(rascil_telescope): + stations_wgs: NDArray[np.float64] = rascil_telescope.get_stations_wgs84() + num_stations = len(stations_wgs) + baseline_length: NDArray[np.float64] = Telescope.get_baseline_lengths(stations_wgs) + assert len(baseline_length) == num_stations * (num_stations - 1) / 2