notebooks/match3d/wradlib_match_workflow.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "nbsphinx": "hidden"
   },
   "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": {},
   "source": [
    "# Match spaceborn SR (GPM/TRMM) with ground radars GR"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The idea is to match ground radar (GR) and space-born radar (SR) measurements in order to create spatially and temporally coicident samples without interpolation. The procedure had been suggested by [Schwaller and Morris (2011)](http://journals.ametsoc.org/doi/pdf/10.1175/2010JTECHA1403.1) and is based on the adaption by [Warren, et. al. (2017)](https://journals.ametsoc.org/doi/pdf/10.1175/JTECH-D-17-0128.1). \n",
    "\n",
    "The basic principle is illustrated in Fig. 2 of the original paper of Schwaller and Morris (2011):\n",
    "\n",
    "![image](fig2_schwaller_morris_2011.png)\n",
    "\n",
    "*Quote Warren, et.al: \"[...] In this approach, intersections between indi vidual SR beams and GR elevation sweeps are identified and the reflectivity values from both  instruments are averaged within a spatial neighborhood around the\n",
    "intersection. Specifically, SR data are averaged in range over the width of the GR beam at the GR range of the intersection, while GR data are averaged in the\n",
    "range–azimuth plane within the footprint of the SR beam. The result is a pair of reflectivity measurements corresponding to approximately the same volume of\n",
    "atmosphere. [...]\".*\n",
    "\n",
    "This becomes clearer in Fig. 3:\n",
    "![image](fig3_schwaller_morris_2011.png)\n",
    "\n",
    "Schwaller, MR, and Morris, KR. 2011. A ground validation network for the Global Precipitation Measurement mission. J. Atmos. Oceanic Technol., 28, 301-319.\n",
    "\n",
    "Warren, R.A., A. Protat, S.T. Siems, H.A. Ramsay, V. Louf, M.J. Manton, and T.A. Kane, 0: Calibrating ground-based radars against TRMM and GPM. J. Atmos. Oceanic Technol., 0,"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Conventions\n",
    "\n",
    "This code is based on the following conventions:\n",
    "\n",
    "- `gr` indicates **g**round **r**adar\n",
    "- `sr` indicates **s**pace-born precipitation **r**adar (TRMM or GPM)\n",
    "\n",
    "The base routines are designed to process one GR sweep at a time. If a full GR volume with `nelev` of sweeps is available, you can iterate over each sweep. In this code, `ee` is an index that points to one of the `nelev` sweeps/elevation angles. Accordingly, a **GR** data set will be organised as an array of shape `(nelev_gr, nray_gr, ngate_gr)`.\n",
    "\n",
    "A **SR** data set is typically organised as arrays with dimensions `(nscan_sr, nray_sr, ngate_sr)`.\n",
    "\n",
    "*Please note*: the GR and SR datasets are managed by using convenience functions from the module file `io_func.py` located in the same directory as this notebook. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import warnings\n",
    "warnings.filterwarnings('ignore')\n",
    "import wradlib as wrl\n",
    "import matplotlib.pyplot as pl\n",
    "import matplotlib as mpl\n",
    "from matplotlib.collections import PatchCollection\n",
    "from matplotlib.colors import from_levels_and_colors\n",
    "from matplotlib.path import Path\n",
    "import matplotlib.patches as patches\n",
    "import matplotlib.cm as cm\n",
    "try:\n",
    "    get_ipython().magic(\"matplotlib inline\")\n",
    "except:\n",
    "    pl.ion()\n",
    "import numpy as np\n",
    "import datetime as dt\n",
    "from osgeo import osr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# define GPM data set\n",
    "gpm_file = wrl.util.get_wradlib_data_file('gpm/2A-CS-151E24S154E30S.GPM.Ku.V7-20170308.20141206-S095002-E095137.004383.V05A.HDF5')\n",
    "# define matching ground radar file\n",
    "gr2gpm_file = wrl.util.get_wradlib_data_file('hdf5/IDR66_20141206_094829.vol.h5')\n",
    "\n",
    "# define TRMM data sets\n",
    "trmm_2a23_file = wrl.util.get_wradlib_data_file('trmm/2A-CS-151E24S154E30S.TRMM.PR.2A23.20100206-S111425-E111526.069662.7.HDF')\n",
    "trmm_2a25_file = wrl.util.get_wradlib_data_file('trmm/2A-CS-151E24S154E30S.TRMM.PR.2A25.20100206-S111425-E111526.069662.7.HDF')\n",
    "\n",
    "# define matching ground radar file\n",
    "gr2trmm_file = wrl.util.get_wradlib_data_file('hdf5/IDR66_20100206_111233.vol.h5')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Space-born precipitation radar parameters\n",
    "sr_pars = {\"trmm\": {\n",
    "    \"zt\": 402500.,  # orbital height of TRMM (post boost)   APPROXIMATION!\n",
    "    \"dr\": 250.,     # gate spacing of TRMM\n",
    "    \"gr_file\": gr2trmm_file,\n",
    "}, \"gpm\": {\n",
    "    \"zt\": 407000.,  # orbital height of GPM                 APPROXIMATION!\n",
    "    \"dr\": 125.,      # gate spacing of GPM\n",
    "    \"gr_file\": gr2gpm_file,\n",
    "}}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set parameters for this procedure\n",
    "bw_sr = 0.71                  # SR beam width\n",
    "platf = \"trmm\"                 # SR platform/product: one out of [\"gpm\", \"trmm\"]\n",
    "zt = sr_pars[platf][\"zt\"]     # SR orbit height (meters)\n",
    "dr_sr = sr_pars[platf][\"dr\"]  # SR gate length (meters)\n",
    "gr_file = sr_pars[platf][\"gr_file\"]\n",
    "ee = 2                        # Index that points to the GR elevation angle to be used"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Data Input"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Read and organize the GR data\n",
    "\n",
    "The following code reads data in ODIM H5 format. If your GR data is in some other format respective adaptions are needed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def _get_tilts(dic):\n",
    "    i = 0\n",
    "    for k in dic.keys():\n",
    "        if 'dataset' in k:\n",
    "            i += 1\n",
    "    return i\n",
    "\n",
    "\n",
    "def read_gr(filename, loaddata=True):\n",
    "\n",
    "    gr_data = wrl.io.read_generic_netcdf(filename)\n",
    "    dat = gr_data['what']['date']\n",
    "    tim = gr_data['what']['time']\n",
    "    date = dt.datetime.strptime(dat + tim, \"%Y%d%m%H%M%S\")\n",
    "    source = gr_data['what']['source']\n",
    "\n",
    "    lon = gr_data['where']['lon']\n",
    "    lat = gr_data['where']['lat']\n",
    "    alt = gr_data['where']['height']\n",
    "\n",
    "    if gr_data['what']['object'] == 'PVOL':\n",
    "        ntilt = _get_tilts(gr_data)\n",
    "    else:\n",
    "        raise ValueError('GR file is no PPI/Volume File')\n",
    "\n",
    "    ngate = np.zeros(ntilt, dtype=np.int16)\n",
    "    nbeam = np.zeros(ntilt)\n",
    "    elang = np.zeros(ntilt)\n",
    "    r0 = np.zeros(ntilt)\n",
    "    dr = np.zeros(ntilt)\n",
    "    a0 = np.zeros(ntilt)\n",
    "\n",
    "    for i in range(0, ntilt):\n",
    "        dset = gr_data['dataset{0}'.format(i+1)]\n",
    "        a0[i] = dset['how']['astart']\n",
    "        elang[i] = dset['where']['elangle']\n",
    "        ngate[i] = dset['where']['nbins']\n",
    "        r0[i] = dset['where']['rstart']\n",
    "        dr[i] = dset['where']['rscale']\n",
    "        nbeam[i] = dset['where']['nrays']\n",
    "\n",
    "    if ((len(np.unique(r0)) != 1) |\n",
    "            (len(np.unique(dr)) != 1) |\n",
    "            (len(np.unique(a0)) != 1) |\n",
    "            (len(np.unique(nbeam)) != 1) |\n",
    "            (nbeam[0] != 360)):\n",
    "        raise ValueError('GroundRadar Data layout dos not match')\n",
    "\n",
    "    gr_dict = {}\n",
    "    gr_dict.update({'source': source, 'date': date, 'lon': lon, 'lat': lat,\n",
    "                    'alt': alt, 'ngate': ngate, 'nbeam': nbeam, 'ntilt': ntilt,\n",
    "                    'r0': r0, 'dr': dr, 'a0': a0, 'elang': elang})\n",
    "    if not loaddata:\n",
    "        return gr_dict\n",
    "\n",
    "    sdate = []\n",
    "    refl = []\n",
    "    for i in range(0, ntilt):\n",
    "        dset = gr_data['dataset{0}'.format(i+1)]\n",
    "        dat = dset['what']['startdate']\n",
    "        tim = dset['what']['starttime']\n",
    "        date = dt.datetime.strptime(dat + tim, \"%Y%d%m%H%M%S\")\n",
    "        sdate.append(date)\n",
    "        data = dset['data1']\n",
    "        quantity = data['what']['quantity']\n",
    "        factor = data['what']['gain']\n",
    "        offset = data['what']['offset']\n",
    "        if quantity == 'DBZH':\n",
    "            dat = data['variables']['data']['data'] * factor + offset\n",
    "            refl.append(dat)\n",
    "\n",
    "    sdate = np.array(sdate)\n",
    "    refl = np.array(refl)\n",
    "\n",
    "    gr_dict.update({'sdate': sdate, 'refl': refl})\n",
    "\n",
    "    return gr_dict"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# read matching GR data\n",
    "gr_data = read_gr(gr_file)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Data handling"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Extract relevant GR data and meta-data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# number of rays in gr sweep\n",
    "nray_gr = gr_data['nbeam'].astype(\"i4\")[ee]\n",
    "# number of gates in gr beam\n",
    "ngate_gr = gr_data['ngate'].astype(\"i4\")[ee]\n",
    "# number of sweeps\n",
    "nelev = gr_data['ntilt']\n",
    "# elevation of sweep (degree)\n",
    "elev_gr = gr_data['elang'][ee]\n",
    "# gate length (meters)\n",
    "dr_gr = gr_data['dr'][ee]\n",
    "# reflectivity array of sweep\n",
    "ref_gr = gr_data['refl'][ee]\n",
    "# sweep datetime stamp\n",
    "date_gr = gr_data['sdate'][ee]\n",
    "# range of first gate\n",
    "r0_gr = gr_data['r0'][ee]                    \n",
    "# azimuth angle of first beam\n",
    "a0_gr = gr_data['a0'][ee]\n",
    "# Longitude of GR\n",
    "lon0_gr = gr_data['lon']\n",
    "# Latitude of GR\n",
    "lat0_gr = gr_data['lat']\n",
    "# Altitude of GR (meters)\n",
    "alt0_gr = gr_data['alt']\n",
    "# Beam width of GR (degree)\n",
    "bw_gr = 1.\n",
    "print(elev_gr, lon0_gr)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Calculate GR Bounding Box"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "coord = wrl.georef.sweep_centroids(nray_gr, dr_gr, ngate_gr, elev_gr)\n",
    "coords = wrl.georef.spherical_to_proj(coord[..., 0], \n",
    "                                      np.degrees(coord[..., 1]),\n",
    "                                      coord[..., 2], \n",
    "                                      (lon0_gr, lat0_gr, alt0_gr))\n",
    "lon = coords[..., 0]\n",
    "lat = coords[..., 1]\n",
    "alt = coords[..., 2]\n",
    "bbox = wrl.zonalstats.get_bbox(lon, lat)\n",
    "print(\"Radar bounding box:\\n\\t%.2f\\n%.2f           %.2f\\n\\t%.2f\" % \n",
    "      (bbox['top'], bbox['left'], bbox['right'], bbox['bottom']))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Read and organize the SR data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# read spaceborn SR data\n",
    "if platf == \"gpm\":\n",
    "    sr_data = wrl.io.read_gpm(gpm_file, bbox)\n",
    "elif platf == \"trmm\":\n",
    "    sr_data = wrl.io.read_trmm(trmm_2a23_file, trmm_2a25_file, bbox)\n",
    "else:\n",
    "    raise(\"Invalid platform\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "refl = sr_data['refl']\n",
    "print(refl)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Extract relevant SR data and meta-data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Longitudes of SR scans\n",
    "sr_lon = sr_data['lon']\n",
    "# Latitudes of SR scans\n",
    "sr_lat = sr_data['lat']\n",
    "# Precip flag\n",
    "pflag = sr_data['pflag']\n",
    "# Number of scans on SR data\n",
    "nscan_sr= sr_data['nscan']\n",
    "# Number of rays in one SR scan\n",
    "nray_sr = sr_data['nray']\n",
    "# Number of gates in one SR ray\n",
    "ngate_sr = sr_data['nbin']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(sr_lon.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Georeferencing"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Set fundamental georeferencing parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Calculate equivalent earth radius\n",
    "wgs84 = wrl.georef.get_default_projection()\n",
    "re1 = wrl.georef.get_earth_radius(lat0_gr, wgs84) \n",
    "print(\"Earth radius 1:\", re1)\n",
    "a = wgs84.GetSemiMajor()\n",
    "b = wgs84.GetSemiMinor()\n",
    "print(\"SemiMajor, SemiMinor:\", a, b)\n",
    "\n",
    "# Set up aeqd-projection gr-centered\n",
    "rad = wrl.georef.proj4_to_osr(('+proj=aeqd +lon_0={lon:f} ' + \n",
    "                                   '+lat_0={lat:f} +a={a:f} ' +\n",
    "                                   '+b={b:f}').format(lon=lon0_gr,\n",
    "                                                      lat=lat0_gr,\n",
    "                                                      a=a, b=b))\n",
    "re2 = wrl.georef.get_earth_radius(lat0_gr, rad)\n",
    "print(\"Earth radius 2:\", re2)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Georeference GR data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# create gr range and azimuth arrays\n",
    "rmax_gr = r0_gr + ngate_gr * dr_gr\n",
    "r_gr = np.arange(0, ngate_gr) * dr_gr + dr_gr/2.\n",
    "az_gr = np.arange(0, nray_gr) - a0_gr\n",
    "print(\"Range/Azi-Shape:\", r_gr.shape, az_gr.shape)\n",
    "\n",
    "# create gr polar grid and calculate aeqd-xyz coordinates\n",
    "gr_polargrid = np.meshgrid(r_gr, az_gr)\n",
    "gr_xyz, rad = wrl.georef.spherical_to_xyz(gr_polargrid[0], gr_polargrid[1], elev_gr, (lon0_gr, lat0_gr, alt0_gr ))\n",
    "print(\"XYZ-Grid-Shape:\", gr_xyz.shape)\n",
    "\n",
    "# create gr poygon array in aeqd-xyz-coordinates\n",
    "gr_poly, rad1 = wrl.georef.spherical_to_polyvert(r_gr, az_gr, elev_gr, (lon0_gr, lat0_gr, alt0_gr))\n",
    "print(gr_poly.shape, 360 * 600)\n",
    "gr_poly.shape = (nray_gr, ngate_gr, 5, 3)\n",
    "\n",
    "# get radar domain (outer ring)\n",
    "gr_domain = gr_xyz[:,-1,0:2]\n",
    "gr_domain = np.vstack((gr_domain, gr_domain[0]))\n",
    "print(\"Domain-Shape:\", gr_domain.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Georeference SR data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sr_x, sr_y = wrl.georef.reproject(sr_lon, sr_lat, \n",
    "                                      projection_source=wgs84, \n",
    "                                      projection_target=rad)\n",
    "sr_xy = np.dstack((sr_x, sr_y))\n",
    "print(\"SR-GRID-Shapes:\", sr_x.shape, sr_y.shape, sr_xy.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Subset relevant SR data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create ZonalData for spatial subsetting (inside GR range domain)\n",
    "\n",
    "# get precip indexes\n",
    "print(sr_xy.shape)\n",
    "print(pflag.shape)\n",
    "\n",
    "precip_mask = (pflag == 2) & wrl.zonalstats.get_clip_mask(sr_xy, gr_domain, rad)\n",
    "\n",
    "# get iscan/iray boolean arrays\n",
    "print(precip_mask.shape, sr_xy.shape)\n",
    "print(pflag.shape, sr_xy.reshape(-1, sr_xy.shape[-1]).shape)\n",
    "pl.imshow(precip_mask)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### SR Parallax Correction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "print(\"NRAY\", nray_sr)\n",
    "print(\"NBIN\", ngate_sr)\n",
    "\n",
    "# use localZenith Angle\n",
    "alpha = sr_data['zenith']\n",
    "beta = abs(-17.04 + np.arange(nray_sr) * bw_sr)\n",
    "\n",
    "# Correct for parallax, get 3D-XYZ-Array\n",
    "#   xyzp_sr: Parallax corrected xyz coordinates\n",
    "#   r_sr_inv: range array from ground to SR platform\n",
    "#   zp: SR bin altitudes\n",
    "xyp_sr, r_sr_inv, z_sr = wrl.georef.correct_parallax(sr_xy, ngate_sr, dr_sr, alpha)\n",
    "print(xyp_sr.shape, r_sr_inv.shape, z_sr.shape)\n",
    "xyzp_sr = np.concatenate((xyp_sr, z_sr[..., np.newaxis]),\n",
    "                   axis=-1)\n",
    "print(sr_xy.shape)\n",
    "print(\"SR_XYP:\", xyp_sr.shape, xyzp_sr.shape, r_sr_inv.shape, z_sr.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Compute spherical coordinates of SR bins with regard to GR"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "r_sr, az_sr, elev_sr = wrl.georef.xyz_to_spherical(xyzp_sr, alt0_gr, proj=rad)    \n",
    "#TODO: hardcoded 1.0\n",
    "mask = (elev_sr > (1.0 - bw_gr/2.)) & (elev_sr < (1.0 + bw_gr/2.))\n",
    "pl.figure()\n",
    "pl.pcolormesh(mask[40,:,:].T)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Compute SR and GR pulse volumes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Calculate distance from orbit $r_s$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rs = wrl.georef.dist_from_orbit(zt, alpha, beta, r_sr_inv, re1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### SR pulse volume"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Small anngle approximation\n",
    "vol_sr2  = np.pi * dr_sr * rs**2 * np.radians(bw_sr / 2.)**2\n",
    "\n",
    "# Or using wradlib's native function\n",
    "vol_sr = wrl.qual.pulse_volume(rs, dr_sr, bw_sr)\n",
    "# Evaluate difference between both approaches\n",
    "print(\"Min. difference (m3):\", (vol_sr - vol_sr2).min())\n",
    "print(\"Max. difference (m3): \", (vol_sr - vol_sr2).max())\n",
    "print(\"Average rel. difference (%):\", round(np.mean(vol_sr-vol_sr2)*100./np.mean(np.mean(vol_sr2)), 4))\n",
    "\n",
    "# Verdict: differences are negligble - use wradlibs's native function!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### GR pulse volume"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# GR pulse volumes\n",
    "#   along one beam\n",
    "vol_gr = wrl.qual.pulse_volume(r_gr, dr_gr, bw_gr)\n",
    "#   with shape (nray_gr, ngate_gr)\n",
    "vol_gr = np.repeat(vol_gr, nray_gr).reshape((nray_gr, ngate_gr), order=\"F\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Calculate horizontal and vertical dimensions $R_s$ and $D_s$ of SR bins"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Rs = 0.5 * (1 +  np.cos(np.radians(alpha)))[:,:,np.newaxis] * rs * np.tan(np.radians(bw_sr/2.))\n",
    "Ds = dr_sr / np.cos(np.radians(alpha))\n",
    "Ds = np.broadcast_to(Ds[..., np.newaxis], Rs.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Median Brightband Width/Height"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(z_sr.shape)\n",
    "print(sr_data['zbb'].shape, sr_data['bbwidth'].shape, sr_data['quality'].shape)\n",
    "ratio, ibb = wrl.qual.get_bb_ratio(sr_data['zbb'], sr_data['bbwidth'], sr_data['quality'], z_sr)\n",
    "print(ratio.shape)\n",
    "zbb = sr_data['zbb'].copy()\n",
    "zbb[~ibb] = np.nan\n",
    "print(np.nanmin(ratio[..., 9]), np.nanmax(ratio[..., 9]))\n",
    "#pl.imshow(ratio, vmin=-10, vmax=30)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Convert SR Ku reflectivities to S-band \n",
    "\n",
    "Based on [Cao et.al (2013)](http://dx.doi.org/10.1002/jgrd.50138)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ref_sr = sr_data['refl'].filled(np.nan)\n",
    "print(ref_sr.shape, ratio.shape)\n",
    "ref_sr_ss = np.zeros_like(ref_sr) * np.nan\n",
    "ref_sr_sh = np.zeros_like(ref_sr) * np.nan\n",
    "\n",
    "a_s, a_h = (wrl.trafo.ku2s.snow, wrl.trafo.ku2s.hail) \n",
    "\n",
    "ia = (ratio >= 1)\n",
    "ref_sr_ss[ia] = ref_sr[ia] + wrl.util.calculate_polynomial(ref_sr[ia], a_s[:,10])\n",
    "ref_sr_sh[ia] = ref_sr[ia] + wrl.util.calculate_polynomial(ref_sr[ia], a_h[:,10])\n",
    "ib = (ratio <= 0)\n",
    "ref_sr_ss[ib] = ref_sr[ib] + wrl.util.calculate_polynomial(ref_sr[ib], a_s[:,0])\n",
    "ref_sr_sh[ib] = ref_sr[ib] + wrl.util.calculate_polynomial(ref_sr[ib], a_h[:,0])\n",
    "im = (ratio > 0) & (ratio < 1)\n",
    "ind = np.round(ratio[im] * 10).astype(np.int)\n",
    "ref_sr_ss[im] = ref_sr[im] + wrl.util.calculate_polynomial(ref_sr[im], a_s[:,ind])\n",
    "ref_sr_sh[im] = ref_sr[im] + wrl.util.calculate_polynomial(ref_sr[im], a_h[:,ind])\n",
    "\n",
    "# Jackson Tan's fix for C-band\n",
    "is_cband = False\n",
    "if (is_cband):\n",
    "    deltas = (ref_sr_ss - ref_sr) * 5.3 / 10.0\n",
    "    ref_sr_ss = outref_sr + deltas\n",
    "    deltah = (ref_sr_sh - ref_sr) * 5.3 / 10.0\n",
    "    ref_sr_sh = ref_sr + deltah\n",
    "\n",
    "ref_sr_ss[ref_sr < 0] = np.nan"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Convert S-band GR reflectivities to Ku-band\n",
    "\n",
    "Using the method of [Liao and Meneghini (2009)](http://dx.doi.org/10.1175/2008JAMC1974.1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Convert S-band GR reflectivities to Ku-band using method of Liao and Meneghini (2009)\n",
    "ref_gr_ku = np.zeros_like(ref_gr) * np.nan\n",
    "\n",
    "# Which zbb value should we take here???\n",
    "#    Q'n'Dirty: just take the mean of all SR profiles\n",
    "#    TODO: Consider zbb for each profile during the matching process\n",
    "\n",
    "# Snow\n",
    "ia = ( gr_xyz[...,2] >= np.nanmean(zbb) )\n",
    "#ref_gr_ku[ia] = wrl.trafo.ku2s.snow[0] + wrl.trafo.ku2s.snow[1]*ref_gr[ia] + wrl.trafo.ku2s.snow[2]*ref_gr[ia]**2\n",
    "ref_gr_ku[ia] = wrl.util.calculate_polynomial(ref_gr[ia], wrl.trafo.s2ku.snow)\n",
    "# Rain\n",
    "ib = ( gr_xyz[...,2] < np.nanmean(zbb) )\n",
    "#ref_gr_ku[ib] = wrl.trafo.ku2s.rain[0] + wrl.trafo.ku2s.rain[1]*ref_gr[ia] + wrl.trafo.ku2s.rain[2]*ref_gr[ia]**2\n",
    "ref_gr_ku[ib] = wrl.util.calculate_polynomial(ref_gr[ib], wrl.trafo.s2ku.rain)\n",
    "\n",
    "# Jackson Tan's fix for C-band\n",
    "is_cband = False\n",
    "if (is_cband):\n",
    "    delta = (ref_gr_ku - ref_gr) * 5.3/10.0\n",
    "    ref_gr_ku = ref_gr + delta"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Matching SR/GR"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Identify which SR rays actually intersect with the GR sweep\n",
    "\n",
    "Based on the above criteria (in radar range, precipitating SR profile) and based on SR elevation angle (with regard to GR)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# First assumption: no valid SR bins (all False)\n",
    "valid = np.asarray(elev_sr, dtype=np.bool)==False\n",
    "print(valid.shape, precip_mask.shape)\n",
    "# SR is inside GR range and is precipitating\n",
    "iscan = precip_mask.nonzero()[0]\n",
    "iray = precip_mask.nonzero()[1]\n",
    "valid[iscan,iray] = True\n",
    "# SR bins intersect with GR sweep\n",
    "valid = valid & (elev_sr >= (elev_gr-bw_gr/2.)) & (elev_sr <= (elev_gr+bw_gr/2.))\n",
    "# Number of matching SR bins per profile\n",
    "nvalids = np.sum(valid, axis=2)\n",
    "# scan and ray indices for profiles with at least one valid bin\n",
    "vscan, vray = np.where(nvalids>0)\n",
    "# number of profiles with at least one valid bin\n",
    "nprof = len(vscan)\n",
    "print(vscan.shape)\n",
    "print(valid.shape)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Averaging SR parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# average coordinates\n",
    "xyz_v1 = xyzp_sr.copy()\n",
    "print(xyz_v1.shape)\n",
    "xyz_v1[~valid] = np.nan\n",
    "xyz_c1 = xyzp_sr.copy()\n",
    "xyz_c1[~valid] = 0\n",
    "c = np.count_nonzero(xyz_c1[..., 0], axis=2)\n",
    "ntotsr = c[vscan, vray]\n",
    "xyz_m1 = np.nanmean(xyz_v1,axis=2)\n",
    "xyz = xyz_m1[vscan, vray]\n",
    "print(xyz.shape, c.shape)\n",
    "\n",
    "# approximate Rs \n",
    "rs_v1 = Rs.copy()\n",
    "rs_v1[~valid] = np.nan\n",
    "rs_m1 = np.nanmax(rs_v1, axis=2)\n",
    "rs_prof = rs_m1[vscan, vray]\n",
    "ds = rs_prof\n",
    "\n",
    "# approximate Ds\n",
    "ds_v1 = Ds.copy()\n",
    "ds_v1[~valid] = np.nan\n",
    "ds_m1 = np.nansum(ds_v1, axis=2)\n",
    "ds_prof = ds_m1[vscan, vray]\n",
    "dz = ds_prof\n",
    "\n",
    "# approximate Vs\n",
    "vs_v1 = vol_sr.copy()\n",
    "vs_v1[~valid] = np.nan\n",
    "vs_m1 = np.nansum(vs_v1, axis=2)\n",
    "vs_prof = vs_m1[vscan, vray]\n",
    "volsr1 = vs_prof\n",
    "\n",
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "fig = pl.figure()\n",
    "ax = fig.add_subplot(111, projection='3d')\n",
    "ax.scatter(xyz_m1[..., 0], xyz_m1[..., 1], xyz_m1[..., 2], c=c.ravel())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Calculate spherical coordinates of SR sample volume with respect to GR"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "r_sr, az_sr, el_rs = wrl.georef.xyz_to_spherical(xyz, alt0_gr, proj=rad)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Average SR reflectivity"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ref_sr_1 = wrl.trafo.idecibel(ref_sr)\n",
    "ref_sr_1[~valid] = np.nan\n",
    "refsr1a = np.nanmean(ref_sr_1, axis=2)[vscan,vray]\n",
    "refsr1a = wrl.trafo.decibel(refsr1a)\n",
    "\n",
    "ref_sr_2 = wrl.trafo.idecibel(ref_sr_ss)\n",
    "ref_sr_2[~valid] = np.nan\n",
    "refsr2a = np.nanmean(ref_sr_2, axis=2)[vscan,vray]\n",
    "refsr2a = wrl.trafo.decibel(refsr2a)\n",
    "\n",
    "ref_sr_3 = wrl.trafo.idecibel(ref_sr_sh)\n",
    "ref_sr_3[~valid] = np.nan\n",
    "refsr3a = np.nanmean(ref_sr_3, axis=2)[vscan,vray]\n",
    "refsr3a = wrl.trafo.decibel(refsr3a)\n",
    "\n",
    "\n",
    "print(refsr1a.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Identify which GR bins actually intersect with the SR footprint"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Create Zonal Data Source of matching SR profiles"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "zds = wrl.zonalstats.DataSource(xyz[:, 0:2].reshape(-1, 2), rad)\n",
    "print(zds.ds.GetLayer().GetFeatureCount())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Create SR beam polygons using $R_s$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tmp_trg_lyr = zds.ds.GetLayer()\n",
    "trg_poly = []\n",
    "for i, feat in enumerate(tmp_trg_lyr):\n",
    "    geom = feat.GetGeometryRef().Buffer(rs_prof[i])\n",
    "    poly = wrl.georef.ogr_to_numpy(geom)\n",
    "    trg_poly.append(poly)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Create ZonalDataPoly for GR and SR profiles"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "print(\"creating\")\n",
    "zdp = wrl.zonalstats.ZonalDataPoly(gr_poly[..., 0:2].reshape(-1, 5, 2), trg_poly, srs=rad)\n",
    "zdp.dump_vector('m3d_zonal_poly_{0}'.format(platf))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "print(\"loading\")\n",
    "obj3 = wrl.zonalstats.GridCellsToPoly('m3d_zonal_poly_{0}'.format(platf))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Average GR volume and reflectivity"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "print(obj3.ix.shape)\n",
    "volgr1 = np.array([np.sum(vol_gr.ravel()[obj3.ix[i]])\n",
    "                   for i in np.arange(len(obj3.ix))[~obj3.check_empty()]])\n",
    "print(volgr1.shape)\n",
    "\n",
    "ref_gr_i = wrl.trafo.idecibel(ref_gr.ravel())\n",
    "ref_gr_ku_i = wrl.trafo.idecibel(ref_gr_ku.ravel())\n",
    "refgr1a = np.array([np.nanmean(ref_gr_i[obj3.ix[i]])\n",
    "             for i in np.arange(len(obj3.ix))[~obj3.check_empty()]])\n",
    "refgr1a = wrl.trafo.decibel(refgr1a)\n",
    "refgr2a = np.array([np.nanmean(ref_gr_ku_i[obj3.ix[i]])\n",
    "             for i in np.arange(len(obj3.ix))[~obj3.check_empty()]])\n",
    "refgr2a = wrl.trafo.decibel(refgr2a)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Diagnostic Plots"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = pl.figure(figsize=(12,5))\n",
    "ax = fig.add_subplot(121, aspect=\"equal\")\n",
    "pl.scatter(refgr1a, refsr1a, marker=\"+\", c=\"black\")\n",
    "pl.plot([0,60],[0,60], linestyle=\"solid\", color=\"black\")\n",
    "pl.xlim(10,50)\n",
    "pl.ylim(10,50)\n",
    "pl.xlabel(\"GR reflectivity (dBZ)\")\n",
    "pl.ylabel(\"SR reflectivity (dBZ)\")\n",
    "ax = fig.add_subplot(122)\n",
    "pl.hist(refgr1a[refsr1a>10], bins=np.arange(-10,50,5), edgecolor=\"None\", label=\"GR\")\n",
    "pl.hist(refsr1a[refsr1a>-10], bins=np.arange(-10,50,5), edgecolor=\"red\", facecolor=\"None\", label=\"SR\")\n",
    "pl.xlabel(\"Reflectivity (dBZ)\")\n",
    "pl.legend()\n",
    "fig.suptitle(\"uncorrected GR vs uncorrected SR\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = pl.figure(figsize=(12,5))\n",
    "ax = fig.add_subplot(121, aspect=\"equal\")\n",
    "pl.scatter(refgr2a, refsr1a, marker=\"+\", c=\"black\")\n",
    "pl.plot([0,60],[0,60], linestyle=\"solid\", color=\"black\")\n",
    "pl.xlim(10,50)\n",
    "pl.ylim(10,50)\n",
    "pl.xlabel(\"GR reflectivity (dBZ)\")\n",
    "pl.ylabel(\"SR reflectivity (dBZ)\")\n",
    "ax = fig.add_subplot(122)\n",
    "pl.hist(refgr1a[refsr1a>10], bins=np.arange(-10,50,5), edgecolor=\"None\", label=\"GR\")\n",
    "pl.hist(refsr1a[refsr1a>-10], bins=np.arange(-10,50,5), edgecolor=\"red\", facecolor=\"None\", label=\"SR\")\n",
    "pl.xlabel(\"Reflectivity (dBZ)\")\n",
    "pl.legend()\n",
    "fig.suptitle(\"corrected GR vs uncorrected SR\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = pl.figure(figsize=(12,5))\n",
    "ax = fig.add_subplot(121, aspect=\"equal\")\n",
    "pl.scatter(refgr2a, refsr2a, marker=\"+\", c=\"black\")\n",
    "pl.plot([0,60],[0,60], linestyle=\"solid\", color=\"black\")\n",
    "pl.xlim(10,50)\n",
    "pl.ylim(10,50)\n",
    "pl.xlabel(\"GR reflectivity (dBZ)\")\n",
    "pl.ylabel(\"SR reflectivity (dBZ)\")\n",
    "ax = fig.add_subplot(122)\n",
    "pl.hist(refgr2a[refsr2a>10], bins=np.arange(-10,50,5), edgecolor=\"None\", label=\"GR\")\n",
    "pl.hist(refsr2a[refsr2a>-10], bins=np.arange(-10,50,5), edgecolor=\"red\", facecolor=\"None\", label=\"SR\")\n",
    "pl.xlabel(\"Reflectivity (dBZ)\")\n",
    "pl.legend()\n",
    "fig.suptitle(\"corrected GR vs corrected SR snow\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = pl.figure(figsize=(12,5))\n",
    "ax = fig.add_subplot(121, aspect=\"equal\")\n",
    "pl.scatter(refgr2a, refsr3a, marker=\"+\", c=\"black\")\n",
    "pl.plot([0,60],[0,60], linestyle=\"solid\", color=\"black\")\n",
    "pl.xlim(10,50)\n",
    "pl.ylim(10,50)\n",
    "pl.xlabel(\"GR reflectivity (dBZ)\")\n",
    "pl.ylabel(\"SR reflectivity (dBZ)\")\n",
    "ax = fig.add_subplot(122)\n",
    "pl.hist(refgr2a[refsr3a>10], bins=np.arange(-10,50,5), edgecolor=\"None\", label=\"GR\")\n",
    "pl.hist(refsr3a[refsr3a>-10], bins=np.arange(-10,50,5), edgecolor=\"red\", facecolor=\"None\", label=\"SR\")\n",
    "pl.xlabel(\"Reflectivity (dBZ)\")\n",
    "pl.legend()\n",
    "fig.suptitle(\"corrected GR vs corrected SR hail\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = pl.figure(figsize=(12,8))\n",
    "ax = fig.add_subplot(121, aspect=\"equal\")\n",
    "pl.scatter(xyz[..., 0], xyz[...,1], c=refsr1a, cmap=pl.cm.jet, vmin=10, vmax=50, edgecolor=\"None\")\n",
    "pl.title(\"SR reflectivity\")\n",
    "pl.xlim(-100000, 150000)\n",
    "pl.ylim(-150000, 150000)\n",
    "pl.grid()\n",
    "ax = fig.add_subplot(122, aspect=\"equal\")\n",
    "pl.scatter(xyz[..., 0], xyz[...,1], c=refgr1a, cmap=pl.cm.jet, vmin=10, vmax=50, edgecolor=\"None\")\n",
    "pl.title(\"GR reflectivity\")\n",
    "pl.xlim(-100000, 150000)\n",
    "pl.ylim(-150000, 150000)\n",
    "pl.grid()\n",
    "fig.suptitle(\"uncorrected GR vs uncorrected SR\")\n",
    "pl.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Former Method using loop"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Lots of containers to store samples (only for one GR sweep angle!)\n",
    "x = np.zeros(nprof)*np.nan        # x coordinate of sample\n",
    "y = np.zeros(nprof)*np.nan        # y coordinate of sample\n",
    "z = np.zeros(nprof)*np.nan        # z coordinate of sample\n",
    "dz = np.zeros(nprof)*np.nan       # depth of sample\n",
    "ds = np.zeros(nprof)*np.nan       # width of sample\n",
    "rs = np.zeros(nprof)*np.nan       # range of sample from GR\n",
    "refsr1 = np.zeros(nprof)*np.nan     # SR reflectivity\n",
    "refsr2 = np.zeros(nprof)*np.nan     # SR reflectivity (S-band, snow)\n",
    "refsr3 = np.zeros(nprof)*np.nan     # SR reflectivity (S-band, hail)  \n",
    "refgr1 = np.zeros(nprof)*np.nan     # GR reflectivity\n",
    "refgr2 = np.zeros(nprof)*np.nan     # GR reflectivity (Ku-band)\n",
    "ntotpr = np.zeros(nprof,dtype=\"i4\")# total number of SR bins in sample\n",
    "nrej1 = np.zeros(nprof,dtype=\"i4\")# number of rejected SR bins in sample\n",
    "ntotgr = np.zeros(nprof,dtype=\"i4\")# total number of GR bins in sample\n",
    "nrej2 = np.zeros(nprof,dtype=\"i4\")# number of rejected GR bins in sample\n",
    "iref1 = np.zeros(nprof)*np.nan    # path-integrated SR reflectivity\n",
    "iref2 = np.zeros(nprof)*np.nan    # path-integrated GR reflectivity\n",
    "stdv1 = np.zeros(nprof)*np.nan    # std. dev. of SR reflectivity in sample\n",
    "stdv2 = np.zeros(nprof)*np.nan    # std. dev. of GR reflectivity in sample\n",
    "volsr = np.zeros(nprof)*np.nan     # total volume of SR bins in sample\n",
    "volgr = np.zeros(nprof)*np.nan     # total volume of GR bins in sample"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "# Loop over relevant SR profiles\n",
    "for ii, (ss, rr)  in enumerate(zip(vscan,vray)):\n",
    "    # Index and count valid bins in each profile\n",
    "    ip = np.where(valid[ss,rr])[0]\n",
    "    numbins = len(ip)\n",
    "    ntotpr[ii]=numbins\n",
    "    if numbins == 0:\n",
    "        continue\n",
    "    # Compute the mean position of these bins\n",
    "    x[ii]=np.mean(xyzp_sr[ss,rr,ip,0])\n",
    "    y[ii]=np.mean(xyzp_sr[ss,rr,ip,1])\n",
    "    z[ii]=np.mean(xyzp_sr[ss,rr,ip,2])\n",
    "    \n",
    "    # Thickness of the layer\n",
    "    dz[ii]=(numbins * dr_sr) / np.cos( np.radians(alpha[ss, rr]) )\n",
    "    \n",
    "    # SR averaging volume\n",
    "    volsr[ii]=np.sum(vol_sr[ss, rr, ip])\n",
    "\n",
    "    # Note mean TRMM beam diameter\n",
    "    ds[ii]=np.radians(bw_sr) * np.mean( ( (zt-z[ii]) / \n",
    "                                         np.cos( np.radians(alpha[ss, rr]) ) ) )\n",
    "    # Note distance from radar\n",
    "    s=np.sqrt(x[ii]**2+y[ii]**2)\n",
    "    rs[ii]=(re2+z[ii]) * np.sin(s/re2) / np.cos(np.radians(elev_gr))\n",
    "    \n",
    "    # This should not be required because we applied ZonalData\n",
    "    ### Check that sample is within radar range\n",
    "    ##if r[ii,jj]+ds[ii,jj]/2. gt rmax then continue\n",
    "\n",
    "    ## THIS IS THE ORIGINAL IDL CODE - IS THIS A BUG???\n",
    "    ##ref1[ii,jj]=MEAN(ref_sr1,/nan)\n",
    "    ##ref3[ii,jj]=MEAN(ref_sr2,/nan)\n",
    "    ##ref4[ii,jj]=MEAN(ref_sr3,/nan)\n",
    "    \n",
    "    # Simple linear average of reflectivity \n",
    "    #   - we can become fancier in the next step\n",
    "    # ATTENTION: NEED TO FLIP ARRAY\n",
    "    \n",
    "    refsr1[ii]=np.nanmean(ref_sr[ss,rr,ip])\n",
    "    refsr2[ii]=np.nanmean(ref_sr_ss[ss,rr,ip])\n",
    "    refsr3[ii]=np.nanmean(ref_sr_sh[ss,rr,ip])\n",
    "    \n",
    "    ## Not sure why we need this...\n",
    "    ### Note the number of rejected bins\n",
    "    ##nrej1[ii,jj]=ROUND(TOTAL(FINITE(ref_sr1,/nan)))\n",
    "    ##if FINITE(stdv1[ii,jj]) eq 0 and np-nrej1[ii,jj] gt 1 then STOP\n",
    "\n",
    "    # SHOULD WE USE ZONALDATA INSTEAD? COULD BE MORE ACCURATE, BUT ALSO SLOWER\n",
    "    # WE COULD BASICALLY START A NEW LOOP HERE AND RUN ZONALDATA BEFORE\n",
    "    \n",
    "    # Compute the horizontal distance to all the GR bins\n",
    "    d = np.sqrt((gr_xyz[...,0]-x[ii])**2 + (gr_xyz[...,1]-y[ii])**2)\n",
    "\n",
    "    # Find all GR bins within the SR beam\n",
    "    aa, bb = np.where(d <= ds[ii]/2.)\n",
    "\n",
    "    # Store the number of bins\n",
    "    ntotgr[ii] = len(aa)\n",
    "\n",
    "    if len(aa) == 0:\n",
    "        continue\n",
    "\n",
    "    # Extract the relevant GR bins\n",
    "\n",
    "    # Compute the GR averaging volume\n",
    "    volgr[ii]=np.sum(vol_gr[aa,bb])\n",
    "\n",
    "    # Average over those bins that exceed the reflectivity threshold \n",
    "    #   IDL code does exponential distance and volume weighting\n",
    "    #   Let's try simple mean first,\n",
    "    #   THEN ZonalStats!\n",
    "    refgr1[ii]=np.nanmean(ref_gr[aa,bb])\n",
    "    refgr2[ii]=np.nanmean(ref_gr_ku[aa,bb])    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = pl.figure(figsize=(12,5))\n",
    "ax = fig.add_subplot(121, aspect=\"equal\")\n",
    "pl.scatter(refgr1, refsr1, marker=\"+\", c=\"black\")\n",
    "pl.plot([0,60],[0,60], linestyle=\"solid\", color=\"black\")\n",
    "pl.xlim(10,50)\n",
    "pl.ylim(10,50)\n",
    "pl.xlabel(\"GR reflectivity (dBZ)\")\n",
    "pl.ylabel(\"SR reflectivity (dBZ)\")\n",
    "ax = fig.add_subplot(122)\n",
    "pl.hist(refgr1[refsr1>10], bins=np.arange(-10,50,5), edgecolor=\"None\", label=\"GR\")\n",
    "pl.hist(refsr1[refsr1>-10], bins=np.arange(-10,50,5), edgecolor=\"red\", facecolor=\"None\", label=\"SR\")\n",
    "pl.xlabel(\"Reflectivity (dBZ)\")\n",
    "pl.legend()\n",
    "fig.suptitle(\"uncorrected GR vs uncorrected SR\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = pl.figure(figsize=(12,5))\n",
    "ax = fig.add_subplot(121, aspect=\"equal\")\n",
    "pl.scatter(refgr2, refsr1, marker=\"+\", c=\"black\")\n",
    "pl.plot([0,60],[0,60], linestyle=\"solid\", color=\"black\")\n",
    "pl.xlim(10,50)\n",
    "pl.ylim(10,50)\n",
    "pl.xlabel(\"GR reflectivity (dBZ)\")\n",
    "pl.ylabel(\"SR reflectivity (dBZ)\")\n",
    "ax = fig.add_subplot(122)\n",
    "pl.hist(refgr1[refsr1>10], bins=np.arange(-10,50,5), edgecolor=\"None\", label=\"GR\")\n",
    "pl.hist(refsr1[refsr1>-10], bins=np.arange(-10,50,5), edgecolor=\"red\", facecolor=\"None\", label=\"SR\")\n",
    "pl.xlabel(\"Reflectivity (dBZ)\")\n",
    "pl.legend()\n",
    "fig.suptitle(\"corrected GR vs uncorrected SR\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = pl.figure(figsize=(12,5))\n",
    "ax = fig.add_subplot(121, aspect=\"equal\")\n",
    "pl.scatter(refgr2, refsr2, marker=\"+\", c=\"black\")\n",
    "pl.plot([0,60],[0,60], linestyle=\"solid\", color=\"black\")\n",
    "pl.xlim(10,50)\n",
    "pl.ylim(10,50)\n",
    "pl.xlabel(\"GR reflectivity (dBZ)\")\n",
    "pl.ylabel(\"SR reflectivity (dBZ)\")\n",
    "ax = fig.add_subplot(122)\n",
    "pl.hist(refgr2[refsr2>10], bins=np.arange(-10,50,5), edgecolor=\"None\", label=\"GR\")\n",
    "pl.hist(refsr2[refsr2>-10], bins=np.arange(-10,50,5), edgecolor=\"red\", facecolor=\"None\", label=\"SR\")\n",
    "pl.xlabel(\"Reflectivity (dBZ)\")\n",
    "pl.legend()\n",
    "fig.suptitle(\"corrected GR vs corrected SR snow\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = pl.figure(figsize=(12,5))\n",
    "ax = fig.add_subplot(121, aspect=\"equal\")\n",
    "pl.scatter(refgr2, refsr3, marker=\"+\", c=\"black\")\n",
    "pl.plot([0,60],[0,60], linestyle=\"solid\", color=\"black\")\n",
    "pl.xlim(10,50)\n",
    "pl.ylim(10,50)\n",
    "pl.xlabel(\"GR reflectivity (dBZ)\")\n",
    "pl.ylabel(\"SR reflectivity (dBZ)\")\n",
    "ax = fig.add_subplot(122)\n",
    "pl.hist(refgr2[refsr3>10], bins=np.arange(-10,50,5), edgecolor=\"None\", label=\"GR\")\n",
    "pl.hist(refsr3[refsr3>-10], bins=np.arange(-10,50,5), edgecolor=\"red\", facecolor=\"None\", label=\"SR\")\n",
    "pl.xlabel(\"Reflectivity (dBZ)\")\n",
    "pl.legend()\n",
    "fig.suptitle(\"corrected GR vs corrected SR hail\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = pl.figure(figsize=(12,8))\n",
    "ax = fig.add_subplot(121, aspect=\"equal\")\n",
    "pl.scatter(x, y, c=refsr1, cmap=pl.cm.jet, vmin=10, vmax=50, edgecolor=\"None\")\n",
    "pl.title(\"SR reflectivity\")\n",
    "pl.xlim(-100000, 150000)\n",
    "pl.ylim(-150000, 150000)\n",
    "pl.grid()\n",
    "ax = fig.add_subplot(122, aspect=\"equal\")\n",
    "pl.scatter(x, y, c=refgr1, cmap=pl.cm.jet, vmin=10, vmax=50, edgecolor=\"None\")\n",
    "pl.title(\"GR reflectivity\")\n",
    "pl.xlim(-100000, 150000)\n",
    "pl.ylim(-150000, 150000)\n",
    "pl.grid()\n",
    "fig.suptitle(\"uncorrected GR vs uncorrected SR\")\n",
    "pl.tight_layout()"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Edit Metadata",
  "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"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}