1.. _example_beam_error_spot_frequencies: 2 3Beam errors at spot frequencies 4=============================== 5 6Summary 7------- 8 9This script demonstrates how the OSKAR Python interface can be used to 10run simulations which compare the performance of a set of telescope models. 11 12The idea was to investigate how errors in the station beams would affect the 13quality of images produced by the SKA-Low telescope at a range of spot 14frequencies across the band, at 50, 70, 110, 137, 160, 230 and 320 MHz. 15A number of "corrupted" visibility data sets were generated at each frequency 16by introducing progressively larger errors in the station beams, and these were 17subtracted from a reference "perfect" data set generated without any 18corruptions. Dirty images were then made of the residual visibilities to 19estimate the amount of noise introduced by the beam errors. 20 21The station beam errors were generated by varying the complex beamforming 22weights randomly within each of the 512 stations over a defined range. 23These errors were kept static throughout each observation to simulate the 24effects of residual calibration errors at the station level. 25Element phase errors were calculated as a function of frequency using a 26Gaussian-distributed cable length tolerance of 15 mm. 27Element gain errors were generated from Gaussian distributions with 28widths varying between 0.005 dB and 0.640 dB. 29 30The sky models and other observation parameters were kept the same for each 31of the data sets. Three different fields were used, with each populated by 32sources from the GLEAM extragalactic catalogue, and the 10 brightest "A-team" 33sources in the sky were also added. 34All observations were set up to span 4 hours and to be symmetric about the 35meridian, with station beams and visibilities evaluated every 60 seconds to 36generate a representative aperture-plane coverage. 37 38These results were written up and presented as part of the SKA System CDR 39documentation pack. 40 41The simulation script 42--------------------- 43 44The script starts (in ``main``) by declaring the axes of the parameter space 45to iterate over: the seven frequencies of interest, the eight gain error 46values, and the three target fields. The interesting parts of the script are 47mainly in three functions: 48 49- The ``run_set`` function uses three nested loops to map out the parameter 50 space, steering the script by creating the :class:`oskar.Sky` model for 51 each field, updating a :class:`oskar.Telescope` model, and setting the 52 relevant parameters in a :class:`oskar.SettingsTree` before simulating the 53 visibilities using :class:`oskar.Interferometer` (in ``make_vis_data``). 54 55- The ``make_diff_image_stats`` function generates the residual visibilities 56 by taking the difference of the data in two :class:`oskar.VisBlock` classes 57 (the "perfect" and "corrupted" data), which are read from two files. 58 Note that a ``ThreadPoolExecutor`` is used to read data from the two files 59 concurrently to improve performance. The residual visibilities are generated 60 using ``block1.cross_correlations()[...] -= block2.cross_correlations()`` 61 and these are passed directly (still inside ``block1``) 62 to :class:`oskar.Imager` via its 63 :meth:`update_from_block() <oskar.Imager.update_from_block()>` method. 64 The images are then returned straight back for analysis using functions 65 from ``numpy``, and the image statistics are stored in a Python dictionary 66 called ``results``. 67 68- Finally, the results dictionary is plotted using ``matplotlib`` in 69 the ``make_plot`` function. 70 71 72.. code-block:: python 73 74 #!/usr/bin/env python3 75 76 """Script to run LOW station beam error simulations at spot frequencies.""" 77 78 from __future__ import print_function, division 79 import concurrent.futures 80 import json 81 import logging 82 import os 83 import sys 84 85 from astropy.io import fits 86 from astropy.time import Time, TimeDelta 87 import matplotlib 88 matplotlib.use('Agg') 89 # pylint: disable=wrong-import-position 90 import matplotlib.pyplot as plt 91 import numpy 92 import oskar 93 94 95 LOG = logging.getLogger() 96 97 98 def bright_sources(): 99 """Returns a list of bright A-team sources.""" 100 # Sgr A: guesstimates only! 101 # For A: data from the Molonglo Southern 4 Jy sample (VizieR). 102 # Others from GLEAM reference paper, Hurley-Walker et al. (2017), Table 2. 103 return numpy.array(( 104 [266.41683, -29.00781, 2000,0,0,0, 0, 0, 0, 3600, 3600, 0], 105 [ 50.67375, -37.20833, 528,0,0,0, 178e6, -0.51, 0, 0, 0, 0], # For 106 [201.36667, -43.01917, 1370,0,0,0, 200e6, -0.50, 0, 0, 0, 0], # Cen 107 [139.52500, -12.09556, 280,0,0,0, 200e6, -0.96, 0, 0, 0, 0], # Hyd 108 [ 79.95833, -45.77889, 390,0,0,0, 200e6, -0.99, 0, 0, 0, 0], # Pic 109 [252.78333, 4.99250, 377,0,0,0, 200e6, -1.07, 0, 0, 0, 0], # Her 110 [187.70417, 12.39111, 861,0,0,0, 200e6, -0.86, 0, 0, 0, 0], # Vir 111 [ 83.63333, 22.01444, 1340,0,0,0, 200e6, -0.22, 0, 0, 0, 0], # Tau 112 [299.86667, 40.73389, 7920,0,0,0, 200e6, -0.78, 0, 0, 0, 0], # Cyg 113 [350.86667, 58.81167, 11900,0,0,0, 200e6, -0.41, 0, 0, 0, 0] # Cas 114 )) 115 116 117 def get_start_time(ra0_deg, length_sec): 118 """Returns optimal start time for field RA and observation length.""" 119 t = Time('2000-01-01 00:00:00', scale='utc', location=('116.764d', '0d')) 120 dt_hours = 24.0 - t.sidereal_time('apparent').hour + (ra0_deg / 15.0) 121 start = t + TimeDelta(dt_hours * 3600.0 - length_sec / 2.0, format='sec') 122 return start.value 123 124 125 def make_vis_data(settings, sky, tel): 126 """Run simulation using supplied settings.""" 127 if os.path.exists(settings['interferometer/oskar_vis_filename']): 128 LOG.info("Skipping simulation, as output data already exist.") 129 return 130 LOG.info("Simulating %s", settings['interferometer/oskar_vis_filename']) 131 sim = oskar.Interferometer(settings=settings) 132 sim.set_sky_model(sky) 133 sim.set_telescope_model(tel) 134 sim.run() 135 136 137 def make_sky_model(sky0, settings, radius_deg, flux_min_outer_jy): 138 """Filter sky model. 139 140 Includes all sources within the given radius, and sources above the 141 specified flux outside this radius. 142 """ 143 # Get pointing centre. 144 ra0_deg = float(settings['observation/phase_centre_ra_deg']) 145 dec0_deg = float(settings['observation/phase_centre_dec_deg']) 146 147 # Create "inner" and "outer" sky models. 148 sky_inner = sky0.create_copy() 149 sky_outer = sky0.create_copy() 150 sky_inner.filter_by_radius(0.0, radius_deg, ra0_deg, dec0_deg) 151 sky_outer.filter_by_radius(radius_deg, 180.0, ra0_deg, dec0_deg) 152 sky_outer.filter_by_flux(flux_min_outer_jy, 1e9) 153 LOG.info("Number of sources in sky0: %d", sky0.num_sources) 154 LOG.info("Number of sources in inner sky model: %d", sky_inner.num_sources) 155 LOG.info("Number of sources in outer sky model above %.3f Jy: %d", 156 flux_min_outer_jy, sky_outer.num_sources) 157 sky_outer.append(sky_inner) 158 LOG.info("Number of sources in output sky model: %d", sky_outer.num_sources) 159 return sky_outer 160 161 162 def make_diff_image_stats(filename1, filename2, use_w_projection, 163 out_image_root=None): 164 """Make an image of the difference between two visibility data sets. 165 166 This function assumes that the observation parameters for both data sets 167 are identical. (It will fail horribly otherwise!) 168 """ 169 # Set up an imager. 170 (hdr1, handle1) = oskar.VisHeader.read(filename1) 171 (hdr2, handle2) = oskar.VisHeader.read(filename2) 172 frequency_hz = hdr1.freq_start_hz 173 fov_ref_frequency_hz = 140e6 174 fov_ref_deg = 5.0 175 fov_deg = fov_ref_deg * (fov_ref_frequency_hz / frequency_hz) 176 imager = oskar.Imager(precision='double') 177 imager.set(fov_deg=fov_deg, image_size=8192, 178 fft_on_gpu=True, grid_on_gpu=True) 179 if out_image_root is not None: 180 imager.output_root = out_image_root 181 182 LOG.info("Imaging differences between '%s' and '%s'", filename1, filename2) 183 block1 = oskar.VisBlock.create_from_header(hdr1) 184 block2 = oskar.VisBlock.create_from_header(hdr2) 185 if hdr1.num_blocks != hdr2.num_blocks: 186 raise RuntimeError("'%s' and '%s' have different dimensions!" % 187 (filename1, filename2)) 188 if use_w_projection: 189 imager.set(algorithm='W-projection') 190 imager.coords_only = True 191 for i_block in range(hdr1.num_blocks): 192 block1.read(hdr1, handle1, i_block) 193 imager.update_from_block(hdr1, block1) 194 imager.coords_only = False 195 imager.check_init() 196 LOG.info("Using %d W-planes", imager.num_w_planes) 197 executor = concurrent.futures.ThreadPoolExecutor(2) 198 for i_block in range(hdr1.num_blocks): 199 tasks_read = [] 200 tasks_read.append(executor.submit(block1.read, hdr1, handle1, i_block)) 201 tasks_read.append(executor.submit(block2.read, hdr2, handle2, i_block)) 202 concurrent.futures.wait(tasks_read) 203 block1.cross_correlations()[...] -= block2.cross_correlations() 204 imager.update_from_block(hdr1, block1) 205 del handle1, handle2, hdr1, hdr2, block1, block2 206 207 # Finalise image and return it to Python. 208 output = imager.finalise(return_images=1) 209 image = output['images'][0] 210 211 LOG.info("Generating image statistics") 212 image_size = imager.image_size 213 box_size = int(0.1 * image_size) 214 centre = image[ 215 (image_size - box_size)//2:(image_size + box_size)//2, 216 (image_size - box_size)//2:(image_size + box_size)//2] 217 del imager 218 return { 219 'image_medianabs': numpy.median(numpy.abs(image)), 220 'image_mean': numpy.mean(image), 221 'image_std': numpy.std(image), 222 'image_rms': numpy.sqrt(numpy.mean(image**2)), 223 'image_centre_mean': numpy.mean(centre), 224 'image_centre_std': numpy.std(centre), 225 'image_centre_rms': numpy.sqrt(numpy.mean(centre**2)) 226 } 227 228 229 def make_plot(prefix, field_name, metric_key, results, 230 axis_freq, axis_gain): 231 """Plot selected results.""" 232 # Get data for contour plot. 233 X, Y = numpy.meshgrid(axis_freq, axis_gain) 234 Z = numpy.zeros(X.shape) 235 for freq, gain, z in numpy.nditer([X, Y, Z], op_flags=['readwrite']): 236 key = '%s_%s_%d_MHz_%.3f_dB' % (prefix, field_name, freq, gain) 237 if key in results: 238 z[...] = numpy.log10(results[key][metric_key]) 239 ax1 = plt.subplot(111) 240 ax1.set_yscale('log') 241 242 # Scatter plot. 243 sp = ax1.scatter(X, Y, c=Z, cmap='plasma') 244 245 # Contour plot. 246 cp = ax1.contour(X, Y, Z, cmap='plasma') 247 levels = cp.levels 248 print(prefix, field_name, len(levels)) 249 if len(levels) > 9: 250 levels = levels[::2] 251 clabels = plt.clabel(cp, levels, inline=False, fontsize=10, fmt='%1.1f') 252 for txt in clabels: 253 txt.set_bbox(dict(facecolor='white', edgecolor='none', pad=1)) 254 255 # Title and axis labels. 256 metric_name = '[ UNKNOWN ]' 257 if metric_key == 'image_centre_rms': 258 metric_name = 'Central RMS [Jy/beam]' 259 elif metric_key == 'image_medianabs': 260 metric_name = 'MEDIAN(ABS(image)) [Jy/beam]' 261 sky_model = 'GLEAM' 262 if 'A-team' in prefix: 263 sky_model = sky_model + ' + A-team' 264 plt.title('%s for %s field (%s)' % (metric_name, field_name, sky_model)) 265 plt.xlabel('Frequency [MHz]') 266 plt.ylabel('Element gain standard deviation [dB]') 267 cbar = plt.colorbar(sp) 268 cbar.set_label('log10(%s)' % metric_name) 269 plt.savefig('%s_%s_%s.png' % (prefix, field_name, metric_key)) 270 plt.close('all') 271 272 273 def run_single(prefix_field, settings, sky, tel, 274 freq_MHz, gain_std_dB, out0_name, results): 275 """Run a single simulation and generate image statistics for it.""" 276 out = '%s_%d_MHz_%.3f_dB' % (prefix_field, freq_MHz, gain_std_dB) 277 if out in results: 278 LOG.info("Using cached results for '%s'", out) 279 return 280 out_name = out + '.vis' 281 gain_std = numpy.power(10.0, gain_std_dB / 20.0) - 1.0 282 tel.override_element_gains(1.0, gain_std) 283 tel.override_element_cable_length_errors(0.015) 284 settings['interferometer/oskar_vis_filename'] = out_name 285 make_vis_data(settings, sky, tel) 286 out_image_root = out 287 use_w_projection = True 288 if str(settings['interferometer/ignore_w_components']).lower() == 'true': 289 use_w_projection = False 290 results[out] = make_diff_image_stats(out0_name, out_name, use_w_projection, 291 out_image_root) 292 293 294 def run_set(prefix, base_settings, fields, axis_freq, axis_gain, plot_only): 295 """Runs a set of simulations.""" 296 if not plot_only: 297 # Load base telescope model. 298 settings = oskar.SettingsTree('oskar_sim_interferometer') 299 settings.from_dict(base_settings) 300 tel = oskar.Telescope(settings=settings) 301 302 # Load base sky model 303 sky0 = oskar.Sky() 304 if 'GLEAM' in prefix: 305 # Load GLEAM catalogue from FITS binary table. 306 hdulist = fits.open('GLEAM_EGC.fits') 307 # pylint: disable=no-member 308 cols = hdulist[1].data[0].array 309 data = numpy.column_stack( 310 (cols['RAJ2000'], cols['DEJ2000'], cols['peak_flux_wide'])) 311 data = data[data[:, 2].argsort()[::-1]] 312 sky_gleam = oskar.Sky.from_array(data) 313 sky0.append(sky_gleam) 314 if 'A-team' in prefix: 315 sky_bright = oskar.Sky.from_array(bright_sources()) 316 sky0.append(sky_bright) 317 318 # Iterate over fields. 319 for field_name, field in fields.items(): 320 # Load result set, if it exists. 321 prefix_field = prefix + '_' + field_name 322 results = {} 323 json_file = prefix_field + '_results.json' 324 if os.path.exists(json_file): 325 with open(json_file, 'r') as input_file: 326 results = json.load(input_file) 327 328 # Iterate over frequencies. 329 if not plot_only: 330 for freq_MHz in axis_freq: 331 # Update settings for field. 332 settings_dict = base_settings.copy() 333 settings_dict.update(field) 334 settings.from_dict(settings_dict) 335 ra_deg = float(settings['observation/phase_centre_ra_deg']) 336 dec_deg = float(settings['observation/phase_centre_dec_deg']) 337 length_sec = float(settings['observation/length']) 338 settings['observation/start_frequency_hz'] = str(freq_MHz * 1e6) 339 settings['observation/start_time_utc'] = get_start_time( 340 ra_deg, length_sec) 341 tel.set_phase_centre(ra_deg, dec_deg) 342 343 # Create the sky model. 344 sky = make_sky_model(sky0, settings, 20.0, 10.0) 345 settings['interferometer/ignore_w_components'] = 'true' 346 if 'A-team' in prefix: 347 settings['interferometer/ignore_w_components'] = 'false' 348 349 # Simulate the 'perfect' case. 350 tel.override_element_gains(1.0, 0.0) 351 tel.override_element_cable_length_errors(0.0) 352 out0_name = '%s_%d_MHz_no_errors.vis' % (prefix_field, freq_MHz) 353 settings['interferometer/oskar_vis_filename'] = out0_name 354 make_vis_data(settings, sky, tel) 355 356 # Simulate the error cases. 357 for gain_std_dB in axis_gain: 358 run_single(prefix_field, settings, sky, tel, 359 freq_MHz, gain_std_dB, out0_name, results) 360 361 # Generate plot for the field. 362 make_plot(prefix, field_name, 'image_centre_rms', 363 results, axis_freq, axis_gain) 364 make_plot(prefix, field_name, 'image_medianabs', 365 results, axis_freq, axis_gain) 366 367 # Save result set. 368 with open(json_file, 'w') as output_file: 369 json.dump(results, output_file, indent=4) 370 371 372 def main(): 373 """Main function.""" 374 handler = logging.StreamHandler(sys.stdout) 375 formatter = logging.Formatter('%(asctime)s - %(levelname)s - %(message)s') 376 handler.setFormatter(formatter) 377 LOG.addHandler(handler) 378 LOG.setLevel(logging.INFO) 379 380 # Define common settings. 381 base_settings = { 382 'simulator': { 383 'double_precision': 'true', 384 'use_gpus': 'true', 385 'max_sources_per_chunk': '23000' 386 }, 387 'observation' : { 388 'frequency_inc_hz': '100e3', 389 'length': '14400.0', 390 'num_time_steps': '240' 391 }, 392 'telescope': { 393 'input_directory': 'SKA1-LOW_SKO-0000422_Rev3_38m_SKALA4_spot_frequencies.tm' 394 }, 395 'interferometer': { 396 'channel_bandwidth_hz': '100e3', 397 'time_average_sec': '1.0', 398 'max_time_samples_per_block': '4' 399 } 400 } 401 402 # Define axes of parameter space. 403 fields = { 404 'EoR0': { 405 'observation/phase_centre_ra_deg': '0.0', 406 'observation/phase_centre_dec_deg': '-27.0' 407 }, 408 'EoR1': { 409 'observation/phase_centre_ra_deg': '60.0', 410 'observation/phase_centre_dec_deg': '-30.0' 411 }, 412 'EoR2': { 413 'observation/phase_centre_ra_deg': '170.0', 414 'observation/phase_centre_dec_deg': '-10.0' 415 } 416 } 417 axis_freq = [50, 70, 110, 137, 160, 230, 320] 418 axis_gain = [0.005, 0.01, 0.02, 0.04, 0.08, 0.16, 0.32, 0.64] 419 420 # GLEAM + A-team sky model simulations. 421 plot_only = False 422 run_set('GLEAM_A-team', base_settings, 423 fields, axis_freq, axis_gain, plot_only) 424 425 426 if __name__ == '__main__': 427 main() 428