pcc_gpm_matching_all_elevation_TEST_onlyGR.ipynb

{
 "cells": [
  {
   "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\n",
    "warnings.filterwarnings('ignore')\n",
    "from satlib import corcor\n",
    "import os\n",
    "import time\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def subplot_cfad(x1, x2, ZP):\n",
    "    \"\"\"\n",
    "    Plot 4 CFAD\n",
    "    \"\"\"\n",
    "    import matplotlib.pyplot as plt\n",
    "    m1 = ~np.isnan(x1) & ~np.isnan(x2)\n",
    "    \n",
    "    x3 = np.array([]); x4 = np.array([]); x5 = np.array([np.nan])\n",
    "\n",
    "    hh = np.array([])\n",
    "    for jj in range(0,35000,100):\n",
    "        hh = np.append(hh,jj)\n",
    "        x3 = np.append(x3,np.nanmean(x1[np.where((x2>jj)&(x2<jj+100))]))\n",
    "        x4 = np.append(x4,np.nanmedian(x1[np.where((x2>jj)&(x2<jj+100))]))\n",
    "        x0 = x1.copy()\n",
    "        if x0[np.where((x2>jj)&(x2<jj+100))].size==0:\n",
    "            x5 = np.append(x5,0)\n",
    "        else:\n",
    "            x5 = np.append(x5,np.nanmax(x0[np.where((x2>jj)&(x2<jj+100))]))\n",
    "    \n",
    "    print hh.shape, x5.shape, x4.shape\n",
    "    plt.hist2d(x1[m1], x2[m1],bins=30, cmap='inferno', vmax=40)\n",
    "    #plt.hexbin(x1[m1], x2[m1],bins=10, cmap='inferno')\n",
    "    plt.colorbar()\n",
    "    plt.xlabel('ZH in DBz')\n",
    "    plt.ylabel('Height in m')\n",
    "    plt.title('CFAD: '+ ZP)\n",
    "    plt.plot(x3,hh, linewidth=2, label='mean')\n",
    "    plt.plot(x4,hh, linewidth=2, label='median')\n",
    "    #plt.plot(x5[1::],hh, linewidth=2, label='max', color='white')\n",
    "    # BB detection with cfad?\n",
    "    #plt.axhline(0, color='white')\n",
    "\n",
    "    plt.legend(loc='upper right', fontsize=10)\n",
    "    plt.grid(color='white')\n",
    "    #plt.close()\n",
    "    plt.ylim(250,7200)\n",
    "    #plt.show()\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from netCDF4 import Dataset\n",
    "#from wradlib import get_clip_mask\n",
    "\n",
    "def read_gpm(filename, bbox):\n",
    "    scan='NS' # NS\n",
    "    pr_data = Dataset(filename, mode=\"r\")\n",
    "    lon = pr_data[scan].variables['Longitude']\n",
    "    lat = pr_data[scan].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 = wrl.zonalstats.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[scan]['ScanTime'].variables['Year'][mask]\n",
    "    month = pr_data[scan]['ScanTime'].variables['Month'][mask]\n",
    "    dayofmonth = pr_data[scan]['ScanTime'].variables['DayOfMonth'][mask]\n",
    "    # dayofyear = pr_data[scan]['ScanTime'].variables['DayOfYear'][mask]\n",
    "    hour = pr_data[scan]['ScanTime'].variables['Hour'][mask]\n",
    "    minute = pr_data[scan]['ScanTime'].variables['Minute'][mask]\n",
    "    second = pr_data[scan]['ScanTime'].variables['Second'][mask]\n",
    "    # secondofday = pr_data[scan]['ScanTime'].variables['SecondOfDay'][mask]\n",
    "    millisecond = pr_data[scan]['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[scan]['PRE'].variables['landSurfaceType'][mask]\n",
    "    pflag = pr_data[scan]['PRE'].variables['flagPrecip'][mask]\n",
    "\n",
    "    # bbflag = pr_data[scan]['CSF'].variables['flagBB'][mask]\n",
    "    zbb = pr_data[scan]['CSF'].variables['heightBB'][mask]\n",
    "    # print(zbb.dtype)\n",
    "    bbwidth = pr_data[scan]['CSF'].variables['widthBB'][mask]\n",
    "    qbb = pr_data[scan]['CSF'].variables['qualityBB'][mask]\n",
    "    qtype = pr_data[scan]['CSF'].variables['qualityTypePrecip'][mask]\n",
    "    ptype = pr_data[scan]['CSF'].variables['typePrecip'][mask]\n",
    "\n",
    "    quality = pr_data[scan]['scanStatus'].variables['dataQuality'][mask]\n",
    "    ## REFL\n",
    "    refl = pr_data[scan]['SLV'].variables['zFactorCorrected'][mask]\n",
    "    #refl = pr_data[scan]['PRE'].variables['zFactorMeasured'][mask]\n",
    "    # print(pr_data[scan]['SLV'].variables['zFactorCorrected']) \n",
    "\n",
    "    zenith = pr_data[scan]['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"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# define GPM data set\n",
    "# ----------------------\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_path = '/automount/ags/velibor/gpmdata/dprV7/2A.GPM.DPR.V7-20170308.20160601-S171836-E185110.012837.V05A.HDF5'\n",
    "#gpm_file_path = '/automount/ags/velibor/gpmdata/dprV7/2A.GPM.DPR.V7-20170308.20150330-S224903-E002135.006166.V05A.HDF5'\n",
    "#gpm_file_path = '/automount/ags/velibor/gpmdata/dprV7/2A.GPM.DPR.V7-20170308.20151216-S014742-E032015.010213.V05A.HDF5'\n",
    "\n",
    "# define matching ground radar file\n",
    "# ----------------------\n",
    "gr_file_path = '/automount/radar-archiv/scans/2014/2014-10/2014-10-07/n_ppi_010deg/2014-10-07--02:37:44,00.mvol'\n",
    "#gr_file_path = '/automount/radar-archiv/scans/2014/2014-10/2014-10-08/n_ppi_010deg/2014-10-08--09:45:00,00.mvol'\n",
    "#gr_file_path = '/automount/radar/scans/2016/2016-06/2016-06-01/n_ppi_010deg/2016-06-01--17:59:50,00.mvol'\n",
    "#gr_file_path = '/automount/radar/scans/2015/2015-03/2015-03-30/n_ppi_010deg/2015-03-30--23:30:03,00.mvol'\n",
    "#gr_file_path = '/automount/radar/scans/2015/2015-12/2015-12-16/n_ppi_010deg/2015-12-16--02:45:01,00.mvol'\n",
    "\n",
    "gpm_file = wrl.util.get_wradlib_data_file(gpm_file_path)\n",
    "gr2gpm_file = wrl.util.get_wradlib_data_file(gr_file_path)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "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": 6,
   "metadata": {
    "collapsed": true
   },
   "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"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def read_gr_boxpol(filename, loaddata=True):\n",
    "    \"\"\"\n",
    "    Read ground radar boxpol scans\n",
    "    \n",
    "    in arbeit....\n",
    "    \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",
    "    \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",
    "    print c_liste\n",
    "    ntilt = len(c_liste)\n",
    "    \n",
    "    ref_data = [[]]*ntilt\n",
    "    rho_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",
    "        \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': dat_, '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",
    "        #dh = grdata['SCAN0']['ZH']['data']\n",
    "        #dv = grdata['SCAN0']['ZV']['data']\n",
    "        #data = (dh + dv)/2.\n",
    "        \n",
    "        data = grdata['SCAN0']['ZH']['data']\n",
    "        \n",
    "        \n",
    "        rho = grdata['SCAN0']['RHOHV']['data']\n",
    "        \n",
    "\n",
    "        #refl.append(data)\n",
    "        ref_data[jj]=data\n",
    "        rho_data[jj]=rho\n",
    "        #ref_data[ref_data<=0]=np.nan\n",
    "\n",
    "        #refl1 = np.array(refl)\n",
    "    gr_dict.update({'refl': ref_data})\n",
    "    gr_dict.update({'rho': rho_data})\n",
    "    \n",
    "    return gr_dict, ntilt\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "/automount/radar-archiv/scans/2014/2014-10/2014-10-07/n_ppi_010deg/2014-10-07--02:37:44,00.mvol\n",
      "['/automount/radar-archiv/scans/2014/2014-10/2014-10-07/n_ppi_010deg/2014-10-07--02:37:44,00.mvol', '/automount/radar-archiv/scans/2014/2014-10/2014-10-07/n_ppi_045deg/2014-10-07--02:37:44,00.mvol', '/automount/radar-archiv/scans/2014/2014-10/2014-10-07/n_ppi_082deg/2014-10-07--02:37:44,00.mvol', '/automount/radar-archiv/scans/2014/2014-10/2014-10-07/n_ppi_110deg/2014-10-07--02:37:44,00.mvol', '/automount/radar-archiv/scans/2014/2014-10/2014-10-07/n_ppi_140deg/2014-10-07--02:37:44,00.mvol', '/automount/radar-archiv/scans/2014/2014-10/2014-10-07/n_ppi_180deg/2014-10-07--02:37:44,00.mvol', '/automount/radar-archiv/scans/2014/2014-10/2014-10-07/n_ppi_280deg/2014-10-07--02:37:44,00.mvol', '/automount/radar-archiv/scans/2014/2014-10/2014-10-07/ppi_1p5deg/2014-10-07--02:37:44,00.mvol', '/automount/radar-archiv/scans/2014/2014-10/2014-10-07/ppi_2p4deg/2014-10-07--02:37:44,00.mvol', '/automount/radar-archiv/scans/2014/2014-10/2014-10-07/ppi_3p4deg/2014-10-07--02:37:44,00.mvol']\n",
      "('ntilit:', 10)\n",
      "/automount/radar-archiv/scans/2014/2014-10/2014-10-07/n_ppi_010deg/2014-10-07--02:37:44,00.mvol\n",
      "/automount/radar-archiv/scans/2014/2014-10/2014-10-07/n_ppi_045deg/2014-10-07--02:37:44,00.mvol\n",
      "/automount/radar-archiv/scans/2014/2014-10/2014-10-07/n_ppi_082deg/2014-10-07--02:37:44,00.mvol\n",
      "/automount/radar-archiv/scans/2014/2014-10/2014-10-07/n_ppi_110deg/2014-10-07--02:37:44,00.mvol\n",
      "/automount/radar-archiv/scans/2014/2014-10/2014-10-07/n_ppi_140deg/2014-10-07--02:37:44,00.mvol\n",
      "/automount/radar-archiv/scans/2014/2014-10/2014-10-07/n_ppi_180deg/2014-10-07--02:37:44,00.mvol\n",
      "/automount/radar-archiv/scans/2014/2014-10/2014-10-07/n_ppi_280deg/2014-10-07--02:37:44,00.mvol\n",
      "/automount/radar-archiv/scans/2014/2014-10/2014-10-07/ppi_1p5deg/2014-10-07--02:37:44,00.mvol\n",
      "/automount/radar-archiv/scans/2014/2014-10/2014-10-07/ppi_2p4deg/2014-10-07--02:37:44,00.mvol\n",
      "/automount/radar-archiv/scans/2014/2014-10/2014-10-07/ppi_3p4deg/2014-10-07--02:37:44,00.mvol\n"
     ]
    }
   ],
   "source": [
    "# read matching GR data\n",
    "gr_data, ntilt = read_gr_boxpol(gr_file)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\"for i in range(len(gr_data.keys())):\\n\\n    print (gr_data.keys()[i])\\n    print(gr_data[gr_data.keys()[i]])\\n    print ('------------------------- \\n \\n \\n')\\n\""
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\"\"\"for i in range(len(gr_data.keys())):\n",
    "\n",
    "    print (gr_data.keys()[i])\n",
    "    print(gr_data[gr_data.keys()[i]])\n",
    "    print ('------------------------- \\n \\n \\n')\n",
    "\"\"\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Beginn der QnD shitty For schleife"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "('----------------Elevationswinkel: ', 0.999755859375)\n",
      "('----------------number of rays in gr sweep: ', 360)\n",
      "('----------------number of gates in gr beam: ', 1000)\n",
      "('----------------number of sweeps: ', 10)\n",
      "('----------------elevation of sweep: ', 0.999755859375)\n",
      "('----------------range of first gate: ', 0.0)\n",
      "('----------------azimuth angle of first beam: ', 0.0)\n",
      "('----------------Longitude of GR: ', 7.071663)\n",
      "('----------------Latitude of GR: ', 50.73052)\n",
      "('----------------Altitude of GR (meters): ', 99.5)\n",
      "('----------------gate length (meters): ', 150.0)\n",
      "________TH______\n",
      "Radar bounding box:\n",
      "\t52.08\n",
      "4.95           9.19\n",
      "\t49.38\n",
      "('----------------Elevationswinkel: ', 4.4989013671875)\n",
      "('----------------number of rays in gr sweep: ', 360)\n",
      "('----------------number of gates in gr beam: ', 1000)\n",
      "('----------------number of sweeps: ', 10)\n",
      "('----------------elevation of sweep: ', 4.4989013671875)\n",
      "('----------------range of first gate: ', 0.0)\n",
      "('----------------azimuth angle of first beam: ', 0.0)\n",
      "('----------------Longitude of GR: ', 7.071663)\n",
      "('----------------Latitude of GR: ', 50.73052)\n",
      "('----------------Altitude of GR (meters): ', 99.5)\n",
      "('----------------gate length (meters): ', 150.0)\n",
      "________TH______\n",
      "Radar bounding box:\n",
      "\t52.07\n",
      "4.96           9.19\n",
      "\t49.39\n",
      "('----------------Elevationswinkel: ', 8.2012939453125)\n",
      "('----------------number of rays in gr sweep: ', 360)\n",
      "('----------------number of gates in gr beam: ', 1100)\n",
      "('----------------number of sweeps: ', 10)\n",
      "('----------------elevation of sweep: ', 8.2012939453125)\n",
      "('----------------range of first gate: ', 0.0)\n",
      "('----------------azimuth angle of first beam: ', 0.0)\n",
      "('----------------Longitude of GR: ', 7.071663)\n",
      "('----------------Latitude of GR: ', 50.73052)\n",
      "('----------------Altitude of GR (meters): ', 99.5)\n",
      "('----------------gate length (meters): ', 100.0)\n",
      "________TH______\n",
      "Radar bounding box:\n",
      "\t51.71\n",
      "5.53           8.61\n",
      "\t49.75\n",
      "('----------------Elevationswinkel: ', 10.997314453125)\n",
      "('----------------number of rays in gr sweep: ', 360)\n",
      "('----------------number of gates in gr beam: ', 1000)\n",
      "('----------------number of sweeps: ', 10)\n",
      "('----------------elevation of sweep: ', 10.997314453125)\n",
      "('----------------range of first gate: ', 0.0)\n",
      "('----------------azimuth angle of first beam: ', 0.0)\n",
      "('----------------Longitude of GR: ', 7.071663)\n",
      "('----------------Latitude of GR: ', 50.73052)\n",
      "('----------------Altitude of GR (meters): ', 99.5)\n",
      "('----------------gate length (meters): ', 100.0)\n",
      "________TH______\n",
      "Radar bounding box:\n",
      "\t51.61\n",
      "5.69           8.46\n",
      "\t49.85\n",
      "('----------------Elevationswinkel: ', 14.0020751953125)\n",
      "('----------------number of rays in gr sweep: ', 360)\n",
      "('----------------number of gates in gr beam: ', 800)\n",
      "('----------------number of sweeps: ', 10)\n",
      "('----------------elevation of sweep: ', 14.0020751953125)\n",
      "('----------------range of first gate: ', 0.0)\n",
      "('----------------azimuth angle of first beam: ', 0.0)\n",
      "('----------------Longitude of GR: ', 7.071663)\n",
      "('----------------Latitude of GR: ', 50.73052)\n",
      "('----------------Altitude of GR (meters): ', 99.5)\n",
      "('----------------gate length (meters): ', 100.0)\n",
      "________TH______\n",
      "Radar bounding box:\n",
      "\t51.43\n",
      "5.98           8.17\n",
      "\t50.03\n",
      "('----------------Elevationswinkel: ', 18.0010986328125)\n",
      "('----------------number of rays in gr sweep: ', 360)\n",
      "('----------------number of gates in gr beam: ', 550)\n",
      "('----------------number of sweeps: ', 10)\n",
      "('----------------elevation of sweep: ', 18.0010986328125)\n",
      "('----------------range of first gate: ', 0.0)\n",
      "('----------------azimuth angle of first beam: ', 0.0)\n",
      "('----------------Longitude of GR: ', 7.071663)\n",
      "('----------------Latitude of GR: ', 50.73052)\n",
      "('----------------Altitude of GR (meters): ', 99.5)\n",
      "('----------------gate length (meters): ', 100.0)\n",
      "________TH______\n",
      "Radar bounding box:\n",
      "\t51.20\n",
      "6.33           7.81\n",
      "\t50.26\n",
      "('----------------Elevationswinkel: ', 27.9986572265625)\n",
      "('----------------number of rays in gr sweep: ', 360)\n",
      "('----------------number of gates in gr beam: ', 350)\n",
      "('----------------number of sweeps: ', 10)\n",
      "('----------------elevation of sweep: ', 27.9986572265625)\n",
      "('----------------range of first gate: ', 0.0)\n",
      "('----------------azimuth angle of first beam: ', 0.0)\n",
      "('----------------Longitude of GR: ', 7.071663)\n",
      "('----------------Latitude of GR: ', 50.73052)\n",
      "('----------------Altitude of GR (meters): ', 99.5)\n",
      "('----------------gate length (meters): ', 100.0)\n",
      "________TH______\n",
      "Radar bounding box:\n",
      "\t51.01\n",
      "6.64           7.51\n",
      "\t50.45\n",
      "('----------------Elevationswinkel: ', 1.4996337890625)\n",
      "('----------------number of rays in gr sweep: ', 360)\n",
      "('----------------number of gates in gr beam: ', 1000)\n",
      "('----------------number of sweeps: ', 10)\n",
      "('----------------elevation of sweep: ', 1.4996337890625)\n",
      "('----------------range of first gate: ', 0.0)\n",
      "('----------------azimuth angle of first beam: ', 0.0)\n",
      "('----------------Longitude of GR: ', 7.071663)\n",
      "('----------------Latitude of GR: ', 50.73052)\n",
      "('----------------Altitude of GR (meters): ', 99.5)\n",
      "('----------------gate length (meters): ', 100.0)\n",
      "________TH______\n",
      "Radar bounding box:\n",
      "\t51.63\n",
      "5.66           8.49\n",
      "\t49.83\n",
      "('----------------Elevationswinkel: ', 2.4005126953125)\n",
      "('----------------number of rays in gr sweep: ', 360)\n",
      "('----------------number of gates in gr beam: ', 1000)\n",
      "('----------------number of sweeps: ', 10)\n",
      "('----------------elevation of sweep: ', 2.4005126953125)\n",
      "('----------------range of first gate: ', 0.0)\n",
      "('----------------azimuth angle of first beam: ', 0.0)\n",
      "('----------------Longitude of GR: ', 7.071663)\n",
      "('----------------Latitude of GR: ', 50.73052)\n",
      "('----------------Altitude of GR (meters): ', 99.5)\n",
      "('----------------gate length (meters): ', 100.0)\n",
      "________TH______\n",
      "Radar bounding box:\n",
      "\t51.63\n",
      "5.66           8.49\n",
      "\t49.83\n",
      "('----------------Elevationswinkel: ', 6.998291015625)\n",
      "('----------------number of rays in gr sweep: ', 360)\n",
      "('----------------number of gates in gr beam: ', 2000)\n",
      "('----------------number of sweeps: ', 10)\n",
      "('----------------elevation of sweep: ', 6.998291015625)\n",
      "('----------------range of first gate: ', 0.0)\n",
      "('----------------azimuth angle of first beam: ', 0.0)\n",
      "('----------------Longitude of GR: ', 7.071663)\n",
      "('----------------Latitude of GR: ', 50.73052)\n",
      "('----------------Altitude of GR (meters): ', 99.5)\n",
      "('----------------gate length (meters): ', 25.0)\n",
      "________TH______\n",
      "Radar bounding box:\n",
      "\t51.18\n",
      "6.37           7.77\n",
      "\t50.28\n"
     ]
    }
   ],
   "source": [
    "SR, GR =  np.array([]), np.array([])\n",
    "X, Y, Z = np.array([]), np.array([]), np.array([])\n",
    "\n",
    "GRnum, GRstd =  np.array([]), np.array([])\n",
    "\n",
    "gr_nr_elev, gr_ele = np.array([]), np.array([])\n",
    "\n",
    "\n",
    "for ijk in range(ntilt):\n",
    "#for ijk in range(1,6,2):\n",
    "\n",
    "    ######### GR data and att\n",
    "    ee = ijk\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",
    "    # 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.###############################einlesen\n",
    "    # reflectivity array of sweep\n",
    "    ref_gr = gr_data['refl'][ee]\n",
    "    # rhohv array of sweep\n",
    "    rho_gr = gr_data['rho'][ee]\n",
    "    \n",
    "    if elev_gr<0.01:\n",
    "        print (\"Elevation zu niedrig\")\n",
    "        pass\n",
    "    else:\n",
    "        ############################################################################## Korrektur der GR Daten\n",
    "        print (\"----------------Elevationswinkel: \", elev_gr)\n",
    "        print (\"----------------number of rays in gr sweep: \", nray_gr)\n",
    "        print (\"----------------number of gates in gr beam: \", ngate_gr)\n",
    "        print (\"----------------number of sweeps: \", nelev)\n",
    "        print (\"----------------elevation of sweep: \", elev_gr)\n",
    "        print (\"----------------range of first gate: \", r0_gr)\n",
    "        print (\"----------------azimuth angle of first beam: \", a0_gr)\n",
    "        print (\"----------------Longitude of GR: \", lon0_gr)\n",
    "        print (\"----------------Latitude of GR: \", lat0_gr)\n",
    "        print (\"----------------Altitude of GR (meters): \",alt0_gr )\n",
    "        print (\"----------------gate length (meters): \",dr_gr )\n",
    "        \n",
    "        \n",
    "        \n",
    "        #print (\"________Offset______\")\n",
    "        #Hardcoded....passt in etwa\n",
    "        #ref_gr = ref_gr + 2\n",
    "\n",
    "\n",
    "        #print (\"________ATTCORR______\")\n",
    "        \n",
    "        \"\"\"        pia_harrison = wrl.atten.correctAttenuationHB(\n",
    "            ref_gr,\n",
    "            coefficients = dict(a=4.57e-5, b=0.731, gate_length=1.0),\n",
    "            mode=\"warn\",\n",
    "            thrs=59.)\n",
    "        pia_harrison[pia_harrison > 4.8] = 4.8\n",
    "\n",
    "        #print (\"________ATTCORR2______\")\n",
    "        ref_gr = ref_gr + pia_harrison\"\"\"\n",
    "\n",
    "\n",
    "        #print (\"________CLUTTER______\")\n",
    "        #rho_th  = 0.85\n",
    "        rho_th = 0.9\n",
    "        ref_gr[rho_gr<= rho_th] = np.nan\n",
    "    \n",
    "        print (\"________TH______\")\n",
    "        #TH = 0\n",
    "        #ref_gr[ref_gr<=TH]=np.nan\n",
    "        \n",
    "        \n",
    "        \n",
    "        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']))\n",
    "    \n",
    "\n",
    "        # 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",
    "\n",
    "        # 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",
    "        \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",
    "\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",
    "\n",
    "        \n",
    "        X = np.append(X,gr_xyz[...,0])\n",
    "        Y = np.append(Y,gr_xyz[...,1])\n",
    "        Z = np.append(Z,gr_xyz[...,2])\n",
    "\n",
    "        GR = np.append(GR, ref_gr)\n",
    "\n",
    "\n",
    "        gr_nr_elev = np.append(gr_nr_elev, ref_gr.shape[0]*ref_gr.shape[1])\n",
    "        gr_ele = np.append(gr_ele, elev_gr)\n",
    "\n",
    "\n",
    "\n",
    "    \n",
    "\n",
    "    \n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[360000. 360000. 396000. 360000. 288000. 198000. 126000. 360000. 360000.\n",
      " 720000.]\n"
     ]
    }
   ],
   "source": [
    "pl.scatter(X[0:360000],Y[0:360000],c=GR[0:360000], edgecolors='none', vmin=0)\n",
    "pl.colorbar()\n",
    "pl.show()\n",
    "print(gr_nr_elev)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def get_my_cmap2():\n",
    "    import matplotlib.cm as cm\n",
    "    my_cmap = cm.get_cmap('jet',40)\n",
    "    my_cmap.set_under('white')\n",
    "    my_cmap.set_over('darkred')\n",
    "    return my_cmap\n",
    "\n",
    "def subplotcfad(Z,GR):\n",
    "    import matplotlib.pyplot as plt\n",
    "    m1 = ~np.isnan(Z) & ~np.isnan(GR)\n",
    "\n",
    "    plt.hist2d(GR[m1], Z[m1],bins=120, cmap=get_my_cmap2(), vmin=0.1)#, vmax=2700)\n",
    "    #plt.hexbin(GR[m1], Z[m1],bins=300, cmap=get_my_cmap2(), vmin=0.1)#, vmax=2700)\n",
    "    plt.colorbar()\n",
    "    plt.xlabel(r'$Z_H$')\n",
    "    plt.ylabel('Height in m')\n",
    "\n",
    "    plt.legend(loc='upper right', fontsize=10)\n",
    "    #plt.grid(color='white')\n",
    "\n",
    "    plt.ylim(0,7200)\n",
    "    plt.xlim(0,50.)\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 3.67 s, sys: 2.1 s, total: 5.77 s\n",
      "Wall time: 5 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "# Elevation unterteilen mit Cumulativer Summe\n",
    "e_pos = gr_nr_elev.cumsum().astype(int)\n",
    "from pcc import get_miub_cmap\n",
    "\n",
    "\n",
    "fig = pl.figure(figsize=(18,6))\n",
    "#plt.suptitle('CFAD 20141007')\n",
    "for eee in range(len(gr_nr_elev)):\n",
    "    ax1 = fig.add_subplot(2,5,eee+1)\n",
    "    \n",
    "    if eee==0: \n",
    "        subplotcfad(Z[0:e_pos[0]],GR[0:e_pos[0]])\n",
    "        #pl.legend(loc='upper left')\n",
    "        pl.title('Elevation: '+str(gr_ele[eee])[0:4])\n",
    "\n",
    "    else:\n",
    "        subplotcfad(Z[e_pos[eee-1]:e_pos[eee]],\n",
    "                   GR[e_pos[eee-1]:e_pos[eee]])\n",
    "        #pl.legend(loc='upper left')\n",
    "        pl.title('Elevation: '+ str(gr_ele[eee])[0:4])\n",
    "\n",
    "\n",
    "pl.tight_layout()\n",
    "#pl.savefig('/automount/ftp/velibor/validation/boxpol_cfad_all_elevation_zdr.pdf')\n",
    "#pl.savefig('/automount/ftp/velibor/validation/boxpol_cfad_all_elevation_zdr.png')\n",
    "\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": []
  },
  {
   "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
}