notebooks/beamblockage/wradlib_beamblock.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "nbsphinx": "hidden",
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "This notebook is part of the $\\omega radlib$ documentation: http://wradlib.org/wradlib-docs.\n",
    "\n",
    "Copyright (c) 2018, $\\omega radlib$ developers.\n",
    "Distributed under the MIT License. See LICENSE.txt for more info."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Beam Blockage Calculation using a DEM"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Here, we derive (**p**artial) **b**eam-**b**lockage (**PBB**) from a **D**igital **E**levation **M**odel (**DEM**). "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "We require\n",
    "- the local radar setup (sitecoords, number of rays, number of bins, antenna elevation, beamwidth, and the range resolution);\n",
    "- a **DEM** with a adequate resolution. \n",
    "\n",
    "Here we use pre-processed data from the [GTOPO30](https://lta.cr.usgs.gov/GTOPO30) and [SRTM](http://www2.jpl.nasa.gov/srtm) missions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "import wradlib as wrl\n",
    "import matplotlib.pyplot as pl\n",
    "import matplotlib as mpl\n",
    "import warnings\n",
    "warnings.filterwarnings('ignore')\n",
    "try:\n",
    "    get_ipython().magic(\"matplotlib inline\")\n",
    "except:\n",
    "    pl.ion()\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Setup for Bonn radar"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "First, we need to define some radar specifications (here: *University of Bonn*)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "sitecoords = (7.071663, 50.73052, 99.5)\n",
    "nrays = 360 # number of rays\n",
    "nbins = 1000 # number of range bins\n",
    "el = 1.0 # vertical antenna pointing angle (deg)\n",
    "bw = 1.0 # half power beam width (deg)\n",
    "range_res = 100. # range resolution (meters)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Create the range, azimuth, and beamradius arrays."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "r = np.arange(nbins) * range_res\n",
    "beamradius = wrl.util.half_power_radius(r, bw)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "We use \n",
    "\n",
    "- [wradlib.georef.sweep_centroids](http://wradlib.org/wradlib-docs/latest/generated/wradlib.georef.sweep_centroids.html)  and \n",
    "- [wradlib.georef.spherical_to_proj](http://wradlib.org/wradlib-docs/latest/generated/wradlib.georef.spherical_to_proj.html) \n",
    "\n",
    "to calculate the spherical coordinates of the bin centroids and their longitude, latitude and altitude."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "coord = wrl.georef.sweep_centroids(nrays, range_res, nbins, el)\n",
    "coords = wrl.georef.spherical_to_proj(coord[..., 0], \n",
    "                                      np.degrees(coord[..., 1]),\n",
    "                                      coord[..., 2], sitecoords)\n",
    "lon = coords[..., 0]\n",
    "lat = coords[..., 1]\n",
    "alt = coords[..., 2]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "polcoords = coords[..., :2]\n",
    "print(\"lon,lat,alt:\", coords.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "rlimits = (lon.min(), lat.min(), lon.max(), lat.max())\n",
    "print(\"Radar bounding box:\\n\\t%.2f\\n%.2f             %.2f\\n\\t%.2f\" % \n",
    "      (lat.max(), lon.min(), lon.max(), lat.min()))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Preprocessing the digitial elevation model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "- Read the DEM from a ``geotiff`` file (in `WRADLIB_DATA`);\n",
    "- clip the region inside the bounding box;\n",
    "- map the DEM values to the polar grid points. \n",
    "\n",
    "*Note*: You can choose between the coarser resolution `bonn_gtopo.tif` (from GTOPO30) and the finer resolution `bonn_new.tif` (from the SRTM mission).\n",
    "\n",
    "The DEM raster data is opened via [wradlib.io.open_raster](http://wradlib.org/wradlib-docs/latest/generated/wradlib.io.open_raster.html) and extracted via [wradlib.georef.extract_raster_dataset](http://wradlib.org/wradlib-docs/latest/generated/wradlib.georef.extract_raster_dataset.html)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#rasterfile = wrl.util.get_wradlib_data_file('geo/bonn_gtopo.tif')\n",
    "rasterfile = wrl.util.get_wradlib_data_file('geo/bonn_new.tif')\n",
    "\n",
    "ds = wrl.io.open_raster(rasterfile)\n",
    "rastervalues, rastercoords, proj = wrl.georef.extract_raster_dataset(ds, nodata=-32768.)\n",
    "\n",
    "# Clip the region inside our bounding box \n",
    "ind = wrl.util.find_bbox_indices(rastercoords, rlimits)\n",
    "rastercoords = rastercoords[ind[1]:ind[3], ind[0]:ind[2], ...]\n",
    "rastervalues = rastervalues[ind[1]:ind[3], ind[0]:ind[2]]\n",
    "\n",
    "# Map rastervalues to polar grid points\n",
    "polarvalues = wrl.ipol.cart2irregular_spline(rastercoords, rastervalues,\n",
    "                                             polcoords, order=3,\n",
    "                                             prefilter=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Calculate Beam-Blockage"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Now we can finally apply the [wradlib.qual.beam_block_frac](http://wradlib.org/wradlib-docs/latest/generated/wradlib.qual.beam_block_frac.html) function to calculate the PBB."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "PBB = wrl.qual.beam_block_frac(polarvalues, alt, beamradius)\n",
    "PBB = np.ma.masked_invalid(PBB)\n",
    "print(PBB.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "So far, we calculated the fraction of beam blockage for each bin.\n",
    "\n",
    "But we need to into account that the radar signal travels along a beam. Cumulative beam blockage (CBB) in one bin along a beam will always be at least as high as the maximum PBB of the preceeding bins (see [wradlib.qual.cum_beam_block_frac](http://wradlib.org/wradlib-docs/latest/generated/wradlib.qual.cum_beam_block_frac.html))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "CBB = wrl.qual.cum_beam_block_frac(PBB)\n",
    "print(CBB.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Visualize Beamblockage"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Now we visualize\n",
    "- the average terrain altitude per radar bin\n",
    "- a beam blockage map\n",
    "- interaction with terrain along a single beam"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "# just a little helper function to style x and y axes of our maps\n",
    "def annotate_map(ax, cm=None, title=\"\"):\n",
    "    ticks = (ax.get_xticks()/1000).astype(np.int)\n",
    "    ax.set_xticklabels(ticks)\n",
    "    ticks = (ax.get_yticks()/1000).astype(np.int)\n",
    "    ax.set_yticklabels(ticks)\n",
    "    ax.set_xlabel(\"Kilometers\")\n",
    "    ax.set_ylabel(\"Kilometers\")\n",
    "    if not cm is None:\n",
    "        pl.colorbar(cm, ax=ax)\n",
    "    if not title==\"\":\n",
    "        ax.set_title(title)\n",
    "    ax.grid()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "fig = pl.figure(figsize=(10, 8))\n",
    "\n",
    "# create subplots\n",
    "ax1 = pl.subplot2grid((2, 2), (0, 0))\n",
    "ax2 = pl.subplot2grid((2, 2), (0, 1))\n",
    "ax3 = pl.subplot2grid((2, 2), (1, 0), colspan=2, rowspan=1)\n",
    "\n",
    "# azimuth angle\n",
    "angle = 225\n",
    "\n",
    "# Plot terrain (on ax1)\n",
    "ax1, dem = wrl.vis.plot_ppi(polarvalues, \n",
    "                            ax=ax1, r=r, \n",
    "                            az=np.degrees(coord[:,0,1]), \n",
    "                            cmap=mpl.cm.terrain, vmin=0.)\n",
    "ax1.plot([0,np.sin(np.radians(angle))*1e5],\n",
    "         [0,np.cos(np.radians(angle))*1e5],\"r-\")\n",
    "ax1.plot(sitecoords[0], sitecoords[1], 'ro')\n",
    "annotate_map(ax1, dem, 'Terrain within {0} km range'.format(np.max(r / 1000.) + 0.1))\n",
    "\n",
    "# Plot CBB (on ax2)\n",
    "ax2, cbb = wrl.vis.plot_ppi(CBB, ax=ax2, r=r, \n",
    "                            az=np.degrees(coord[:,0,1]),\n",
    "                            cmap=mpl.cm.PuRd, vmin=0, vmax=1)\n",
    "annotate_map(ax2, cbb, 'Beam-Blockage Fraction')\n",
    "\n",
    "# Plot single ray terrain profile on ax3\n",
    "bc, = ax3.plot(r / 1000., alt[angle, :], '-b',\n",
    "               linewidth=3, label='Beam Center')\n",
    "b3db, = ax3.plot(r / 1000., (alt[angle, :] + beamradius), ':b',\n",
    "                 linewidth=1.5, label='3 dB Beam width')\n",
    "ax3.plot(r / 1000., (alt[angle, :] - beamradius), ':b')\n",
    "ax3.fill_between(r / 1000., 0.,\n",
    "                 polarvalues[angle, :],\n",
    "                 color='0.75')\n",
    "ax3.set_xlim(0., np.max(r / 1000.) + 0.1)\n",
    "ax3.set_ylim(0., 3000)\n",
    "ax3.set_xlabel('Range (km)')\n",
    "ax3.set_ylabel('Altitude (m)')\n",
    "ax3.grid()\n",
    "\n",
    "axb = ax3.twinx()\n",
    "bbf, = axb.plot(r / 1000., CBB[angle, :], '-k',\n",
    "                label='BBF')\n",
    "axb.set_ylabel('Beam-blockage fraction')\n",
    "axb.set_ylim(0., 1.)\n",
    "axb.set_xlim(0., np.max(r / 1000.) + 0.1)\n",
    "\n",
    "\n",
    "legend = ax3.legend((bc, b3db, bbf), \n",
    "                    ('Beam Center', '3 dB Beam width', 'BBF'),\n",
    "                    loc='upper left', fontsize=10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Visualize Beam Propagation showing earth curvature"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Now we visualize\n",
    "- interaction with terrain along a single beam\n",
    "\n",
    "In this representation the earth curvature is shown. For this we assume the earth a sphere with exactly 6370000 m radius. This is needed to get the height ticks at nice position."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "def height_formatter(x, pos):\n",
    "    x = (x - 6370000) / 1000\n",
    "    fmt_str = '{:g}'.format(x)\n",
    "    return fmt_str\n",
    "    \n",
    "def range_formatter(x, pos):\n",
    "    x = x / 1000.\n",
    "    fmt_str = '{:g}'.format(x)\n",
    "    return fmt_str"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "- The [wradlib.vis.create_cg](http://wradlib.org/wradlib-docs/latest/generated/wradlib.vis.create_cg.html)-function is facilitated to create the curved geometries. \n",
    "- The actual data is plottet as (theta, range) on the parasite axis. \n",
    "- Some tweaking is needed to get the final plot look nice."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "fig = pl.figure(figsize=(10, 6))\n",
    "\n",
    "cgax, caax, paax = wrl.vis.create_cg('RHI', fig, 111)\n",
    "\n",
    "# azimuth angle\n",
    "angle = 225\n",
    "\n",
    "# fix grid_helper\n",
    "er = 6370000\n",
    "gh = cgax.get_grid_helper()\n",
    "gh.grid_finder.grid_locator2._nbins=80\n",
    "gh.grid_finder.grid_locator2._steps=[1,2,4,5,10]\n",
    "\n",
    "# calculate beam_height and arc_distance for ke=1\n",
    "# means line of sight\n",
    "bhe = wrl.georef.bin_altitude(r, 0, sitecoords[2], re=er, ke=1.)\n",
    "ade = wrl.georef.bin_distance(r, 0, sitecoords[2], re=er, ke=1.)\n",
    "nn0 = np.zeros_like(r)\n",
    "# for nice plotting we assume earth_radius = 6370000 m\n",
    "ecp = nn0 + er\n",
    "# theta (arc_distance sector angle)\n",
    "thetap = - np.degrees(ade/er) + 90.0\n",
    "\n",
    "# zero degree elevation with standard refraction\n",
    "bh0 = wrl.georef.bin_altitude(r, 0, sitecoords[2], re=er)\n",
    "\n",
    "# plot (ecp is earth surface normal null)\n",
    "bes, = paax.plot(thetap, ecp, '-k', linewidth=3, label='Earth Surface NN')\n",
    "bc, = paax.plot(thetap, ecp + alt[angle, :], '-b', linewidth=3, label='Beam Center')\n",
    "bc0r, = paax.plot(thetap, ecp + bh0 + alt[angle, 0] , '-g', label='0 deg Refraction')\n",
    "bc0n, = paax.plot(thetap, ecp + bhe + alt[angle, 0], '-r', label='0 deg line of sight')\n",
    "b3db, = paax.plot(thetap, ecp + alt[angle, :]  + beamradius, ':b', label='+3 dB Beam width')\n",
    "paax.plot(thetap, ecp + alt[angle, :] - beamradius, ':b', label='-3 dB Beam width')\n",
    "\n",
    "# orography\n",
    "paax.fill_between(thetap, ecp,\n",
    "                  ecp + polarvalues[angle, :],\n",
    "                  color='0.75')\n",
    "\n",
    "# shape axes\n",
    "cgax.set_xlim(0, np.max(ade))\n",
    "cgax.set_ylim([ecp.min()-1000, ecp.max()+2500])\n",
    "caax.grid(True, axis='x')\n",
    "cgax.grid(True, axis='y')\n",
    "cgax.axis['top'].toggle(all=False)\n",
    "caax.yaxis.set_major_locator(mpl.ticker.MaxNLocator(steps=[1,2,4,5,10], nbins=20, prune='both'))\n",
    "caax.xaxis.set_major_locator(mpl.ticker.MaxNLocator())\n",
    "caax.yaxis.set_major_formatter(mpl.ticker.FuncFormatter(height_formatter))\n",
    "caax.xaxis.set_major_formatter(mpl.ticker.FuncFormatter(range_formatter))\n",
    "\n",
    "caax.set_xlabel('Range (km)')\n",
    "caax.set_ylabel('Altitude (km)')\n",
    "    \n",
    "legend = paax.legend((bes, bc0n, bc0r, bc, b3db), \n",
    "                    ('Earth Surface NN', '0 deg line of sight', '0 deg std refraction', 'Beam Center', '3 dB Beam width'),\n",
    "                    loc='upper left', fontsize=10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Go back to [Read DEM Raster Data](#Read-DEM-Raster-Data), change the rasterfile to use the other resolution DEM and process again."
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Slideshow",
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "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.6.3"
  },
  "livereveal": {
   "scroll": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}