diff --git a/lofarimaging.ipynb b/lofarimaging.ipynb index 8b145a7..fd8f3e0 100644 --- a/lofarimaging.ipynb +++ b/lofarimaging.ipynb @@ -61,7 +61,8 @@ "import warnings\n", "from mpl_toolkits.axes_grid1 import make_axes_locatable\n", "import matplotlib.axes as maxes\n", - "from matplotlib.ticker import FormatStrFormatter" + "from matplotlib.ticker import FormatStrFormatter\n", + "import re" ] }, { @@ -118,8 +119,7 @@ "outputs": [], "source": [ "data_dir = \"./test\"\n", - "caltable_dir = \"./test/CalTables\" # Root directory under which station information is stored in subdirectories DE602C/etc/, RS106/etc/, ...\n", - "station_name = 'DE603'" + "caltable_dir = \"./test/CalTables\" # Root directory under which station information is stored in subdirectories DE602C/etc/, RS106/etc/, ..." ] }, { @@ -139,7 +139,7 @@ "metadata": {}, "outputs": [], "source": [ - "station_type = get_station_type(station_name)" + "os.makedirs('results', exist_ok=True)" ] }, { @@ -147,22 +147,13 @@ "execution_count": 11, "metadata": {}, "outputs": [], - "source": [ - "os.makedirs('results', exist_ok=True)" - ] - }, - { - "cell_type": "code", - "execution_count": 12, - "metadata": {}, - "outputs": [], "source": [ "start1 = time.time()" ] }, { "cell_type": "code", - "execution_count": 13, + "execution_count": 12, "metadata": {}, "outputs": [ { @@ -170,7 +161,7 @@ "output_type": "stream", "text": [ "Searching for available files in ./test\n", - "0: ./test/20200525_084652_mode_5_xst_sb100.dat\n", + "0: ./test/20170621_072634_mode_sparseeven_xst_sb350.dat\n", "1: ./test/20170720_095816_mode_3_xst_sb297.dat\n" ] } @@ -185,27 +176,58 @@ }, { "cell_type": "code", - "execution_count": 58, + "execution_count": 33, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ - "File selected: ./test/20200525_084652_mode_5_xst_sb100.dat\n" + "File selected: ./test/20170720_095816_mode_3_xst_sb297.dat\n" ] } ], "source": [ "# Select a file\n", - "xst_filename = files[0]\n", + "xst_filename = \"./test/20170720_095816_mode_3_xst_sb297.dat\" # files[0]\n", "\n", "print(\"File selected:\", xst_filename)" ] }, { "cell_type": "code", - "execution_count": 15, + "execution_count": 34, + "metadata": {}, + "outputs": [], + "source": [ + "station_name = \"DE603\"\n", + "\n", + "rcu_mode = re.search(\"mode_([^_]*)\", xst_filename).groups(0)[0]\n", + "subband = int(re.search(\"sb([0-9]*)\", xst_filename).groups(0)[0])\n", + "obsdatestr, obstimestr, *_ = os.path.basename(xst_filename).rstrip(\".dat\").split(\"_\")" + ] + }, + { + "cell_type": "code", + "execution_count": 35, + "metadata": {}, + "outputs": [], + "source": [ + "rcu_mode = rcu_mode.replace(\"sparseeven\", \"sparse_even\").replace(\"sparseodd\", \"sparse_odd\")" + ] + }, + { + "cell_type": "code", + "execution_count": 36, + "metadata": {}, + "outputs": [], + "source": [ + "station_type = get_station_type(station_name)" + ] + }, + { + "cell_type": "code", + "execution_count": 37, "metadata": {}, "outputs": [ { @@ -225,10 +247,6 @@ } ], "source": [ - "# Distill metadata from filename\n", - "cubename = os.path.basename(xst_filename)\n", - "obsdatestr, obstimestr, _, rcu_mode, _, subbandname = cubename.rstrip(\".dat\").split(\"_\")\n", - "subband = int(subbandname[2:])\n", "station_name = get_full_station_name(station_name, rcu_mode)\n", "\n", "# Get the data\n", @@ -246,8 +264,7 @@ "obstime = datetime.datetime.strptime(obsdatestr + \":\" + obstimestr, '%Y%m%d:%H%M%S')\n", "\n", "# Confirm the data has been read correctly\n", - "print(f\"\"\"Filename: {cubename}\n", - "Station: {station_name}\n", + "print(f\"\"\"Station: {station_name}\n", "Station type: {station_type}\n", "Subband: {subband}\n", "Timestep: {timestep}\n", @@ -260,27 +277,7 @@ }, { "cell_type": "code", - "execution_count": 16, - "metadata": {}, - "outputs": [ - { - "data": { - "text/plain": [ - "'20170720_095816_DE603LBA_SB297'" - ] - }, - "execution_count": 16, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "fname" - ] - }, - { - "cell_type": "code", - "execution_count": 17, + "execution_count": 40, "metadata": {}, "outputs": [], "source": [ @@ -290,7 +287,7 @@ }, { "cell_type": "code", - "execution_count": 18, + "execution_count": 41, "metadata": {}, "outputs": [], "source": [ @@ -299,7 +296,7 @@ }, { "cell_type": "code", - "execution_count": 19, + "execution_count": 42, "metadata": {}, "outputs": [], "source": [ @@ -315,7 +312,7 @@ }, { "cell_type": "code", - "execution_count": 20, + "execution_count": 43, "metadata": {}, "outputs": [ { @@ -324,7 +321,7 @@ "((29, 96, 96), (29, 96, 96), (96, 96), dtype('complex128'))" ] }, - "execution_count": 20, + "execution_count": 43, "metadata": {}, "output_type": "execute_result" } @@ -408,7 +405,7 @@ }, { "cell_type": "code", - "execution_count": 23, + "execution_count": 46, "metadata": {}, "outputs": [], "source": [ @@ -418,7 +415,7 @@ }, { "cell_type": "code", - "execution_count": 24, + "execution_count": 47, "metadata": {}, "outputs": [], "source": [ @@ -452,7 +449,7 @@ }, { "cell_type": "code", - "execution_count": 27, + "execution_count": 49, "metadata": {}, "outputs": [], "source": [ @@ -487,7 +484,7 @@ }, { "cell_type": "code", - "execution_count": 29, + "execution_count": 51, "metadata": {}, "outputs": [ { diff --git a/lofarimaging/hdf5util.py b/lofarimaging/hdf5util.py index fa29cbd..d168637 100644 --- a/lofarimaging/hdf5util.py +++ b/lofarimaging/hdf5util.py @@ -101,6 +101,9 @@ def write_hdf5(filename: str, xst_data: np.ndarray, visibilities: np.ndarray, sk """ short_station_name = station_name[:5] + if subtracted is None: + subtracted = [] + with h5py.File(filename, 'a') as h5file: new_obsname = get_new_obsname(h5file) obs_group = h5file.create_group(new_obsname) diff --git a/lofarimaging/singlestationutil.py b/lofarimaging/singlestationutil.py index 5685fc8..4c8700d 100644 --- a/lofarimaging/singlestationutil.py +++ b/lofarimaging/singlestationutil.py @@ -100,16 +100,18 @@ def freq_from_sb(sb: int, rcu_mode: Union[str, int] = 1): 58007812.5 """ clock = 200e6 - if int(rcu_mode) == 6: - clock = 160e6 - freq_offset = 0 - if int(rcu_mode) == 5: - freq_offset = 100e6 - elif int(rcu_mode) == 6: - freq_offset = 160e6 - elif int(rcu_mode) == 7: - freq_offset = 200e6 + + if 'sparse' not in str(rcu_mode): + if int(rcu_mode) == 6: + clock = 160e6 + + if int(rcu_mode) == 5: + freq_offset = 100e6 + elif int(rcu_mode) == 6: + freq_offset = 160e6 + elif int(rcu_mode) == 7: + freq_offset = 200e6 sb_bandwidth = 0.5 * clock / 512. freq = (sb * sb_bandwidth) + freq_offset @@ -150,6 +152,10 @@ def find_caltable(field_name: str, rcu_mode: Union[str, int], caltable_dir='calt filename += "-HBA-170_230.dat" elif str(rcu_mode) == '7': filename += "-HBA-210_250.dat" + elif str(rcu_mode) == 'sparse_even': + filename += "-LBA_SPARSE_EVEN-10_90.dat" + elif str(rcu_mode) == 'sparse_odd': + filename += "-LBA_SPARSE_ODD-10_90.dat" else: raise RuntimeError("Unexpected mode: " + str(rcu_mode) + " for field_name " + str(field_name)) @@ -317,15 +323,24 @@ def get_station_pqr(station_name: str, rcu_mode: Union[str, int], db): full_station_name = get_full_station_name(station_name, rcu_mode) station_type = get_station_type(full_station_name) - if 'LBA' in station_name or str(rcu_mode) in ('1', '2', '3', '4', 'inner', 'outer'): - # Get the PQR positions for an individual station - station_pqr = db.antenna_pqr(full_station_name) - - # Exception: for Dutch stations (sparse not yet accommodated) - if (station_type == 'core' or station_type == 'remote') and int(rcu_mode) in (3, 4): - station_pqr = station_pqr[0:48, :] - elif (station_type == 'core' or station_type == 'remote') and int(rcu_mode) in (1, 2): - station_pqr = station_pqr[48:, :] + if 'LBA' in station_name or str(rcu_mode) in ('1', '2', '3', '4', 'inner', 'outer', 'sparse_even', 'sparse_odd', 'sparse'): + if (station_type == 'core' or station_type == 'remote'): + if str(rcu_mode) in ('3', '4', 'inner'): + station_pqr = db.antenna_pqr(full_station_name)[0:48, :] + elif str(rcu_mode) in ('1', '2', 'outer'): + station_pqr = db.antenna_pqr(full_station_name)[48:, :] + elif rcu_mode in ('sparse_even', 'sparse'): + all_pqr = db.antenna_pqr(full_station_name) + # Indices 0, 49, 2, 51, 4, 53, ... + station_pqr = np.ravel(np.column_stack((all_pqr[:48:2], all_pqr[49::2]))).reshape(48, 3) + elif rcu_mode == 'sparse_odd': + all_pqr = db.antenna_pqr(full_station_name) + # Indices 1, 48, 3, 50, 5, 52, ... + station_pqr = np.ravel(np.column_stack((all_pqr[1:48:2], all_pqr[48::2]))).reshape(48, 3) + else: + raise RuntimeError("Cannot select subset of LBA antennas for mode " + rcu_mode) + else: + station_pqr = db.antenna_pqr(full_station_name) elif 'HBA' in station_name or str(rcu_mode) in ('5', '6', '7', '8'): selected_dipole_config = { 'intl': GENERIC_INT_201512, 'remote': GENERIC_REMOTE_201512, 'core': GENERIC_CORE_201512 @@ -549,14 +564,13 @@ def get_full_station_name(station_name: str, rcu_mode: Union[str, int]) -> str: return station_name if str(rcu_mode) in ('1', '2', 'outer'): - if len(station_name) == 5: - station_name += "LBA" + station_name += "LBA" elif str(rcu_mode) in ('3', '4', 'inner'): - if len(station_name) == 5: - station_name += "LBA" + station_name += "LBA" + elif 'sparse' in str(rcu_mode): + station_name += "LBA" elif str(rcu_mode) in ('5', '6', '7'): - if len(station_name) == 5: - station_name += "HBA" + station_name += "HBA" else: raise Exception("Unexpected rcu_mode: ", rcu_mode) @@ -645,6 +659,14 @@ def make_xst_plots(xst_data: np.ndarray, >>> type(leafletmap) + + >>> xst_data = read_acm_cube("test/20170621_072634_sb350_xst.dat", "remote")[0] + >>> obstime = datetime.datetime(2017, 6, 21, 7, 26, 34) + >>> sky_fig, ground_fig, leafletmap = make_xst_plots(xst_data, "RS509", obstime, 350, \ + 'sparse_even', \ + caltable_dir="test/CalTables", \ + hdf5_filename="test/test.h5") + Maximum at 2m east, -2m north of station center (lat/long 53.40884, 6.78531) """ if extent is None: extent = [-150, 150, -150, 150] @@ -853,7 +875,7 @@ def reimage_sky(h5: h5py.File, obsnum: str, db: lofarantpos.db.LofarAntennaDatab sky_data = h5[obsnum]["sky_img"] freq = h5[obsnum].attrs['frequency'] marked_bodies_lmn = dict(zip(h5[obsnum].attrs["source_names"], h5[obsnum].attrs["source_lmn"])) - visibilities = h5[obsnum]['calibrated_data'].value + visibilities = h5[obsnum]['calibrated_data'][:] visibilities_xx = visibilities[0::2, 0::2] visibilities_yy = visibilities[1::2, 1::2] # Stokes I @@ -902,7 +924,7 @@ def reimage_nearfield(h5: h5py.File, obsnum: str, db: lofarantpos.db.LofarAntenn rcu_mode = h5[obsnum].attrs['rcu_mode'] freq = h5[obsnum].attrs['frequency'] marked_bodies_lmn = dict(zip(h5[obsnum].attrs["source_names"], h5[obsnum].attrs["source_lmn"])) - visibilities = h5[obsnum]['calibrated_data'].value + visibilities = h5[obsnum]['calibrated_data'][:] visibilities_xx = visibilities[0::2, 0::2] visibilities_yy = visibilities[1::2, 1::2] # Stokes I diff --git a/test/20170621_072634_sb350_xst.dat b/test/20170621_072634_sb350_xst.dat new file mode 100644 index 0000000..52a36b7 Binary files /dev/null and b/test/20170621_072634_sb350_xst.dat differ diff --git a/test/CalTables/RS509/CalTable-509-LBA_SPARSE_EVEN-10_90.dat b/test/CalTables/RS509/CalTable-509-LBA_SPARSE_EVEN-10_90.dat new file mode 100644 index 0000000..f7839cf Binary files /dev/null and b/test/CalTables/RS509/CalTable-509-LBA_SPARSE_EVEN-10_90.dat differ