From 3fb4a9f56f5585c1e92ecd48f690b9aa4f928c35 Mon Sep 17 00:00:00 2001 From: Tammo Jan Dijkema Date: Fri, 28 Feb 2020 13:54:45 +0100 Subject: [PATCH] Move out single station util --- LOFAR LBA imaging tutorial v4.ipynb | 12 +- lofarimaging/__init__.py | 1 + lofarimaging/lofarimaging.py | 561 +--------------------------- lofarimaging/singlestationutil.py | 526 ++++++++++++++++++++++++++ 4 files changed, 551 insertions(+), 549 deletions(-) create mode 100644 lofarimaging/singlestationutil.py diff --git a/LOFAR LBA imaging tutorial v4.ipynb b/LOFAR LBA imaging tutorial v4.ipynb index 7dde0a6..73ff050 100644 --- a/LOFAR LBA imaging tutorial v4.ipynb +++ b/LOFAR LBA imaging tutorial v4.ipynb @@ -85,7 +85,7 @@ }, { "cell_type": "code", - "execution_count": 52, + "execution_count": 6, "metadata": {}, "outputs": [], "source": [ @@ -829,7 +829,7 @@ "name": "stdout", "output_type": "stream", "text": [ - "Time elapsed: 50.90 s\n" + "Time elapsed: 36.33 s\n" ] } ], @@ -890,19 +890,19 @@ }, { "cell_type": "code", - "execution_count": 56, + "execution_count": 50, "metadata": {}, "outputs": [ { "data": { "text/html": [ - "
" + "
" ], "text/plain": [ - "" + "" ] }, - "execution_count": 56, + "execution_count": 50, "metadata": {}, "output_type": "execute_result" } diff --git a/lofarimaging/__init__.py b/lofarimaging/__init__.py index 437c486..e4f0ca5 100644 --- a/lofarimaging/__init__.py +++ b/lofarimaging/__init__.py @@ -1,2 +1,3 @@ from .lofarimaging import * from .maputil import * +from .singlestationutil import * diff --git a/lofarimaging/lofarimaging.py b/lofarimaging/lofarimaging.py index 82b4645..33f8566 100644 --- a/lofarimaging/lofarimaging.py +++ b/lofarimaging/lofarimaging.py @@ -2,102 +2,13 @@ import numpy as np import numexpr as ne -import os -from matplotlib.pyplot import imread -import folium -import datetime -import lofargeotiff +from astropy.coordinates import SkyCoord, SkyOffsetFrame, CartesianRepresentation -from scipy import ndimage -from lofarantpos.db import LofarAntennaDatabase - -import matplotlib.pyplot as plt -from matplotlib.ticker import FormatStrFormatter -from matplotlib import cm -from matplotlib.colors import ListedColormap, Normalize -import warnings -from mpl_toolkits.axes_grid1 import make_axes_locatable -import matplotlib.axes as maxes - -from astropy.coordinates import SkyCoord, GCRS, EarthLocation, AltAz, SkyOffsetFrame, CartesianRepresentation, get_sun -import astropy.units as u -from astropy.time import Time - -import lofarantpos -from packaging import version - -from .maputil import get_map, make_leaflet_map - -__all__ = ["sb_from_freq", "freq_from_sb", "find_caltable", "read_caltable", - "rcus_in_station", "read_acm_cube", "nearfield_imager", - "sky_imager", "ground_imager", "get_station_pqr", "skycoord_to_lmn", "make_ground_image"] +__all__ = ["nearfield_imager", "sky_imager", "ground_imager", "skycoord_to_lmn"] __version__ = "1.5.0" - -# Configurations for HBA observations with a single dipole activated per tile. -GENERIC_INT_201512 = [0, 5, 3, 1, 8, 3, 12, 15, 10, 13, 11, 5, 12, 12, 5, 2, 10, 8, 0, 3, 5, 1, 4, 0, 11, 6, 2, 4, 9, - 14, 15, 3, 7, 5, 13, 15, 5, 6, 5, 12, 15, 7, 1, 1, 14, 9, 4, 9, 3, 9, 3, 13, 7, 14, 7, 14, 2, 8, - 8, 0, 1, 4, 2, 2, 12, 15, 5, 7, 6, 10, 12, 3, 3, 12, 7, 4, 6, 0, 5, 9, 1, 10, 10, 11, 5, 11, 7, 9, - 7, 6, 4, 4, 15, 4, 1, 15] -GENERIC_CORE_201512 = [0, 10, 4, 3, 14, 0, 5, 5, 3, 13, 10, 3, 12, 2, 7, 15, 6, 14, 7, 5, 7, 9, 0, 15, 0, 10, 4, 3, 14, - 0, 5, 5, 3, 13, 10, 3, 12, 2, 7, 15, 6, 14, 7, 5, 7, 9, 0, 15] -GENERIC_REMOTE_201512 = [0, 13, 12, 4, 11, 11, 7, 8, 2, 7, 11, 2, 10, 2, 6, 3, 8, 3, 1, 7, 1, 15, 13, 1, 11, 1, 12, 7, - 10, 15, 8, 2, 12, 13, 9, 13, 4, 5, 5, 12, 5, 5, 9, 11, 15, 12, 2, 15] - -assert (version.parse(lofarantpos.__version__) >= version.parse("0.4.0")) - - -def sb_from_freq(freq: float, rcu_mode='1'): - """ - Convert subband number to central frequency - - Args: - rcu_mode: rcu mode - freq: frequency in Hz - - Returns: - int: subband number - """ - clock = 200e6 - if rcu_mode == '6': - clock = 160e6 - sb_bandwidth = 0.5 * clock / 512. - freq_offset = 0 - if rcu_mode == '5': - freq_offset = 100e6 - elif rcu_mode == '6': - freq_offset = 160e6 - elif rcu_mode == '7': - freq_offset = 200e6 - sb = round((freq - freq_offset) / sb_bandwidth) - return int(sb) - - -def freq_from_sb(sb: int, rcu_mode='1'): - """ - Convert central frequency to subband number - - Args: - rcu_mode: rcu mode - sb: subband number - - Returns: - float: frequency in Hz - """ - clock = 200e6 - if rcu_mode == '6': - clock = 160e6 - freq_offset = 0 - if rcu_mode == '5': - freq_offset = 100e6 - elif rcu_mode == '6': - freq_offset = 160e6 - elif rcu_mode == '7': - freq_offset = 200e6 - sb_bandwidth = 0.5 * clock / 512. - freq = (sb * sb_bandwidth) + freq_offset - return freq +SPEED_OF_LIGHT = 299792458.0 def skycoord_to_lmn(pos: SkyCoord, phasecentre: SkyCoord): @@ -123,174 +34,29 @@ def skycoord_to_lmn(pos: SkyCoord, phasecentre: SkyCoord): return dc.y.value, dc.z.value, dc.x.value - 1 -def find_caltable(field_name: str, rcu_mode: str, config_dir='caltables'): - """ - Find the file of a caltable. - - Args: - field_name: Name of the antenna field, e.g. 'DE602LBA' - rcu_mode: Receiver mode for which the calibration table is requested. - Probably should be 'inner' or 'outer' - config_dir: Root directory under which station information is stored in - subdirectories DE602C/etc/, RS106/etc/, ... - - Returns: - str: filename if it exists, None if nothing found - """ - station, field = field_name[0:5].upper(), field_name[5:].upper() - station_number = station[2:5] - - # Map to the correct file depending on the RCU mode - if rcu_mode == 'outer' and 'LBA' in field_name: - filename = os.path.join(config_dir, f"CalTable-{station_number}-LBA_OUTER-10_90.dat") - elif rcu_mode == 'inner' and 'LBA' in field_name: - filename = os.path.join(config_dir, f"CalTable-{station_number}-LBA_INNER-10_90.dat") - else: - filename = os.path.join(config_dir, f"CalTable_{station_number}_mode{rcu_mode}.dat") - - if os.path.exists(filename): - return filename - - # If the original folder structure is kept - if rcu_mode == 'outer' and 'LBA' in field_name: - filename = os.path.join(config_dir, f"{station}/CalTable-{station_number}-LBA_OUTER-10_90.dat") - elif rcu_mode == 'inner' and 'LBA' in field_name: - filename = os.path.join(config_dir, f"{station}/CalTable-{station_number}-LBA_INNER-10_90.dat") - else: - filename = os.path.join(config_dir, f"{station}/CalTable_{station_number}_mode{rcu_mode}.dat") - - if os.path.exists(filename): - return filename - else: - return None - - -def read_caltable(filename: str, num_subbands=512): - """ - Read a station's calibration table. - - Args: - filename: Filename with the caltable - num_subbands: Number of subbands - - Returns: - Tuple(Dict[str, str], np.array): A tuple containing a dict with - the header lines, and a 2D numpy.array of complex numbers - representing the station gain coefficients. - """ - infile = open(filename, 'rb') - - header_lines = [] - - try: - while True: - header_lines.append(infile.readline().decode('utf8').strip()) - if 'HeaderStop' in header_lines[-1]: - break - except UnicodeDecodeError: - # No header; close and open again - infile.close() - infile = open(filename, 'rb') - - caldata = np.fromfile(infile, dtype=np.complex128) - num_rcus = len(caldata) // num_subbands - - infile.close() - - header_dict = {key: val for key, val in [line.split(" = ") - for line in header_lines[1:-1]]} - - return header_dict, caldata.reshape((num_subbands, num_rcus)) - - -def rcus_in_station(station_type: str): - """ - Give the number of RCUs in a station, given its type. - - Args: - station_type: Kind of station that produced the correlation. One of - 'core', 'remote', 'intl'. - """ - return {'core': 96, 'remote': 96, 'intl': 192}[station_type] - - -def read_acm_cube(filename: str, station_type: str): - """ - Read an ACM binary data cube (function from Michiel) - - Args: - filename: File containing the array correlation matrix. - station_type: Kind of station that produced the correlation. One of - 'core', 'remote', 'intl'. - - Returns: - np.array: 3D cube of complex numbers, with indices [time slots, rcu, rcu]. - - Examples: - >>> cube = read_acm_cube('20170720_095816_xst.dat', 'intl') - >>> cube.shape - (29, 192, 192) - """ - num_rcu = rcus_in_station(station_type) - data = np.fromfile(filename, dtype=np.complex128) - time_slots = int(len(data) / num_rcu / num_rcu) - return data.reshape((time_slots, num_rcu, num_rcu)) - - -SPEED_OF_LIGHT = 299792458.0 - - -def get_station_pqr(station_name: str, station_type: str, array_type: str, db): - if 'LBA' in station_name: - # Get the PQR positions for an individual station - station_pqr = db.antenna_pqr(station_name) - - # Exception: for Dutch stations (sparse not yet accommodated) - if (station_type == 'core' or station_type == 'remote') and array_type == 'inner': - station_pqr = station_pqr[0:48, :] - elif (station_type == 'core' or station_type == 'remote') and array_type == 'outer': - station_pqr = station_pqr[48:, :] - elif 'HBA' in station_name: - selected_dipole_config = { - 'intl': GENERIC_INT_201512, 'remote': GENERIC_REMOTE_201512, 'core': GENERIC_CORE_201512 - } - selected_dipoles = selected_dipole_config[station_type] + \ - np.arange(len(selected_dipole_config[station_type])) * 16 - station_pqr = db.hba_dipole_pqr(station_name)[selected_dipoles] - else: - raise RuntimeError("Station name did not contain LBA or HBA, could not load antenna positions") - - return station_pqr.astype('float32') - - def sky_imager(visibilities, baselines, freq, npix_l, npix_m): """Do a Fourier transform for sky imaging""" - img = np.zeros([npix_m, npix_l], dtype=np.float32) + img = np.zeros([npix_m, npix_l], dtype=np.complex128) - with warnings.catch_warnings(): - warnings.filterwarnings("ignore", message="Casting complex values to real discards the imaginary part") - for m_ix, m in enumerate(np.linspace(-1, 1, npix_l)): - for l_ix, l in enumerate(np.linspace(1, -1, npix_m)): - img[m_ix, l_ix] = np.mean(visibilities * - np.exp(-2j * np.pi * freq * - (baselines[:, :, 0] * l + baselines[:, :, 1] * m) / SPEED_OF_LIGHT)) - return img + for m_ix, m in enumerate(np.linspace(-1, 1, npix_l)): + for l_ix, l in enumerate(np.linspace(1, -1, npix_m)): + img[m_ix, l_ix] = np.mean(visibilities * + np.exp(-2j * np.pi * freq * + (baselines[:, :, 0] * l + baselines[:, :, 1] * m) / SPEED_OF_LIGHT)) + return np.abs(img) def ground_imager(visibilities, freq, npix_p, npix_q, dims, station_pqr, height=1.5): """Do a Fourier transform for ground imaging""" - img = np.zeros([npix_q, npix_p], dtype=np.complex) + img = np.zeros([npix_q, npix_p], dtype=np.complex128) - with warnings.catch_warnings(): - warnings.filterwarnings("ignore", message="Casting complex values to real discards the imaginary part") - for q_ix, q in enumerate(np.linspace(dims[2], dims[3], npix_q)): - for p_ix, p in enumerate(np.linspace(dims[0], dims[1], npix_p)): - r = height - pqr = np.array([p, q, r], dtype=np.float32) - antdist = np.linalg.norm(station_pqr - pqr[np.newaxis, :], axis=1) - groundbase = antdist[:, np.newaxis] - antdist[np.newaxis, :] - img[q_ix, p_ix] = np.mean(visibilities * np.exp(-2j * np.pi * freq * (-groundbase) / SPEED_OF_LIGHT)) - return np.abs(img) + for q_ix, q in enumerate(np.linspace(dims[2], dims[3], npix_q)): + for p_ix, p in enumerate(np.linspace(dims[0], dims[1], npix_p)): + r = height + pqr = np.array([p, q, r], dtype=np.float32) + antdist = np.linalg.norm(station_pqr - pqr[np.newaxis, :], axis=1) + groundbase = antdist[:, np.newaxis] - antdist[np.newaxis, :] + img[q_ix, p_ix] = np.mean(visibilities * np.exp(-2j * np.pi * freq * (-groundbase) / SPEED_OF_LIGHT)) def nearfield_imager(visibilities, baseline_indices, freqs, npix_p, npix_q, extent, station_pqr, height=1.5, @@ -343,294 +109,3 @@ def nearfield_imager(visibilities, baseline_indices, freqs, npix_p, npix_q, exte img /= len(freqs) * len(baseline_indices) return img - - -def make_ground_image(xst_filename, - station_name, - caltable_dir, - extent=None, - pixels_per_metre=0.5, - sky_vmin=None, - sky_vmax=None, - ground_vmin=None, - ground_vmax=None, - height=1.5, - map_zoom=19, - sky_only=False, - opacity=0.6): - """Make a ground image""" - cubename = os.path.basename(xst_filename) - - if extent is None: - extent = [-150, 150, -150, 150] - - if station_name[0] == "C": - station_type = "core" - elif station_name[0] == "R" or station_name[:5] == "PL611": - station_type = "remote" - else: - station_type = "intl" - - try: - os.mkdir('results') - except FileExistsError: - pass - - # Distill metadata from filename - obsdatestr, obstimestr, _, rcu_mode, _, subbandname = cubename.rstrip(".dat").split("_") - subband = int(subbandname[2:]) - - # Needed for NL stations: inner (rcu_mode 3/4), outer (rcu_mode 1/2), (sparse tbd) - # Should be set to 'inner' if station type = 'intl' - array_type = None - if rcu_mode in ('1', '2'): - array_type = 'outer' - elif rcu_mode in ('3', '4'): - array_type = 'inner' - elif rcu_mode in ('5', '6', '7'): - array_type = rcu_mode - else: - raise Exception("Unexpected rcu_mode: ", rcu_mode) - - # Get the data - fname = f"{obsdatestr}_{obstimestr}_{station_name}_SB{subband}" - - npix_l, npix_m = 131, 131 - freq = freq_from_sb(subband, rcu_mode=rcu_mode) - - # Which slice in time to visualise - timestep = 0 - - # For ground imaging - ground_resolution = pixels_per_metre # pixels per metre for ground_imaging, default is 0.5 pixel/metre - - obstime = datetime.datetime.strptime(obsdatestr + ":" + obstimestr, '%Y%m%d:%H%M%S') - - cube = read_acm_cube(xst_filename, station_type) - - # Apply calibration - - caltable_filename = find_caltable(station_name, rcu_mode=array_type, - config_dir=caltable_dir) - - cal_header = None - if caltable_filename is None: - print('No calibration table found... cube remains uncalibrated!') - else: - cal_header, cal_data = read_caltable(caltable_filename) - - rcu_gains = cal_data[subband, :] - rcu_gains = np.array(rcu_gains, dtype=np.complex64) - gain_matrix = rcu_gains[np.newaxis, :] * np.conj(rcu_gains[:, np.newaxis]) - cube = cube / gain_matrix - - # Split into the XX and YY polarisations (RCUs) - # This needs to be modified in future for LBA sparse - cube_xx = cube[:, 0::2, 0::2] - cube_yy = cube[:, 1::2, 1::2] - visibilities_all = cube_xx + cube_yy - - # Stokes I for specified timestep - visibilities = visibilities_all[timestep] - - # Setup the database - db = LofarAntennaDatabase() - - station_pqr = get_station_pqr(station_name, station_type, array_type, db) - - # Rotate station_pqr to a north-oriented xyz frame, where y points North, in a plane through the station. - rotation = db.rotation_from_north(station_name) - - pqr_to_xyz = np.array([[np.cos(-rotation), -np.sin(-rotation), 0], - [np.sin(-rotation), np.cos(-rotation), 0], - [0, 0, 1]]) - - station_xyz = (pqr_to_xyz @ station_pqr.T).T - - baselines = station_xyz[:, np.newaxis, :] - station_xyz[np.newaxis, :, :] - - # Make a sky image, by numerically Fourier-transforming from visibilities to image plane - from matplotlib.patches import Circle - - # Fourier transform, and account for the rotation (rotation is positive in this space) - # visibilities = cube_xx[2,:,:] - img = sky_imager(visibilities, baselines, freq, npix_l, npix_m) - img = ndimage.interpolation.rotate(img, -rotation, reshape=False, mode='nearest') - - obstime_astropy = Time(obstime) - # Determine positions of Cas A and Cyg A - station_earthlocation = EarthLocation.from_geocentric(*(db.phase_centres[station_name] * u.m)) - zenith = AltAz(az=0 * u.deg, alt=90 * u.deg, obstime=obstime_astropy, - location=station_earthlocation).transform_to(GCRS) - - marked_bodies = { - 'Cas A': SkyCoord(ra=350.85 * u.deg, dec=58.815 * u.deg), - 'Cyg A': SkyCoord(ra=299.868 * u.deg, dec=40.734 * u.deg), - # 'Per A': SkyCoord.from_name("Perseus A"), - # 'Her A': SkyCoord.from_name("Hercules A"), - # 'Cen A': SkyCoord.from_name("Centaurus A"), - # '?': SkyCoord.from_name("J101415.9+105106"), - # '3C295': SkyCoord.from_name("3C295"), - # 'Moon': get_moon(obstime_astropy, location=station_earthlocation).transform_to(GCRS), - 'Sun': get_sun(obstime_astropy) - # '3C196': SkyCoord.from_name("3C196") - } - - marked_bodies_lmn = {} - for body_name, body_coord in marked_bodies.items(): - # print(body_name, body_coord.separation(zenith), body_coord.separation(zenith)) - if body_coord.transform_to(AltAz(location=station_earthlocation, obstime=obstime_astropy)).alt > 0: - marked_bodies_lmn[body_name] = skycoord_to_lmn(marked_bodies[body_name], zenith) - - # Plot the resulting sky image - fig, ax = plt.subplots(figsize=(10, 10)) - - circle1 = Circle((0, 0), 1.0, edgecolor='k', fill=False, facecolor='none', alpha=0.3) - ax.add_artist(circle1) - - cimg = ax.imshow(img, origin='lower', cmap=cm.Spectral_r, extent=(1, -1, -1, 1), - clip_path=circle1, clip_on=True, vmin=sky_vmin, vmax=sky_vmax) - divider = make_axes_locatable(ax) - cax = divider.append_axes("right", size="5%", pad=0.2, axes_class=maxes.Axes) - fig.colorbar(cimg, cax=cax, orientation="vertical", format="%.1e") - - ax.set_xlim(1, -1) - - ax.set_xticks(np.arange(-1, 1.1, 0.5)) - ax.xaxis.set_major_formatter(FormatStrFormatter('%.1f')) - ax.set_yticks(np.arange(-1, 1.1, 0.5)) - ax.yaxis.set_major_formatter(FormatStrFormatter('%.1f')) - - # Labels - ax.set_xlabel('$ℓ$', fontsize=14) - ax.set_ylabel('$m$', fontsize=14) - - for body_name, lmn in marked_bodies_lmn.items(): - ax.plot([lmn[0]], [lmn[1]], marker='x', color='white', mew=0.5) - ax.annotate(body_name, (lmn[0], lmn[1])) - - ax.text(0.5, 1.05, f"Sky image for {station_name}", - fontsize=17, ha='center', va='bottom', transform=ax.transAxes) - ax.text(0.5, 1.02, f"SB {subband} ({freq / 1e6:.1f} MHz), {str(obstime)[:16]}", - fontsize=12, ha='center', va='bottom', transform=ax.transAxes) - - # Plot the compass directions - ax.text(0.9, 0, 'W', horizontalalignment='center', verticalalignment='center', color='w', fontsize=17) - ax.text(-0.9, 0, 'E', horizontalalignment='center', verticalalignment='center', color='w', fontsize=17) - ax.text(0, 0.9, 'N', horizontalalignment='center', verticalalignment='center', color='w', fontsize=17) - ax.text(0, -0.9, 'S', horizontalalignment='center', verticalalignment='center', color='w', fontsize=17) - - plt.savefig(f'results/{fname}_sky_calibrated.png', bbox_inches='tight', dpi=200) - plt.close(fig) - - if sky_only: - return img - - npix_x, npix_y = int(ground_resolution * (extent[1] - extent[0])), int(ground_resolution * (extent[3] - extent[2])) - - os.environ["NUMEXPR_NUM_THREADS"] = "3" - - # Select a subset of visibilities, only the lower triangular part - baseline_indices = np.tril_indices(visibilities.shape[0]) - - visibilities_selection = visibilities[baseline_indices] - - img = nearfield_imager(visibilities_selection.flatten()[:, np.newaxis], - np.array(baseline_indices).T, - [freq], npix_x, npix_y, extent, station_xyz, height=height) - - # Correct for taking only lower triangular part - img = np.real(2 * img) - - # Convert bottom left and upper right to PQR just for lofargeo - pmin, qmin, _ = pqr_to_xyz.T @ (np.array([extent[0], extent[2], 0])) - pmax, qmax, _ = pqr_to_xyz.T @ (np.array([extent[1], extent[3], 0])) - lon_center, lat_center, _ = lofargeotiff.pqr_to_longlatheight([0, 0, 0], station_name) - lon_min, lat_min, _ = lofargeotiff.pqr_to_longlatheight([pmin, qmin, 0], station_name) - lon_max, lat_max, _ = lofargeotiff.pqr_to_longlatheight([pmax, qmax, 0], station_name) - - background_map = get_map(lon_min, lon_max, lat_min, lat_max, zoom=map_zoom) - - # Make colors semi-transparent in the lower 3/4 of the scale - cmap = cm.Spectral_r - cmap_with_alpha = cmap(np.arange(cmap.N)) - cmap_with_alpha[:, -1] = np.clip(np.linspace(0, 1.5, cmap.N), 0., 1.) - cmap_with_alpha = ListedColormap(cmap_with_alpha) - - # Plot the resulting image - fig = plt.figure(figsize=(10, 10), constrained_layout=True) - ax = fig.add_subplot(111, ymargin=-0.4) - ax.imshow(background_map, extent=extent) - cimg = ax.imshow(img, origin='lower', cmap=cmap_with_alpha, extent=extent, - alpha=0.7, vmin=ground_vmin, vmax=ground_vmax) - - divider = make_axes_locatable(ax) - cax = divider.append_axes("right", size="5%", pad=0.2, axes_class=maxes.Axes) - cbar = fig.colorbar(cimg, cax=cax, orientation="vertical", format="%.1e") - cbar.set_alpha(1.0) - cbar.draw_all() - # cbar.set_ticks([]) - - ax.set_xlabel('$W-E$ (metres)', fontsize=14) - ax.set_ylabel('$S-N$ (metres)', fontsize=14) - - ax.text(0.5, 1.05, f"Near field image for {station_name}", - fontsize=17, ha='center', va='bottom', transform=ax.transAxes) - ax.text(0.5, 1.02, f"SB {subband} ({freq / 1e6:.1f} MHz), {str(obstime)[:16]}", - fontsize=12, ha='center', va='bottom', transform=ax.transAxes) - - # Change limits to match the original specified extent in the localnorth frame - ax.set_xlim(extent[0], extent[1]) - ax.set_ylim(extent[2], extent[3]) - ax.tick_params(axis='both', which='both', length=0) - - # Place the NSEW coordinate directions - ax.text(0.95, 0.5, 'E', color='w', fontsize=18, transform=ax.transAxes, ha='center', va='center') - ax.text(0.05, 0.5, 'W', color='w', fontsize=18, transform=ax.transAxes, ha='center', va='center') - ax.text(0.5, 0.95, 'N', color='w', fontsize=18, transform=ax.transAxes, ha='center', va='center') - ax.text(0.5, 0.05, 'S', color='w', fontsize=18, transform=ax.transAxes, ha='center', va='center') - - ground_vmin_img, ground_vmax_img = cimg.get_clim() - ax.contour(img, np.linspace(ground_vmin_img, ground_vmax_img, 15), origin='lower', cmap=cm.Greys, - extent=extent, linewidths=0.5, alpha=opacity) - ax.grid(True, alpha=0.3) - plt.savefig(f"results/{fname}_nearfield_calibrated.png", bbox_inches='tight', dpi=200) - plt.close(fig) - - vmin, vmax = cimg.get_clim() - folium_overlay = cmap_with_alpha(Normalize(vmin=vmin, vmax=vmax)(img))[::-1, :] - - maxpixel_ypix, maxpixel_xpix = np.unravel_index(np.argmax(img), img.shape) - maxpixel_x = np.interp(maxpixel_xpix, [0, npix_x], [extent[0], extent[1]]) - maxpixel_y = np.interp(maxpixel_ypix, [0, npix_y], [extent[2], extent[3]]) - [maxpixel_p, maxpixel_q, _] = pqr_to_xyz.T @ np.array([maxpixel_x, maxpixel_y, height]) - maxpixel_lon, maxpixel_lat, _ = lofargeotiff.pqr_to_longlatheight([maxpixel_p, maxpixel_q], station_name) - - # Show location of maximum if not at the image border - if 2 < maxpixel_xpix < npix_x - 2 and 2 < maxpixel_ypix < npix_y - 2: - print(f"Maximum at {maxpixel_x:.0f}m east, {maxpixel_y:.0f}m north of station center " + - f"(lat/long {maxpixel_lat:.5f}, {maxpixel_lon:.5f})") - - obstime = datetime.datetime.strptime(obsdatestr + ":" + obstimestr, '%Y%m%d:%H%M%S') - - tags = {"datafile": xst_filename, - "generated_with": f"lofarimaging v{__version__}", - "caltable": caltable_filename, - "subband": subband, - "frequency": freq, - "extent_xyz": extent, - "height": height, - "station": station_name, - "pixels_per_metre": pixels_per_metre} - if cal_header is not None: - if "CalTableHeader.Observation.Date" in cal_header: - tags["calibration_obsdate"] = cal_header["CalTableHeader.Observation.Date"] - if "CalTableHeader.Calibration.Date" in cal_header: - tags["calibration_date"] = cal_header["CalTableHeader.Calibration.Date"] - if "CalTableHeader.Comment" in cal_header: - tags["calibration_comment"] = cal_header["CalTableHeader.Comment"] - lofargeotiff.write_geotiff(img, f"results/{fname}_nearfield_calibrated.tiff", - (pmin, qmin), (pmax, qmax), stationname=station_name, - obsdate=obstime, tags=tags) - - return make_leaflet_map(folium_overlay, lon_center, lat_center, lon_min, lat_min, lon_max, lat_max) diff --git a/lofarimaging/singlestationutil.py b/lofarimaging/singlestationutil.py new file mode 100644 index 0000000..1de4b67 --- /dev/null +++ b/lofarimaging/singlestationutil.py @@ -0,0 +1,526 @@ +"""Functions for working with LOFAR single station data""" + +import numpy as np +import os +import datetime +import lofargeotiff + +from scipy import ndimage + +from lofarantpos.db import LofarAntennaDatabase + +import matplotlib.pyplot as plt +from matplotlib.ticker import FormatStrFormatter +from matplotlib import cm +from matplotlib.colors import ListedColormap, Normalize +from mpl_toolkits.axes_grid1 import make_axes_locatable +import matplotlib.axes as maxes + +from astropy.coordinates import SkyCoord, GCRS, EarthLocation, AltAz, get_sun +import astropy.units as u +from astropy.time import Time + +import lofarantpos +from packaging import version + +from .maputil import get_map, make_leaflet_map +from .lofarimaging import nearfield_imager, sky_imager, skycoord_to_lmn + + +__all__ = ["sb_from_freq", "freq_from_sb", "find_caltable", "read_caltable", + "rcus_in_station", "read_acm_cube", "get_station_pqr", + "make_ground_image"] + +__version__ = "1.5.0" + +# Configurations for HBA observations with a single dipole activated per tile. +GENERIC_INT_201512 = [0, 5, 3, 1, 8, 3, 12, 15, 10, 13, 11, 5, 12, 12, 5, 2, 10, 8, 0, 3, 5, 1, 4, 0, 11, 6, 2, 4, 9, + 14, 15, 3, 7, 5, 13, 15, 5, 6, 5, 12, 15, 7, 1, 1, 14, 9, 4, 9, 3, 9, 3, 13, 7, 14, 7, 14, 2, 8, + 8, 0, 1, 4, 2, 2, 12, 15, 5, 7, 6, 10, 12, 3, 3, 12, 7, 4, 6, 0, 5, 9, 1, 10, 10, 11, 5, 11, 7, 9, + 7, 6, 4, 4, 15, 4, 1, 15] +GENERIC_CORE_201512 = [0, 10, 4, 3, 14, 0, 5, 5, 3, 13, 10, 3, 12, 2, 7, 15, 6, 14, 7, 5, 7, 9, 0, 15, 0, 10, 4, 3, 14, + 0, 5, 5, 3, 13, 10, 3, 12, 2, 7, 15, 6, 14, 7, 5, 7, 9, 0, 15] +GENERIC_REMOTE_201512 = [0, 13, 12, 4, 11, 11, 7, 8, 2, 7, 11, 2, 10, 2, 6, 3, 8, 3, 1, 7, 1, 15, 13, 1, 11, 1, 12, 7, + 10, 15, 8, 2, 12, 13, 9, 13, 4, 5, 5, 12, 5, 5, 9, 11, 15, 12, 2, 15] + +assert (version.parse(lofarantpos.__version__) >= version.parse("0.4.0")) + + +def sb_from_freq(freq: float, rcu_mode='1'): + """ + Convert subband number to central frequency + + Args: + rcu_mode: rcu mode + freq: frequency in Hz + + Returns: + int: subband number + """ + clock = 200e6 + if rcu_mode == '6': + clock = 160e6 + sb_bandwidth = 0.5 * clock / 512. + freq_offset = 0 + if rcu_mode == '5': + freq_offset = 100e6 + elif rcu_mode == '6': + freq_offset = 160e6 + elif rcu_mode == '7': + freq_offset = 200e6 + sb = round((freq - freq_offset) / sb_bandwidth) + return int(sb) + + +def freq_from_sb(sb: int, rcu_mode='1'): + """ + Convert central frequency to subband number + + Args: + rcu_mode: rcu mode + sb: subband number + + Returns: + float: frequency in Hz + """ + clock = 200e6 + if rcu_mode == '6': + clock = 160e6 + freq_offset = 0 + if rcu_mode == '5': + freq_offset = 100e6 + elif rcu_mode == '6': + freq_offset = 160e6 + elif rcu_mode == '7': + freq_offset = 200e6 + sb_bandwidth = 0.5 * clock / 512. + freq = (sb * sb_bandwidth) + freq_offset + return freq + + +def find_caltable(field_name: str, rcu_mode: str, config_dir='caltables'): + """ + Find the file of a caltable. + + Args: + field_name: Name of the antenna field, e.g. 'DE602LBA' + rcu_mode: Receiver mode for which the calibration table is requested. + Probably should be 'inner' or 'outer' + config_dir: Root directory under which station information is stored in + subdirectories DE602C/etc/, RS106/etc/, ... + + Returns: + str: filename if it exists, None if nothing found + """ + station, field = field_name[0:5].upper(), field_name[5:].upper() + station_number = station[2:5] + + # Map to the correct file depending on the RCU mode + if rcu_mode == 'outer' and 'LBA' in field_name: + filename = os.path.join(config_dir, f"CalTable-{station_number}-LBA_OUTER-10_90.dat") + elif rcu_mode == 'inner' and 'LBA' in field_name: + filename = os.path.join(config_dir, f"CalTable-{station_number}-LBA_INNER-10_90.dat") + else: + filename = os.path.join(config_dir, f"CalTable_{station_number}_mode{rcu_mode}.dat") + + if os.path.exists(filename): + return filename + + # If the original folder structure is kept + if rcu_mode == 'outer' and 'LBA' in field_name: + filename = os.path.join(config_dir, f"{station}/CalTable-{station_number}-LBA_OUTER-10_90.dat") + elif rcu_mode == 'inner' and 'LBA' in field_name: + filename = os.path.join(config_dir, f"{station}/CalTable-{station_number}-LBA_INNER-10_90.dat") + else: + filename = os.path.join(config_dir, f"{station}/CalTable_{station_number}_mode{rcu_mode}.dat") + + if os.path.exists(filename): + return filename + else: + return None + + +def read_caltable(filename: str, num_subbands=512): + """ + Read a station's calibration table. + + Args: + filename: Filename with the caltable + num_subbands: Number of subbands + + Returns: + Tuple(Dict[str, str], np.array): A tuple containing a dict with + the header lines, and a 2D numpy.array of complex numbers + representing the station gain coefficients. + """ + infile = open(filename, 'rb') + + header_lines = [] + + try: + while True: + header_lines.append(infile.readline().decode('utf8').strip()) + if 'HeaderStop' in header_lines[-1]: + break + except UnicodeDecodeError: + # No header; close and open again + infile.close() + infile = open(filename, 'rb') + + caldata = np.fromfile(infile, dtype=np.complex128) + num_rcus = len(caldata) // num_subbands + + infile.close() + + header_dict = {key: val for key, val in [line.split(" = ") + for line in header_lines[1:-1]]} + + return header_dict, caldata.reshape((num_subbands, num_rcus)) + + +def rcus_in_station(station_type: str): + """ + Give the number of RCUs in a station, given its type. + + Args: + station_type: Kind of station that produced the correlation. One of + 'core', 'remote', 'intl'. + """ + return {'core': 96, 'remote': 96, 'intl': 192}[station_type] + + +def read_acm_cube(filename: str, station_type: str): + """ + Read an ACM binary data cube (function from Michiel) + + Args: + filename: File containing the array correlation matrix. + station_type: Kind of station that produced the correlation. One of + 'core', 'remote', 'intl'. + + Returns: + np.array: 3D cube of complex numbers, with indices [time slots, rcu, rcu]. + + Examples: + >>> cube = read_acm_cube('20170720_095816_xst.dat', 'intl') + >>> cube.shape + (29, 192, 192) + """ + num_rcu = rcus_in_station(station_type) + data = np.fromfile(filename, dtype=np.complex128) + time_slots = int(len(data) / num_rcu / num_rcu) + return data.reshape((time_slots, num_rcu, num_rcu)) + + +def get_station_pqr(station_name: str, station_type: str, array_type: str, db): + if 'LBA' in station_name: + # Get the PQR positions for an individual station + station_pqr = db.antenna_pqr(station_name) + + # Exception: for Dutch stations (sparse not yet accommodated) + if (station_type == 'core' or station_type == 'remote') and array_type == 'inner': + station_pqr = station_pqr[0:48, :] + elif (station_type == 'core' or station_type == 'remote') and array_type == 'outer': + station_pqr = station_pqr[48:, :] + elif 'HBA' in station_name: + selected_dipole_config = { + 'intl': GENERIC_INT_201512, 'remote': GENERIC_REMOTE_201512, 'core': GENERIC_CORE_201512 + } + selected_dipoles = selected_dipole_config[station_type] + \ + np.arange(len(selected_dipole_config[station_type])) * 16 + station_pqr = db.hba_dipole_pqr(station_name)[selected_dipoles] + else: + raise RuntimeError("Station name did not contain LBA or HBA, could not load antenna positions") + + return station_pqr.astype('float32') + + +def make_ground_image(xst_filename, + station_name, + caltable_dir, + extent=None, + pixels_per_metre=0.5, + sky_vmin=None, + sky_vmax=None, + ground_vmin=None, + ground_vmax=None, + height=1.5, + map_zoom=19, + sky_only=False, + opacity=0.6): + """Make a ground image""" + cubename = os.path.basename(xst_filename) + + if extent is None: + extent = [-150, 150, -150, 150] + + if station_name[0] == "C": + station_type = "core" + elif station_name[0] == "R" or station_name[:5] == "PL611": + station_type = "remote" + else: + station_type = "intl" + + try: + os.mkdir('results') + except FileExistsError: + pass + + # Distill metadata from filename + obsdatestr, obstimestr, _, rcu_mode, _, subbandname = cubename.rstrip(".dat").split("_") + subband = int(subbandname[2:]) + + # Needed for NL stations: inner (rcu_mode 3/4), outer (rcu_mode 1/2), (sparse tbd) + # Should be set to 'inner' if station type = 'intl' + array_type = None + if rcu_mode in ('1', '2'): + array_type = 'outer' + elif rcu_mode in ('3', '4'): + array_type = 'inner' + elif rcu_mode in ('5', '6', '7'): + array_type = rcu_mode + else: + raise Exception("Unexpected rcu_mode: ", rcu_mode) + + # Get the data + fname = f"{obsdatestr}_{obstimestr}_{station_name}_SB{subband}" + + npix_l, npix_m = 131, 131 + freq = freq_from_sb(subband, rcu_mode=rcu_mode) + + # Which slice in time to visualise + timestep = 0 + + # For ground imaging + ground_resolution = pixels_per_metre # pixels per metre for ground_imaging, default is 0.5 pixel/metre + + obstime = datetime.datetime.strptime(obsdatestr + ":" + obstimestr, '%Y%m%d:%H%M%S') + + cube = read_acm_cube(xst_filename, station_type) + + # Apply calibration + + caltable_filename = find_caltable(station_name, rcu_mode=array_type, + config_dir=caltable_dir) + + cal_header = None + if caltable_filename is None: + print('No calibration table found... cube remains uncalibrated!') + else: + cal_header, cal_data = read_caltable(caltable_filename) + + rcu_gains = cal_data[subband, :] + rcu_gains = np.array(rcu_gains, dtype=np.complex64) + gain_matrix = rcu_gains[np.newaxis, :] * np.conj(rcu_gains[:, np.newaxis]) + cube = cube / gain_matrix + + # Split into the XX and YY polarisations (RCUs) + # This needs to be modified in future for LBA sparse + cube_xx = cube[:, 0::2, 0::2] + cube_yy = cube[:, 1::2, 1::2] + visibilities_all = cube_xx + cube_yy + + # Stokes I for specified timestep + visibilities = visibilities_all[timestep] + + # Setup the database + db = LofarAntennaDatabase() + + station_pqr = get_station_pqr(station_name, station_type, array_type, db) + + # Rotate station_pqr to a north-oriented xyz frame, where y points North, in a plane through the station. + rotation = db.rotation_from_north(station_name) + + pqr_to_xyz = np.array([[np.cos(-rotation), -np.sin(-rotation), 0], + [np.sin(-rotation), np.cos(-rotation), 0], + [0, 0, 1]]) + + station_xyz = (pqr_to_xyz @ station_pqr.T).T + + baselines = station_xyz[:, np.newaxis, :] - station_xyz[np.newaxis, :, :] + + # Make a sky image, by numerically Fourier-transforming from visibilities to image plane + from matplotlib.patches import Circle + + # Fourier transform, and account for the rotation (rotation is positive in this space) + # visibilities = cube_xx[2,:,:] + img = sky_imager(visibilities, baselines, freq, npix_l, npix_m) + img = ndimage.interpolation.rotate(img, -rotation, reshape=False, mode='nearest') + + obstime_astropy = Time(obstime) + # Determine positions of Cas A and Cyg A + station_earthlocation = EarthLocation.from_geocentric(*(db.phase_centres[station_name] * u.m)) + zenith = AltAz(az=0 * u.deg, alt=90 * u.deg, obstime=obstime_astropy, + location=station_earthlocation).transform_to(GCRS) + + marked_bodies = { + 'Cas A': SkyCoord(ra=350.85 * u.deg, dec=58.815 * u.deg), + 'Cyg A': SkyCoord(ra=299.868 * u.deg, dec=40.734 * u.deg), + # 'Per A': SkyCoord.from_name("Perseus A"), + # 'Her A': SkyCoord.from_name("Hercules A"), + # 'Cen A': SkyCoord.from_name("Centaurus A"), + # '?': SkyCoord.from_name("J101415.9+105106"), + # '3C295': SkyCoord.from_name("3C295"), + # 'Moon': get_moon(obstime_astropy, location=station_earthlocation).transform_to(GCRS), + 'Sun': get_sun(obstime_astropy) + # '3C196': SkyCoord.from_name("3C196") + } + + marked_bodies_lmn = {} + for body_name, body_coord in marked_bodies.items(): + # print(body_name, body_coord.separation(zenith), body_coord.separation(zenith)) + if body_coord.transform_to(AltAz(location=station_earthlocation, obstime=obstime_astropy)).alt > 0: + marked_bodies_lmn[body_name] = skycoord_to_lmn(marked_bodies[body_name], zenith) + + # Plot the resulting sky image + fig, ax = plt.subplots(figsize=(10, 10)) + + circle1 = Circle((0, 0), 1.0, edgecolor='k', fill=False, facecolor='none', alpha=0.3) + ax.add_artist(circle1) + + cimg = ax.imshow(img, origin='lower', cmap=cm.Spectral_r, extent=(1, -1, -1, 1), + clip_path=circle1, clip_on=True, vmin=sky_vmin, vmax=sky_vmax) + divider = make_axes_locatable(ax) + cax = divider.append_axes("right", size="5%", pad=0.2, axes_class=maxes.Axes) + fig.colorbar(cimg, cax=cax, orientation="vertical", format="%.1e") + + ax.set_xlim(1, -1) + + ax.set_xticks(np.arange(-1, 1.1, 0.5)) + ax.xaxis.set_major_formatter(FormatStrFormatter('%.1f')) + ax.set_yticks(np.arange(-1, 1.1, 0.5)) + ax.yaxis.set_major_formatter(FormatStrFormatter('%.1f')) + + # Labels + ax.set_xlabel('$ℓ$', fontsize=14) + ax.set_ylabel('$m$', fontsize=14) + + for body_name, lmn in marked_bodies_lmn.items(): + ax.plot([lmn[0]], [lmn[1]], marker='x', color='white', mew=0.5) + ax.annotate(body_name, (lmn[0], lmn[1])) + + ax.text(0.5, 1.05, f"Sky image for {station_name}", + fontsize=17, ha='center', va='bottom', transform=ax.transAxes) + ax.text(0.5, 1.02, f"SB {subband} ({freq / 1e6:.1f} MHz), {str(obstime)[:16]}", + fontsize=12, ha='center', va='bottom', transform=ax.transAxes) + + # Plot the compass directions + ax.text(0.9, 0, 'W', horizontalalignment='center', verticalalignment='center', color='w', fontsize=17) + ax.text(-0.9, 0, 'E', horizontalalignment='center', verticalalignment='center', color='w', fontsize=17) + ax.text(0, 0.9, 'N', horizontalalignment='center', verticalalignment='center', color='w', fontsize=17) + ax.text(0, -0.9, 'S', horizontalalignment='center', verticalalignment='center', color='w', fontsize=17) + + plt.savefig(f'results/{fname}_sky_calibrated.png', bbox_inches='tight', dpi=200) + plt.close(fig) + + if sky_only: + return img + + npix_x, npix_y = int(ground_resolution * (extent[1] - extent[0])), int(ground_resolution * (extent[3] - extent[2])) + + os.environ["NUMEXPR_NUM_THREADS"] = "3" + + # Select a subset of visibilities, only the lower triangular part + baseline_indices = np.tril_indices(visibilities.shape[0]) + + visibilities_selection = visibilities[baseline_indices] + + img = nearfield_imager(visibilities_selection.flatten()[:, np.newaxis], + np.array(baseline_indices).T, + [freq], npix_x, npix_y, extent, station_xyz, height=height) + + # Correct for taking only lower triangular part + img = np.real(2 * img) + + # Convert bottom left and upper right to PQR just for lofargeo + pmin, qmin, _ = pqr_to_xyz.T @ (np.array([extent[0], extent[2], 0])) + pmax, qmax, _ = pqr_to_xyz.T @ (np.array([extent[1], extent[3], 0])) + lon_center, lat_center, _ = lofargeotiff.pqr_to_longlatheight([0, 0, 0], station_name) + lon_min, lat_min, _ = lofargeotiff.pqr_to_longlatheight([pmin, qmin, 0], station_name) + lon_max, lat_max, _ = lofargeotiff.pqr_to_longlatheight([pmax, qmax, 0], station_name) + + background_map = get_map(lon_min, lon_max, lat_min, lat_max, zoom=map_zoom) + + # Make colors semi-transparent in the lower 3/4 of the scale + cmap = cm.Spectral_r + cmap_with_alpha = cmap(np.arange(cmap.N)) + cmap_with_alpha[:, -1] = np.clip(np.linspace(0, 1.5, cmap.N), 0., 1.) + cmap_with_alpha = ListedColormap(cmap_with_alpha) + + # Plot the resulting image + fig = plt.figure(figsize=(10, 10), constrained_layout=True) + ax = fig.add_subplot(111, ymargin=-0.4) + ax.imshow(background_map, extent=extent) + cimg = ax.imshow(img, origin='lower', cmap=cmap_with_alpha, extent=extent, + alpha=0.7, vmin=ground_vmin, vmax=ground_vmax) + + divider = make_axes_locatable(ax) + cax = divider.append_axes("right", size="5%", pad=0.2, axes_class=maxes.Axes) + cbar = fig.colorbar(cimg, cax=cax, orientation="vertical", format="%.1e") + cbar.set_alpha(1.0) + cbar.draw_all() + # cbar.set_ticks([]) + + ax.set_xlabel('$W-E$ (metres)', fontsize=14) + ax.set_ylabel('$S-N$ (metres)', fontsize=14) + + ax.text(0.5, 1.05, f"Near field image for {station_name}", + fontsize=17, ha='center', va='bottom', transform=ax.transAxes) + ax.text(0.5, 1.02, f"SB {subband} ({freq / 1e6:.1f} MHz), {str(obstime)[:16]}", + fontsize=12, ha='center', va='bottom', transform=ax.transAxes) + + # Change limits to match the original specified extent in the localnorth frame + ax.set_xlim(extent[0], extent[1]) + ax.set_ylim(extent[2], extent[3]) + ax.tick_params(axis='both', which='both', length=0) + + # Place the NSEW coordinate directions + ax.text(0.95, 0.5, 'E', color='w', fontsize=18, transform=ax.transAxes, ha='center', va='center') + ax.text(0.05, 0.5, 'W', color='w', fontsize=18, transform=ax.transAxes, ha='center', va='center') + ax.text(0.5, 0.95, 'N', color='w', fontsize=18, transform=ax.transAxes, ha='center', va='center') + ax.text(0.5, 0.05, 'S', color='w', fontsize=18, transform=ax.transAxes, ha='center', va='center') + + ground_vmin_img, ground_vmax_img = cimg.get_clim() + ax.contour(img, np.linspace(ground_vmin_img, ground_vmax_img, 15), origin='lower', cmap=cm.Greys, + extent=extent, linewidths=0.5, alpha=opacity) + ax.grid(True, alpha=0.3) + plt.savefig(f"results/{fname}_nearfield_calibrated.png", bbox_inches='tight', dpi=200) + plt.close(fig) + + vmin, vmax = cimg.get_clim() + folium_overlay = cmap_with_alpha(Normalize(vmin=vmin, vmax=vmax)(img))[::-1, :] + + maxpixel_ypix, maxpixel_xpix = np.unravel_index(np.argmax(img), img.shape) + maxpixel_x = np.interp(maxpixel_xpix, [0, npix_x], [extent[0], extent[1]]) + maxpixel_y = np.interp(maxpixel_ypix, [0, npix_y], [extent[2], extent[3]]) + [maxpixel_p, maxpixel_q, _] = pqr_to_xyz.T @ np.array([maxpixel_x, maxpixel_y, height]) + maxpixel_lon, maxpixel_lat, _ = lofargeotiff.pqr_to_longlatheight([maxpixel_p, maxpixel_q], station_name) + + # Show location of maximum if not at the image border + if 2 < maxpixel_xpix < npix_x - 2 and 2 < maxpixel_ypix < npix_y - 2: + print(f"Maximum at {maxpixel_x:.0f}m east, {maxpixel_y:.0f}m north of station center " + + f"(lat/long {maxpixel_lat:.5f}, {maxpixel_lon:.5f})") + + obstime = datetime.datetime.strptime(obsdatestr + ":" + obstimestr, '%Y%m%d:%H%M%S') + + tags = {"datafile": xst_filename, + "generated_with": f"lofarimaging v{__version__}", + "caltable": caltable_filename, + "subband": subband, + "frequency": freq, + "extent_xyz": extent, + "height": height, + "station": station_name, + "pixels_per_metre": pixels_per_metre} + if cal_header is not None: + if "CalTableHeader.Observation.Date" in cal_header: + tags["calibration_obsdate"] = cal_header["CalTableHeader.Observation.Date"] + if "CalTableHeader.Calibration.Date" in cal_header: + tags["calibration_date"] = cal_header["CalTableHeader.Calibration.Date"] + if "CalTableHeader.Comment" in cal_header: + tags["calibration_comment"] = cal_header["CalTableHeader.Comment"] + lofargeotiff.write_geotiff(img, f"results/{fname}_nearfield_calibrated.tiff", + (pmin, qmin), (pmax, qmax), stationname=station_name, + obsdate=obstime, tags=tags) + + return make_leaflet_map(folium_overlay, lon_center, lat_center, lon_min, lat_min, lon_max, lat_max)