From c96bf22bb3db6a5a8a93762323e50152ef314afb Mon Sep 17 00:00:00 2001 From: gigjozsa Date: Wed, 3 Nov 2021 17:37:35 +0100 Subject: [PATCH] Added notebook to calculate circularly averaged attenuation --- notebooks/03h_radial_attenuation.ipynb | 1636 ++++++++++++++++++++++++ 1 file changed, 1636 insertions(+) create mode 100644 notebooks/03h_radial_attenuation.ipynb diff --git a/notebooks/03h_radial_attenuation.ipynb b/notebooks/03h_radial_attenuation.ipynb new file mode 100644 index 000000000..e33f87db7 --- /dev/null +++ b/notebooks/03h_radial_attenuation.ipynb @@ -0,0 +1,1636 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Calculate average path propagation/attenuation according to ITU-R P.452 (16)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## License\n", + "\n", + "```\n", + "Calculate path propagation/attenuation according to ITU-R P.452 (16)\n", + "in presence of geography and under flat-earth assumption and comparison.\n", + "Copyright (C) 2021+ Gyula Jozsa (gjozsa@mipfr.de) \n", + "& Benjamin Winkel (bwinkel@mpifr.de)\n", + "\n", + "This program is free software; you can redistribute it and/or\n", + "modify it under the terms of the GNU General Public License\n", + "as published by the Free Software Foundation; either version 2\n", + "of the License, or (at your option) any later version.\n", + "\n", + "This program is distributed in the hope that it will be useful,\n", + "but WITHOUT ANY WARRANTY; without even the implied warranty of\n", + "MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the\n", + "GNU General Public License for more details.\n", + "\n", + "You should have received a copy of the GNU General Public License\n", + "along with this program; if not, write to the Free Software\n", + "Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.\n", + "```" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false, + "jupyter": { + "outputs_hidden": false + } + }, + "outputs": [], + "source": [ + "%matplotlib inline" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false, + "jupyter": { + "outputs_hidden": false + } + }, + "outputs": [], + "source": [ + "from collections import namedtuple, OrderedDict\n", + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "from pycraf import pathprof\n", + "from pycraf import conversions as cnv\n", + "from astropy import units as u\n", + "from astropy.utils.console import ProgressBar\n", + "import os.path\n", + "import copy\n", + "import time" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "With pycraf it is very easy to query SRTM data for arbitrary positions. However, you will need to have to download such SRTM data, either in advance or during your session. Please see the [documentation](https://bwinkel.github.io/pycraf/pathprof/working_with_srtm.html) for more details." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Use the following, if you already have the SRTMDATA environment variable set, but not downloaded any \".hgt\" files so far." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "# pathprof.SrtmConf.set(download='missing', server='viewpano')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "If you neither have the environment variable set nor \".hgt\" files downloaded (you can use any target directory, of course):" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "# pathprof.SrtmConf.set(download='missing', server='viewpano', srtm_dir='.')" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "tags": [] + }, + "source": [ + "If you do not have the environment variable set but you have all .hgt files downloaded:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "pathprof.SrtmConf.set(srtm_dir='/Users/jozsa/software/srtm/all')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Averaging path propagation loss maps\n", + "\n", + "Often it is desired to generate maps of path attenuation values, e.g.,\n", + "to quickly determine regions where the necessary separations between a\n", + "potential interferer and the victim terminal would be too small,\n", + "potentially leading to radio frequency interference (RFI). In this\n", + "example we use attenuation maps to average the attenuation at a\n", + "certain distance from the over all directions, for a set of telescopes\n", + "and a set of frequencies. We will use our knowledge from notebook 03c\n", + "to do so. We first encapsule the functionality in some functions and\n", + "objects." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### The backbone information\n", + "We create a lookup table of telescopes (i.e. a dictionary), in which\n", + "we encapsule their names, geographical position, height above the\n", + "ground, and their operating frequencies. We also create a dictionary\n", + "of landmarks, for orientation." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "Telescope = namedtuple('telescope',\n", + " ['name', 'coordinates', 'height', 'frequencies'])\n", + "telinfo = {\n", + " 'Wb': Telescope('RT Westerbork', (6.60334 * u.deg, 52.91474 * u.deg),\n", + " 22.5 * u.m, (610. * u.MHz, 6650. * u.MHz)),\n", + " 'Eb': Telescope('RT Effelsberg', (6.88361 * u.deg, 50.52483 * u.deg),\n", + " 50. * u.m, (610. * u.MHz, 6650. * u.MHz))\n", + "}\n", + "\n", + "Site = namedtuple('site', ['name', 'coord', 'pixoff', 'color'])\n", + "sites = OrderedDict([\n", + " # ID, tuple(Name,\n", + " # (lon, lat),\n", + " # (xoff, yoff),\n", + " # color\n", + " # )\n", + " ('RT Effelsberg (Eb)',\n", + " Site('RT Effelsberg (Eb)',\n", + " (6.88361, 50.52483),\n", + " (20, +30),\n", + " 'k')\n", + " ),\n", + " ('Cologne',\n", + " Site('Cologne',\n", + " (6.9602, 50.9377),\n", + " (20, +30),\n", + " 'k')\n", + " ),\n", + " ('Bonn',\n", + " Site('Bonn',\n", + " (7.0982, 50.7375),\n", + " (20, +30),\n", + " 'k')\n", + " ),\n", + " ('RT Westerbork (Wb)',\n", + " Site('RT Westerbork (Wb)',\n", + " (6.60334, 52.91474),\n", + " (20, +30),\n", + " 'k')\n", + " ),\n", + " ('Assen',\n", + " Site('Assen',\n", + " (6.5625, 52.99667),\n", + " (20, +30),\n", + " 'k')\n", + " ),\n", + " ('Groningen',\n", + " Site('Groningen',\n", + " (6.56667, 53.21917),\n", + " (20, +30),\n", + " 'k')\n", + " ),\n", + " ('Istanbul',\n", + " Site('Istanbul',\n", + " (28.955, 41.013611),\n", + " (20, +30),\n", + " 'k')\n", + " ),\n", + " ('London',\n", + " Site('London',\n", + " (0.1275, 51.507222),\n", + " (20, +30),\n", + " 'k')\n", + " ),\n", + " ('Saint Petersburg',\n", + " Site('Saint Petersburg',\n", + " (30.3, 59.95),\n", + " (20, +30),\n", + " 'k')\n", + " ),\n", + " ('Berlin',\n", + " Site('Berlin',\n", + " (13.383333, 52.516667),\n", + " (20, +30),\n", + " 'k')\n", + " ),\n", + " ('Madrid',\n", + " Site('Madrid',\n", + " (3.716667, 40.383333),\n", + " (20, +30),\n", + " 'k')\n", + " ),\n", + " ('Kiev',\n", + " Site('Kiev',\n", + " (30.523333, 50.45),\n", + " (20, +30),\n", + " 'k')\n", + " ),\n", + " ('Rome',\n", + " Site('Rome',\n", + " (12.5, 41.9),\n", + " (20, +30),\n", + " 'k')\n", + " ),\n", + " ('Paris',\n", + " Site('Paris',\n", + " (2.3508, 48.8567),\n", + " (20, +30),\n", + " 'k')\n", + " ),\n", + " ('Bucharest',\n", + " Site('Bucharest',\n", + " (26.103889, 44.4325),\n", + " (20, +30),\n", + " 'k')\n", + " ),\n", + " ('Minsk',\n", + " Site('Minsk',\n", + " (27.566667, 53.9),\n", + " (20, +30),\n", + " 'k')\n", + " ),\n", + " ('Vienna',\n", + " Site('Vienna',\n", + " (16.366667, 48.2),\n", + " (20, +30),\n", + " 'k')\n", + " ),\n", + " ('Hamburg',\n", + " Site('Hamburg',\n", + " (10.001389, 53.565278),\n", + " (20, +30),\n", + " 'k')\n", + " ),\n", + " ('Budapest',\n", + " Site('Budapest',\n", + " (19.051389, 47.4925),\n", + " (20, +30),\n", + " 'k')\n", + " ),\n", + " ('Warsaw',\n", + " Site('Warsaw',\n", + " (21.016667, 52.233333),\n", + " (20, +30),\n", + " 'k')\n", + " ),\n", + " ('Barcelona',\n", + " Site('Barcelona',\n", + " (2.183333, 41.383333),\n", + " (20, +30),\n", + " 'k')\n", + " ),\n", + " ('Munich',\n", + " Site('Munich',\n", + " (11.566667, 48.133333),\n", + " (20, +30),\n", + " 'k')\n", + " ),\n", + " ('Kharkiv',\n", + " Site('Kharkiv',\n", + " (36.231389, 50.004444),\n", + " (20, +30),\n", + " 'k')\n", + " ),\n", + " ('Milan',\n", + " Site('Milan',\n", + " (9.183333, 45.466667),\n", + " (20, +30),\n", + " 'k')\n", + " ),\n", + " ('Prague',\n", + " Site('Prague',\n", + " (14.416667, 50.083333),\n", + " (20, +30),\n", + " 'k')\n", + " ),\n", + " ('Nizhny Novgorod',\n", + " Site('Nizhny Novgorod',\n", + " (44.0075, 56.326944),\n", + " (20, +30),\n", + " 'k')\n", + " ),\n", + " ('Sofia',\n", + " Site('Sofia',\n", + " (23.33, 42.7),\n", + " (20, +30),\n", + " 'k')\n", + " ),\n", + " ('Kazan',\n", + " Site('Kazan',\n", + " (49.134722, 55.790278),\n", + " (20, +30),\n", + " 'k')\n", + " ),\n", + " ('Brussels',\n", + " Site('Brussels',\n", + " (4.35, 50.85),\n", + " (20, +30),\n", + " 'k')\n", + " ),\n", + " ('Samara',\n", + " Site('Samara',\n", + " (50.140833, 53.202778),\n", + " (20, +30),\n", + " 'k')\n", + " ),\n", + " ('Belgrade',\n", + " Site('Belgrade',\n", + " (20.466667, 44.816667),\n", + " (20, +30),\n", + " 'k')\n", + " ),\n", + " ('Rostov-on-Don',\n", + " Site('Rostov-on-Don',\n", + " (39.7, 47.233333),\n", + " (20, +30),\n", + " 'k')\n", + " ),\n", + " ('Birmingham',\n", + " Site('Birmingham',\n", + " (1.893611, 52.483056),\n", + " (20, +30),\n", + " 'k')\n", + " ),\n", + " ('Ufa',\n", + " Site('Ufa',\n", + " (55.966667, 54.75),\n", + " (20, +30),\n", + " 'k')\n", + " ),\n", + " ('Cologne',\n", + " Site('Cologne',\n", + " (6.952778, 50.936389),\n", + " (20, +30),\n", + " 'k')\n", + " ),\n", + " ('Perm',\n", + " Site('Perm',\n", + " (56.316667, 58),\n", + " (20, +30),\n", + " 'k')\n", + " ),\n", + " ('Voronezh',\n", + " Site('Voronezh',\n", + " (39.210556, 51.671667),\n", + " (20, +30),\n", + " 'k')\n", + " ),\n", + " ('Volgograd',\n", + " Site('Volgograd',\n", + " (44.516667, 48.7),\n", + " (20, +30),\n", + " 'k')\n", + " ),\n", + " ('Odessa',\n", + " Site('Odessa',\n", + " (30.733333, 46.466667),\n", + " (20, +30),\n", + " 'k')\n", + " ),\n", + " ('Naples',\n", + " Site('Naples',\n", + " (14.25, 40.833333),\n", + " (20, +30),\n", + " 'k')\n", + " ),\n", + " ('Dnipropetrovsk',\n", + " Site('Dnipropetrovsk',\n", + " (34.983333, 48.45),\n", + " (20, +30),\n", + " 'k')\n", + " ),\n", + " ('Mariehamn',\n", + " Site('Mariehamn',\n", + " (19.9, 60.116667),\n", + " (20, +30),\n", + " 'k')\n", + " ),\n", + " ('Tirana',\n", + " Site('Tirana',\n", + " (19.816667, 41.31666667),\n", + " (20, +30),\n", + " 'k')\n", + " ),\n", + " ('Andorra la Vella',\n", + " Site('Andorra la Vella',\n", + " (1.516667, 42.5),\n", + " (20, +30),\n", + " 'k')\n", + " ),\n", + " ('Yerevan',\n", + " Site('Yerevan',\n", + " (44.5, 40.16666667),\n", + " (20, +30),\n", + " 'k')\n", + " ),\n", + " ('Baku',\n", + " Site('Baku',\n", + " (49.866667, 40.38333333),\n", + " (20, +30),\n", + " 'k')\n", + " ),\n", + " ('Sarajevo',\n", + " Site('Sarajevo',\n", + " (18.416667, 43.86666667),\n", + " (20, +30),\n", + " 'k')\n", + " ),\n", + " ('Zagreb',\n", + " Site('Zagreb',\n", + " (16, 45.8),\n", + " (20, +30),\n", + " 'k')\n", + " ),\n", + " ('Nicosia',\n", + " Site('Nicosia',\n", + " (33.366667, 35.16666667),\n", + " (20, +30),\n", + " 'k')\n", + " ),\n", + " ('Copenhagen',\n", + " Site('Copenhagen',\n", + " (12.583333, 55.66666667),\n", + " (20, +30),\n", + " 'k')\n", + " ),\n", + " ('Tallinn',\n", + " Site('Tallinn',\n", + " (24.716667, 59.43333333),\n", + " (20, +30),\n", + " 'k')\n", + " ),\n", + " ('Torshavn',\n", + " Site('Torshavn',\n", + " (-6.766667, 62),\n", + " (20, +30),\n", + " 'k')\n", + " ),\n", + " ('Helsinki',\n", + " Site('Helsinki',\n", + " (24.933333, 60.16666667),\n", + " (20, +30),\n", + " 'k')\n", + " ),\n", + " ('Tbilisi',\n", + " Site('Tbilisi',\n", + " (44.833333, 41.68333333),\n", + " (20, +30),\n", + " 'k')\n", + " ),\n", + " ('Gibraltar',\n", + " Site('Gibraltar',\n", + " (-5.35, 36.13333333),\n", + " (20, +30),\n", + " 'k')\n", + " ),\n", + " ('Athens',\n", + " Site('Athens',\n", + " (23.733333, 37.98333333),\n", + " (20, +30),\n", + " 'k')\n", + " ),\n", + " ('Saint Peter Port',\n", + " Site('Saint Peter Port',\n", + " (-2.533333, 49.45),\n", + " (20, +30),\n", + " 'k')\n", + " ),\n", + " ('Reykjavik',\n", + " Site('Reykjavik',\n", + " (-21.95, 64.15),\n", + " (20, +30),\n", + " 'k')\n", + " ),\n", + " ('Dublin',\n", + " Site('Dublin',\n", + " (-6.233333, 53.31666667),\n", + " (20, +30),\n", + " 'k')\n", + " ),\n", + " ('Douglas',\n", + " Site('Douglas',\n", + " (-4.483333, 54.15),\n", + " (20, +30),\n", + " 'k')\n", + " ),\n", + " ('Saint Helier',\n", + " Site('Saint Helier',\n", + " (-2.1, 49.18333333),\n", + " (20, +30),\n", + " 'k')\n", + " ),\n", + " ('Pristina',\n", + " Site('Pristina',\n", + " (21.166667, 42.66666667),\n", + " (20, +30),\n", + " 'k')\n", + " ),\n", + " ('Riga',\n", + " Site('Riga',\n", + " (24.1, 56.95),\n", + " (20, +30),\n", + " 'k')\n", + " ),\n", + " ('Vaduz',\n", + " Site('Vaduz',\n", + " (9.516667, 47.13333333),\n", + " (20, +30),\n", + " 'k')\n", + " ),\n", + " ('Vilnius',\n", + " Site('Vilnius',\n", + " (25.316667, 54.68333333),\n", + " (20, +30),\n", + " 'k')\n", + " ),\n", + " ('Luxembourg',\n", + " Site('Luxembourg',\n", + " (6.116667, 49.6),\n", + " (20, +30),\n", + " 'k')\n", + " ),\n", + " ('Skopje',\n", + " Site('Skopje',\n", + " (21.433333, 42),\n", + " (20, +30),\n", + " 'k')\n", + " ),\n", + " ('Valletta',\n", + " Site('Valletta',\n", + " (14.5, 35.88333333),\n", + " (20, +30),\n", + " 'k')\n", + " ),\n", + " ('Chisinau',\n", + " Site('Chisinau',\n", + " (28.85, 47),\n", + " (20, +30),\n", + " 'k')\n", + " ),\n", + " ('Monaco',\n", + " Site('Monaco',\n", + " (7.416667, 43.73333333),\n", + " (20, +30),\n", + " 'k')\n", + " ),\n", + " ('Podgorica',\n", + " Site('Podgorica',\n", + " (19.266667, 42.43333333),\n", + " (20, +30),\n", + " 'k')\n", + " ),\n", + " ('Amsterdam',\n", + " Site('Amsterdam',\n", + " (4.916667, 52.35),\n", + " (20, +30),\n", + " 'k')\n", + " ),\n", + " ('North Nicosia',\n", + " Site('North Nicosia',\n", + " (33.366667, 35.183333),\n", + " (20, +30),\n", + " 'k')\n", + " ),\n", + " ('Oslo',\n", + " Site('Oslo',\n", + " (10.75, 59.91666667),\n", + " (20, +30),\n", + " 'k')\n", + " ),\n", + " ('Lisbon',\n", + " Site('Lisbon',\n", + " (-9.133333, 38.71666667),\n", + " (20, +30),\n", + " 'k')\n", + " ),\n", + " ('Moscow',\n", + " Site('Moscow',\n", + " (37.6, 55.75),\n", + " (20, +30),\n", + " 'k')\n", + " ),\n", + " ('San Marino',\n", + " Site('San Marino',\n", + " (12.416667, 43.93333333),\n", + " (20, +30),\n", + " 'k')\n", + " ),\n", + " ('Bratislava',\n", + " Site('Bratislava',\n", + " (17.116667, 48.15),\n", + " (20, +30),\n", + " 'k')\n", + " ),\n", + " ('Ljubljana',\n", + " Site('Ljubljana',\n", + " (14.516667, 46.05),\n", + " (20, +30),\n", + " 'k')\n", + " ),\n", + " ('Longyearbyen',\n", + " Site('Longyearbyen',\n", + " (15.633333, 78.21666667),\n", + " (20, +30),\n", + " 'k')\n", + " ),\n", + " ('Stockholm',\n", + " Site('Stockholm',\n", + " (18.05, 59.33333333),\n", + " (20, +30),\n", + " 'k')\n", + " ),\n", + " ('Bern',\n", + " Site('Bern',\n", + " (7.466667, 46.91666667),\n", + " (20, +30),\n", + " 'k')\n", + " ),\n", + " ('Ankara',\n", + " Site('Ankara',\n", + " (32.866667, 39.93333333),\n", + " (20, +30),\n", + " 'k')\n", + " ),\n", + " ('Kyiv',\n", + " Site('Kyiv',\n", + " (30.516667, 50.43333333),\n", + " (20, +30),\n", + " 'k')\n", + " ),\n", + " ('London',\n", + " Site('London',\n", + " (-0.083333, 51.5),\n", + " (20, +30),\n", + " 'k')\n", + " ),\n", + " ('Vatican City',\n", + " Site('Vatican City',\n", + " (12.45, 41.9),\n", + " (20, +30),\n", + " 'k')\n", + " )\n", + " ])\n" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "jupyter": { + "outputs_hidden": false + } + }, + "source": [ + "### Topographical information\n", + "We create a function that can generate a topographical map and a\n", + "hardcopy" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false, + "jupyter": { + "outputs_hidden": false + } + }, + "outputs": [], + "source": [ + "def plotheightmap(heightmap, sites, lons, lats, outname=None,\n", + " overwrite=True, title=''):\n", + " \"\"\"Plot map showing height above geoid (topographic map)\n", + "\n", + " Takes a numpy array heightmap, including an array of longitudes\n", + " and latitudes and generates a topographic map. An additional input\n", + " is an ordered dict of named tuples:\n", + "\n", + " ID: Site(Name,\n", + " (lon, lat),\n", + " (xoff, yoff),\n", + " color\n", + " )\n", + "\n", + " indicating landmarks (cities), which are pointed at by arrows and\n", + " whose names are displayed on the plot. Name is the name of the\n", + " site, (lon, lat), the longitude and latitude pair, (xoff, yoff)\n", + " the offset of the name on the map, and color the (matplotlib)\n", + " color of the name.\n", + "\n", + " The map can be exported to a file outname, overwriting any\n", + " existing file if overwrite is True. Optionally, the map gets a\n", + " title specified with the keyword title, placed on the top of the\n", + " map.\n", + "\n", + " Parameters:\n", + " heightmap (ndarray): Height map im m above geoid,\n", + " as provided by function\n", + " pycraf.pathprof.srtm_height_map\n", + " sites (OrderedDict): Dictionary with site information\n", + " lons (ndarray): Longitudes,\n", + " as provided by function\n", + " pycraf.pathprof.srtm_height_map\n", + " lats (ndarray): Latitudes,\n", + " as provided by function\n", + " pycraf.pathprof.srtm_height_map\n", + " outname (str or None): Name of output file\n", + " overwrite (bool): Overwrite existing output file? (True:yes)\n", + " title (str): Optional title of map\n", + "\n", + " Returns:\n", + " None\n", + " \"\"\"\n", + " _lons = lons.to(u.deg).value\n", + " _lats = lats.to(u.deg).value\n", + " _heightmap = heightmap.to(u.m).value\n", + "\n", + " vmin, vmax = -20, 950\n", + " terrain_cmap, terrain_norm = pathprof.terrain_cmap_factory(\n", + " sealevel=0.5, vmax=vmax)\n", + " _heightmap[_heightmap < 0] = 0.51 # fix for coastal region\n", + "\n", + " # plt.close()\n", + "\n", + " fig = plt.figure(figsize=(10, 10))\n", + " ax = fig.add_axes((0.1, 0.12, 0.88, 0.83))\n", + " cbax = fig.add_axes((0.15, 0.145, 0.78, .02))\n", + " # ax = fig.add_axes((0., 0., 1.0, 1.0))\n", + " # cbax = fig.add_axes((0., 0., 1.0, .02))\n", + "\n", + " cim = ax.imshow(\n", + " _heightmap,\n", + " origin='lower', interpolation='nearest',\n", + " cmap=terrain_cmap, norm=terrain_norm,\n", + " # vmin=vmin, vmax=vmax, # deprecated in newer matplotlib versions\n", + " extent=(_lons[0], _lons[-1], _lats[0], _lats[-1]),\n", + " )\n", + " cbar = fig.colorbar(\n", + " cim, cax=cbax, orientation='horizontal'\n", + " )\n", + " aspect = abs(_lons[-1] - _lons[0]) / abs(_lats[-1] - _lats[0])\n", + " ax.set_aspect(aspect)\n", + " cbar.set_label(r'Height (amsl)', color='k')\n", + " for t in cbax.xaxis.get_major_ticks():\n", + " t.tick1line.set_visible(True)\n", + " t.tick2line.set_visible(True)\n", + " t.label1.set_visible(False)\n", + " t.label2.set_visible(True)\n", + " ctics = np.arange(0, 1350, 200)\n", + " cbar.set_ticks(ctics)\n", + " ctics = cbar.get_ticks()\n", + " cbar.ax.set_xticklabels(map('{:.0f} m'.format, ctics), color='k')\n", + " ax.set_xlabel('Longitude [deg]')\n", + " ax.set_ylabel('Latitude [deg]')\n", + "\n", + " ax.xaxis.tick_top()\n", + " ax.xaxis.set_label_position('top')\n", + " ax.set_title(title, y=0.95, color='k', size=16)\n", + "\n", + " for site_id, site in sites.items():\n", + " color = site.color\n", + " # color = 'k'\n", + " ax.annotate(\n", + " site.name, xy=site.coord, xytext=site.pixoff,\n", + " textcoords='offset points', color=color,\n", + " arrowprops=dict(arrowstyle=\"->\", color=color)\n", + " )\n", + "\n", + " ax.xaxis.tick_top()\n", + " ax.xaxis.set_label_position('top')\n", + " # plt.show()\n", + " if isinstance(outname, type('')):\n", + " if overwrite:\n", + " if os.path.isfile(outname):\n", + " os.remove(outname)\n", + " plt.savefig(outname)\n", + "\n", + " return" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Plotting an attenuation map\n", + "We saw how to generate and plot the attenuation map in notebook 3c.\n", + "Here, we create a plotting funtion to do so." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false, + "jupyter": { + "outputs_hidden": false + } + }, + "outputs": [], + "source": [ + "def plot_attenmap(total_atten, sites, lons, lats, outname=None,\n", + " overwrite=True, vmin=-5, vmax=175, title='',\n", + " ctics=20):\n", + " \"\"\"Plot map showing attenuation map\n", + "\n", + " Takes a pycraf array total_atten (attenuation in dB), an array of\n", + " longitudes and an array of latitudes and generates a topographic\n", + " map. An additional input is an ordered dict of named tuples:\n", + "\n", + " ID: Site(Name,\n", + " (lon, lat),\n", + " (xoff, yoff),\n", + " color\n", + " )\n", + "\n", + " indicating landmarks (cities), which are pointed at by arrows and\n", + " whose names are displayed on the plot. Name is the name of the\n", + " site, (lon, lat), the longitude and latitude pair, (xoff, yoff)\n", + " the offset of the name on the map, and color the (matplotlib)\n", + " color of the name.\n", + "\n", + " The map can be exported to a file outname, overwriting any\n", + " existing file if overwrite is True. vmin and vmax indicate maximum\n", + " and minimum of the color scheme. Optionally, the map gets a title\n", + " specified with the keyword title, placed on the top of the\n", + " map. The parameter ctics determines the difference (in dB) between\n", + " tics on the color bar.\n", + "\n", + " Parameters:\n", + " heightmap (ndarray): Height map im m above geoid,\n", + " as provided by function\n", + " pycraf.pathprof.srtm_height_map\n", + " sites (OrderedDict): Dictionary with site information\n", + " lons (ndarray): Longitudes,\n", + " as provided by function\n", + " pycraf.pathprof.srtm_height_map\n", + " lats (ndarray): Latitudes,\n", + " as provided by function\n", + " pycraf.pathprof.srtm_height_map\n", + " outname (str or None): Name of output file\n", + " overwrite (bool): Overwrite existing output file (True:yes)?\n", + " vmin (float): Minimum of colour scheme\n", + " vmax (float): Maximum of colour scheme\n", + " title (str): Optional title of map\n", + " ctics (float): Difference between tickmarks on colour bar\n", + "\n", + " Returns:\n", + " None\n", + " \"\"\"\n", + " _lons = lons.to(u.deg).value\n", + " _lats = lats.to(u.deg).value\n", + "\n", + " # plt.close()\n", + "\n", + " fig = plt.figure(figsize=(10, 10))\n", + " ax = fig.add_axes((0.1, 0.08, 0.88, 0.87))\n", + " cbax = fig.add_axes((0.15, 0.1, 0.78, .02))\n", + " cim = ax.imshow(\n", + " total_atten.to(cnv.dB).value,\n", + " origin='lower', interpolation='nearest', cmap='inferno_r',\n", + " vmin=vmin, vmax=vmax,\n", + " extent=(_lons[0], _lons[-1], _lats[0], _lats[-1]),\n", + " )\n", + " cbar = fig.colorbar(\n", + " cim, cax=cbax, orientation='horizontal'\n", + " )\n", + " ax.set_aspect(abs(_lons[-1] - _lons[0]) / abs(_lats[-1] - _lats[0]))\n", + " cbar.set_label(r'Path propagation loss', color='w')\n", + " cbax.xaxis.set_label_position('top')\n", + " for t in cbax.xaxis.get_major_ticks():\n", + " t.tick1line.set_visible(True)\n", + " t.tick2line.set_visible(True)\n", + " t.tick2line.set_color('w')\n", + " t.label1.set_visible(False)\n", + " t.label2.set_visible(True)\n", + "\n", + " ctics = np.arange(0, 300, ctics)\n", + " cbar.set_ticks(ctics)\n", + " ctics = cbar.get_ticks()\n", + " cbar.ax.set_xticklabels(map('{:.0f} dB'.format, ctics), color='w')\n", + " ax.set_xlabel('Longitude [deg]')\n", + " ax.set_ylabel('Latitude [deg]')\n", + "\n", + " for site_id, site in sites.items():\n", + " color = site.color\n", + " color = 'g'\n", + " ax.annotate(\n", + " site.name, xy=site.coord, xytext=site.pixoff,\n", + " textcoords='offset points', color=color,\n", + " arrowprops=dict(arrowstyle=\"->\", color=color)\n", + " )\n", + "\n", + " ax.xaxis.tick_top()\n", + " ax.xaxis.set_label_position('top')\n", + " ax.set_title(title, y=0.95, color='w', size=16)\n", + " if isinstance(outname, type('')):\n", + " if overwrite:\n", + " if os.path.isfile(outname):\n", + " os.remove(outname)\n", + " plt.savefig(outname)" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "tags": [] + }, + "source": [ + "### Statistics along a circle\n", + "Given the attenuation map, we are now able to create statistics about\n", + "the attenuation. Often, we do not need to predict the attenuation at a\n", + "specific geographical point, but we want to know, what the statistics\n", + "of the attenuation is at a certain distance from the receiving\n", + "antenna, minimum, maximum, percentiles, median. Using pycraf structs,\n", + "we generate this type of statistics with the following function." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "def averatt(atten_map, hprof_cache, lo_perc=5., hi_perc=95., rmax=10*u.km,\n", + " bins=100):\n", + " \"\"\"Return circular statistics of attenuation\n", + "\n", + " Given an attenuation map atten_map in addition with either a map\n", + " providing the physical distances between the reference point (at\n", + " which attenuation = 0) and the single pixels or the output of a\n", + " run of the pycraf.pathprof.height_map_data method, the function\n", + " creates statistics along concentric rings (with a constant geoid\n", + " distance to the centre), which are returned as a list of numpy\n", + " arrays. First, the distance between 0 and rmax is binned into bins\n", + " equally-sized distance intervals. For the pixels in each distance\n", + " interval, the following statistics is calculated for the\n", + " attenuation map:\n", + "\n", + " - median\n", + " - percentile with percentage given by lo_perc\n", + " - percentile with percentage given by hi_perc\n", + " - maximum\n", + " - minimum\n", + "\n", + " In addition the average distance for each bin is calculated. Each\n", + " statistics is converted to a numpy array. The function returns a\n", + " dictionary with the following arrays for the given keys:\n", + "\n", + " 'radius': average radii\n", + " 'median': median at each interval\n", + " 'perc_lo': percentile for lo_perc percentage at each interval\n", + " 'perc_hi': percentile for hi_perc percentage at each interval\n", + " 'minimum': minimum for each interval\n", + " 'maximum': maximum for each interval\n", + "\n", + " To e.g. study the median of the attenuation with radius, one can\n", + " use the 'radius' array and the 'median' array etc.\n", + "\n", + " Arguments:\n", + " atten_map (pycraf array): Attenuation map\n", + " hprof_cache (numpy array or pycraf height profile object): map\n", + " with physical distances between reference point and pixel\n", + " lo_perc (float): Lower percentage for the calculation of lower\n", + " percentile\n", + " hi_perc (float): Higher percentage for the calculation of upper\n", + " percentile\n", + " rmax (float, int, or astropy length): maximum radius (in km for\n", + " int or float)\n", + " bins (int): number of bins\n", + "\n", + " Return:\n", + " dictionary with statistics vs. radius.\n", + " \"\"\"\n", + "\n", + " # Distance map in km\n", + " if (isinstance(hprof_cache, type(np.zeros((1, 1))))):\n", + " dis = hprof_cache\n", + " else:\n", + " dis = hprof_cache['dist_map']\n", + "\n", + " # Flatten maps\n", + " # allatt = _total_atten.reshape(dis.shape[0]*dis.shape[1]).to(cnv.dimless)\n", + " # allatt = np.power(10., _total_atten.reshape(\n", + " # dis.shape[0]*dis.shape[1]).value/10.)\n", + " allatt = atten_map.reshape(dis.shape[0]*dis.shape[1]).value\n", + " alldis = dis.reshape(dis.shape[0]*dis.shape[1])\n", + "\n", + " # Create bins\n", + " if isinstance(rmax, type(float)) or isinstance(rmax, type(int)):\n", + " rmaxi = rmax\n", + " else:\n", + " rmaxi = rmax.to(u.km).value\n", + " binsar = np.linspace(0, rmaxi, bins+1)\n", + "\n", + " # Get stats\n", + " digitized = np.digitize(alldis, binsar)\n", + "\n", + " averatt = {}\n", + " averatt['bin_dists'] = \\\n", + " np.array([alldis[digitized == i].mean()\n", + " for i in range(1, len(binsar))])\n", + " averatt['bin_lo'] = \\\n", + " np.array([np.percentile(allatt[digitized == i], lo_perc)\n", + " for i in range(1, len(binsar))])\n", + " averatt['bin_med'] = \\\n", + " np.array([np.percentile(allatt[digitized == i], 50)\n", + " for i in range(1, len(binsar))])\n", + " averatt['bin_hi'] = \\\n", + " np.array([np.percentile(allatt[digitized == i], hi_perc)\n", + " for i in range(1, len(binsar))])\n", + " averatt['bin_min'] = \\\n", + " np.array([allatt[digitized == i].min()\n", + " for i in range(1, len(binsar))])\n", + " averatt['bin_max'] = \\\n", + " np.array([allatt[digitized == i].max()\n", + " for i in range(1, len(binsar))])\n", + "\n", + " return averatt" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "tags": [] + }, + "source": [ + "### Plotting statistical profiles\n", + "Now that we generated the statistics of the attenuation as a function\n", + "of radius but averaged over equidistant rings, we can plot the\n", + "results. The following function is able to generate an average\n", + "attenuation profile, and, using an inlay, another plot showing the\n", + "difference to an approximation of a flat earth." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "def plotaveratt(averatts, title=None, outname=None,\n", + " overwrite=True, ymin=100, ymax=200, averatts2=None,\n", + " inlay=(0.3, 0.13, 0.65, 0.25), avmd=2, grid=10):\n", + " \"\"\"Plot circularly averaged attenuations with radius\n", + "\n", + " Given the output of function averatt (attenuation statistics along\n", + " concentric rings), plots median as a black line with green area\n", + " for upper and lower percentile, and cyan area between min and max\n", + " for each radius. If outname is not None, a file with the viewgraph\n", + " will be generated. If overwrite is True any existing file with the\n", + " same name will be overwritten. ymin and ymax are the minimimum and\n", + " the maximum of the viewgraph. If averatts2 is not None, produces\n", + " red line showing the median in averatts2, and an inlay with\n", + " dimensions specified by the inlay parameter, a list or tuple with\n", + " the position of the lower left corner given by the first two\n", + " members, and width and height relative to the plot widths and\n", + " height. The inlay shows above quantities of averatts minus the\n", + " median in averatts2. Instead of the average, inlay contains an\n", + " vertical line indicating the average difference between the\n", + " medians in averatts and averatts2, starting at radius avmd (in\n", + " km). Grid denotes the vertical grid line separation in the main\n", + " viewgraph.\n", + "\n", + " Arguments:\n", + " averatts (dict): Output of averatt function\n", + " title (str or None): Title of the plot\n", + " outname (str or None): Name of output file\n", + " overwrite (bool): Overwrite existing output? (True: yes)\n", + " ymin (float): Minimum on ordinate\n", + " ymax (float): Maximum on ordinate\n", + " averatts (dict): Output of averatt function (second dict)\n", + " avmd (float): Minimum radius for averaging inlay curve\n", + " inlay (tuple): Dimensions of inlay\n", + " grid (float): Grid line separation in the main viewgraph\n", + "\n", + " Returns:\n", + " None\n", + " \"\"\"\n", + "\n", + " # bin_means_dB = (cnv.dimless*bin_means).to(cnv.dB).value\n", + " # bin_upper_dB = (cnv.dimless*(bin_means+bin_stdev)).to(cnv.dB).value\n", + " # bin_lower_dB = (cnv.dimless*(bin_means-bin_stdev)).to(cnv.dB).value\n", + " bin_dists = averatts['bin_dists']\n", + " bin_med_dB = averatts['bin_med']\n", + " bin_hi_dB = averatts['bin_hi']\n", + " bin_lo_dB = averatts['bin_lo']\n", + " bin_max_dB = averatts['bin_max']\n", + " bin_min_dB = averatts['bin_min']\n", + "\n", + " # plt.close()\n", + " fig = plt.figure(figsize=(10, 5))\n", + " ax = fig.add_axes((0.1, 0.08, 0.88, 0.87))\n", + " ax.plot(bin_dists, bin_med_dB, 'k-')\n", + " ax.fill_between(bin_dists, bin_lo_dB, bin_hi_dB,\n", + " fc='g', ec='r', linewidth=0, linestyle='None')\n", + " ax.fill_between(bin_dists, bin_lo_dB, bin_min_dB,\n", + " fc='cyan', linewidth=0, linestyle='None', alpha=0.5)\n", + " ax.fill_between(bin_dists, bin_hi_dB, bin_max_dB,\n", + " fc='cyan', linewidth=0, linestyle='None', alpha=0.5)\n", + " ax.set_xlabel('Distance [km]', size=24)\n", + " ax.set_ylabel('Average attenuation [dB]', size=24)\n", + " ax.xaxis.set_tick_params(labelsize=16)\n", + " ax.yaxis.set_tick_params(labelsize=16)\n", + " ax.set_ylim(ymin, ymax)\n", + " if grid:\n", + " ax.yaxis.grid()\n", + "\n", + " if isinstance(title, type('')):\n", + " ax.set_title(title, size=24)\n", + "\n", + " if not isinstance(averatts2, type(None)):\n", + " bin_med_dB = averatts2['bin_med']\n", + " ax.plot(bin_dists, bin_med_dB, 'r-')\n", + " if not isinstance(inlay, type(None)):\n", + " ax = fig.add_axes(inlay)\n", + " bin_diff_dB = averatts['bin_med']-bin_med_dB\n", + " ad = np.average(bin_diff_dB[bin_dists > avmd])\n", + " bin_adiff_dB = bin_diff_dB[bin_dists > avmd]*0+ad\n", + " bin_hi_dBd = averatts['bin_hi']-bin_med_dB\n", + " bin_lo_dBd = averatts['bin_lo']-bin_med_dB\n", + " bin_max_dBd = averatts['bin_max']-bin_med_dB\n", + " bin_min_dBd = averatts['bin_min']-bin_med_dB\n", + " ax.plot(bin_dists, bin_diff_dB, 'k-')\n", + " ax.plot(bin_dists[bin_dists > avmd], bin_adiff_dB, 'r-')\n", + " ax.fill_between(bin_dists, bin_lo_dBd, bin_hi_dBd, fc='g',\n", + " ec='r', linewidth=0, linestyle='None')\n", + " ax.fill_between(bin_dists, bin_lo_dBd, bin_min_dBd, fc='cyan',\n", + " linewidth=0, linestyle='None', alpha=0.5)\n", + " ax.fill_between(bin_dists, bin_hi_dBd, bin_max_dBd, fc='cyan',\n", + " linewidth=0, linestyle='None', alpha=0.5)\n", + " ax.set_ylabel('Difference [dB]', size=10)\n", + " ax.xaxis.set_tick_params(labelsize=10)\n", + " ax.yaxis.set_tick_params(labelsize=10)\n", + " shift = (ax.get_ylim()[1]-ax.get_ylim()[0])*0.05\n", + " ax.annotate('{:.0f} dB'.format(ad),\n", + " xy=(bin_dists[-1]*0.9, ad+shift), color='r')\n", + " if grid:\n", + " ax.yaxis.grid()\n", + "\n", + " if isinstance(outname, type('')):\n", + " if overwrite:\n", + " if os.path.isfile(outname):\n", + " os.remove(outname)\n", + " plt.savefig(outname)\n", + " return" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Generating multiple maps and profiles\n", + "As laid out above, we want to create attenuation maps and average\n", + "profiles for several telescopes and multiple frequencies. The\n", + "following function generates those maps for a single telescope, but\n", + "multiple frequencies. It will also not only generate a map/profile\n", + "assuming normal geography, but also assuming a flat earth." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "def getattdata(telname, tellist, params, genparams, transmitters):\n", + " \"\"\"Batch production of topographic and attenuation maps and statistics\n", + "\n", + " For a telescope with telname the function constructs a\n", + " topographical map, a set of attenuation maps, and radial\n", + " statistics based on the attenuation maps. tellist is a\n", + " dictionary, of which one of the keys must match telname. The\n", + " values are a named tuple containing the telescope name, a tuple\n", + " with its coordinates (longitude and latitude in astropy units),\n", + " its height above the ground in astropy units, and an iterable\n", + " containg the frequencies at which the telescope operates in\n", + " astropy units. The function identifies the telescope in question\n", + " from tellist. params is a dictionary of parameters for the\n", + " calculation of the attenuation map, see description of functions\n", + " pycraf.pathprof.srtm_height_map, pycraf.pathprof.height_map_data,\n", + " and pycraf.pathprof.atten_map_fast:\n", + "\n", + " key description\n", + " map_size_lon: Size of map in geographical longitude in astropy units\n", + " map_size_lat: Size of map in geographical latitude in astropy units\n", + " omega: Fraction of path over sea (astropy units)\n", + " temperature: Temperature in astropy units\n", + " pressure: Air pressure in astropy units\n", + " timepercent: Time percentage in astropy units, see P.452 for explanation\n", + " zone_t: Clutter type for transmitter terminal\n", + " zone_r: Clutter type for receiver terminal\n", + " map_resolution: Resolution of map in astropy units\n", + "\n", + " genparms is a dictionary with parameters specific to the\n", + " generation of radial statistics, see function averatt in this\n", + " module:\n", + " key description\n", + " lo_perc: Lower percentage for the calculation of lower percentile\n", + " hi_perc: Higher percentage for the calculation of upper percentile\n", + " rmax: Maximum radius (in km for int or float)\n", + " bins: Number of bins\n", + "\n", + " transmitters is a dictionary of transmitters, in which the keys\n", + " are the transmitted frequencies in astropy units and the value the\n", + " corresponding height above the ground of the transmitter.\n", + "\n", + " For the given telescope the function produces a topographical map,\n", + " then for each frequency an attenuation map (see function\n", + " pycraf.pathprof.atten_map_fast) and a radial statistics object\n", + " (see function averatt in this module), to then create the\n", + " attenuation map and corresponding statistics under the assumption\n", + " of a flat geoid geometry (\"flat earth\"). The output is a\n", + " dictionary:\n", + "\n", + " key description\n", + " telname Telescope name\n", + " frequencies Frequencies as listed for telescope\n", + " lons Longitudes along longitude axis\n", + " lats Latitudes along latitude axis\n", + " heightmap A topographical map (ndarray)\n", + " distmap A map showing the physical distance w.r.t. the central pixel\n", + " attenmaps List of attenuation maps, one per frequency\n", + " attenmaps_0 List of attenuation maps under flat-earth approximation,\n", + " one per frequency\n", + " averatts List of statistics objects as returned by averatt, one per\n", + " frequency\n", + " averatts_0 List of statistics objects as returned by averatt, one per\n", + " frequency, flat earth approximation\n", + "\n", + " Arguments:\n", + " telname (str): Telescope name\n", + " tellist (dict): List of available telescopes\n", + " params (dict): Parameter list for pycraf functions\n", + " genparams (dict): Parameter list for averatt function (this module)\n", + " transmitters(dict): Dictionary with frequencies and heights\n", + "\n", + " Returns:\n", + " A dictionary (see description)\n", + " \"\"\"\n", + " start_time_tot = time.time()\n", + " output = {}\n", + " output['telname'] = telname\n", + " telpars = tellist[telname]\n", + " output['frequencies'] = telpars.frequencies\n", + " print('{}: Creating height map'.format(telname))\n", + " start_time_loc = time.time()\n", + " lons, lats, heightmap = pathprof.srtm_height_map(\n", + " telpars.coordinates[0], telpars.coordinates[1],\n", + " params['map_size_lon'], params['map_size_lat'],\n", + " map_resolution=params['map_resolution']\n", + " )\n", + " print('Time after start: {} minutes'.format((time.time() -\n", + " start_time_tot)/60.))\n", + " print('Time for this task: {} minutes'.format((time.time() -\n", + " start_time_loc)/60.))\n", + " output['lons'] = lons\n", + " output['lats'] = lats\n", + " output['heightmap'] = heightmap\n", + " print('{}: Creating terrain information map'.format(telname))\n", + " start_time_loc = time.time()\n", + " hprof_cache = pathprof.height_map_data(\n", + " telpars.coordinates[0], telpars.coordinates[1],\n", + " params['map_size_lon'], params['map_size_lat'],\n", + " map_resolution=params['map_resolution'],\n", + " zone_t=params['zone_t'], zone_r=params['zone_r'],\n", + " )\n", + " print('Time after start: {} minutes'.format((time.time() -\n", + " start_time_tot)/60.))\n", + " print('Time for this task: {} minutes'.format((time.time() -\n", + " start_time_loc)/60.))\n", + " output['distmap'] = hprof_cache['dist_map']\n", + " attenmaps = []\n", + " attenmaps_0 = []\n", + " averatts = []\n", + " averatts_0 = []\n", + " for frequency in output['frequencies']:\n", + " print('{}: Creating attenuation map at'.format(telname),\n", + " ' {:.0f} MHz'.format(frequency.to(u.MHz).value))\n", + " start_time_loc = time.time()\n", + "\n", + " h_rg = transmitters[frequency.to(u.MHz)]\n", + " attenmaps += [pathprof.atten_map_fast(\n", + " frequency,\n", + " params['temperature'],\n", + " params['pressure'],\n", + " telpars.height,\n", + " h_rg,\n", + " params['timepercent'],\n", + " hprof_cache\n", + " )['L_b']]\n", + " print('Time after start: {} minutes'.format((time.time() -\n", + " start_time_tot)/60.))\n", + " print('Time for this task: {} minutes'.format((time.time() -\n", + " start_time_loc)/60.))\n", + "\n", + " print('{}: Creating attenuation map '.format(telname),\n", + " 'at {:.0f} MHz, flat earth assumption'.format(frequency.to(u.MHz).value))\n", + " start_time_loc = time.time()\n", + " height_profs_resc = copy.deepcopy(hprof_cache['height_profs'])\n", + " hprof_cache['height_profs'][:] = 0.\n", + " attenmaps_0 += [pathprof.atten_map_fast(\n", + " frequency,\n", + " params['temperature'],\n", + " params['pressure'],\n", + " telpars.height,\n", + " h_rg,\n", + " params['timepercent'],\n", + " hprof_cache\n", + " )['L_b']]\n", + " hprof_cache['height_profs'] = height_profs_resc\n", + " print('Time after start: {} minutes'.format((time.time() -\n", + " start_time_tot)/60.))\n", + " print('Time for this task: {} minutes'.format((time.time() -\n", + " start_time_loc)/60.))\n", + "\n", + " print('{}: Creating stats for attenuation map'.format(telname),\n", + " ' at {:.0f} MHz'.format(frequency.to(u.MHz).value))\n", + " start_time_loc = time.time()\n", + " a = attenmaps[-1]\n", + " averatts += [averatt(attenmaps[-1], hprof_cache, **genparams)]\n", + " print('Time after start: {} minutes'.format((time.time() -\n", + " start_time_tot)/60.))\n", + " print('Time for this task: {} minutes'.format((time.time() -\n", + " start_time_loc)/60.))\n", + "\n", + " print('{}: Creating stats for attenuation map'.format(telname),\n", + " ' at {:.0f} MHz,'.format(frequency.to(u.MHz).value),\n", + " ' flat earth assumption')\n", + " start_time_loc = time.time()\n", + " averatts_0 += [averatt(attenmaps_0[-1], hprof_cache, **genparams)]\n", + " print('Time after start: {} minutes'.format((time.time() -\n", + " start_time_tot)/60.))\n", + " print('Time for this task: {} minutes'.format((time.time() -\n", + " start_time_loc)/60.))\n", + "\n", + " output['attenmaps'] = attenmaps\n", + " output['attenmaps_0'] = attenmaps_0\n", + " output['averatts'] = averatts\n", + " output['averatts_0'] = averatts_0\n", + " return output" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Batch plotting\n", + "The following function will plot the data that we generated in the\n", + "previous step, a topographical map for any frequency, and an\n", + "attenuation map and a statistical profile for a set of frequencies." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "def plotattdata(attdata, tellist, sites, att_lo=2, att_hi=90,\n", + " ctics=10, ave_lo=0.0, ave_hi=100, grid=10,\n", + " inlay=(0.3, 0.13, 0.65, 0.25), avmd=2):\n", + " \"\"\"Generate topographical and attenuation maps and attenuation profiles\n", + "\n", + " Use a dictionary of maps and other data as generated by function\n", + " getattdata (this module) to create a set of maps. The dict attdata\n", + " is described as output in function getattdata, tellist is a list\n", + " of telescope properties as described in function getattdata, sites\n", + " is a list of landmarks as described in function plot_attenmap\n", + " (this module). The function will first plot a topographic map and\n", + " generate a hardcopy with the name xxx_hei.pdf, where xxx is the\n", + " short id of the telescope. For each frequency, the attenuation\n", + " maps will be plotted, and hardcopies xxx_att_yyyy.pdf will\n", + " generated, where xxx is the short id of the telescope, and yyyy\n", + " the frequency in MHz. The maximum and minimum of the colour scheme\n", + " are determined by the percentiles at percentages att_hi (maximum)\n", + " and att_lo (minimum). The colour bar tick mark separation is given\n", + " by the ctics parameter. Finally, circularly averaged attenuation\n", + " profiles are plotted with a hardcopy xxx_rat_yyyy.pdf (xxx and\n", + " yyyy same as above), see function plotaveratt (this module) for\n", + " details. The maximum and minimum of the ordinate are determined by\n", + " the percentiles at percentages ave_hi (maximum) and ave_lo\n", + " (minimum), the grid separation along the ordinate by the parameter\n", + " grid (in dB), the size of the inlay by the parameter inlay (see\n", + " plotaveratt), and the minimum range for averaging the difference\n", + " between conventional profiles and flat-earth profiles is given by\n", + " the parameter avmd.\n", + "\n", + " Parameters:\n", + " attdata (dict): Dictionary of attenuation maps and profiles\n", + " tellist (dict): List of available telescopes\n", + " sites (OrderedDict): Dictionary with site information\n", + " att_lo (float): Percentage of lower percentile to calculate\n", + " attenuation map colour scheme\n", + " att_hi (float): Percentage of higher percentile to calculate\n", + " attenuation map colour scheme\n", + " ctics (float): Difference between tickmarks on colour bar\n", + " ave_lo (float): Percentage of lower percentile to calculate\n", + " attenuation profile ordinate limits\n", + " ave_hi (float): Percentage of higher percentile to calculate\n", + " attenuation profile ordinate limits\n", + " grid (float): Grid line separation in the main viewgraph\n", + " inlay (tuple): Dimensions of inlay on viewgraph\n", + " avmd (float): Minimum radius for averaging inlay curve\n", + "\n", + " Returns:\n", + " void\n", + " \"\"\"\n", + " heightmap = attdata['heightmap']\n", + " lons = attdata['lons']\n", + " lats = attdata['lats']\n", + " telname = attdata['telname']\n", + " telongname = tellist[telname].name\n", + " frequencies = attdata['frequencies']\n", + " attenmaps = attdata['attenmaps']\n", + " attenmaps_0 = attdata['attenmaps_0']\n", + " averatts = attdata['averatts']\n", + " averatts_0 = attdata['averatts_0']\n", + "\n", + " print('{}: Plotting height map'.format(telname))\n", + "\n", + " plotheightmap(heightmap, sites, lons, lats,\n", + " outname='{}_hei.pdf'.format(telname),\n", + " title='Area around {}'.format(telongname))\n", + " for freqnum in range(len(frequencies)):\n", + " attmap = attenmaps[freqnum]\n", + " reshaped = attmap.reshape(attmap.shape[0]*attmap.shape[1]).value\n", + " reshdist = attdata['distmap'].reshape(attmap.shape[0]*attmap.shape[1])\n", + " reshaped = reshaped[reshdist < averatts[freqnum]['bin_dists'].max()]\n", + " vmin = np.percentile(reshaped, att_lo)\n", + " vmax = np.percentile(reshaped, att_hi)\n", + " ymin = np.percentile(reshaped, ave_lo)\n", + " ymax = np.percentile(reshaped, ave_hi)+5.\n", + " freqn = frequencies[freqnum].to(u.MHz).value\n", + " print('{}: Plotting attenuation map'.format(telname),\n", + " 'at {:.0f} MHz'.format(freqn))\n", + " plot_attenmap(attmap, sites, lons, lats,\n", + " outname='{}_att_{:04.0f}.pdf'.format(telname, freqn),\n", + " vmin=vmin, vmax=vmax,\n", + " title='{} {:.0f} MHz'.format(telongname, freqn),\n", + " ctics=ctics)\n", + " print('{}: Plotting average attenuation '.format(telname),\n", + " 'at {:.0f} MHz'.format(freqn))\n", + " plotaveratt(averatts[freqnum], averatts2=averatts_0[freqnum],\n", + " title='{} {:.0f} MHz'.format(telongname, freqn),\n", + " outname='{}_rat_{:04.0f}.pdf'.format(telname, freqn),\n", + " ymin=ymin, ymax=ymax, grid=grid,\n", + " inlay=inlay, avmd=avmd)\n", + " return" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Putting things together, a study of attenuation profiles for different telescopes and frequencies\n", + "We are now ready to generate a set of plots. Our main question is how predicted attenuation profiles differ between frequencies, depending on the location, and depending on the flat earth assumption. We choose two frequencies, for which we consider these questions, $610\\,\\mathrm{MHz}$, and $6650\\,\\mathrm{MHz}$. The chosen sites are the WSRT in The Netherlands in a rather flat terrain, and the Effelsberg radio telescope in a middle mountain range. We create very large maps, such that we can generate attenuation profiles of up to $500\\,\\mathrm{km}$. The time percentage (the fraction of time in which the attenuation is below the value calculated according to [ITU-R Recommendation P.452](https://www.itu.int/rec/R-REC-P.452-16-201507-I/en)) is chosen to be 2%, the statistically allowable loss for the Radio Astronomy Services. We choose the high map resolution of $3^{\\prime\\prime}$. We then generate the profile statistics. For each telescope for each frequency we generate and attenuation map taking geography into account and using the flat earth approximation. The profile viewgraphs show:\n", + "\n", + "Black line: attenuation, 50% percentile\n", + "Green area: attenuation, between 5% perzentile and 95% perzentile\n", + "Cyan area: attenuation, enclosing minimum and maximum\n", + "Red line: flat-earth approximation\n", + "\n", + "The inlay shows the difference of the main quantities and the red line (flat earth) in the same colours, while the red line shows the average of the difference between the median and the flat earth line beyond a distance of $10\\,\\mathrm{km}$. This is done in two steps, the time consuming calculation of the maps (it takes a long time), and the faster generation of the plots." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Characteristics of transmitters, currently height above ground\n", + "transmitters = {\n", + " 610. * u.MHz : 30 * u.m,\n", + " 6650. * u.MHz: 20 * u.m\n", + "}\n", + "\n", + "# Parameters for pycraf attenuation map calculations\n", + "pars = {\n", + " 'map_size_lon': 3.0 * u.deg, \n", + " 'map_size_lat': 3.0 * u.deg,\n", + " 'omega': 0. * u.percent, # fraction of path over sea\n", + " 'temperature': 290. * u.K,\n", + " 'pressure': 1013. * u.hPa,\n", + " 'timepercent': 2 * u.percent, # see P.452 for explanation\n", + " 'zone_t': pathprof.CLUTTER.UNKNOWN,\n", + " 'zone_r': pathprof.CLUTTER.UNKNOWN,\n", + " 'map_resolution': 30. * u.arcsec\n", + "}\n", + "\n", + "# Parameters for circular statistics on attenuation maps\n", + "genpars = {\n", + " 'lo_perc': 5,\n", + " 'hi_perc' : 95,\n", + "# 'rmax': 500 * u.km,\n", + "# 'bins': 500\n", + " 'rmax': 200 * u.km,\n", + " 'bins': 50\n", + "}\n", + "\n", + "telprocess = {\n", + " # att_lo, att_hi, ave_lo, ave_hi, grid, ctics, inlay, avmd\n", + " # input to plotdata\n", + " 'Wb': (2, 90, 0.0, 100, 10, 10, (0.3, 0.13, 0.65, 0.25), 10),\n", + " 'Eb': (2, 90, 0.0, 100, 10, 20, (0.3, 0.13, 0.65, 0.25), 10)\n", + "}\n", + "\n", + "plt.close()\n", + "attdata = {}\n", + "for telescope in telprocess.keys():\n", + " attdata[telescope] = getattdata(telescope, telinfo, pars, genpars, transmitters)\n", + " " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "telprocess = {\n", + " # att_lo, att_hi, ave_lo, ave_hi, grid, ctics, inlay, avmd\n", + " # input to plotdata\n", + " 'Wb': (2, 90, 0.0, 100, 10, 10, (0.3, 0.13, 0.65, 0.25), 80),\n", + " 'Eb': (2, 90, 0.0, 100, 10, 20, (0.3, 0.13, 0.65, 0.25), 80)\n", + "}\n", + "for telescope in telprocess.keys():\n", + " plotattdata(attdata[telescope], telinfo, sites,\n", + " att_lo=telprocess[telescope][0], \n", + " att_hi=telprocess[telescope][1], \n", + " ave_lo=telprocess[telescope][2], \n", + " ave_hi=telprocess[telescope][3], \n", + " grid=telprocess[telescope][4],\n", + " ctics=telprocess[telescope][5],\n", + " inlay=telprocess[telescope][6], \n", + " avmd=telprocess[telescope][7] \n", + " )\n", + " " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Results\n", + "As one can see, the flat earth is a better approximation to the Westerbork telescope than Effelsberg, which is not a surprise. There is nearly no difference between a flat earth assumption and the prediction using the geographical information for the Westerbork telescope. In the far field the difference between Effelsberg using geographical data and the flat earth is about 30 dB, or a factor of 1000." + ] + } + ], + "metadata": { + "anaconda-cloud": {}, + "kernelspec": { + "display_name": "Python [conda env:pycraf3.7_02]", + "language": "python", + "name": "conda-env-pycraf3.7_02-py" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.7.11" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +}