gpm_matching-test-australien.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "# Step By Step Matching verstehen!!!!\n",
    "# Alle Steps Plotten!\n",
    "# Warum entfallen bei steileren Elevationswinkel die Messungen ????!!!!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "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": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# define GPM data set\n",
    "#gpm_file_path = '/automount/ags/velibor/gpmdata/dprV7/2A.GPM.DPR.V7-20170308.20141007-S015721-E032951.003445.V05A.HDF5'\n",
    "#gpm_file_path = '/automount/ags/velibor/gpmdata/dprV7/2A.GPM.DPR.V7-20170308.20141008-S084747-E102017.003465.V05A.HDF5'\n",
    "#gpm_file = wrl.util.get_wradlib_data_file(gpm_file_path)\n",
    "\n",
    "\n",
    "\n",
    "#gpm_file = wrl.util.get_wradlib_data_file('/user/velibor/SHKGPM/Australien-gpm/volmatch/GPM_example/2A-RW-BRS.GPM.Ku.V6-20160118.20141206-S095002-E095137.004383.V04A.HDF5')\n",
    "gpm_file = wrl.util.get_wradlib_data_file('/automount/ags/velibor/gpmdata/dprV7/2A.GPM.DPR.V7-20170308.20141206-S083332-E100603.004383.V05A.HDF5')\n",
    "# define matching ground radar file\n",
    "#gr_file_path = '/automount/radar-archiv/scans/2014/2014-10/2014-10-08/n_ppi_280deg/2014-10-08--09:45:00,00.mvol'\n",
    "#gr2gpm_file = wrl.util.get_wradlib_data_file(gr_file_path)\n",
    "\n",
    "gr2gpm_file = wrl.util.get_wradlib_data_file('/user/velibor/SHKGPM/Australien-gpm/volmatch/GPM_example/IDR66_20141206_094829.vol.h5')\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "/user/velibor/SHKGPM/Australien-gpm/volmatch/GPM_example/IDR66_20141206_094829.vol.h5\n",
      "/automount/ags/velibor/gpmdata/dprV7/2A.GPM.DPR.V7-20170308.20141206-S083332-E100603.004383.V05A.HDF5\n"
     ]
    }
   ],
   "source": [
    "print (gr2gpm_file)\n",
    "print (gpm_file)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Space-born precipitation radar parameters\n",
    "sr_pars = {\"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": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Set parameters for this procedure\n",
    "bw_sr = 0.71                  # SR beam width\n",
    "platf = \"gpm\"                 # 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",
    "\n",
    "\n",
    "ee = 2     ###########?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/velibor/anaconda2/envs/wradlib/lib/python2.7/site-packages/IPython/core/formatters.py:92: DeprecationWarning: DisplayFormatter._ipython_display_formatter_default is deprecated: use @default decorator instead.\n",
      "  def _ipython_display_formatter_default(self):\n",
      "/home/velibor/anaconda2/envs/wradlib/lib/python2.7/site-packages/IPython/core/formatters.py:98: DeprecationWarning: DisplayFormatter._formatters_default is deprecated: use @default decorator instead.\n",
      "  def _formatters_default(self):\n",
      "/home/velibor/anaconda2/envs/wradlib/lib/python2.7/site-packages/IPython/core/formatters.py:677: DeprecationWarning: PlainTextFormatter._deferred_printers_default is deprecated: use @default decorator instead.\n",
      "  def _deferred_printers_default(self):\n",
      "/home/velibor/anaconda2/envs/wradlib/lib/python2.7/site-packages/IPython/core/formatters.py:669: DeprecationWarning: PlainTextFormatter._singleton_printers_default is deprecated: use @default decorator instead.\n",
      "  def _singleton_printers_default(self):\n",
      "/home/velibor/anaconda2/envs/wradlib/lib/python2.7/site-packages/IPython/core/formatters.py:672: DeprecationWarning: PlainTextFormatter._type_printers_default is deprecated: use @default decorator instead.\n",
      "  def _type_printers_default(self):\n",
      "/home/velibor/anaconda2/envs/wradlib/lib/python2.7/site-packages/IPython/core/formatters.py:669: DeprecationWarning: PlainTextFormatter._singleton_printers_default is deprecated: use @default decorator instead.\n",
      "  def _singleton_printers_default(self):\n",
      "/home/velibor/anaconda2/envs/wradlib/lib/python2.7/site-packages/IPython/core/formatters.py:672: DeprecationWarning: PlainTextFormatter._type_printers_default is deprecated: use @default decorator instead.\n",
      "  def _type_printers_default(self):\n",
      "/home/velibor/anaconda2/envs/wradlib/lib/python2.7/site-packages/IPython/core/formatters.py:677: DeprecationWarning: PlainTextFormatter._deferred_printers_default is deprecated: use @default decorator instead.\n",
      "  def _deferred_printers_default(self):\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "'/user/velibor/SHKGPM/Australien-gpm/volmatch/GPM_example/IDR66_20141206_094829.vol.h5'"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sr_pars['gpm']['gr_file']#['zt', 'dr', 'gr_file']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true
   },
   "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\n",
    "\n",
    "def read_gpm(filename, bbox):\n",
    "    pr_data = Dataset(filename, mode=\"r\")\n",
    "    lon = pr_data['NS'].variables['Longitude']\n",
    "    lat = pr_data['NS'].variables['Latitude']\n",
    "\n",
    "    poly = [[bbox['left'], bbox['bottom']],\n",
    "            [bbox['left'], bbox['top']],\n",
    "            [bbox['right'], bbox['top']],\n",
    "            [bbox['right'], bbox['bottom']],\n",
    "            [bbox['left'], bbox['bottom']]]\n",
    "    mask = get_clip_mask(np.dstack((lon[:], lat[:])), poly)\n",
    "    mask = np.nonzero(np.count_nonzero(mask, axis=1))\n",
    "    lon = lon[mask]\n",
    "    lat = lat[mask]\n",
    "\n",
    "    year = pr_data['NS']['ScanTime'].variables['Year'][mask]\n",
    "    month = pr_data['NS']['ScanTime'].variables['Month'][mask]\n",
    "    dayofmonth = pr_data['NS']['ScanTime'].variables['DayOfMonth'][mask]\n",
    "    # dayofyear = pr_data['NS']['ScanTime'].variables['DayOfYear'][mask]\n",
    "    hour = pr_data['NS']['ScanTime'].variables['Hour'][mask]\n",
    "    minute = pr_data['NS']['ScanTime'].variables['Minute'][mask]\n",
    "    second = pr_data['NS']['ScanTime'].variables['Second'][mask]\n",
    "    # secondofday = pr_data['NS']['ScanTime'].variables['SecondOfDay'][mask]\n",
    "    millisecond = pr_data['NS']['ScanTime'].variables['MilliSecond'][mask]\n",
    "    date_array = zip(year, month, dayofmonth,\n",
    "                     hour, minute, second,\n",
    "                     millisecond.astype(np.int32) * 1000)\n",
    "    pr_time = np.array(\n",
    "        [dt.datetime(d[0], d[1], d[2], d[3], d[4], d[5], d[6]) for d in\n",
    "         date_array])\n",
    "\n",
    "    sfc = pr_data['NS']['PRE'].variables['landSurfaceType'][mask]\n",
    "    pflag = pr_data['NS']['PRE'].variables['flagPrecip'][mask]\n",
    "\n",
    "    # bbflag = pr_data['NS']['CSF'].variables['flagBB'][mask]\n",
    "    zbb = pr_data['NS']['CSF'].variables['heightBB'][mask]\n",
    "    # print(zbb.dtype)\n",
    "    bbwidth = pr_data['NS']['CSF'].variables['widthBB'][mask]\n",
    "    qbb = pr_data['NS']['CSF'].variables['qualityBB'][mask]\n",
    "    qtype = pr_data['NS']['CSF'].variables['qualityTypePrecip'][mask]\n",
    "    ptype = pr_data['NS']['CSF'].variables['typePrecip'][mask]\n",
    "\n",
    "    quality = pr_data['NS']['scanStatus'].variables['dataQuality'][mask]\n",
    "    refl = pr_data['NS']['SLV'].variables['zFactorCorrected'][mask]\n",
    "    # print(pr_data['NS']['SLV'].variables['zFactorCorrected'])\n",
    "\n",
    "    zenith = pr_data['NS']['PRE'].variables['localZenithAngle'][mask]\n",
    "\n",
    "    pr_data.close()\n",
    "\n",
    "    # Check for bad data\n",
    "    if max(quality) != 0:\n",
    "        raise ValueError('GPM contains Bad Data')\n",
    "\n",
    "    pflag = pflag.astype(np.int8)\n",
    "\n",
    "    # Determine the dimensions\n",
    "    ndim = refl.ndim\n",
    "    if ndim != 3:\n",
    "        raise ValueError('GPM Dimensions do not match! '\n",
    "                         'Needed 3, given {0}'.format(ndim))\n",
    "\n",
    "    tmp = refl.shape\n",
    "    nscan = tmp[0]\n",
    "    nray = tmp[1]\n",
    "    nbin = tmp[2]\n",
    "\n",
    "    # Reverse direction along the beam\n",
    "    refl = np.flip(refl, axis=-1)\n",
    "\n",
    "    # Change pflag=1 to pflag=2 to be consistent with 'Rain certain' in TRMM\n",
    "    pflag[pflag == 1] = 2\n",
    "\n",
    "    # Simplify the precipitation types\n",
    "    ptype = (ptype / 1e7).astype(np.int16)\n",
    "\n",
    "    # Simplify the surface types\n",
    "    imiss = (sfc == -9999)\n",
    "    sfc = (sfc / 1e2).astype(np.int16) + 1\n",
    "    sfc[imiss] = 0\n",
    "\n",
    "    # Set a quality indicator for the BB and precip type data\n",
    "    # TODO: Why is the `quality` variable overwritten?\n",
    "\n",
    "    quality = np.zeros((nscan, nray), dtype=np.uint8)\n",
    "\n",
    "    i1 = ((qbb == 0) | (qbb == 1)) & (qtype == 1)\n",
    "    quality[i1] = 1\n",
    "\n",
    "    i2 = ((qbb > 1) | (qtype > 2))\n",
    "    quality[i2] = 2\n",
    "\n",
    "    gpm_data = {}\n",
    "    gpm_data.update({'nscan': nscan, 'nray': nray, 'nbin': nbin,\n",
    "                     'date': pr_time, 'lon': lon, 'lat': lat,\n",
    "                     'pflag': pflag, 'ptype': ptype, 'zbb': zbb,\n",
    "                     'bbwidth': bbwidth, 'sfc': sfc, 'quality': quality,\n",
    "                     'refl': refl, 'zenith': zenith})\n",
    "\n",
    "    return gpm_data\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def read_gr_boxpol(filename, loaddata=True):\n",
    "    \"\"\"\n",
    "    Read ground radar boxpol scans\n",
    "    \n",
    "    in arbeit....\n",
    "    zur zeit nur ein scan (1.5)\n",
    "    \"\"\"\n",
    "    sdate = []\n",
    "    refl = []\n",
    "    \n",
    "    gr_data = wrl.io.read_generic_netcdf(gr_file)\n",
    "    grdata, grattrs = wrl.io.read_GAMIC_hdf5(gr_file)\n",
    "    print gr_file\n",
    "    \n",
    "    dat_ = gr_data['what']['date']\n",
    "    print dat_\n",
    "    dat = dat_[0:10]\n",
    "    tim = dat_[11:-1]\n",
    "    print (dat, tim)\n",
    "    # '2014-10-07T02:42:02Z'\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",
    "    import glob\n",
    "    c_liste = glob.glob(gr_file.replace('n_ppi_010deg','*ppi*'))\n",
    "    ntilt = len(c_liste)\n",
    "    \n",
    "    ref_data = [[]]*ntilt\n",
    "\n",
    "\n",
    "    print ('ntilit:',ntilt)    \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 jj in range(ntilt):\n",
    "        print c_liste[jj]\n",
    "        \n",
    "        gr_data = wrl.io.read_generic_netcdf(c_liste[jj])\n",
    "        grdata, grattrs = wrl.io.read_GAMIC_hdf5(c_liste[jj])\n",
    "\n",
    "        #for i in range(0, ntilt):\n",
    "        dset = gr_data['scan0']\n",
    "        # Azimuth Start\n",
    "        a0[jj] = dset['how']['azi_start'] \n",
    "        # Elevations Winkel\n",
    "        elang[jj] = dset['how']['elevation']\n",
    "        # Anzahl bins\n",
    "        ngate[jj] = dset['how']['bin_count']\n",
    "        # Start und Step der range\n",
    "        r0[jj] = dset['how']['range_start']\n",
    "        dr[jj] = dset['how']['range_step'] * dset['how']['range_samples']\n",
    "        # Anzahl Rays\n",
    "        nbeam[jj] = dset['how']['ray_count']\n",
    "\n",
    "        print ('azi start: ', a0)    \n",
    "        print ('elevation: ', elang)    \n",
    "        print ('bin count: ', ngate)    \n",
    "        print ('range start: ', r0)    \n",
    "        print ('range step: ', dr)    \n",
    "        print ('ray count: ', nbeam)    \n",
    "\n",
    "\n",
    "        gr_dict = {}\n",
    "        gr_dict.update({'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",
    "\n",
    "\n",
    "        #for i in range(0, ntilt):\n",
    "            #dset = gr_data['scan0']\n",
    "            #data = dset['variables']['moment_10']['data']\n",
    "\n",
    "        data = grdata['SCAN0']['ZH']['data']\n",
    "        \n",
    "\n",
    "        #refl.append(data)\n",
    "        ref_data[jj]=data\n",
    "\n",
    "        #refl1 = np.array(refl)\n",
    "    gr_dict.update({'refl': ref_data})\n",
    "    \n",
    "    return gr_dict\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# read matching GR data\n",
    "gr_data = read_gr(gr_file)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "elang\n",
      "[ 0.5         0.89999998  1.29999995  1.79999995  2.4000001   3.0999999\n",
      "  4.19999981  5.5999999   7.4000001  10.         13.30000019 17.89999962\n",
      " 23.89999962 32.        ]\n",
      "------------------------- \n",
      " \n",
      " \n",
      "\n",
      "sdate\n",
      "[datetime.datetime(2014, 6, 12, 9, 48, 29)\n",
      " datetime.datetime(2014, 6, 12, 9, 49, 2)\n",
      " datetime.datetime(2014, 6, 12, 9, 49, 31)\n",
      " datetime.datetime(2014, 6, 12, 9, 49, 58)\n",
      " datetime.datetime(2014, 6, 12, 9, 50, 20)\n",
      " datetime.datetime(2014, 6, 12, 9, 50, 37)\n",
      " datetime.datetime(2014, 6, 12, 9, 50, 54)\n",
      " datetime.datetime(2014, 6, 12, 9, 51, 11)\n",
      " datetime.datetime(2014, 6, 12, 9, 51, 28)\n",
      " datetime.datetime(2014, 6, 12, 9, 51, 45)\n",
      " datetime.datetime(2014, 6, 12, 9, 52, 2)\n",
      " datetime.datetime(2014, 6, 12, 9, 52, 20)\n",
      " datetime.datetime(2014, 6, 12, 9, 52, 38)\n",
      " datetime.datetime(2014, 6, 12, 9, 52, 56)]\n",
      "------------------------- \n",
      " \n",
      " \n",
      "\n",
      "r0\n",
      "[0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]\n",
      "------------------------- \n",
      " \n",
      " \n",
      "\n",
      "source\n",
      "RAD:AU66,PLC:MtStapl\n",
      "------------------------- \n",
      " \n",
      " \n",
      "\n",
      "nbeam\n",
      "[360. 360. 360. 360. 360. 360. 360. 360. 360. 360. 360. 360. 360. 360.]\n",
      "------------------------- \n",
      " \n",
      " \n",
      "\n",
      "lon\n",
      "153.24000549316406\n",
      "------------------------- \n",
      " \n",
      " \n",
      "\n",
      "ntilt\n",
      "14\n",
      "------------------------- \n",
      " \n",
      " \n",
      "\n",
      "a0\n",
      "[-0.5 -0.5 -0.5 -0.5 -0.5 -0.5 -0.5 -0.5 -0.5 -0.5 -0.5 -0.5 -0.5 -0.5]\n",
      "------------------------- \n",
      " \n",
      " \n",
      "\n",
      "ngate\n",
      "[600 600 600 600 600 600 600 600 600 600 600 600 600 600]\n",
      "------------------------- \n",
      " \n",
      " \n",
      "\n",
      "refl\n",
      "[[[-32.  -32.  -32.  ... -32.  -32.  -32. ]\n",
      "  [-32.  -32.  -30.  ... -32.  -32.  -32. ]\n",
      "  [-32.  -32.   -9.  ... -32.  -32.  -32. ]\n",
      "  ...\n",
      "  [-32.  -32.  -32.  ...   3.    2.    2. ]\n",
      "  [-32.  -32.  -32.  ...  -3.  -32.  -32. ]\n",
      "  [-32.  -32.  -32.  ... -32.  -32.  -32. ]]\n",
      "\n",
      " [[-32.  -32.  -30.  ...  -1.5   0.   -1.5]\n",
      "  [-32.  -32.  -30.  ...  -1.5  -1.5  -3. ]\n",
      "  [-32.  -32.   -9.  ...  -1.5 -32.  -32. ]\n",
      "  ...\n",
      "  [-32.  -32.  -32.  ...   6.    6.    5. ]\n",
      "  [-32.  -32.  -32.  ...   3.   -1.5   2. ]\n",
      "  [-32.  -32.  -32.  ...   1.    0.    1. ]]\n",
      "\n",
      " [[-32.  -32.   -7.5 ...   5.    0.    5. ]\n",
      "  [-32.  -32.   -6.  ...   4.    1.    4. ]\n",
      "  [-32.  -32.   -6.  ...   2.    1.    2. ]\n",
      "  ...\n",
      "  [-32.  -32.   -9.  ...   6.    7.    8. ]\n",
      "  [-32.  -32.   -9.  ...   7.    6.    2. ]\n",
      "  [-32.  -32.  -30.  ...   6.    4.    4. ]]\n",
      "\n",
      " ...\n",
      "\n",
      " [[ 32.5 -30.    4.  ... -32.  -32.  -32. ]\n",
      "  [-32.  -30.    5.  ... -32.  -32.  -32. ]\n",
      "  [ 37.5 -30.    1.  ... -32.  -32.  -32. ]\n",
      "  ...\n",
      "  [-32.  -30.    0.  ... -32.  -32.  -32. ]\n",
      "  [-32.    7.    3.  ... -32.  -32.  -32. ]\n",
      "  [-32.  -30.   -3.  ... -32.  -32.  -32. ]]\n",
      "\n",
      " [[ 33.5 -30.    3.  ... -32.  -32.  -32. ]\n",
      "  [ 34.5 -30.    5.  ... -32.  -32.  -32. ]\n",
      "  [ 34.5 -30.   -3.  ... -32.  -32.  -32. ]\n",
      "  ...\n",
      "  [-32.  -30.    3.  ... -32.  -32.  -32. ]\n",
      "  [ 36.5 -30.    1.  ... -32.  -32.  -32. ]\n",
      "  [ 35.  -30.    3.  ... -32.  -32.  -32. ]]\n",
      "\n",
      " [[ 36.  -30.    7.  ... -32.  -32.  -32. ]\n",
      "  [-32.  -30.    1.  ... -32.  -32.  -32. ]\n",
      "  [ 36.5 -30.    4.  ... -32.  -32.  -32. ]\n",
      "  ...\n",
      "  [ 37.  -30.    5.  ... -32.  -32.  -32. ]\n",
      "  [ 37.  -30.    5.  ... -32.  -32.  -32. ]\n",
      "  [ 37.  -30.   11.  ... -32.  -32.  -32. ]]]\n",
      "------------------------- \n",
      " \n",
      " \n",
      "\n",
      "date\n",
      "2014-06-12 09:48:29\n",
      "------------------------- \n",
      " \n",
      " \n",
      "\n",
      "lat\n",
      "-27.71809959411621\n",
      "------------------------- \n",
      " \n",
      " \n",
      "\n",
      "alt\n",
      "174.99999701976776\n",
      "------------------------- \n",
      " \n",
      " \n",
      "\n",
      "dr\n",
      "[250. 250. 250. 250. 250. 250. 250. 250. 250. 250. 250. 250. 250. 250.]\n",
      "------------------------- \n",
      " \n",
      " \n",
      "\n"
     ]
    }
   ],
   "source": [
    "\n",
    "for i in range(len(gr_data.keys())):\n",
    "    print (gr_data.keys()[i])\n",
    "    print(gr_data[gr_data.keys()[i]])\n",
    "    print ('------------------------- \\n \\n \\n')\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "((360, 600),)\n"
     ]
    }
   ],
   "source": [
    "######### GR data and att\n",
    "ee = 11\n",
    "# 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(ref_gr.shape, )\n",
    "\n",
    "\n",
    "##############\n",
    "##############\n",
    "### synthetic#\n",
    "##############\n",
    "##############\n",
    "#ref_gr = np.ones(ref_gr.shape)*40"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Radar bounding box:\n",
      "\t-26.44\n",
      "151.80           154.68\n",
      "\t-29.00\n"
     ]
    }
   ],
   "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": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "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\")\n",
    "sr_data = wrl.io.read_gpm(gpm_file, bbox)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(77, 49, 176)\n"
     ]
    }
   ],
   "source": [
    "pl.subplot(1,2,1)\n",
    "refl = sr_data['refl']\n",
    "print refl.shape\n",
    "pl.imshow(refl[:,:,12])\n",
    "\n",
    "pl.subplot(1,2,2)\n",
    "#print (gr_data['refl'][0,:,:].shape)\n",
    "pl.imshow(gr_data['refl'][0])\n",
    "pl.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "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",
    "\n",
    "##############\n",
    "##############\n",
    "### synthetic#\n",
    "##############\n",
    "##############\n",
    "pflag = np.ones(pflag.shape)*2\n",
    "\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": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "((77, 49), (77, 49), (77, 49))\n",
      "(77, 49, 176)\n",
      "['nscan', 'pflag', 'zbb', 'zenith', 'nray', 'lon', 'ptype', 'lat', 'refl', 'date', 'nbin', 'quality', 'sfc', 'bbwidth']\n"
     ]
    }
   ],
   "source": [
    "\n",
    "print (sr_lon.shape, sr_lat.shape, pflag.shape, )\n",
    "print (nscan_sr, nray_sr, ngate_sr)\n",
    "print (sr_data.keys())\n",
    "#print (sr_data['refl'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "('Earth radius 1:', 6373541.039814419)\n",
      "('SemiMajor, SemiMinor:', 6378137.0, 6356752.314245179)\n",
      "('Earth radius 2:', 6373541.039814381)\n"
     ]
    }
   ],
   "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": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "('Range/Azi-Shape:', (600,), (360,))\n",
      "('XYZ-Grid-Shape:', (360, 600, 3))\n",
      "((216000, 5, 3), 216000)\n",
      "('Domain-Shape:', (361, 2))\n"
     ]
    }
   ],
   "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": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "('SR-GRID-Shapes:', (77, 49), (77, 49), (77, 49, 2))\n"
     ]
    }
   ],
   "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": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[2.]\n"
     ]
    }
   ],
   "source": [
    "pl.imshow(pflag)\n",
    "pl.colorbar()\n",
    "print (np.unique(pflag))\n",
    "#pl.plot(pflag)\n",
    "pl.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(77, 49, 2)\n",
      "(77, 49)\n",
      "((77, 49), (77, 49, 2))\n",
      "((77, 49), (3773, 2))\n",
      "[False  True]\n"
     ]
    }
   ],
   "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",
    "####################!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n",
    "#########!!!!!!!!!!!!!!!!!!!! eigentlich == 2 ===?????\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)\n",
    "\n",
    "#pl.colorbar()\n",
    "print (np.unique(precip_mask))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Parallax Correction\n",
    "---------------------------"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "('NRAY', 49)\n",
      "('NBIN', 176)\n",
      "((77, 49, 176, 2), (176,), (77, 49, 176))\n",
      "(77, 49, 2)\n",
      "('SR_XYP:', (77, 49, 176, 2), (77, 49, 176, 3), (176,), (77, 49, 176))\n"
     ]
    }
   ],
   "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": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "pl.subplot(1,2,1)\n",
    "pl.imshow(xyp_sr[:,:,160,0])\n",
    "pl.subplot(1,2,2)\n",
    "pl.imshow(xyzp_sr[:,:,160,2])\n",
    "pl.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Compute spherical coordinates of SR bins with regard to GR\n",
    "-------------------------------"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "collapsed": false
   },
   "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[9,:,:].T))\n",
    "mask.shape\n",
    "pl.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Compute SR and GR pulse volumes\n",
    "------------------------------------\n",
    "\n",
    "## Calculate distance from orbit $r_s$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "rs = wrl.georef.dist_from_orbit(zt, alpha, beta, r_sr_inv, re1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# SR pulse volume\n",
    "--------------------------------"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "('Min. difference (m3):', 57223.422050476074)\n",
      "('Max. difference (m3): ', 70339.55933856964)\n",
      "('Average rel. difference (%):', 0.0026)\n",
      "(77, 49, 176)\n"
     ]
    }
   ],
   "source": [
    "\n",
    "# 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",
    "\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!\n",
    "\n",
    "print (vol_sr.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# GR Pulse Volume\n",
    "----------------------------------"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(360, 600)\n"
     ]
    }
   ],
   "source": [
    "\n",
    "\n",
    "# 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\")\n",
    "\n",
    "print (vol_gr.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Calculate horizontal and vertical dimensions $R_S$ and $D_S$ of SR bin\n",
    "---------------------------------"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "((77, 49, 176), (77, 49, 176))\n"
     ]
    }
   ],
   "source": [
    "\n",
    "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)\n",
    "\n",
    "print (Rs.shape, Ds.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Median Brightband Width/Height\n",
    "----------------------------"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(77, 49, 176)\n",
      "((77, 49), (77, 49), (77, 49))\n",
      "(77, 49, 176)\n",
      "(-4.195068728630307, -4.102380785848437)\n"
     ]
    }
   ],
   "source": [
    "\n",
    "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",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Convert SR Ku in S-Band Ref and vice versa\n",
    "----------------------------------"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "((77, 49, 176), (77, 49, 176))\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "'\\nHier ein eigenes Convert zwischen KU und X BAND\\n'"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    " ##############\n",
    "##############\n",
    "### synthetic#\n",
    "##############\n",
    "##############\n",
    "#ref_sr = sr_data['refl'].filled(np.nan)\n",
    "ref_sr1 = np.ma.getdata(sr_data['refl'].filled(1))\n",
    "ref_sr = np.ones(ref_sr1.shape)*30\n",
    "\n",
    "\n",
    "print(ref_sr.shape, ratio.shape)\n",
    "\n",
    "'''\n",
    "Hier ein eigenes Convert zwischen KU und X BAND\n",
    "'''"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Matching SR / GR\n",
    "------------------"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "((77, 49, 176), (77, 49))\n",
      "(609,)\n",
      "(77, 49, 176)\n",
      "(77, 49)\n"
     ]
    }
   ],
   "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",
    "\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",
    "\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",
    "\n",
    "# Number of matching SR bins per profile\n",
    "nvalids = np.sum(valid, axis=2)\n",
    "\n",
    "# scan and ray indices for profiles with at least one valid bin\n",
    "vscan, vray = np.where(nvalids>0)\n",
    "\n",
    "# number of profiles with at least one valid bin\n",
    "nprof = len(vscan)\n",
    "print(vscan.shape)\n",
    "print(valid.shape)\n",
    "print(nvalids.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(77, 49, 176, 3)\n",
      "((609, 3), (77, 49))\n"
     ]
    }
   ],
   "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",
    "pppl = ax.scatter(xyz_m1[..., 0], xyz_m1[..., 1], xyz_m1[..., 2], c=c.ravel())\n",
    "fig.colorbar(pppl, shrink=0.5, aspect=5)\n",
    "ax.set_zlim(0, 200000)\n",
    "pl.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Calculate spherical coordinates of SR sample volume with respect to GR\n",
    "--------------------------------------------"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {
    "collapsed": false
   },
   "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\n",
    "------------------------------------------"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(609,)\n"
     ]
    }
   ],
   "source": [
    "\n",
    "\n",
    "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)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Identify which GR bins actually intersect with the SR footprint\n",
    "---------------------\n",
    "## Create Zonal Data Source of matching SR profile"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "609\n"
     ]
    }
   ],
   "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\n",
    "--------------------------"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {
    "collapsed": false
   },
   "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\n",
    "------------------------"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "creating\n",
      "CPU times: user 50.6 s, sys: 371 ms, total: 50.9 s\n",
      "Wall time: 50.8 s\n"
     ]
    }
   ],
   "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": 38,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "loading\n",
      "CPU times: user 1min 19s, sys: 441 ms, total: 1min 20s\n",
      "Wall time: 1min 20s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "print(\"loading\")\n",
    "#obj3 = wrl.zonalstats.GridCellsToPoly('m3d_zonal_poly_{0}'.format(platf))\n",
    "obj3 = wrl.zonalstats.GridCellsToPoly(zdp)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Average GR volume and reflectivity\n",
    "---------------------------"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(609,)\n",
      "(609,)\n",
      "CPU times: user 546 ms, sys: 11.8 ms, total: 558 ms\n",
      "Wall time: 554 ms\n"
     ]
    }
   ],
   "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": "code",
   "execution_count": 40,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "##### betrachten von noederschlags relevanten Reflek und Threshold DPR Ku = 15 DBZ\n",
    "from satlib import corcor"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Threshold festlegen\n",
    "---------------------------------"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "TH = 15\n",
    "refsr1a[refsr1a<=TH]=np.nan\n",
    "refgr1a[refgr1a<=TH]=np.nan\n",
    "ref_gr[ref_gr<=TH]=np.nan\n",
    "ref_sr[ref_sr<=TH]=np.nan"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# PLOT\n",
    "______________"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "fig = pl.figure(figsize=(14,14))\n",
    "ax = fig.add_subplot(221, aspect=\"equal\")\n",
    "pl.scatter(refgr1a, refsr1a, marker=\"+\", c=\"black\")\n",
    "pl.plot([10,50],[10,50], linestyle=\"solid\", color=\"black\")\n",
    "#pl.xlim(10,50)\n",
    "#pl.ylim(10,50)\n",
    "pl.grid()\n",
    "pl.xlabel(\"GR reflectivity (dBZ)\")\n",
    "pl.ylabel(\"SR reflectivity (dBZ)\")\n",
    "pl.title(corcor(refgr1a, refsr1a))\n",
    "\n",
    "ax = fig.add_subplot(222)\n",
    "m1 = ~np.isnan(refgr1a)& ~np.isnan(refsr1a)\n",
    "\n",
    "pl.hist(refgr1a[m1], bins=20, edgecolor=\"blue\", facecolor=\"grey\", label=\"GR\")\n",
    "pl.hist(refsr1a[m1], bins=20, edgecolor=\"red\", facecolor=\"None\", label=\"SR\", lw=3)\n",
    "pl.xlabel(\"Reflectivity (dBZ)\")\n",
    "pl.legend()\n",
    "pl.grid()\n",
    "pl.title('Matched SR and GR Ref')\n",
    "pl.xlim(10,50)\n",
    "\n",
    "ax = fig.add_subplot(224)\n",
    "r_gr = ref_gr.reshape(ref_gr.shape[0]*ref_gr.shape[1])\n",
    "r_sr = ref_sr.reshape(ref_sr.shape[0]*ref_sr.shape[1]*ref_sr.shape[2])\n",
    "m1, m2 = ~np.isnan(r_gr), ~np.isnan(r_sr)\n",
    "\n",
    "pl.hist(r_gr[m1], bins=20, edgecolor=\"blue\", facecolor=\"grey\", label=\"GR\", normed=1)\n",
    "pl.hist(r_sr[m2], bins=20, edgecolor=\"red\", facecolor=\"None\", label=\"SR\", lw=3, normed=1)\n",
    "pl.xlabel(\"Reflectivity (dBZ)\")\n",
    "pl.legend()\n",
    "pl.grid()\n",
    "pl.title('Normed Raw SR and GR Ref')\n",
    "pl.xlim(10,50)\n",
    "#fig.suptitle(\"uncorrected GR vs uncorrected SR\")\n",
    "\n",
    "pl.tight_layout()\n",
    "pl.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'fig = pl.figure(figsize=(12,8))\\nax = fig.add_subplot(131, aspect=\"equal\")\\npl.scatter(xyz[..., 0], xyz[...,1], c=refsr1a, cmap=pl.cm.jet, vmin=TH, vmax=50, edgecolor=\"None\")\\npl.title(\"SR reflectivity\")\\npl.xlim(-50000, 50000)\\npl.ylim(-50000, 50000)\\npl.grid()\\npl.colorbar(shrink=0.5)\\n\\nax = fig.add_subplot(132, aspect=\"equal\")\\npl.scatter(xyz[..., 0], xyz[...,1], c=refgr1a, cmap=pl.cm.jet, vmin=TH, vmax=50, edgecolor=\"None\")\\npl.title(\"GR reflectivity\")\\npl.xlim(-50000, 50000)\\npl.ylim(-50000, 50000)\\npl.grid()\\npl.colorbar(shrink=0.5)\\nfig.suptitle(\"uncorrected GR vs uncorrected SR\")\\npl.tight_layout()\\n\\nax = fig.add_subplot(133, aspect=\"equal\")\\npl.scatter(xyz[..., 0], xyz[...,1], c=refsr1a-refgr1a,vmin=-5,vmax=5 ,edgecolor=\"None\")\\npl.title(\"GR reflectivity\")\\npl.xlim(-50000, 50000)\\npl.ylim(-50000, 50000)\\npl.grid()\\npl.colorbar(shrink=0.5)\\nfig.suptitle(\"uncorrected SR - GR\")\\npl.tight_layout()\\npl.show()'"
      ]
     },
     "execution_count": 43,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\"\"\"fig = pl.figure(figsize=(12,8))\n",
    "ax = fig.add_subplot(131, aspect=\"equal\")\n",
    "pl.scatter(xyz[..., 0], xyz[...,1], c=refsr1a, cmap=pl.cm.jet, vmin=TH, vmax=50, edgecolor=\"None\")\n",
    "pl.title(\"SR reflectivity\")\n",
    "pl.xlim(-50000, 50000)\n",
    "pl.ylim(-50000, 50000)\n",
    "pl.grid()\n",
    "pl.colorbar(shrink=0.5)\n",
    "\n",
    "ax = fig.add_subplot(132, aspect=\"equal\")\n",
    "pl.scatter(xyz[..., 0], xyz[...,1], c=refgr1a, cmap=pl.cm.jet, vmin=TH, vmax=50, edgecolor=\"None\")\n",
    "pl.title(\"GR reflectivity\")\n",
    "pl.xlim(-50000, 50000)\n",
    "pl.ylim(-50000, 50000)\n",
    "pl.grid()\n",
    "pl.colorbar(shrink=0.5)\n",
    "fig.suptitle(\"uncorrected GR vs uncorrected SR\")\n",
    "pl.tight_layout()\n",
    "\n",
    "ax = fig.add_subplot(133, aspect=\"equal\")\n",
    "pl.scatter(xyz[..., 0], xyz[...,1], c=refsr1a-refgr1a,vmin=-5,vmax=5 ,edgecolor=\"None\")\n",
    "pl.title(\"GR reflectivity\")\n",
    "pl.xlim(-50000, 50000)\n",
    "pl.ylim(-50000, 50000)\n",
    "pl.grid()\n",
    "pl.colorbar(shrink=0.5)\n",
    "fig.suptitle(\"uncorrected SR - GR\")\n",
    "pl.tight_layout()\n",
    "pl.show()\"\"\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "fig = pl.figure(figsize=(12,8))\n",
    "ax = fig.add_subplot(131, aspect=\"auto\")\n",
    "pl.scatter(xyz[..., 1], xyz[...,2], c=refsr1a, cmap=pl.cm.jet, vmin=TH, vmax=50, edgecolor=\"None\")\n",
    "pl.title(\"SR reflectivity\")\n",
    "pl.xlim(np.nanmin(xyz[...,1]),np.nanmax(xyz[...,1]))\n",
    "pl.ylim(0,np.nanmax(xyz[...,2]))\n",
    "pl.grid()\n",
    "pl.colorbar(shrink=0.5)\n",
    "\n",
    "ax = fig.add_subplot(132, aspect=\"auto\")\n",
    "pl.scatter(xyz[..., 1], xyz[...,2], c=refgr1a, cmap=pl.cm.jet, vmin=TH, vmax=50, edgecolor=\"None\")\n",
    "pl.title(\"GR reflectivity\")\n",
    "pl.xlim(np.nanmin(xyz[...,1]),np.nanmax(xyz[...,1]))\n",
    "pl.ylim(0,np.nanmax(xyz[...,2]))\n",
    "pl.grid()\n",
    "pl.colorbar(shrink=0.5)\n",
    "#fig.suptitle(\"uncorrected GR vs uncorrected SR\")\n",
    "pl.tight_layout()\n",
    "\n",
    "ax = fig.add_subplot(133, aspect=\"auto\")\n",
    "pl.scatter(xyz[..., 1], xyz[...,2], c=refsr1a-refgr1a, cmap='seismic',vmin=-10,vmax=10 ,edgecolor=\"None\")\n",
    "pl.title(\"SR - GR reflectivity\")\n",
    "pl.xlim(np.nanmin(xyz[...,1]),np.nanmax(xyz[...,1]))\n",
    "pl.ylim(0,np.nanmax(xyz[...,2]))\n",
    "pl.grid()\n",
    "pl.colorbar(shrink=0.5)\n",
    "#fig.suptitle(\"uncorrected SR - GR\")\n",
    "pl.tight_layout()\n",
    "pl.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "# ????? Wie groß ist der Fehler in unseren Breiten zwischen 3d Korrigier und nicht korrigiert ???"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "('Corellation: ', 'nan')\n"
     ]
    }
   ],
   "source": [
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "fig = pl.figure(figsize=(18,6))\n",
    "ax1 = fig.add_subplot(131, projection='3d')\n",
    "pppl = ax1.scatter(xyz[..., 0], xyz[..., 1], xyz[..., 2], c=refsr1a.ravel(), vmin=TH, vmax=50)\n",
    "fig.colorbar(pppl, shrink=0.5, aspect=5)\n",
    "ax1.set_title('Ref SR')\n",
    "\n",
    "ax2 = fig.add_subplot(132, projection='3d')\n",
    "pppl = ax2.scatter(xyz[..., 0], xyz[..., 1], xyz[..., 2], c=refgr1a.ravel(), vmin=TH, vmax=50)\n",
    "fig.colorbar(pppl, shrink=0.5, aspect=5)\n",
    "ax2.set_title('Ref GR')\n",
    "\n",
    "ax3 = fig.add_subplot(133, projection='3d')\n",
    "pppl = ax3.scatter(xyz[..., 0], xyz[..., 1], xyz[..., 2], c=refsr1a.ravel()-refgr1a.ravel(),\n",
    "                   cmap='seismic',vmin=-10,vmax=10)\n",
    "fig.colorbar(pppl, shrink=0.5, aspect=5)\n",
    "ax3.set_title('diffRef SR-GR')\n",
    "\n",
    "from satlib import corcor\n",
    "print ('Corellation: ',corcor( refgr1a, refsr1a))\n",
    "pl.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'ax = fig.add_subplot(111)\\nGRi, SRi = refgr1a, refsr1a\\nmmi = ~np.isnan(GRi)& ~np.isnan(SRi)\\n\\npl.hist(SRi[mmi]-GRi[mmi],bins=50, facecolor=\"grey\", alpha=0.6)\\npl.axvline(np.mean(SRi[mmi]-GRi[mmi]), color=\\'blue\\', label=\\'Mean diff\\')\\npl.axvline(np.median(SRi[mmi]-GRi[mmi]), color=\\'green\\', label=\\'Median diff\\')\\npl.axvline(0, color=\\'black\\',ls=\\'--\\',lw=3)\\n\\nxmin, xmax = np.nanmin(SRi[mmi]-GRi[mmi]), np.nanmax(SRi[mmi]-GRi[mmi])\\nlimits = np.nanmax([abs(xmin),abs(xmax)])\\npl.xlim(-limits,limits)\\n\\npl.title(\\'Diff SR-GR, bias: \\'+str( np.sum(SRi[mmi]-GRi[mmi])/len(SRi[mmi]-GRi[mmi]))+\\' dbz\\')\\npl.grid()\\npl.legend()\\npl.xlabel(\\'GR overestimate----------------------------Diff in dbz----------------------------SR overestimate\\')\\npl.ylabel(\\'number\\')'"
      ]
     },
     "execution_count": 46,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\"\"\"ax = fig.add_subplot(111)\n",
    "GRi, SRi = refgr1a, refsr1a\n",
    "mmi = ~np.isnan(GRi)& ~np.isnan(SRi)\n",
    "\n",
    "pl.hist(SRi[mmi]-GRi[mmi],bins=50, facecolor=\"grey\", alpha=0.6)\n",
    "pl.axvline(np.mean(SRi[mmi]-GRi[mmi]), color='blue', label='Mean diff')\n",
    "pl.axvline(np.median(SRi[mmi]-GRi[mmi]), color='green', label='Median diff')\n",
    "pl.axvline(0, color='black',ls='--',lw=3)\n",
    "\n",
    "xmin, xmax = np.nanmin(SRi[mmi]-GRi[mmi]), np.nanmax(SRi[mmi]-GRi[mmi])\n",
    "limits = np.nanmax([abs(xmin),abs(xmax)])\n",
    "pl.xlim(-limits,limits)\n",
    "\n",
    "pl.title('Diff SR-GR, bias: '+str( np.sum(SRi[mmi]-GRi[mmi])/len(SRi[mmi]-GRi[mmi]))+' dbz')\n",
    "pl.grid()\n",
    "pl.legend()\n",
    "pl.xlabel('GR overestimate----------------------------Diff in dbz----------------------------SR overestimate')\n",
    "pl.ylabel('number')\"\"\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.text.Text at 0x7f63cc934b50>"
      ]
     },
     "execution_count": 47,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pl.figure(figsize=(14,14))\n",
    "pl.scatter(refgr1a, refsr1a, c=xyz[...,2],marker=\"+\",s=40 )\n",
    "pl.plot([10,50],[10,50], linestyle=\"solid\", color=\"black\")\n",
    "pl.colorbar()\n",
    "#pl.xlim(10,50)\n",
    "#pl.ylim(10,50)\n",
    "pl.grid()\n",
    "pl.xlabel(\"GR reflectivity (dBZ)\")\n",
    "pl.ylabel(\"SR reflectivity (dBZ)\")\n",
    "pl.title(corcor(refgr1a, refsr1a))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(609,)"
      ]
     },
     "execution_count": 48,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "refgr1a.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(609, 3)"
      ]
     },
     "execution_count": 49,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "xyz.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "112.0"
      ]
     },
     "execution_count": 50,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "2.8*40"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.0016666666666666668"
      ]
     },
     "execution_count": 51,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "1./600."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'pl.pcolormesh(sr_data[\"lon\"], sr_data[\"lat\"], ref_sr[:,:,10], vmin=0, vmax=40, cmap=\\'jet\\')\\npl.colorbar()\\npl.figure()\\n\\nax, pm = wrl.vis.plot_ppi(ref_gr, vmin=0, vmax=40, cmap=\\'jet\\')\\ncbar = pl.colorbar(pm, shrink=0.75)\\ncbar.set_label(\"Reflectivity (dBZ)\")\\n\\npl.figure(figsize=(10,5))\\npl.subplot(1,2,1)\\npl.scatter(xyz[...,0], xyz[...,1], c=refsr1a, label=str(elev_gr),vmin=0, vmax=40, s=40)\\npl.colorbar();pl.xlim(-150000,150000);pl.ylim(-150000,150000)\\npl.subplot(1,2,2)\\npl.scatter(xyz[...,0], xyz[...,1], c=refgr1a, label=str(elev_gr),vmin=0, vmax=40, s=40)\\npl.colorbar();pl.xlim(-150000,150000);pl.ylim(-150000,150000)\\npl.show()\\n        '"
      ]
     },
     "execution_count": 52,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\"\"\"pl.pcolormesh(sr_data[\"lon\"], sr_data[\"lat\"], ref_sr[:,:,10], vmin=0, vmax=40, cmap='jet')\n",
    "pl.colorbar()\n",
    "pl.figure()\n",
    "\n",
    "ax, pm = wrl.vis.plot_ppi(ref_gr, vmin=0, vmax=40, cmap='jet')\n",
    "cbar = pl.colorbar(pm, shrink=0.75)\n",
    "cbar.set_label(\"Reflectivity (dBZ)\")\n",
    "\n",
    "pl.figure(figsize=(10,5))\n",
    "pl.subplot(1,2,1)\n",
    "pl.scatter(xyz[...,0], xyz[...,1], c=refsr1a, label=str(elev_gr),vmin=0, vmax=40, s=40)\n",
    "pl.colorbar();pl.xlim(-150000,150000);pl.ylim(-150000,150000)\n",
    "pl.subplot(1,2,2)\n",
    "pl.scatter(xyz[...,0], xyz[...,1], c=refgr1a, label=str(elev_gr),vmin=0, vmax=40, s=40)\n",
    "pl.colorbar();pl.xlim(-150000,150000);pl.ylim(-150000,150000)\n",
    "pl.show()\n",
    "        \"\"\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'pl.xlim(0,40000)\\npl.ylim(0,400)'"
      ]
     },
     "execution_count": 53,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pl.scatter(gr_xyz[:,:,0], gr_xyz[:,:,1], color='green')\n",
    "pl.scatter(sr_x,sr_y)\n",
    "\n",
    "pl.scatter(xyz[...,0], xyz[...,1],color='red')\n",
    "\n",
    "#pl.scatter(gr_polargrid[0],gr_polargrid[1])\n",
    "#pl.scatter(sr_x,sr_y, color='green')\n",
    "pl.xlim(-10000,10000)\n",
    "pl.ylim(-10000,10000)\n",
    "\n",
    "\n",
    "\"\"\"pl.xlim(0,40000)\n",
    "pl.ylim(0,400)\"\"\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "fig = pl.figure()\n",
    "ax = fig.add_subplot(111, projection='3d')\n",
    "pppl = ax.scatter(xyz_m1[..., 0], xyz_m1[..., 1], xyz_m1[..., 2])\n",
    "#pl.savefig('/automount/ftp/velibor/validation/MATCHINGPROB1.pdf')\n",
    "pl.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 70,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "fig = pl.figure(figsize=(10,10))\n",
    "ax = fig.add_subplot(111, projection='3d')\n",
    "pppl = ax.scatter(gr_xyz[:,:,0], gr_xyz[:,:,1], gr_xyz[:,:,2],c=ref_gr, edgecolors='none')\n",
    "pl.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 69,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "((609, 3), (609,))\n"
     ]
    }
   ],
   "source": [
    "fig = pl.figure(figsize=(10,10))\n",
    "ax = fig.add_subplot(111, projection='3d')\n",
    "print(xyz.shape, refsr1a.shape)\n",
    "pppl = ax.scatter(xyz[...,0],xyz[...,1],xyz[...,2], edgecolors='none', s=100, color='black')\n",
    "\n",
    "pl.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "((609, 3), (609,))\n"
     ]
    }
   ],
   "source": [
    "fig = pl.figure(figsize=(10,10))\n",
    "ax = fig.add_subplot(111, projection='3d')\n",
    "print(xyz.shape, refsr1a.shape)\n",
    "pppl = ax.scatter(gr_xyz[:,:,0], gr_xyz[:,:,1], gr_xyz[:,:,2],c=ref_gr, edgecolors='none')\n",
    "pppl = ax.scatter(xyz[...,0],xyz[...,1],xyz[...,2], edgecolors='none', s=100, c='black')\n",
    "pppl = ax.plot_wireframe(xyz[...,0],xyz[...,1],xyz[...,2])\n",
    "\n",
    "ax.set_xlim(-30000, 30000)\n",
    "ax.set_ylim(-30000, 30000)\n",
    "ax.set_zlim(0, 7000)\n",
    "pl.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "pl.figure(figsize=(10,5))\n",
    "pl.subplot(1,2,1)\n",
    "pl.scatter(xyz[...,0], xyz[...,2], c=refsr1a, label=str(elev_gr),vmin=0, vmax=40, s=40)\n",
    "pl.colorbar();pl.xlim(-150000,150000);pl.ylim(0,7000)\n",
    "pl.subplot(1,2,2)\n",
    "pl.scatter(xyz[...,0], xyz[...,2], c=refgr1a, label=str(elev_gr),vmin=0, vmax=40, s=40)\n",
    "pl.colorbar();pl.xlim(-150000,150000);pl.ylim(0,7000)\n",
    "pl.show()\n",
    "        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 76,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "fig = pl.figure(figsize=(10,10))\n",
    "ax = fig.add_subplot(111, projection='3d')\n",
    "pppl = ax.plot_surface(gr_xyz[:,:,0], gr_xyz[:,:,1], gr_xyz[:,:,2])\n",
    "pl.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}