pcc_gpm_matching_all_elevation.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "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",
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "import glob\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "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,20000,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": null,
   "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, MS, HS\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",
    "    #DPR Altitude in m \n",
    "    dpr_alt = pr_data[scan]['navigation'].variables['dprAlt'][mask]\n",
    "    print ('dpr_alt: ', dpr_alt)\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",
    "    \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",
    "    #### DFR\n",
    "    #refl_ns = pr_data['NS']['PRE'].variables['zFactorMeasured'][mask]\n",
    "    ## Auf Ka zuschneiden\n",
    "    #refl_ns = refl_ns[:,12:37]\n",
    "    \n",
    "    #refl_ms = pr_data['MS']['PRE'].variables['zFactorMeasured'][mask]\n",
    "    #refl = refl_ns - refl_ms\n",
    "    \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, 'dpr_alt':dpr_alt})\n",
    "\n",
    "    return gpm_data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "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",
    "\n",
    "\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",
    "#gpm_file_path = '/automount/ags/velibor/gpmdata/dprV7/2A.GPM.DPR.V7-20170308.20141217-S050923-E064153.004552.V05A.HDF5'\n",
    "#gpm_file_path = '/automount/ags/velibor/gpmdata/dprV7/2A.GPM.DPR.V7-20170308.20150128-S164142-E181415.005213.V05A.HDF5'\n",
    "#gpm_file_path = '/automount/ags/velibor/gpmdata/dprV7/2A.GPM.DPR.V7-20170308.20151211-S195806-E213039.010147.V05A.HDF5'\n",
    "###gpm_file_path = '/automount/ags/velibor/gpmdata/dprV7/2A.GPM.DPR.V7-20170308.20170519-S102835-E120109.018308.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",
    "\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",
    "#gr_file_path = '/automount/radar-archiv/scans/2014/2014-12/2014-12-17/n_ppi_010deg/2014-12-17--05:45:00,00.mvol'\n",
    "#gr_file_path = '/automount/radar/scans/2015/2015-01/2015-01-28/n_ppi_010deg/2015-01-28--17:22:08,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-12/2015-12-11/n_ppi_010deg/2015-12-11--22:35:00,00.mvol'\n",
    "###gr_file_path = '/automount/radar/scans/2017/2017-05/2017-05-19/n_ppi_010deg/n_ppi_010deg_12345_20170519110333_00.h5'\n",
    "\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": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "gr_file_path"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Space-born precipitation radar parameters\n",
    "sr_pars = {\"gpm\": {\n",
    "    \"zt\": 409000.,  # orbital height of GPM                 ##################APPROXIMATION!\n",
    "    \"dr\": 125.,      # gate spacing of GPM\n",
    "    \"gr_file\": gr2gpm_file,\n",
    "}}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "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": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def read_gr_boxpol(gr_file, 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",
    "    # Bei ALTEN einstellungen\n",
    "    c_liste = glob.glob(gr_file.replace('n_ppi_010deg','*ppi*'))\n",
    "    print c_liste\n",
    "    \n",
    "    #bei NEUEN -.- \n",
    "    #new_gr_file = gr_file.replace('n_ppi_010deg','*ppi*')\n",
    "    #c_liste = sorted(glob.glob(new_gr_file[0:-9]+'*'))[::2]\n",
    "\n",
    "    print c_liste\n",
    "    ntilt = len(c_liste)\n",
    "    \n",
    "    ref_data = [[]]*ntilt\n",
    "    rho_data = [[]]*ntilt\n",
    "    kdp_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",
    "\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",
    "        data = grdata['SCAN0']['ZH']['data']\n",
    "        \n",
    "        \n",
    "        rho = grdata['SCAN0']['RHOHV']['data']\n",
    "        kdp = grdata['SCAN0']['KDP']['data']\n",
    "        \n",
    "\n",
    "        #refl.append(data)\n",
    "        ref_data[jj]=data\n",
    "        rho_data[jj]=rho\n",
    "        kdp_data[jj]=kdp\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",
    "    gr_dict.update({'kdp': kdp_data})\n",
    "    \n",
    "    return gr_dict, ntilt\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "\"\"\"glob.glob(gr_file.replace('n_ppi_010deg','*ppi*'))\n",
    "new_gr_file = gr_file.replace('n_ppi_010deg','*ppi*')\n",
    "print 'new file', new_gr_file[0:-9]+'*'\n",
    "sorted(glob.glob(new_gr_file[0:-9]+'*'))\n",
    "\"\"\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# read matching GR data\n",
    "gr_data, ntilt = read_gr_boxpol(gr_file_path)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "\"\"\"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')\"\"\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Beginn der QnD shitty For loop"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "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 np.array([0, 3]):\n",
    "\n",
    "for ijk in np.array([0,7, 8, 1, 9, 2, 3, 4, 5, 6]):\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",
    "    # kdp array of sweep\n",
    "    kdp_gr = gr_data['kdp'][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",
    "        \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 (\"________Offset______\")\n",
    "        #Hardcoded....passt in etwa\n",
    "        \"\"\"ref_gr = ref_gr + 2\"\"\"\n",
    "        \n",
    "        # print (\"________ATTCORR______\")\n",
    "        # http://docs.wradlib.org/en/latest/notebooks/attenuation/wradlib_attenuation.html\n",
    "        #alpha_xband = 0.28\n",
    "        #ref_gr = ref_gr + (alpha_xband * kdp_gr)\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",
    "        #ref_gr = pia_harrison\"\"\"\n",
    "        \n",
    "        print (\"________CLUTTER______\")\n",
    "        \"\"\"rho_th = 0.98\n",
    "        ref_gr[rho_gr<= rho_th] = np.nan\"\"\"\n",
    "    \n",
    "        #print (\"________TH______\")\n",
    "        \"\"\"TH = 15\n",
    "        ref_gr[ref_gr<=TH]=np.nan\"\"\"\n",
    "        \n",
    "        \n",
    "        print (\"________Spatial Thresholds_________\")\n",
    "        \n",
    "        ##############\n",
    "        ##############\n",
    "        ### synthetic#\n",
    "        ##############\n",
    "        ##############\n",
    "        #ref_gr = np.ones(ref_gr.shape)*25\n",
    "        \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",
    "        \n",
    "        #################################################################################\n",
    "        # Offset\n",
    "        #ref_gr = ref_gr + 2\n",
    "        \n",
    "        #AttCorr HB\n",
    "        \"\"\"pia_hb= 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_hb[pia_hb > 4.8] = 4.8\n",
    "        \n",
    "        ref_gr = ref_gr + pia_hb\"\"\"\n",
    "        \n",
    "        # Clutter\n",
    "        rho_th = 0.95\n",
    "        #ref_gr[(rho_gr<= rho_th) & (alt<=2000)] = np.nan\n",
    "        ref_gr[(rho_gr<= rho_th)] = np.nan\n",
    "        print (rho_gr.shape, alt.shape)\n",
    "        \n",
    "        #GPM Threshold\n",
    "        TH = 15\n",
    "        ref_gr[ref_gr<=TH]=np.nan\n",
    "        ####################################################################################\n",
    "        \n",
    "        \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",
    "        ######################################################################################Read Gpm file\n",
    "        sr_data = read_gpm(gpm_file, bbox)\n",
    "        #sr_data = wrl.io.read_gpm(gpm_file, bbox)\n",
    "        print ('...SR_Data loaded...')\n",
    "        # 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",
    "        \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']\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",
    "        \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",
    "        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",
    "\n",
    "        precip_mask = (pflag >= 2) & wrl.zonalstats.get_clip_mask(sr_xy, gr_domain, rad)\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",
    "        \n",
    "        \n",
    "        xyzp_sr = np.concatenate((xyp_sr, z_sr[..., np.newaxis]),\n",
    "                           axis=-1)\n",
    "        r_sr, az_sr, elev_sr = wrl.georef.xyz_to_spherical(xyzp_sr, alt0_gr, proj=rad)\n",
    "        \n",
    "        #elev_gr\n",
    "        mask = (elev_sr > (elev_gr - bw_gr/2.)) & (elev_sr < (elev_gr + bw_gr/2.))\n",
    "        \n",
    "        rs = wrl.georef.dist_from_orbit(zt, alpha, beta, r_sr_inv, re1)\n",
    "\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",
    "        # 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",
    "        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",
    "        ratio, ibb = wrl.qual.get_bb_ratio(sr_data['zbb'], sr_data['bbwidth'], sr_data['quality'], z_sr)\n",
    "\n",
    "        zbb = sr_data['zbb'].copy()\n",
    "        zbb[~ibb] = np.nan\n",
    "        \n",
    "        ####################################################################BB entfernen\n",
    "        #ref_gr_ku[gr_xyz[...,2] >= np.nanmean(zbb)] = np.nan\n",
    "        \n",
    "\n",
    "        \n",
    "        \n",
    "\n",
    "        ##############\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",
    "        # 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",
    "\n",
    "        # average coordinates\n",
    "        xyz_v1 = xyzp_sr.copy()\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",
    "\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",
    "        r_sr, az_sr, el_rs = wrl.georef.xyz_to_spherical(xyz, alt0_gr, proj=rad)\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",
    "        print ('... start zonal stats')\n",
    "\n",
    "        zds = wrl.zonalstats.DataSource(xyz[:, 0:2].reshape(-1, 2), rad)\n",
    "\n",
    "\n",
    "        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)    \n",
    "\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))\n",
    "\n",
    "        #obj3 = wrl.zonalstats.GridCellsToPoly('m3d_zonal_poly_{0}'.format(platf))\n",
    "        obj3 = wrl.zonalstats.GridCellsToPoly(zdp)\n",
    "\n",
    "        print ('... end zonal stats')\n",
    "        #print(obj3.ix.shape)\n",
    "        \n",
    "        print ('... averaging GR bins...')\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",
    "        ####################### Mitteln von GR bins\n",
    "        ref_gr_i = wrl.trafo.idecibel(ref_gr.ravel())\n",
    "\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",
    "        \n",
    "        print ('... STD of GR bins...')\n",
    "        ############################# Berechnung von STD\n",
    "        refgr1std = np.array([np.nanstd(ref_gr.ravel()[obj3.ix[i]])\n",
    "                     for i in np.arange(len(obj3.ix))[~obj3.check_empty()]])\n",
    "        #refgr1std = wrl.trafo.decibel(refgr1std)\n",
    "        \n",
    "        print ('... NUM of GR bins...')\n",
    "        ########################## Berechnung von Num\n",
    "        refgr1num = np.array([np.count_nonzero(~np.isnan(ref_gr.ravel()[obj3.ix[i]]))\n",
    "                     for i in np.arange(len(obj3.ix))[~obj3.check_empty()]])\n",
    "        \n",
    "\n",
    "        \n",
    "        \n",
    "        print ('Durchgang: ', ijk)\n",
    "\n",
    "        \n",
    "        X = np.append(X,xyz[...,0])\n",
    "        Y = np.append(Y,xyz[...,1])\n",
    "        Z = np.append(Z,xyz[...,2])\n",
    "\n",
    "        SR = np.append(SR, refsr1a)\n",
    "        GR = np.append(GR, refgr1a)\n",
    "        GRstd = np.append(GRstd, refgr1std)\n",
    "        GRnum = np.append(GRnum, refgr1num)\n",
    "\n",
    "        gr_nr_elev = np.append(gr_nr_elev, refgr1a.shape)\n",
    "        gr_ele = np.append(gr_ele, elev_gr)\n",
    "        \n",
    "        \n",
    "\n",
    "\n",
    "\n",
    "    \n",
    "\n",
    "    \n",
    "  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "0 - No precipitation by both Ku and Ka\n",
    "1 - Precipitation by Ka, no rain by Ku\n",
    "10 - Precipitation by Ku, no rain by Ka\n",
    "11 - Precipitation by both Ku and Ka\n",
    "-9999 -  Missing value\n",
    "\"\"\"\n",
    "print(np.unique(pflag))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "np.mean(sr_data['dpr_alt'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from satlib import corcor\n",
    "from pcc import get_miub_cmap\n",
    "TH = 15\n",
    "SR[SR<=TH]=np.nan\n",
    "GR[GR<=TH]=np.nan\n",
    "\n",
    "\n",
    "fig = pl.figure(figsize=(14,14))\n",
    "ax = fig.add_subplot(221, aspect=\"equal\")\n",
    "pl.scatter(SR, GR, 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.ylabel(\"GR reflectivity (dBZ)\")\n",
    "pl.xlabel(\"SR reflectivity (dBZ)\")\n",
    "pl.title(corcor(GR, SR))\n",
    "\n",
    "ax = fig.add_subplot(222)\n",
    "m0 = ~np.isnan(GR)& ~np.isnan(SR)\n",
    "\n",
    "pl.hist(GR[m0], bins=20, edgecolor=\"blue\", facecolor=\"grey\", label=\"GR\")\n",
    "pl.hist(SR[m0], 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(223)\n",
    "pl.hist2d(SR[m0], GR[m0],bins=20)\n",
    "\n",
    "pl.ylabel(\"GR reflectivity (dBZ)\")\n",
    "pl.xlabel(\"SR reflectivity (dBZ)\")\n",
    "pl.legend()\n",
    "pl.grid()\n",
    "pl.colorbar()\n",
    "#pl.xlim(10,50)\n",
    "\n",
    "\n",
    "\n",
    "pl.tight_layout()\n",
    "pl.savefig('/automount/ftp/velibor/validation/matching1.pdf')\n",
    "pl.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "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(X, Y, Z, c=SR.ravel(), vmin=0, vmax=40, cmap=get_miub_cmap())\n",
    "fig.colorbar(pppl, shrink=0.5, aspect=5)\n",
    "ax1.set_title('Ref SR')\n",
    "#ax1.set_xlim3d(-50000,50000)\n",
    "#ax1.set_ylim3d(-50000,50000)\n",
    "#ax1.set_zlim3d(0,50000)\n",
    "ax1.set_xlabel(\"x\")\n",
    "ax1.set_ylabel(\"y\")\n",
    "ax1.set_zlabel(\"z\")\n",
    "\n",
    "ax2 = fig.add_subplot(132, projection='3d')\n",
    "pppl = ax2.scatter(X, Y, Z, c=GR.ravel(), vmin=0, vmax=40, cmap=get_miub_cmap())\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(X, Y, Z, c=SR.ravel()-GR.ravel(),\n",
    "                   cmap='coolwarm',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( SR, GR))\n",
    "pl.savefig('/automount/ftp/velibor/validation/matching2.pdf')\n",
    "pl.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "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', alpha=0.2)\n",
    "pppl = ax.scatter(xyz[...,0],xyz[...,1],xyz[...,2], edgecolors='none', s=300, c='black')\n",
    "ax.set_xlim(-10000, 10000)\n",
    "ax.set_ylim(-10000, 10000)\n",
    "ax.set_zlim(0, 2000)\n",
    "pl.show()\"\"\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "fig = pl.figure(figsize=(12,8))\n",
    "ax = fig.add_subplot(131, aspect=\"auto\")\n",
    "pl.scatter(Y,Z, c=SR,  vmin=0, vmax=40, edgecolor=\"None\", cmap=get_miub_cmap())\n",
    "pl.title(\"SR reflectivity\")\n",
    "pl.xlim(np.nanmin(Y),np.nanmax(Y))\n",
    "pl.ylim(0,8000)\n",
    "pl.grid()\n",
    "pl.colorbar(shrink=0.5)\n",
    "\n",
    "ax = fig.add_subplot(132, aspect=\"auto\")\n",
    "pl.scatter(Y,Z, c=GR,  vmin=0, vmax=40, edgecolor=\"None\", cmap=get_miub_cmap())\n",
    "pl.title(\"GR reflectivity\")\n",
    "pl.xlim(np.nanmin(Y),np.nanmax(Y))\n",
    "pl.ylim(0,8000)\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(Y,Z, c=SR-GR, cmap='jet',vmin=-10,vmax=10 ,edgecolor=\"None\")\n",
    "pl.title(\"SR - GR reflectivity\")\n",
    "pl.xlim(np.nanmin(Y),np.nanmax(Y))\n",
    "pl.ylim(0,8000)\n",
    "pl.grid()\n",
    "pl.colorbar(shrink=0.5)\n",
    "#fig.suptitle(\"uncorrected SR - GR\")\n",
    "pl.tight_layout()\n",
    "pl.savefig('/automount/ftp/velibor/validation/matching3.png')\n",
    "pl.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "fig = pl.figure(figsize=(12,8))\n",
    "ax = fig.add_subplot(131, aspect=\"auto\")\n",
    "pl.scatter(X,Z, c=SR, cmap=get_miub_cmap(), vmin=0, vmax=40, edgecolor=\"None\")\n",
    "pl.title(\"SR reflectivity\")\n",
    "pl.xlim(np.nanmin(X),np.nanmax(X))\n",
    "pl.ylim(np.nanmin(Z),8000)\n",
    "pl.grid()\n",
    "pl.colorbar(shrink=0.5)\n",
    "\n",
    "ax = fig.add_subplot(132, aspect=\"auto\")\n",
    "pl.scatter(X,Z, c=GR, cmap=get_miub_cmap(), vmin=0, vmax=40, edgecolor=\"None\")\n",
    "pl.title(\"GR reflectivity\")\n",
    "pl.xlim(np.nanmin(X),np.nanmax(X))\n",
    "pl.ylim(np.nanmin(Z),8000)\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(X,Z, c=SR-GR, cmap='jet',vmin=-10,vmax=10 ,edgecolor=\"None\")\n",
    "pl.title(\"SR - GR reflectivity\")\n",
    "pl.xlim(np.nanmin(X),np.nanmax(X))\n",
    "pl.ylim(np.nanmin(Z),8000)\n",
    "pl.grid()\n",
    "pl.colorbar(shrink=0.5)\n",
    "#fig.suptitle(\"uncorrected SR - GR\")\n",
    "pl.tight_layout()\n",
    "pl.savefig('/automount/ftp/velibor/validation/matching4.pdf')\n",
    "pl.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "fig = pl.figure(figsize=(12,8))\n",
    "ax = fig.add_subplot(131, aspect=\"equal\")\n",
    "pl.scatter(X,Y, c=SR, cmap=get_miub_cmap(), vmin=0, vmax=40, edgecolor=\"None\",s=50)\n",
    "pl.title(\"SR reflectivity\")\n",
    "pl.xlim(np.nanmin(X),np.nanmax(X))\n",
    "pl.ylim(np.nanmin(Y),np.nanmax(Y))\n",
    "pl.grid()\n",
    "pl.colorbar(shrink=0.5)\n",
    "\n",
    "ax = fig.add_subplot(132, aspect=\"equal\")\n",
    "pl.scatter(X,Y, c=GR, cmap=get_miub_cmap(), vmin=0, vmax=40, edgecolor=\"None\",s=50)\n",
    "pl.title(\"GR reflectivity\")\n",
    "pl.xlim(np.nanmin(X),np.nanmax(X))\n",
    "pl.ylim(np.nanmin(Y),np.nanmax(Y))\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(X,Y, c=SR-GR, cmap='jet',vmin=-5,vmax=5 ,edgecolor=\"None\",s=50)\n",
    "pl.title(\"SR - GR reflectivity\")\n",
    "pl.xlim(np.nanmin(X),np.nanmax(X))\n",
    "pl.ylim(np.nanmin(Y),np.nanmax(Y))\n",
    "pl.grid()\n",
    "pl.colorbar(shrink=0.5)\n",
    "#fig.suptitle(\"uncorrected SR - GR\")\n",
    "pl.tight_layout()\n",
    "pl.savefig('/automount/ftp/velibor/validation/matching5.png')\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": [
    "def get_my_cmap():\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 plot_cfad(x1, x2,vmaxi=None, bini=100, tit=None):\n",
    "    \"\"\"\n",
    "    Plot 4 CFAD\n",
    "    \"\"\"\n",
    "    import matplotlib.pyplot as plt\n",
    "    \n",
    "    m1 = ~np.isnan(x1) & ~np.isnan(x2)\n",
    "    \n",
    "    x3 = np.array([]); x4 = np.array([]); x5 = np.array([])\n",
    "\n",
    "    hh = np.array([])\n",
    "    for jj in range(0,20000,50):\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",
    "    \n",
    "    fig = plt.figure(figsize=(10,12))\n",
    "    ax = fig.add_subplot(1,1,1)\n",
    "    plt.hist2d(x1[m1], x2[m1],bins=bini, cmap=get_my_cmap(), vmin=0.1,vmax=vmaxi)\n",
    "    #plt.hist2d(x1[m1], x2[m1],bins=120, cmap='inferno')\n",
    "    #plt.hist2d(x1[m1], x2[m1], bins=111, cmap='inferno')\n",
    "    \n",
    "    cbar = plt.colorbar()\n",
    "    cbar.ax.set_ylabel('', fontsize=20, rotation=0)\n",
    "    cbar.ax.tick_params(labelsize=20) \n",
    "    plt.xlabel(r'$Z_H$ in DBz', fontsize=20)\n",
    "    plt.ylabel('Height in m', fontsize=20)\n",
    "    plt.title(tit, fontsize=20)\n",
    "    plt.xticks(fontsize=20)\n",
    "    plt.yticks(fontsize=20)\n",
    "    #plt.plot(x3,hh, linewidth=2, label='mean')\n",
    "    plt.plot(x4,hh, linewidth=3,color='black', label='median')\n",
    "    plt.ylim(0,7000)\n",
    "    plt.xlim(0,50)\n",
    "    \n",
    "    # Major ticks every 20, minor ticks every 5\n",
    "    major_xticks = np.arange(0, 50, 10)\n",
    "    minor_xticks = np.arange(0, 50, 5)\n",
    "    major_yticks = np.arange(0, 7000, 1000)\n",
    "    minor_yticks = np.arange(0, 7000, 200)\n",
    "\n",
    "    ax.set_xticks(major_xticks)\n",
    "    ax.set_xticks(minor_xticks, minor=True)\n",
    "    ax.set_yticks(major_yticks)\n",
    "    ax.set_yticks(minor_yticks, minor=True)\n",
    "\n",
    "    # And a corresponding grid\n",
    "    ax.grid(which='both')\n",
    "\n",
    "    # Or if you want different settings for the grids:\n",
    "    ax.grid(which='minor', alpha=0.5)\n",
    "    ax.grid(which='major', alpha=0.8)\n",
    "\n",
    "    \n",
    "    \n",
    "    #plt.legend(loc='upper right', fontsize=10)\n",
    "    #plt.grid(color='white')\n",
    "    #plt.savefig('/automount/ags/velibor/plot/GPM_boxpol/topevents/boxpol_cfad_'+ZP+'.png')\n",
    "    #plt.close()\n",
    "    #plt.savefig('/automount/ftp/velibor/validation/__irs_cfad.png')\n",
    "    \n",
    "    #plt.show()\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "ref_sr.shape\n",
    "z = 1000 * (np.arange(0,176,1)*0.125)\n",
    "h = np.ones([45,49,176])\n",
    "z_wm = h*z"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "cmap = pl.cm.jet\n",
    "cmaplist = [cmap(i) for i in range(cmap.N)]\n",
    "cmaplist[0] = cmap.from_list('Costom cmap', cmaplist, cmap.N)\n",
    "import matplotlib as mpl\n",
    "bounds = np.arange(-300,0,10)\n",
    "norm = mpl.colors.BoundaryNorm(bounds, cmap.N)\n",
    "\n",
    "zdiff = xyzp_sr[20,:,:,2]- z_wm[20,:,:]\n",
    "print (zdiff.shape)\n",
    "xxx, zzz = np.arange(-25,25,1), np.arange(0,176,1)\n",
    "print(xxx.shape, zzz.shape)\n",
    "fig =  pl.figure(figsize=(15,8))\n",
    "ax = fig.add_subplot(1,1,1)\n",
    "pl.pcolormesh(xxx,zzz,zdiff.T, cmap=cmap, norm=norm)#, vmin=-300, vmax=0)\n",
    "cbar = pl.colorbar()\n",
    "cbar.set_label('Difference in meters')\n",
    "pl.xlabel('Scan Angle')\n",
    "pl.ylabel('R bin for Height')\n",
    "pl.title('How big is the Error if we using the GPM (Ka) bins as 125m step array')\n",
    "\n",
    "\n",
    "\"\"\"# Major ticks every 20, minor ticks every 5\n",
    "ma_xticks = np.arange(-12, 12, 5)\n",
    "mi_xticks = np.arange(-12, 12, 1)\n",
    "\n",
    "ma_yticks = np.arange(0, 176, 10)\n",
    "mi_yticks = np.arange(0, 176, 1)\n",
    "\n",
    "\n",
    "ax.set_xticks(ma_xticks)\n",
    "ax.set_xticks(mi_xticks, minor=True)\n",
    "ax.set_yticks(ma_yticks)\n",
    "ax.set_yticks(mi_yticks, minor=True)\n",
    "\n",
    "# And a corresponding grid\n",
    "ax.grid(which='both')\n",
    "\n",
    "# Or if you want different settings for the grids:\n",
    "ax.grid(which='minor', alpha=0.5)\n",
    "ax.grid(which='major', alpha=1)\"\"\"\n",
    "\n",
    "pl.xlim(-25,24)\n",
    "pl.ylim(0,175)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plot_cfad(ref_sr,z_wm,tit='SR not matched', bini=30, vmaxi=220)\n",
    "pl.savefig('/automount/ftp/velibor/validation/cfadSR_notmatched_hight_without_corretion.png')\n",
    "\n",
    "pl.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plot_cfad(ref_sr,xyzp_sr[:,:,:,2],tit='SR not matched', bini=30, vmaxi=300)\n",
    "pl.savefig('/automount/ftp/velibor/validation/cfadSR_notmatched.png')\n",
    "\n",
    "pl.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plot_cfad(SR,Z,tit='SR', bini=30, vmaxi=30)\n",
    "pl.savefig('/automount/ftp/velibor/validation/cfadSR.png')\n",
    "\n",
    "pl.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plot_cfad(GR,Z,tit='GR', bini=30, vmaxi=30)\n",
    "pl.savefig('/automount/ftp/velibor/validation/CFADGR.png')\n",
    "pl.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "def plot_cfad_diff(x1, x2,vmaxi=None, bini=100, tit=None):\n",
    "    \"\"\"\n",
    "    Plot 4 CFAD\n",
    "    \"\"\"\n",
    "    import matplotlib.pyplot as plt\n",
    "    \n",
    "    m1 = ~np.isnan(x1) & ~np.isnan(x2)\n",
    "    \n",
    "    x3 = np.array([]); x4 = np.array([]); x5 = np.array([])\n",
    "\n",
    "    hh = np.array([])\n",
    "    for jj in range(0,20000,50):\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",
    "    \n",
    "    fig = plt.figure(figsize=(10,12))\n",
    "    ax = fig.add_subplot(1,1,1)\n",
    "    plt.hist2d(x1[m1], x2[m1],bins=bini, cmap=get_my_cmap(), vmin=0.1,vmax=vmaxi)\n",
    "    #plt.hist2d(x1[m1], x2[m1],bins=120, cmap='inferno')\n",
    "    #plt.hist2d(x1[m1], x2[m1], bins=111, cmap='inferno')\n",
    "    \n",
    "    cbar = plt.colorbar()\n",
    "    cbar.ax.set_ylabel('', fontsize=20, rotation=0)\n",
    "    cbar.ax.tick_params(labelsize=20) \n",
    "    plt.xlabel(r'$Z_H$ in DBz', fontsize=20)\n",
    "    plt.ylabel('Height in m', fontsize=20)\n",
    "    plt.title(tit, fontsize=20)\n",
    "    plt.xticks(fontsize=20)\n",
    "    plt.yticks(fontsize=20)\n",
    "    #plt.plot(x3,hh, linewidth=2, label='mean')\n",
    "    plt.plot(x4,hh, linewidth=3,color='black', label='median')\n",
    "    plt.ylim(0,7000)\n",
    "    plt.xlim(-20,20)\n",
    "    plt.axvline(0, color='black', lw=2, linestyle='--')\n",
    "    \n",
    "    # Major ticks every 20, minor ticks every 5\n",
    "    major_xticks = np.arange(-20,20, 10)\n",
    "    minor_xticks = np.arange(-20,20, 5)\n",
    "    major_yticks = np.arange(0, 7000, 1000)\n",
    "    minor_yticks = np.arange(0, 7000, 200)\n",
    "\n",
    "    ax.set_xticks(major_xticks)\n",
    "    ax.set_xticks(minor_xticks, minor=True)\n",
    "    ax.set_yticks(major_yticks)\n",
    "    ax.set_yticks(minor_yticks, minor=True)\n",
    "\n",
    "    # And a corresponding grid\n",
    "    ax.grid(which='both')\n",
    "\n",
    "    # Or if you want different settings for the grids:\n",
    "    ax.grid(which='minor', alpha=0.5)\n",
    "    ax.grid(which='major', alpha=0.8)\n",
    "\n",
    "    \n",
    "    \n",
    "    #plt.legend(loc='upper right', fontsize=10)\n",
    "    #plt.grid(color='white')\n",
    "    #plt.savefig('/automount/ags/velibor/plot/GPM_boxpol/topevents/boxpol_cfad_'+ZP+'.png')\n",
    "    #plt.close()\n",
    "    #plt.savefig('/automount/ftp/velibor/validation/__irs_cfad.png')\n",
    "    \n",
    "    #plt.show()\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "diffi = GR-SR\n",
    "maske = ~np.isnan(diffi) & ~np.isnan(Z)\n",
    "plot_cfad_diff(diffi[maske],Z[maske],tit='GR-SR', bini=30, vmaxi=30)\n",
    "\n",
    "pl.savefig('/automount/ftp/velibor/validation/CFADGR_SR.png')\n",
    "pl.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "pl.figure(figsize=(14,10))\n",
    "\n",
    "SRcal, GRcal = SR.copy(), GR.copy()\n",
    "\n",
    "SRcal[SRcal>30.]=np.nan\n",
    "GRcal[GRcal>30.]=np.nan\n",
    "SRcal[Z>2000.]=np.nan\n",
    "GRcal[Z>2000.]=np.nan\n",
    "\n",
    "mm = ~np.isnan(GRcal)& ~np.isnan(SRcal)\n",
    "\n",
    "\n",
    "pl.hist(GRcal[mm]-SRcal[mm],bins=100, facecolor=\"grey\", alpha=0.6)\n",
    "pl.axvline(np.mean(GRcal[mm]-SRcal[mm]), color='blue', label='Mean diff')\n",
    "pl.axvline(np.median(GRcal[mm]-SRcal[mm]), color='green', label='Median diff')\n",
    "pl.axvline(0, color='black',ls='--',lw=3)\n",
    "\n",
    "xmin, xmax = np.nanmin(GRcal[mm]-SRcal[mm]), np.nanmax(GRcal[mm]-SRcal[mm])\n",
    "limits = np.nanmax([abs(xmin),abs(xmax)])\n",
    "pl.xlim(-limits,limits)\n",
    "\n",
    "pl.title('Diff GR-SR, bias: '+str( np.sum(GRcal[mm]-SRcal[mm])/len(GRcal[mm]-SRcal[mm]))+' dbz')\n",
    "pl.grid()\n",
    "pl.legend()\n",
    "pl.xlabel('SR overestimate----------------------------Diff in dbz----------------------------GR overestimate')\n",
    "pl.ylabel('number')\n",
    "pl.savefig('/automount/ftp/velibor/validation/matching_bias.pdf')\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": [
    "\n",
    "mm = ~np.isnan(GR) & ~np.isnan(SR)\n",
    "#m2 = \n",
    "\n",
    "pl.figure(figsize=(14,5))\n",
    "pl.hist(SR[mm],bins=50, facecolor=\"None\", alpha=1, edgecolor='black', histtype='step', label='SR', normed=1)\n",
    "pl.xlim(left=0)\n",
    "\n",
    "\n",
    "pl.hist(GR[mm],bins=50, facecolor=\"None\", alpha=1, edgecolor='blue', histtype='step', label='GR', normed=1)\n",
    "pl.axvline(0, color='black',ls='--',lw=3)\n",
    "\n",
    "\n",
    "pl.grid()\n",
    "pl.legend(fontsize=20)\n",
    "#pl.ylim(top=90)\n",
    "pl.ylabel('P(Z)', fontsize=20)\n",
    "pl.xlabel(r'$Z_H$ in dbz', fontsize=20)\n",
    "pl.xticks(fontsize=20); pl.yticks(fontsize=20)\n",
    "pl.xlim(0,80); pl.ylim(0,0.1)\n",
    "pl.savefig('/automount/ftp/velibor/validation/matching_hist.pdf')\n",
    "pl.savefig('/automount/ftp/velibor/validation/matching_hist.png')\n",
    "pl.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "aa=0.9\n",
    "# Elevation unterteilen mit Cumulativer Summe\n",
    "e_pos = gr_nr_elev.cumsum().astype(int)\n",
    "\n",
    "pl.figure(figsize=(16,10))\n",
    "\n",
    "for eee in range(len(gr_nr_elev)):\n",
    "    pl.subplot(2,5,eee+1)\n",
    "    \n",
    "    if eee==0: \n",
    "        pl.scatter(SR[0:e_pos[0]],GR[0:e_pos[0]],alpha=aa,label='elevation: '+str(gr_ele[0])[0:5],color='black',s=10)\n",
    "        pl.title('r: ' + corcor(SR[0:e_pos[0]],GR[0:e_pos[0]])[0:4])\n",
    "        pl.legend(loc='upper left')\n",
    "    else:\n",
    "        pl.scatter(SR[e_pos[eee-1]:e_pos[eee]],GR[e_pos[eee-1]:e_pos[eee]],alpha=aa,\n",
    "                   label='elevation: ' + str(gr_ele[eee])[0:5],color='black',s=10)\n",
    "        pl.title('r: ' + corcor(SR[e_pos[eee-1]:e_pos[eee]],GR[e_pos[eee-1]:e_pos[eee]])[0:6])\n",
    "        pl.legend(loc='upper left')\n",
    "    pl.xlabel('SR in dbz')\n",
    "    pl.ylabel('GR in dbz')\n",
    "    pl.grid()\n",
    "    pl.xlim(0,45)\n",
    "    pl.ylim(0,45)\n",
    "    pl.plot([0,45],[0,45], color='black', ls='--')\n",
    "\n",
    "pl.tight_layout()\n",
    "pl.savefig('/automount/ftp/velibor/validation/corr_elev.pdf')\n",
    "pl.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Elevation unterteilen mit Cumulativer Summe\n",
    "e_pos = gr_nr_elev.cumsum().astype(int)\n",
    "from pcc import get_miub_cmap\n",
    "\n",
    "pl.figure(figsize=(30,10))\n",
    "\n",
    "for eee in range(len(gr_nr_elev)):\n",
    "    pl.subplot(2,5,eee+1)\n",
    "    \n",
    "    if eee==0: \n",
    "        pl.scatter(X[0:e_pos[0]],Y[0:e_pos[0]],c=GR[0:e_pos[0]],label=str(gr_ele[0]),\n",
    "                   cmap=get_miub_cmap(), vmin=0, vmax=40)\n",
    "        pl.legend(loc='upper left')\n",
    "    else:\n",
    "        pl.scatter(X[e_pos[eee-1]:e_pos[eee]],Y[e_pos[eee-1]:e_pos[eee]],c=GR[e_pos[eee-1]:e_pos[eee]],\n",
    "                   label=str(gr_ele[eee]), cmap=get_miub_cmap(), vmin=0, vmax=40)\n",
    "        pl.legend(loc='upper left')\n",
    "    pl.xlabel('X in meters')\n",
    "    pl.ylabel('Y in meters')\n",
    "    pl.grid()\n",
    "    pl.xlim(-150000,150000)\n",
    "    pl.ylim(-150000,150000)\n",
    "    pl.colorbar()\n",
    "\n",
    "\n",
    "pl.tight_layout()\n",
    "pl.savefig('/automount/ftp/velibor/validation/match_ppiGR.pdf')\n",
    "pl.savefig('/automount/ftp/velibor/validation/match_ppiGR.png')\n",
    "\n",
    "\n",
    "pl.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Elevation unterteilen mit Cumulativer Summe\n",
    "e_pos = gr_nr_elev.cumsum().astype(int)\n",
    "from pcc import get_miub_cmap\n",
    "\n",
    "pl.figure(figsize=(30,10))\n",
    "\n",
    "for eee in range(len(gr_nr_elev)):\n",
    "    pl.subplot(2,5,eee+1)\n",
    "    \n",
    "    if eee==0: \n",
    "        pl.scatter(X[0:e_pos[0]],Y[0:e_pos[0]],c=SR[0:e_pos[0]],label=str(gr_ele[0]),\n",
    "                   cmap=get_miub_cmap(), vmin=0, vmax=40)\n",
    "        pl.legend(loc='upper left')\n",
    "    else:\n",
    "        pl.scatter(X[e_pos[eee-1]:e_pos[eee]],Y[e_pos[eee-1]:e_pos[eee]],c=SR[e_pos[eee-1]:e_pos[eee]],\n",
    "                   label=str(gr_ele[eee]), cmap=get_miub_cmap(), vmin=0, vmax=40)\n",
    "        pl.legend(loc='upper left')\n",
    "    pl.xlabel('X in meters')\n",
    "    pl.ylabel('Y in meters')\n",
    "    pl.grid()\n",
    "    pl.xlim(-150000,150000)\n",
    "    pl.ylim(-150000,150000)\n",
    "    pl.colorbar()\n",
    "\n",
    "\n",
    "pl.tight_layout()\n",
    "pl.savefig('/automount/ftp/velibor/validation/match_ppiSR.pdf')\n",
    "pl.savefig('/automount/ftp/velibor/validation/match_ppiSR.png')\n",
    "\n",
    "\n",
    "pl.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Elevation unterteilen mit Cumulativer Summe\n",
    "e_pos = gr_nr_elev.cumsum().astype(int)\n",
    "from pcc import get_miub_cmap\n",
    "\n",
    "pl.figure(figsize=(30,10))\n",
    "\n",
    "for eee in range(len(gr_nr_elev)):\n",
    "    pl.subplot(2,5,eee+1)\n",
    "    \n",
    "    if eee==0: \n",
    "        pl.scatter(X[0:e_pos[0]],Z[0:e_pos[0]],c=SR[0:e_pos[0]],label=str(gr_ele[0]),\n",
    "                   cmap=get_miub_cmap(), vmin=0, vmax=40)\n",
    "        pl.legend(loc='upper left')\n",
    "    else:\n",
    "        pl.scatter(X[e_pos[eee-1]:e_pos[eee]],Z[e_pos[eee-1]:e_pos[eee]],c=SR[e_pos[eee-1]:e_pos[eee]],\n",
    "                   label=str(gr_ele[eee]), cmap=get_miub_cmap(), vmin=0, vmax=40)\n",
    "        pl.legend(loc='upper left')\n",
    "    pl.xlabel('SR in dbz')\n",
    "    pl.ylabel('GR in dbz')\n",
    "    pl.grid()\n",
    "    pl.xlim(-150000,150000)\n",
    "    pl.ylim(0,7000)\n",
    "    pl.colorbar()\n",
    "\n",
    "\n",
    "pl.tight_layout()\n",
    "pl.savefig('/automount/ftp/velibor/validation/match_XppiSR.png')\n",
    "\n",
    "pl.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Elevation unterteilen mit Cumulativer Summe\n",
    "e_pos = gr_nr_elev.cumsum().astype(int)\n",
    "from pcc import get_miub_cmap\n",
    "\n",
    "pl.figure(figsize=(30,10))\n",
    "\n",
    "for eee in range(len(gr_nr_elev)):\n",
    "    pl.subplot(2,5,eee+1)\n",
    "    \n",
    "    if eee==0: \n",
    "        pl.scatter(X[0:e_pos[0]],Z[0:e_pos[0]],c=GR[0:e_pos[0]],label=str(gr_ele[0]),\n",
    "                   cmap=get_miub_cmap(), vmin=0, vmax=40)\n",
    "        pl.legend(loc='upper left')\n",
    "    else:\n",
    "        pl.scatter(X[e_pos[eee-1]:e_pos[eee]],Z[e_pos[eee-1]:e_pos[eee]],c=GR[e_pos[eee-1]:e_pos[eee]],\n",
    "                   label=str(gr_ele[eee]), cmap=get_miub_cmap(), vmin=0, vmax=40)\n",
    "        pl.legend(loc='upper left')\n",
    "    pl.xlabel('SR in dbz')\n",
    "    pl.ylabel('GR in dbz')\n",
    "    pl.grid()\n",
    "    pl.xlim(-150000,150000)\n",
    "    pl.ylim(0,20000)\n",
    "    pl.colorbar()\n",
    "\n",
    "\n",
    "pl.tight_layout()\n",
    "pl.savefig('/automount/ftp/velibor/validation/match_XppiGR.pdf')\n",
    "\n",
    "pl.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Elevation unterteilen mit Cumulativer Summe\n",
    "e_pos = gr_nr_elev.cumsum().astype(int)\n",
    "from pcc import get_miub_cmap\n",
    "\n",
    "pl.figure(figsize=(30,10))\n",
    "\n",
    "for eee in range(len(gr_nr_elev)):\n",
    "    pl.subplot(2,5,eee+1)\n",
    "    \n",
    "    if eee==0: \n",
    "        pl.scatter(X[0:e_pos[0]],Y[0:e_pos[0]],c=SR[0:e_pos[0]]-GR[0:e_pos[0]],label=str(gr_ele[0]),\n",
    "                   cmap='jet', vmin=-10, vmax=10)\n",
    "        pl.legend(loc='upper left')\n",
    "    else:\n",
    "        pl.scatter(X[e_pos[eee-1]:e_pos[eee]],Y[e_pos[eee-1]:e_pos[eee]],c=SR[e_pos[eee-1]:e_pos[eee]]-GR[e_pos[eee-1]:e_pos[eee]],\n",
    "                   label=str(gr_ele[eee]), cmap='jet', vmin=-10, vmax=10)\n",
    "        pl.legend(loc='upper left')\n",
    "    pl.xlabel('X')\n",
    "    pl.ylabel('Y')\n",
    "    pl.title('SR-GR')\n",
    "    pl.grid()\n",
    "    pl.xlim(-150000,150000)\n",
    "    pl.ylim(-150000,150000)\n",
    "    pl.colorbar()\n",
    "\n",
    "\n",
    "pl.tight_layout()\n",
    "pl.savefig('/automount/ftp/velibor/validation/match_ppiSR.png')\n",
    "\n",
    "pl.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Elevation unterteilen mit Cumulativer Summe\n",
    "e_pos = gr_nr_elev.cumsum().astype(int)\n",
    "from pcc import get_miub_cmap\n",
    "\n",
    "pl.figure(figsize=(30,10))\n",
    "\n",
    "for eee in range(len(gr_nr_elev)):\n",
    "    pl.subplot(2,5,eee+1)\n",
    "    \n",
    "    if eee==0: \n",
    "        pl.scatter(X[0:e_pos[0]],Z[0:e_pos[0]],c=SR[0:e_pos[0]]-GR[0:e_pos[0]],label=str(gr_ele[0]),\n",
    "                   cmap='jet', vmin=-10, vmax=10)\n",
    "        pl.legend(loc='upper left')\n",
    "    else:\n",
    "        pl.scatter(X[e_pos[eee-1]:e_pos[eee]],Z[e_pos[eee-1]:e_pos[eee]],c=SR[e_pos[eee-1]:e_pos[eee]]-GR[e_pos[eee-1]:e_pos[eee]],\n",
    "                   label=str(gr_ele[eee]), cmap='jet', vmin=-10, vmax=10)\n",
    "        pl.legend(loc='upper left')\n",
    "    pl.xlabel('X')\n",
    "    pl.ylabel('Z')\n",
    "    pl.grid()\n",
    "    pl.xlim(-150000,150000)\n",
    "    #pl.ylim(0,7000)\n",
    "    pl.colorbar()\n",
    "\n",
    "\n",
    "pl.tight_layout()\n",
    "pl.savefig('/automount/ftp/velibor/validation/match_XppiSR.pdf')\n",
    "pl.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "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([])\n",
    "\n",
    "    hh = np.array([])\n",
    "    for jj in range(0,20000,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",
    "    \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",
    "\n",
    "    plt.legend(loc='upper right', fontsize=10)\n",
    "    plt.grid(color='white')\n",
    "    plt.show()\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "diff_SRGR = SR - GR\n",
    "pl.figure(figsize=(12,10)) \n",
    "\"\"\"pl.subplot(1,3,1)\n",
    "subplot_cfad(SR,Z,'SR')\n",
    "pl.ylim(0,12000)\n",
    "pl.subplot(1,3,2)\n",
    "subplot_cfad(GR,Z,'GR')\n",
    "pl.ylim(0,12000)\n",
    "pl.subplot(1,3,3)\"\"\"\n",
    "subplot_cfad(diff_SRGR,Z,'SR - GR')\n",
    "pl.ylim(0,12000)\n",
    "\n",
    "pl.tight_layout()\n",
    "\"\"\"pl.savefig('/automount/ftp/velibor/validation/cfad_SR_GR.pdf')\n",
    "pl.savefig('/automount/ftp/velibor/validation/cfad_SR_GR.png')\"\"\"\n",
    "\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": [
    "x1, x2 = 0, X\n",
    "y1, y2 = 0, Y\n",
    "\n",
    "dist = np.hypot(x2 - x1, y2 - y1)\n",
    "\n",
    "dist = dist/1000.\n",
    "\n",
    "aa=0.9\n",
    "# Elevation unterteilen mit Cumulativer Summe\n",
    "e_pos = gr_nr_elev.cumsum().astype(int)\n",
    "\n",
    "pl.figure(figsize=(20,10))\n",
    "\n",
    "for eee in range(len(gr_nr_elev)):\n",
    "    pl.subplot(2,5,eee+1)\n",
    "    \n",
    "    if eee==0: \n",
    "        sc = pl.scatter(SR[0:e_pos[0]],\n",
    "                   GR[0:e_pos[0]],\n",
    "                   c=dist[0:e_pos[0]],\n",
    "                   alpha=aa,label='elevation: '+str(gr_ele[0])[0:5],\n",
    "                   color='black',\n",
    "                   s=10,\n",
    "                   vmin=0.,\n",
    "                   vmax=90,\n",
    "                   marker='o')\n",
    "        sc.set_facecolor(\"none\")\n",
    "        pl.title('r: ' + corcor(SR[0:e_pos[0]],GR[0:e_pos[0]])[0:4])\n",
    "        pl.legend(loc='upper left')\n",
    "        cbar = pl.colorbar(shrink=0.7)\n",
    "        cbar.ax.set_ylabel('Distance from radar center in km', fontsize=10)\n",
    "        \n",
    "    else:\n",
    "        sc = pl.scatter(SR[e_pos[eee-1]:e_pos[eee]],\n",
    "                   GR[e_pos[eee-1]:e_pos[eee]],\n",
    "                   c=dist[e_pos[eee-1]:e_pos[eee]],\n",
    "                   vmin = 0.,\n",
    "                   vmax = 90.,\n",
    "                   alpha=aa,\n",
    "                   label='elevation: ' + str(gr_ele[eee])[0:5],\n",
    "                   color='black',\n",
    "                   s=10,\n",
    "                   marker='o')\n",
    "        sc.set_facecolor(\"none\")\n",
    "        pl.title('r: ' + corcor(SR[e_pos[eee-1]:e_pos[eee]],GR[e_pos[eee-1]:e_pos[eee]])[0:6])\n",
    "        pl.legend(loc='upper left')\n",
    "        cbar = pl.colorbar(shrink=0.7)\n",
    "        cbar.ax.set_ylabel('Distance from radar center in km', fontsize=10)\n",
    "        \n",
    "    pl.xlabel('SR in dbz')\n",
    "    pl.ylabel('GR in dbz')\n",
    "    pl.grid()\n",
    "    pl.xlim(15,45)\n",
    "    pl.ylim(15,45)\n",
    "    pl.plot([0,45],[0,45], color='black', ls='--')\n",
    "\n",
    "pl.tight_layout()\n",
    "pl.savefig('/automount/ftp/velibor/validation/corr_elev_dist.pdf')\n",
    "pl.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "x1, x2 = 0, X\n",
    "y1, y2 = 0, Y\n",
    "\n",
    "dist = np.hypot(x2 - x1, y2 - y1)\n",
    "\n",
    "dist = dist/1000.\n",
    "\n",
    "aa=0.9\n",
    "# Elevation unterteilen mit Cumulativer Summe\n",
    "e_pos = gr_nr_elev.cumsum().astype(int)\n",
    "\n",
    "pl.figure(figsize=(20,10))\n",
    "\n",
    "for eee in range(len(gr_nr_elev)):\n",
    "    pl.subplot(2,5,eee+1)\n",
    "    \n",
    "    if eee==0: \n",
    "        sc = pl.scatter(SR[0:e_pos[0]],\n",
    "                   GR[0:e_pos[0]],\n",
    "                   c=Z[0:e_pos[0]],\n",
    "                   alpha=aa,label='elevation: '+str(gr_ele[0])[0:5],\n",
    "                   color='black',\n",
    "                   s=10,\n",
    "                   vmin=0.,\n",
    "                   vmax=7000,\n",
    "                   marker='o')\n",
    "        sc.set_facecolor(\"none\")\n",
    "        pl.title('r: ' + corcor(SR[0:e_pos[0]],GR[0:e_pos[0]])[0:4])\n",
    "        pl.legend(loc='upper left')\n",
    "        cbar = pl.colorbar(shrink=0.7)\n",
    "        cbar.ax.set_ylabel('Hight in m', fontsize=10)\n",
    "        \n",
    "    else:\n",
    "        sc = pl.scatter(SR[e_pos[eee-1]:e_pos[eee]],\n",
    "                   GR[e_pos[eee-1]:e_pos[eee]],\n",
    "                   c=Z[e_pos[eee-1]:e_pos[eee]],\n",
    "                   vmin = 0.,\n",
    "                   vmax = 7000.,\n",
    "                   alpha=aa,\n",
    "                   label='elevation: ' + str(gr_ele[eee])[0:5],\n",
    "                   color='black',\n",
    "                   s=10,\n",
    "                   marker='o')\n",
    "        sc.set_facecolor(\"none\")\n",
    "        pl.title('r: ' + corcor(SR[e_pos[eee-1]:e_pos[eee]],GR[e_pos[eee-1]:e_pos[eee]])[0:6])\n",
    "        pl.legend(loc='upper left')\n",
    "        cbar = pl.colorbar(shrink=0.7)\n",
    "        cbar.ax.set_ylabel('Hight in m', fontsize=10)\n",
    "        \n",
    "    pl.xlabel('SR in dbz')\n",
    "    pl.ylabel('GR in dbz')\n",
    "    pl.grid()\n",
    "    pl.xlim(15,45)\n",
    "    pl.ylim(15,45)\n",
    "    pl.plot([0,45],[0,45], color='black', ls='--')\n",
    "\n",
    "pl.tight_layout()\n",
    "pl.savefig('/automount/ftp/velibor/validation/corr_elev_hight.pdf')\n",
    "pl.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "x1, x2 = 0, X\n",
    "y1, y2 = 0, Y\n",
    "\n",
    "dist = np.hypot(x2 - x1, y2 - y1)\n",
    "\n",
    "dist = dist/10000.\n",
    "\n",
    "aa=0.9\n",
    "# Elevation unterteilen mit Cumulativer Summe\n",
    "e_pos = gr_nr_elev.cumsum().astype(int)\n",
    "\n",
    "pl.figure(figsize=(20,10))\n",
    "\n",
    "for eee in range(len(gr_nr_elev)):\n",
    "    pl.subplot(2,5,eee+1)\n",
    "    \n",
    "    if eee==0: \n",
    "        sc = pl.scatter(SR[0:e_pos[0]],\n",
    "                   GR[0:e_pos[0]],\n",
    "                   c=Z[0:e_pos[0]],\n",
    "                   alpha=aa,label='elevation: '+str(gr_ele[0])[0:5],\n",
    "                   color='black',\n",
    "                   s=10*GRstd[0:e_pos[0]]/10,\n",
    "                   vmin=1.,\n",
    "                   #vmax=25,\n",
    "                   marker='o')\n",
    "        sc.set_facecolor(\"none\")\n",
    "        pl.title('r: ' + corcor(SR[0:e_pos[0]],GR[0:e_pos[0]])[0:4])\n",
    "        pl.legend(loc='upper left')\n",
    "        cbar = pl.colorbar(shrink=0.7)\n",
    "        cbar.ax.set_ylabel('Std in dbz', fontsize=10)\n",
    "        \n",
    "    else:\n",
    "        sc = pl.scatter(SR[e_pos[eee-1]:e_pos[eee]],\n",
    "                   GR[e_pos[eee-1]:e_pos[eee]],\n",
    "                   c=Z[e_pos[eee-1]:e_pos[eee]],\n",
    "                   vmin = 1.,\n",
    "                   #vmax = 25.,\n",
    "                   alpha=aa,\n",
    "                   label='elevation: ' + str(gr_ele[eee])[0:5],\n",
    "                   color='black',\n",
    "                   s=10*GRstd[e_pos[eee-1]:e_pos[eee]],\n",
    "                   marker='o')\n",
    "        sc.set_facecolor(\"none\")\n",
    "        pl.title('r: ' + corcor(SR[e_pos[eee-1]:e_pos[eee]],GR[e_pos[eee-1]:e_pos[eee]])[0:6])\n",
    "        pl.legend(loc='upper left')\n",
    "        cbar = pl.colorbar(shrink=0.7)\n",
    "        cbar.ax.set_ylabel('Std in dbz', fontsize=10)\n",
    "        \n",
    "    pl.xlabel('SR in dbz')\n",
    "    pl.ylabel('GR in dbz')\n",
    "    pl.grid()\n",
    "    pl.xlim(15,45)\n",
    "    pl.ylim(15,45)\n",
    "    pl.plot([0,45],[0,45], color='black', ls='--')\n",
    "\n",
    "pl.tight_layout()\n",
    "pl.savefig('/automount/ftp/velibor/validation/corr_elev_hight_std.pdf')\n",
    "pl.savefig('/automount/ftp/velibor/validation/corr_elev_hight_std.png')\n",
    "pl.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "x1, x2 = 0, X\n",
    "y1, y2 = 0, Y\n",
    "\n",
    "dist = np.hypot(x2 - x1, y2 - y1)\n",
    "\n",
    "dist = dist/10000.\n",
    "\n",
    "aa=0.9\n",
    "# Elevation unterteilen mit Cumulativer Summe\n",
    "e_pos = gr_nr_elev.cumsum().astype(int)\n",
    "\n",
    "pl.figure(figsize=(20,10))\n",
    "\n",
    "for eee in range(len(gr_nr_elev)):\n",
    "    pl.subplot(2,5,eee+1)\n",
    "    \n",
    "    if eee==0: \n",
    "        sc = pl.scatter(SR[0:e_pos[0]],\n",
    "                   GR[0:e_pos[0]],\n",
    "                   c=Z[0:e_pos[0]],\n",
    "                   alpha=aa,label='elevation: '+str(gr_ele[0])[0:5],\n",
    "                   color='black',\n",
    "                   s=GRnum[0:e_pos[0]]/10,\n",
    "                   vmin=0.,\n",
    "                   vmax=7000,\n",
    "                   marker='o')\n",
    "        sc.set_facecolor(\"none\")\n",
    "        pl.title('r: ' + corcor(SR[0:e_pos[0]],GR[0:e_pos[0]])[0:4])\n",
    "        pl.legend(loc='upper left')\n",
    "        cbar = pl.colorbar(shrink=0.7)\n",
    "        cbar.ax.set_ylabel('Hight in m', fontsize=10)\n",
    "        \n",
    "    else:\n",
    "        sc = pl.scatter(SR[e_pos[eee-1]:e_pos[eee]],\n",
    "                   GR[e_pos[eee-1]:e_pos[eee]],\n",
    "                   c=Z[e_pos[eee-1]:e_pos[eee]],\n",
    "                   vmin = 0.,\n",
    "                   vmax = 7000.,\n",
    "                   alpha=aa,\n",
    "                   label='elevation: ' + str(gr_ele[eee])[0:5],\n",
    "                   color='black',\n",
    "                   s=GRnum[e_pos[eee-1]:e_pos[eee]]/10,\n",
    "                   marker='o')\n",
    "        sc.set_facecolor(\"none\")\n",
    "        pl.title('r: ' + corcor(SR[e_pos[eee-1]:e_pos[eee]],GR[e_pos[eee-1]:e_pos[eee]])[0:6])\n",
    "        pl.legend(loc='upper left')\n",
    "        cbar = pl.colorbar(shrink=0.7)\n",
    "        cbar.ax.set_ylabel('Hight in m', fontsize=10)\n",
    "        \n",
    "    pl.xlabel('SR in dbz')\n",
    "    pl.ylabel('GR in dbz')\n",
    "    pl.grid()\n",
    "    pl.xlim(15,45)\n",
    "    pl.ylim(15,45)\n",
    "    pl.plot([0,45],[0,45], color='black', ls='--')\n",
    "\n",
    "pl.tight_layout()\n",
    "pl.savefig('/automount/ftp/velibor/validation/corr_elev_hight_num.pdf')\n",
    "pl.savefig('/automount/ftp/velibor/validation/corr_elev_hight_num.png')\n",
    "pl.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "x1, x2 = 0, X\n",
    "y1, y2 = 0, Y\n",
    "\n",
    "dist = np.hypot(x2 - x1, y2 - y1)\n",
    "\n",
    "dist = dist/10000.\n",
    "\n",
    "aa=0.9\n",
    "# Elevation unterteilen mit Cumulativer Summe\n",
    "e_pos = gr_nr_elev.cumsum().astype(int)\n",
    "\n",
    "pl.figure(figsize=(20,10))\n",
    "\n",
    "for eee in range(len(gr_nr_elev)):\n",
    "    pl.subplot(2,5,eee+1)\n",
    "    \n",
    "    if eee==0: \n",
    "        sc = pl.scatter(SR[0:e_pos[0]],\n",
    "                   GR[0:e_pos[0]],\n",
    "                   c=GRnum[0:e_pos[0]],\n",
    "                   alpha=aa,label='elevation: '+str(gr_ele[0])[0:5],\n",
    "                   color='black',\n",
    "                   s=10,\n",
    "                   vmin=0.1,\n",
    "                   vmax=500,\n",
    "                   marker='o')\n",
    "        sc.set_facecolor(\"none\")\n",
    "        pl.title('r: ' + corcor(SR[0:e_pos[0]],GR[0:e_pos[0]])[0:4])\n",
    "        pl.legend(loc='upper left')\n",
    "        cbar = pl.colorbar(shrink=0.7)\n",
    "        cbar.ax.set_ylabel('log(#)', fontsize=10)\n",
    "        \n",
    "    else:\n",
    "        sc = pl.scatter(SR[e_pos[eee-1]:e_pos[eee]],\n",
    "                   GR[e_pos[eee-1]:e_pos[eee]],\n",
    "                   c=GRnum[e_pos[eee-1]:e_pos[eee]],\n",
    "                   vmin = 0.1,\n",
    "                   vmax = 500.,\n",
    "                   alpha=aa,\n",
    "                   label='elevation: ' + str(gr_ele[eee])[0:5],\n",
    "                   color='black',\n",
    "                   s=10,\n",
    "                   marker='o')\n",
    "        sc.set_facecolor(\"none\")\n",
    "        pl.title('r: ' + corcor(SR[e_pos[eee-1]:e_pos[eee]],GR[e_pos[eee-1]:e_pos[eee]])[0:6])\n",
    "        pl.legend(loc='upper left')\n",
    "        cbar = pl.colorbar(shrink=0.7)\n",
    "        cbar.ax.set_ylabel('#', fontsize=10)\n",
    "        \n",
    "    pl.xlabel('SR in dbz')\n",
    "    pl.ylabel('GR in dbz')\n",
    "    pl.grid()\n",
    "    pl.xlim(15,45)\n",
    "    pl.ylim(15,45)\n",
    "    pl.plot([0,45],[0,45], color='black', ls='--')\n",
    "\n",
    "pl.tight_layout()\n",
    "pl.savefig('/automount/ftp/velibor/validation/corr_elev_num.pdf')\n",
    "pl.savefig('/automount/ftp/velibor/validation/corr_elev_num.png')\n",
    "pl.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "x1, x2 = 0, X\n",
    "y1, y2 = 0, Y\n",
    "\n",
    "dist = np.hypot(x2 - x1, y2 - y1)\n",
    "\n",
    "dist = dist/10000.\n",
    "\n",
    "aa=0.9\n",
    "# Elevation unterteilen mit Cumulativer Summe\n",
    "e_pos = gr_nr_elev.cumsum().astype(int)\n",
    "\n",
    "pl.figure(figsize=(20,10))\n",
    "\n",
    "for eee in range(len(gr_nr_elev)):\n",
    "    pl.subplot(2,5,eee+1)\n",
    "    \n",
    "    if eee==0: \n",
    "        sc = pl.scatter(SR[0:e_pos[0]],\n",
    "                   GR[0:e_pos[0]],\n",
    "                   c=GRstd[0:e_pos[0]],\n",
    "                   alpha=aa,label='elevation: '+str(gr_ele[0])[0:5],\n",
    "                   color='black',\n",
    "                   s=10,\n",
    "                   vmin=1.,\n",
    "                   vmax=5,\n",
    "                   marker='o')\n",
    "        sc.set_facecolor(\"none\")\n",
    "        pl.title('r: ' + corcor(SR[0:e_pos[0]],GR[0:e_pos[0]])[0:4])\n",
    "        pl.legend(loc='upper left')\n",
    "        cbar = pl.colorbar(shrink=0.7)\n",
    "        cbar.ax.set_ylabel('Std in dbz', fontsize=10)\n",
    "        \n",
    "    else:\n",
    "        sc = pl.scatter(SR[e_pos[eee-1]:e_pos[eee]],\n",
    "                   GR[e_pos[eee-1]:e_pos[eee]],\n",
    "                   c=GRstd[e_pos[eee-1]:e_pos[eee]],\n",
    "                   vmin = 1.,\n",
    "                   vmax = 5.,\n",
    "                   alpha=aa,\n",
    "                   label='elevation: ' + str(gr_ele[eee])[0:5],\n",
    "                   color='black',\n",
    "                   s=10,\n",
    "                   marker='o')\n",
    "        sc.set_facecolor(\"none\")\n",
    "        pl.title('r: ' + corcor(SR[e_pos[eee-1]:e_pos[eee]],GR[e_pos[eee-1]:e_pos[eee]])[0:6])\n",
    "        pl.legend(loc='upper left')\n",
    "        cbar = pl.colorbar(shrink=0.7)\n",
    "        cbar.ax.set_ylabel('Std in dbz', fontsize=10)\n",
    "        \n",
    "    pl.xlabel('SR in dbz')\n",
    "    pl.ylabel('GR in dbz')\n",
    "    pl.grid()\n",
    "    pl.xlim(15,45)\n",
    "    pl.ylim(15,45)\n",
    "    pl.plot([0,45],[0,45], color='black', ls='--')\n",
    "\n",
    "pl.tight_layout()\n",
    "pl.savefig('/automount/ftp/velibor/validation/corr_elev_std.pdf')\n",
    "pl.savefig('/automount/ftp/velibor/validation/corr_elev_std.png')\n",
    "pl.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Elevation unterteilen mit Cumulativer Summe\n",
    "e_pos = gr_nr_elev.cumsum().astype(int)\n",
    "from pcc import get_miub_cmap\n",
    "\n",
    "pl.figure(figsize=(30,10))\n",
    "\n",
    "for eee in range(len(gr_nr_elev)):\n",
    "    pl.subplot(2,5,eee+1)\n",
    "    \n",
    "    if eee==0: \n",
    "        pl.scatter(X[0:e_pos[0]],Z[0:e_pos[0]],c=GRnum[0:e_pos[0]],label=str(gr_ele[0])[0:4],\n",
    "                   cmap='jet', vmin=0, vmax=500)\n",
    "        pl.legend(loc='upper left')\n",
    "    else:\n",
    "        pl.scatter(X[e_pos[eee-1]:e_pos[eee]],Z[e_pos[eee-1]:e_pos[eee]],c=GRnum[e_pos[eee-1]:e_pos[eee]],\n",
    "                   label=str(gr_ele[eee])[0:4], cmap='jet', vmin=0, vmax=500)\n",
    "        pl.legend(loc='upper left')\n",
    "    pl.xlabel('X')\n",
    "    pl.ylabel('Z')\n",
    "    pl.grid()\n",
    "    pl.xlim(-150000,150000)\n",
    "    #pl.ylim(0,8000)\n",
    "    pl.colorbar()\n",
    "\n",
    "\n",
    "pl.tight_layout()\n",
    "pl.savefig('/automount/ftp/velibor/validation/GRnum.pdf')\n",
    "pl.savefig('/automount/ftp/velibor/validation/GRnum.png')\n",
    "pl.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "\"\"\"time_now = time.strftime(\"%Y%m%d--%H%M%S\")\n",
    "match_run = 'match_' + time_now\n",
    "savepfad = '/automount/ags/velibor/data/matching/'\n",
    "os.mkdir(savepfad+match_run)\n",
    "\n",
    "names = np.array(['X','Y','Z','SR','GR'])\n",
    "np.save(savepfad + match_run+ '/match' + sr_data['date'][1].strftime(\"%Y%m%d--%H%M\")+'.npy',\n",
    "       [X,\n",
    "       Y,\n",
    "       Z,\n",
    "       SR,\n",
    "       GR,\n",
    "       names])\n",
    "\n",
    "!jupyter nbconvert --output-dir={savepfad+match_run} --to html pcc_gpm_matching_all_elevation.ipynb\"\"\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "gr_mask = ~np.isnan(GRstd)\n",
    "pl.hist(GRstd[gr_mask])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import matplotlib.colors\n",
    "\n",
    "n = 10 # how many lines to draw or number of discrete color levels\n",
    "cmap = pl.get_cmap(\"jet\", n)\n",
    "\n",
    "pl.scatter(SR,GR,c=Z)\n",
    "\n",
    "\n",
    "norm= matplotlib.colors.BoundaryNorm(np.arange(0,n+1)-0.5, n)\n",
    "sm = pl.cm.ScalarMappable(cmap=cmap, norm=norm)\n",
    "sm.set_array([])\n",
    "pl.colorbar(sm, ticks=np.arange(0,n))\n",
    "\n",
    "pl.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# setup the plot\n",
    "fig, ax = pl.subplots(1,1, figsize=(12,6))\n",
    "\n",
    "# define the data\n",
    "x = SR\n",
    "y = GR\n",
    "tag = Z\n",
    "#tag[10:12] = 0 # make sure there are some 0 values to showup as grey\n",
    "\n",
    "# define the colormap\n",
    "cmap = pl.cm.jet\n",
    "# extract all colors from the .jet map\n",
    "cmaplist = [cmap(i) for i in range(cmap.N)]\n",
    "# force the first color entry to be grey\n",
    "cmaplist[0] = (.5,.5,.5,1.0)\n",
    "# create the new map\n",
    "cmap = cmap.from_list('Custom cmap', cmaplist, cmap.N)\n",
    "\n",
    "# define the bins and normalize\n",
    "#bounds = np.linspace(0,5000,10)\n",
    "bounds = np.array([0,1000,2000,3000,20000])\n",
    "norm = mpl.colors.BoundaryNorm(bounds, cmap.N)\n",
    "\n",
    "# make the scatter\n",
    "scat = ax.scatter(x,y,c=tag,cmap=cmap, norm=norm)\n",
    "# create a second axes for the colorbar\n",
    "ax2 = fig.add_axes([0.95, 0.1, 0.03, 0.8])\n",
    "cb = mpl.colorbar.ColorbarBase(ax2, cmap=cmap, norm=norm, ticks=bounds, boundaries=bounds, format='%1i')\n",
    "\n",
    "ax.set_title(' ')\n",
    "ax2.set_ylabel('z in m', size=12)\n",
    "ax.grid()\n",
    "pl.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "ZZ = Z.copy()\n",
    "SRR = SR.copy()\n",
    "GRR = GR.copy()\n",
    "GN = GRnum.copy()\n",
    "\n",
    "_num = 150\n",
    "_h1, _h2 = 2000, 3000\n",
    "\n",
    "\n",
    "# BB entfernen\n",
    "\"\"\"ZZ[(Z>=_h1) & (Z<=_h2)] = np.nan\n",
    "SRR[(Z>=_h1) & (Z<=_h2)] = np.nan\n",
    "GRR[(Z>=_h1) & (Z<=_h2)] = np.nan\n",
    "GN[(Z>=_h1) & (Z<=_h2)] = np.nan\"\"\"\n",
    "\n",
    "# NUM entfernen\n",
    "#ZZ[(GRnum < _num) ] = np.nan\n",
    "#SRR[(GRnum < _num) ] = np.nan\n",
    "#GRR[(GRnum < _num)] = np.nan\n",
    "#GN[(GRnum < _num) ] = np.nan\n",
    "\n",
    "cmap = pl.cm.Set1\n",
    "# extract all colors from the .jet map\n",
    "cmaplist = [cmap(i) for i in range(cmap.N)]\n",
    "# force the first color entry to be grey\n",
    "cmaplist[0] = (.5,.5,.5,1.0)\n",
    "# create the new map\n",
    "cmap = cmap.from_list('Custom cmap', cmaplist, cmap.N)\n",
    "\n",
    "# define the bins and normalize\n",
    "#bounds = np.linspace(0,5000,10)\n",
    "bounds = np.arange(0,7000,1000)\n",
    "norm = mpl.colors.BoundaryNorm(bounds, cmap.N)\n",
    "\n",
    "aa=0.9\n",
    "# Elevation unterteilen mit Cumulativer Summe\n",
    "e_pos = gr_nr_elev.cumsum().astype(int)\n",
    "\n",
    "pl.figure(figsize=(20,10))\n",
    "\n",
    "for eee in range(len(gr_nr_elev)):\n",
    "    pl.subplot(2,5,eee+1)\n",
    "    \n",
    "    if eee==0: \n",
    "        sc = pl.scatter(SRR[0:e_pos[0]],\n",
    "                   GRR[0:e_pos[0]],\n",
    "                   c=ZZ[0:e_pos[0]],\n",
    "                   alpha=aa,label='elevation: '+str(gr_ele[0])[0:5],\n",
    "                   color='black',\n",
    "                   s=GN[0:e_pos[0]]/10,\n",
    "                   vmin=0.,\n",
    "                   vmax=7000,\n",
    "                   marker='o',\n",
    "                   cmap=cmap, norm=norm)\n",
    "        sc.set_facecolor(\"none\")\n",
    "        pl.title('r: ' + corcor(SRR[0:e_pos[0]],GRR[0:e_pos[0]])[0:4])\n",
    "        pl.legend(loc='upper left')\n",
    "        cbar = pl.colorbar(shrink=0.7, cmap=cmap, norm=norm, ticks=bounds, boundaries=bounds)\n",
    "        cbar.ax.set_ylabel('Hight in m', fontsize=10)\n",
    "        \n",
    "    else:\n",
    "        sc = pl.scatter(SRR[e_pos[eee-1]:e_pos[eee]],\n",
    "                   GRR[e_pos[eee-1]:e_pos[eee]],\n",
    "                   c=ZZ[e_pos[eee-1]:e_pos[eee]],\n",
    "                   vmin = 0.,\n",
    "                   vmax = 7000.,\n",
    "                   alpha=aa,\n",
    "                   label='elevation: ' + str(gr_ele[eee])[0:5],\n",
    "                   color='black',\n",
    "                   s=GN[e_pos[eee-1]:e_pos[eee]]/10,\n",
    "                   marker='o', \n",
    "                   cmap=cmap, norm=norm)\n",
    "        sc.set_facecolor(\"none\")\n",
    "        pl.title('r: ' + corcor(SRR[e_pos[eee-1]:e_pos[eee]],GRR[e_pos[eee-1]:e_pos[eee]])[0:6])\n",
    "        pl.legend(loc='upper left')\n",
    "        cbar = pl.colorbar(shrink=0.7, cmap=cmap, norm=norm, ticks=bounds, boundaries=bounds)\n",
    "        cbar.ax.set_ylabel('Hight in m', fontsize=10)\n",
    "        \n",
    "    pl.xlabel('SR in dbz')\n",
    "    pl.ylabel('GR in dbz')\n",
    "    pl.grid()\n",
    "    pl.xlim(15,45)\n",
    "    pl.ylim(15,45)\n",
    "    pl.plot([0,45],[0,45], color='black', ls='--')\n",
    "\n",
    "pl.tight_layout()\n",
    "pl.savefig('/automount/ftp/velibor/validation/new_corr_elev_hight_num_withBB.pdf')\n",
    "pl.savefig('/automount/ftp/velibor/validation/new_corr_elev_hight_num_withBB.png')\n",
    "pl.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%%time\n",
    "ZZ = Z.copy()\n",
    "SRR = SR.copy()\n",
    "GRR = GR.copy()\n",
    "GN = GRnum.copy()\n",
    "\n",
    "_num = 150\n",
    "_h1, _h2 = 2000, 3000\n",
    "\n",
    "\n",
    "# BB entfernen\n",
    "ZZ[(Z>=_h1) & (Z<=_h2)] = np.nan\n",
    "SRR[(Z>=_h1) & (Z<=_h2)] = np.nan\n",
    "GRR[(Z>=_h1) & (Z<=_h2)] = np.nan\n",
    "GN[(Z>=_h1) & (Z<=_h2)] = np.nan\n",
    "\n",
    "# NUM entfernen\n",
    "#ZZ[(GRnum < _num) ] = np.nan\n",
    "#SRR[(GRnum < _num) ] = np.nan\n",
    "#GRR[(GRnum < _num)] = np.nan\n",
    "#GN[(GRnum < _num) ] = np.nan\n",
    "\n",
    "cmap = pl.cm.jet\n",
    "# extract all colors from the .jet map\n",
    "cmaplist = [cmap(i) for i in range(cmap.N)]\n",
    "# force the first color entry to be grey\n",
    "cmaplist[0] = (.5,.5,.5,1.0)\n",
    "# create the new map\n",
    "cmap = cmap.from_list('Custom cmap', cmaplist, cmap.N)\n",
    "\n",
    "# define the bins and normalize\n",
    "#bounds = np.linspace(0,5000,10)\n",
    "bounds = np.arange(0,7000,1000)\n",
    "norm = mpl.colors.BoundaryNorm(bounds, cmap.N)\n",
    "\n",
    "aa=0.9\n",
    "# Elevation unterteilen mit Cumulativer Summe\n",
    "e_pos = gr_nr_elev.cumsum().astype(int)\n",
    "\n",
    "pl.figure(figsize=(20,10))\n",
    "\n",
    "for eee in range(len(gr_nr_elev)):\n",
    "    pl.subplot(2,5,eee+1)\n",
    "    \n",
    "    if eee==0: \n",
    "        sc = pl.scatter(SRR[0:e_pos[0]],\n",
    "                   GRR[0:e_pos[0]],\n",
    "                   c=ZZ[0:e_pos[0]],\n",
    "                   alpha=aa,label='elevation: '+str(gr_ele[0])[0:5],\n",
    "                   color='black',\n",
    "                   s=GN[0:e_pos[0]]/10,\n",
    "                   vmin=0.,\n",
    "                   vmax=7000,\n",
    "                   marker='o',\n",
    "                   cmap=cmap, norm=norm)\n",
    "        sc.set_facecolor(\"none\")\n",
    "        pl.title('r: ' + corcor(SRR[0:e_pos[0]],GRR[0:e_pos[0]])[0:4])\n",
    "        pl.legend(loc='upper left')\n",
    "        cbar = pl.colorbar(shrink=0.7, cmap=cmap, norm=norm, ticks=bounds, boundaries=bounds)\n",
    "        cbar.ax.set_ylabel('Hight in m', fontsize=10)\n",
    "        \n",
    "    else:\n",
    "        sc = pl.scatter(SRR[e_pos[eee-1]:e_pos[eee]],\n",
    "                   GRR[e_pos[eee-1]:e_pos[eee]],\n",
    "                   c=ZZ[e_pos[eee-1]:e_pos[eee]],\n",
    "                   vmin = 0.,\n",
    "                   vmax = 7000.,\n",
    "                   alpha=aa,\n",
    "                   label='elevation: ' + str(gr_ele[eee])[0:5],\n",
    "                   color='black',\n",
    "                   s=GN[e_pos[eee-1]:e_pos[eee]]/10,\n",
    "                   marker='o', \n",
    "                   cmap=cmap, norm=norm)\n",
    "        sc.set_facecolor(\"none\")\n",
    "        pl.title('r: ' + corcor(SRR[e_pos[eee-1]:e_pos[eee]],GRR[e_pos[eee-1]:e_pos[eee]])[0:6])\n",
    "        pl.legend(loc='upper left')\n",
    "        cbar = pl.colorbar(shrink=0.7, cmap=cmap, norm=norm, ticks=bounds, boundaries=bounds)\n",
    "        cbar.ax.set_ylabel('Hight in m', fontsize=10)\n",
    "        \n",
    "    pl.xlabel('SR in dbz')\n",
    "    pl.ylabel('GR in dbz')\n",
    "    pl.grid()\n",
    "    pl.xlim(15,45)\n",
    "    pl.ylim(15,45)\n",
    "    pl.plot([0,45],[0,45], color='black', ls='--')\n",
    "\n",
    "pl.tight_layout()\n",
    "pl.savefig('/automount/ftp/velibor/validation/new_corr_elev_hight_num_withoutBB.pdf')\n",
    "pl.savefig('/automount/ftp/velibor/validation/new_corr_elev_hight_num_withoutBB.png')\n",
    "pl.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "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",
    "pl.figure(figsize=(30,10))\n",
    "\n",
    "for eee in range(len(gr_nr_elev)):\n",
    "    pl.subplot(2,5,eee+1)\n",
    "    \n",
    "    if eee==0: \n",
    "        pl.scatter(X[0:e_pos[0]],Y[0:e_pos[0]],c=Z[0:e_pos[0]],label=str(gr_ele[0])[0:4],\n",
    "                   cmap='jet', \n",
    "                   vmin=0,\n",
    "                   vmax=1000)\n",
    "        pl.legend(loc='upper left')\n",
    "    else:\n",
    "        pl.scatter(X[e_pos[eee-1]:e_pos[eee]],Y[e_pos[eee-1]:e_pos[eee]],c=Z[e_pos[eee-1]:e_pos[eee]],\n",
    "                   label=str(gr_ele[eee])[0:4], cmap='jet',\n",
    "                   vmin=0,\n",
    "                   vmax=1000)\n",
    "        print(np.nanmin(Z[e_pos[eee-1]:e_pos[eee]]))\n",
    "        pl.legend(loc='upper left')\n",
    "    pl.xlabel('X')\n",
    "    pl.ylabel('Y')\n",
    "    pl.grid()\n",
    "    pl.xlim(-150000,150000)\n",
    "    pl.ylim(-150000,150000)\n",
    "    #pl.ylim(0,8000)\n",
    "    pl.colorbar()\n",
    "\n",
    "\n",
    "pl.tight_layout()\n",
    "pl.savefig('/automount/ftp/velibor/validation/SR_hights.pdf')\n",
    "pl.savefig('/automount/ftp/velibor/validation/SR_hights.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": false
   },
   "outputs": [],
   "source": [
    "Zh = Z.copy()\n",
    "SRh = SR.copy()\n",
    "GRh = GR.copy()\n",
    "Gh = GRnum.copy()\n",
    "\n",
    "SR_top = SRh[Zh>2800]\n",
    "SR_bb = SRh[(Zh<2800) & (Zh>2200)]\n",
    "SR_bot = SRh[Zh<2200]\n",
    "\n",
    "GR_top = GRh[Zh>2800]\n",
    "GR_bb = GRh[(Zh<2800) & (Zh>2200)]\n",
    "GR_bot = GRh[Zh<2200]\n",
    "\n",
    "mm_top = ~np.isnan(GR_top) & ~np.isnan(SR_top)\n",
    "mm_bb = ~np.isnan(GR_bb) & ~np.isnan(SR_bb)\n",
    "mm_bot = ~np.isnan(GR_bot) & ~np.isnan(SR_bot)\n",
    "\n",
    "\n",
    "pl.figure(figsize=(16,8))\n",
    "pl.subplot(3,1,1)\n",
    "pl.hist(SR_top[mm_top],bins=50, facecolor=\"None\", alpha=1, edgecolor='black', histtype='step', label='SR above BB', normed=1,ls='-')\n",
    "pl.hist(GR_top[mm_top],bins=50, facecolor=\"None\", alpha=1, edgecolor='blue', histtype='step', label='GR above BB', normed=1,ls='-')\n",
    "pl.grid()\n",
    "pl.legend(fontsize=20)\n",
    "#pl.ylim(top=90)\n",
    "pl.ylabel('P(Z)', fontsize=20)\n",
    "pl.xlabel(r'$Z_H$ in dbz', fontsize=20)\n",
    "pl.xticks(fontsize=20); pl.yticks(fontsize=20)\n",
    "pl.xlim(0,60); pl.ylim(0,0.2)\n",
    "\n",
    "pl.subplot(3,1,2)\n",
    "pl.hist(SR_bb[mm_bb],bins=50, facecolor=\"None\", alpha=1, edgecolor='black', histtype='step', label='SR in BB', normed=1,ls='-')\n",
    "pl.hist(GR_bb[mm_bb],bins=50, facecolor=\"None\", alpha=1, edgecolor='blue', histtype='step', label='GR in BB', normed=1,ls='-')\n",
    "pl.grid()\n",
    "pl.legend(fontsize=20)\n",
    "#pl.ylim(top=90)\n",
    "pl.ylabel('P(Z)', fontsize=20)\n",
    "pl.xlabel(r'$Z_H$ in dbz', fontsize=20)\n",
    "pl.xticks(fontsize=20); pl.yticks(fontsize=20)\n",
    "pl.xlim(0,60); pl.ylim(0,0.2)\n",
    "\n",
    "pl.subplot(3,1,3)\n",
    "pl.hist(SR_bot[mm_bot],bins=50, facecolor=\"None\", alpha=1, edgecolor='black', histtype='step', label='SR below BB', normed=1,ls='-')\n",
    "pl.hist(GR_bot[mm_bot],bins=50, facecolor=\"None\", alpha=1, edgecolor='blue', histtype='step', label='GR below BB', normed=1,ls='-')\n",
    "pl.grid()\n",
    "pl.legend(fontsize=20)\n",
    "#pl.ylim(top=90)\n",
    "pl.ylabel('P(Z)', fontsize=20)\n",
    "pl.xlabel(r'$Z_H$ in dbz', fontsize=20)\n",
    "pl.xticks(fontsize=20); pl.yticks(fontsize=20)\n",
    "pl.xlim(0,60); pl.ylim(0,0.2)\n",
    "\n",
    "\n",
    "pl.tight_layout()\n",
    "pl.savefig('/automount/ftp/velibor/validation/matching_hist.pdf')\n",
    "pl.savefig('/automount/ftp/velibor/validation/matching_hist.png')\n",
    "pl.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "len(GR_bot[mm_bot])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def stats(a, b):\n",
    "    \"\"\"beschriften und Dict!\"\"\"\n",
    "    mask = ~np.isnan(a) & ~np.isnan(b)\n",
    "    H = len(a[mask])\n",
    "    bias = np.nansum(a - b)/H\n",
    "    rmse = np.sqrt(np.nansum(((a -b )**2.0)/H))\n",
    "    \n",
    "    return bias, rmse\n",
    "    \n",
    "    \n",
    "    \n",
    "\n",
    "ZZ = Z.copy()\n",
    "SRR = SR.copy()\n",
    "GRR = GR.copy()\n",
    "GN = GRnum.copy()\n",
    "\n",
    "_num = 150\n",
    "_h1, _h2 = 2200, 2800\n",
    "\n",
    "\n",
    "# BB entfernen\n",
    "\"\"\"ZZ[(Z>=_h1) & (Z<=_h2)] = np.nan\n",
    "SRR[(Z>=_h1) & (Z<=_h2)] = np.nan\n",
    "GRR[(Z>=_h1) & (Z<=_h2)] = np.nan\n",
    "GN[(Z>=_h1) & (Z<=_h2)] = np.nan\"\"\"\n",
    "\n",
    "# NUM entfernen\n",
    "#ZZ[(GRnum < _num) ] = np.nan\n",
    "#SRR[(GRnum < _num) ] = np.nan\n",
    "#GRR[(GRnum < _num)] = np.nan\n",
    "#GN[(GRnum < _num) ] = np.nan\n",
    "\n",
    "cmap = pl.cm.Set1\n",
    "# extract all colors from the .jet map\n",
    "cmaplist = [cmap(i) for i in range(cmap.N)]\n",
    "# force the first color entry to be grey\n",
    "cmaplist[0] = (.5,.5,.5,1.0)\n",
    "# create the new map\n",
    "cmap = cmap.from_list('Custom cmap', cmaplist, cmap.N)\n",
    "\n",
    "# define the bins and normalize\n",
    "#bounds = np.linspace(0,5000,10)\n",
    "bounds = np.arange(0,7000,1000)\n",
    "norm = mpl.colors.BoundaryNorm(bounds, cmap.N)\n",
    "\n",
    "aa=0.9\n",
    "# Elevation unterteilen mit Cumulativer Summe\n",
    "e_pos = gr_nr_elev.cumsum().astype(int)\n",
    "\n",
    "pl.figure(figsize=(26,8))\n",
    "\n",
    "for eee in range(len(gr_nr_elev)):\n",
    "    pl.subplot(2,5,eee+1)\n",
    "    \n",
    "    if eee==0: \n",
    "        sc = pl.scatter(SRR[0:e_pos[0]],\n",
    "                   GRR[0:e_pos[0]],\n",
    "                   c=ZZ[0:e_pos[0]],\n",
    "                   alpha=aa,label='elevation: '+str(gr_ele[0])[0:5],\n",
    "                   color='black',\n",
    "                   s=GN[0:e_pos[0]]/10,\n",
    "                   vmin=0.,\n",
    "                   vmax=7000,\n",
    "                   marker='o',\n",
    "                   cmap=cmap, norm=norm)\n",
    "        sc.set_facecolor(\"none\")\n",
    "        pl.title('r: ' + corcor(SRR[0:e_pos[0]],GRR[0:e_pos[0]])[0:4]+\n",
    "                '\\n bias: ' +  str(stats(SRR[0:e_pos[0]],GRR[0:e_pos[0]])[0])[0:4]+\n",
    "                 '\\n rmse: ' + str(stats(SRR[0:e_pos[0]],GRR[0:e_pos[0]])[1])[0:4]\n",
    "                )\n",
    "        pl.legend(loc='upper left',markerscale=0)\n",
    "        cbar = pl.colorbar(shrink=0.7, cmap=cmap, norm=norm, ticks=bounds, boundaries=bounds)\n",
    "        cbar.ax.set_ylabel('Hight in m', fontsize=10)\n",
    "        \n",
    "    else:\n",
    "        sc = pl.scatter(SRR[e_pos[eee-1]:e_pos[eee]],\n",
    "                   GRR[e_pos[eee-1]:e_pos[eee]],\n",
    "                   c=ZZ[e_pos[eee-1]:e_pos[eee]],\n",
    "                   vmin = 0.,\n",
    "                   vmax = 7000.,\n",
    "                   alpha=aa,\n",
    "                   label='elevation: ' + str(gr_ele[eee])[0:5],\n",
    "                   color='black',\n",
    "                   s=GN[e_pos[eee-1]:e_pos[eee]]/10,\n",
    "                   marker='o', \n",
    "                   cmap=cmap, norm=norm)\n",
    "        sc.set_facecolor(\"none\")\n",
    "        pl.title('r: ' + corcor(SRR[e_pos[eee-1]:e_pos[eee]],GRR[e_pos[eee-1]:e_pos[eee]])[0:4]+\n",
    "                '\\n bias: ' +  str(stats(SRR[e_pos[eee-1]:e_pos[eee]],GRR[e_pos[eee-1]:e_pos[eee]])[0])[0:4] +\n",
    "                 '\\n rmse: ' + str(stats(SRR[e_pos[eee-1]:e_pos[eee]],GRR[e_pos[eee-1]:e_pos[eee]])[1])[0:4]\n",
    "                )\n",
    "        pl.legend(loc='upper left',markerscale=0)\n",
    "        cbar = pl.colorbar(shrink=0.7, cmap=cmap, norm=norm, ticks=bounds, boundaries=bounds)\n",
    "        cbar.ax.set_ylabel('Hight in m', fontsize=10)\n",
    "        \n",
    "    pl.xlabel('SR in dbz')\n",
    "    pl.ylabel('GR in dbz')\n",
    "    pl.grid()\n",
    "    pl.xlim(15,45)\n",
    "    pl.ylim(15,45)\n",
    "    pl.plot([0,45],[0,45], color='black', ls='--')\n",
    "\n",
    "pl.tight_layout()\n",
    "pl.savefig('/automount/ftp/velibor/validation/new_corr_elev_hight_num_BB_stats.pdf')\n",
    "pl.savefig('/automount/ftp/velibor/validation/new_corr_elev_hight_num_BB_stats.png')\n",
    "pl.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def stats(a, b):\n",
    "    \"\"\"beschriften und Dict!\"\"\"\n",
    "    mask = ~np.isnan(a) & ~np.isnan(b)\n",
    "    H = len(a[mask])\n",
    "    bias = np.nansum(a - b)/H\n",
    "    rmse = np.sqrt(np.nansum(((a -b )**2.0)/H))\n",
    "    \n",
    "    return bias, rmse\n",
    "    \n",
    "ff = 15    \n",
    "ff2 = 20    \n",
    "\n",
    "ZZ = Z.copy()\n",
    "SRR = SR.copy()\n",
    "GRR = GR.copy()\n",
    "GN = GRnum.copy()\n",
    "\n",
    "_num = 150\n",
    "_h1, _h2 = 2200, 2800\n",
    "\n",
    "\n",
    "# BB entfernen\n",
    "\"\"\"ZZ[(Z>=_h1) & (Z<=_h2)] = np.nan\n",
    "SRR[(Z>=_h1) & (Z<=_h2)] = np.nan\n",
    "GRR[(Z>=_h1) & (Z<=_h2)] = np.nan\n",
    "GN[(Z>=_h1) & (Z<=_h2)] = np.nan\"\"\"\n",
    "\n",
    "# NUM entfernen\n",
    "#ZZ[(GRnum < _num) ] = np.nan\n",
    "#SRR[(GRnum < _num) ] = np.nan\n",
    "#GRR[(GRnum < _num)] = np.nan\n",
    "#GN[(GRnum < _num) ] = np.nan\n",
    "\n",
    "cmap = pl.cm.jet\n",
    "# extract all colors from the .jet map\n",
    "cmaplist = [cmap(i) for i in range(cmap.N)]\n",
    "# force the first color entry to be grey\n",
    "cmaplist[0] = (.5,.5,.5,1.0)\n",
    "# create the new map\n",
    "cmap = cmap.from_list('Custom cmap', cmaplist, cmap.N)\n",
    "\n",
    "# define the bins and normalize\n",
    "#bounds = np.linspace(0,5000,10)\n",
    "bounds = np.arange(0,7000,500)\n",
    "norm = mpl.colors.BoundaryNorm(bounds, cmap.N)\n",
    "\n",
    "aa=0.9\n",
    "# Elevation unterteilen mit Cumulativer Summe\n",
    "e_pos = gr_nr_elev.cumsum().astype(int)\n",
    "\n",
    "pl.figure(figsize=(26,8))\n",
    "\n",
    "for eee in range(len(gr_nr_elev)):\n",
    "    pl.subplot(2,5,eee+1)\n",
    "    \n",
    "    if eee==0: \n",
    "        sc = pl.scatter(SRR[0:e_pos[0]],\n",
    "                   GRR[0:e_pos[0]],\n",
    "                   c=ZZ[0:e_pos[0]],\n",
    "                   alpha=aa,label=r'$\\Phi$: '+str(gr_ele[0])[0:4],\n",
    "                   color='black',\n",
    "                   s=40,#GN[0:e_pos[0]]/10,\n",
    "                   vmin=0.,\n",
    "                   vmax=7000,\n",
    "                   marker='o',\n",
    "                   cmap=cmap, norm=norm)\n",
    "        sc.set_facecolor(\"none\")\n",
    "        pl.title('r: ' + corcor(SRR[0:e_pos[0]],GRR[0:e_pos[0]])[0:4]+\n",
    "                '\\n bias: ' +  str(stats(SRR[0:e_pos[0]],GRR[0:e_pos[0]])[0])[0:4]#+\n",
    "                 #'\\n rmse: ' + str(stats(SRR[0:e_pos[0]],GRR[0:e_pos[0]])[1])[0:4]\n",
    "                ,fontsize=ff)\n",
    "        pl.legend(loc='upper left',markerscale=0, fontsize=ff, scatterpoints=1)\n",
    "        cbar = pl.colorbar(shrink=0.7, cmap=cmap, norm=norm, ticks=bounds, boundaries=bounds)\n",
    "        cbar.ax.set_ylabel('Hight in m', fontsize=ff)\n",
    "        cbar.ax.tick_params(labelsize=ff)\n",
    "        \n",
    "    else:\n",
    "        sc = pl.scatter(SRR[e_pos[eee-1]:e_pos[eee]],\n",
    "                   GRR[e_pos[eee-1]:e_pos[eee]],\n",
    "                   c=ZZ[e_pos[eee-1]:e_pos[eee]],\n",
    "                   vmin = 0.,\n",
    "                   vmax = 7000.,\n",
    "                   alpha=aa,\n",
    "                   label=r'$\\Phi$: ' + str(gr_ele[eee])[0:4],\n",
    "                   color='black',\n",
    "                   s=40,#GN[e_pos[eee-1]:e_pos[eee]]/10,\n",
    "                   marker='o', \n",
    "                   cmap=cmap, norm=norm)\n",
    "        sc.set_facecolor(\"none\")\n",
    "        pl.title('r: ' + corcor(SRR[e_pos[eee-1]:e_pos[eee]],GRR[e_pos[eee-1]:e_pos[eee]])[0:4]+\n",
    "                '\\n bias: ' +  str(stats(SRR[e_pos[eee-1]:e_pos[eee]],GRR[e_pos[eee-1]:e_pos[eee]])[0])[0:4] #+\n",
    "                 #'\\n rmse: ' + str(stats(SRR[e_pos[eee-1]:e_pos[eee]],GRR[e_pos[eee-1]:e_pos[eee]])[1])[0:4]\n",
    "                ,fontsize=ff)\n",
    "        pl.legend(loc='upper left',markerscale=0, fontsize=ff, scatterpoints=1)\n",
    "        cbar = pl.colorbar(shrink=0.7, cmap=cmap, norm=norm, ticks=bounds, boundaries=bounds)\n",
    "        cbar.ax.set_ylabel('Hight in m', fontsize=ff)\n",
    "        cbar.ax.tick_params(labelsize=ff)\n",
    "        \n",
    "    pl.xlabel('SR in dbz',fontsize=ff)\n",
    "    pl.ylabel('GR in dbz',fontsize=ff)\n",
    "    pl.grid()\n",
    "    pl.xlim(15,45)\n",
    "    pl.ylim(15,45)\n",
    "    pl.plot([0,45],[0,45], color='black', ls='--')\n",
    "    pl.xticks(fontsize=ff), pl.yticks(fontsize=ff)\n",
    "\n",
    "pl.tight_layout()\n",
    "pl.savefig('/automount/ftp/velibor/validation/new_corr_elev_hight_num_BB_stats2.pdf')\n",
    "pl.savefig('/automount/ftp/velibor/validation/new_corr_elev_hight_num_BB_stats2.png')\n",
    "pl.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def stats(a, b):\n",
    "    \"\"\"beschriften und Dict!\"\"\"\n",
    "    mask = ~np.isnan(a) & ~np.isnan(b)\n",
    "    H = len(a[mask])\n",
    "    bias = np.nansum(a - b)/H\n",
    "    rmse = np.sqrt(np.nansum(((a -b )**2.0)/H))\n",
    "    \n",
    "    return bias, rmse\n",
    "    \n",
    "    \n",
    "    \n",
    "\n",
    "ZZ = Z.copy()\n",
    "SRR = SR.copy()\n",
    "GRR = GR.copy()\n",
    "GN = GRnum.copy()\n",
    "\n",
    "_num = 150\n",
    "_h1, _h2 = 2000, 3000\n",
    "\n",
    "\n",
    "# BB entfernen\n",
    "ZZ[(Z>=_h1) & (Z<=_h2)] = np.nan\n",
    "SRR[(Z>=_h1) & (Z<=_h2)] = np.nan\n",
    "GRR[(Z>=_h1) & (Z<=_h2)] = np.nan\n",
    "GN[(Z>=_h1) & (Z<=_h2)] = np.nan\n",
    "\n",
    "# NUM entfernen\n",
    "#ZZ[(GRnum < _num) ] = np.nan\n",
    "#SRR[(GRnum < _num) ] = np.nan\n",
    "#GRR[(GRnum < _num)] = np.nan\n",
    "#GN[(GRnum < _num) ] = np.nan\n",
    "\n",
    "cmap = pl.cm.jet\n",
    "# extract all colors from the .jet map\n",
    "cmaplist = [cmap(i) for i in range(cmap.N)]\n",
    "# force the first color entry to be grey\n",
    "cmaplist[0] = (.5,.5,.5,1.0)\n",
    "# create the new map\n",
    "cmap = cmap.from_list('Custom cmap', cmaplist, cmap.N)\n",
    "\n",
    "# define the bins and normalize\n",
    "#bounds = np.linspace(0,5000,10)\n",
    "bounds = np.arange(0,7000,500)\n",
    "norm = mpl.colors.BoundaryNorm(bounds, cmap.N)\n",
    "\n",
    "aa=0.9\n",
    "# Elevation unterteilen mit Cumulativer Summe\n",
    "e_pos = gr_nr_elev.cumsum().astype(int)\n",
    "\n",
    "pl.figure(figsize=(26,8))\n",
    "\n",
    "for eee in range(len(gr_nr_elev)):\n",
    "    pl.subplot(2,5,eee+1)\n",
    "    \n",
    "    if eee==0: \n",
    "        sc = pl.scatter(SRR[0:e_pos[0]],\n",
    "                   GRR[0:e_pos[0]],\n",
    "                   c=ZZ[0:e_pos[0]],\n",
    "                   alpha=aa,label='elevation: '+str(gr_ele[0])[0:4],\n",
    "                   color='black',\n",
    "                   s=40,#GN[0:e_pos[0]]/10,\n",
    "                   vmin=0.,\n",
    "                   vmax=7000,\n",
    "                   marker='o',\n",
    "                   cmap=cmap, norm=norm)\n",
    "        sc.set_facecolor(\"none\")\n",
    "        pl.title('r: ' + corcor(SRR[0:e_pos[0]],GRR[0:e_pos[0]])[0:4]+\n",
    "                '\\n bias: ' +  str(stats(SRR[0:e_pos[0]],GRR[0:e_pos[0]])[0])[0:4]#+\n",
    "                 #'\\n rmse: ' + str(stats(SRR[0:e_pos[0]],GRR[0:e_pos[0]])[1])[0:4]\n",
    "                ,fontsize=20)\n",
    "        pl.legend(loc='upper left',markerscale=0)\n",
    "        cbar = pl.colorbar(shrink=0.7, cmap=cmap, norm=norm, ticks=bounds, boundaries=bounds)\n",
    "        cbar.ax.set_ylabel('Hight in m', fontsize=10)\n",
    "        \n",
    "    else:\n",
    "        sc = pl.scatter(SRR[e_pos[eee-1]:e_pos[eee]],\n",
    "                   GRR[e_pos[eee-1]:e_pos[eee]],\n",
    "                   c=ZZ[e_pos[eee-1]:e_pos[eee]],\n",
    "                   vmin = 0.,\n",
    "                   vmax = 7000.,\n",
    "                   alpha=aa,\n",
    "                   label='elevation: ' + str(gr_ele[eee])[0:4],\n",
    "                   color='black',\n",
    "                   s=40,#GN[e_pos[eee-1]:e_pos[eee]]/10,\n",
    "                   marker='o', \n",
    "                   cmap=cmap, norm=norm)\n",
    "        sc.set_facecolor(\"none\")\n",
    "        pl.title('r: ' + corcor(SRR[e_pos[eee-1]:e_pos[eee]],GRR[e_pos[eee-1]:e_pos[eee]])[0:4]+\n",
    "                '\\n bias: ' +  str(stats(SRR[e_pos[eee-1]:e_pos[eee]],GRR[e_pos[eee-1]:e_pos[eee]])[0])[0:4] #+\n",
    "                 #'\\n rmse: ' + str(stats(SRR[e_pos[eee-1]:e_pos[eee]],GRR[e_pos[eee-1]:e_pos[eee]])[1])[0:4]\n",
    "                ,fontsize=20)\n",
    "        pl.legend(loc='upper left',markerscale=0)\n",
    "        cbar = pl.colorbar(shrink=0.7, cmap=cmap, norm=norm, ticks=bounds, boundaries=bounds)\n",
    "        cbar.ax.set_ylabel('Hight in m', fontsize=10)\n",
    "        \n",
    "    pl.xlabel('SR in dbz')\n",
    "    pl.ylabel('GR in dbz')\n",
    "    pl.grid()\n",
    "    pl.xlim(0,45)\n",
    "    pl.ylim(0,45)\n",
    "    pl.plot([0,45],[0,45], color='black', ls='--')\n",
    "\n",
    "pl.tight_layout()\n",
    "pl.savefig('/automount/ftp/velibor/validation/new_corr_elev_hight_num_BB_stats2_withoutBB.pdf')\n",
    "pl.savefig('/automount/ftp/velibor/validation/new_corr_elev_hight_num_BB_stats2_withoutBB.png')\n",
    "pl.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def stats(a, b):\n",
    "    \"\"\"beschriften und Dict!\"\"\"\n",
    "    mask = ~np.isnan(a) & ~np.isnan(b)\n",
    "    H = len(a[mask])\n",
    "    bias = np.nansum(a - b)/H\n",
    "    rmse = np.sqrt(np.nansum(((a -b )**2.0)/H))\n",
    "    \n",
    "    return bias, rmse\n",
    "    \n",
    "    \n",
    "    \n",
    "\n",
    "ZZ = Z.copy()\n",
    "SRR = SR.copy()\n",
    "GRR = GR.copy()\n",
    "GN = GRnum.copy()\n",
    "\n",
    "_num = 150\n",
    "_h1, _h2 = 2200, 2800\n",
    "\n",
    "\n",
    "# BB entfernen\n",
    "\"\"\"ZZ[(Z>=_h1) & (Z<=_h2)] = np.nan\n",
    "SRR[(Z>=_h1) & (Z<=_h2)] = np.nan\n",
    "GRR[(Z>=_h1) & (Z<=_h2)] = np.nan\n",
    "GN[(Z>=_h1) & (Z<=_h2)] = np.nan\"\"\"\n",
    "\n",
    "# NUM entfernen\n",
    "#ZZ[(GRnum < _num) ] = np.nan\n",
    "#SRR[(GRnum < _num) ] = np.nan\n",
    "#GRR[(GRnum < _num)] = np.nan\n",
    "#GN[(GRnum < _num) ] = np.nan\n",
    "\n",
    "cmap = pl.cm.jet\n",
    "# extract all colors from the .jet map\n",
    "cmaplist = [cmap(i) for i in range(cmap.N)]\n",
    "# force the first color entry to be grey\n",
    "cmaplist[0] = (.5,.5,.5,1.0)\n",
    "# create the new map\n",
    "cmap = cmap.from_list('Custom cmap', cmaplist, cmap.N)\n",
    "\n",
    "# define the bins and normalize\n",
    "#bounds = np.linspace(0,5000,10)\n",
    "bounds = np.arange(0,7000,500)\n",
    "norm = mpl.colors.BoundaryNorm(bounds, cmap.N)\n",
    "\n",
    "aa=0.9\n",
    "# Elevation unterteilen mit Cumulativer Summe\n",
    "e_pos = gr_nr_elev.cumsum().astype(int)\n",
    "\n",
    "pl.figure(figsize=(26,8))\n",
    "\n",
    "for eee in range(len(gr_nr_elev)):\n",
    "    pl.subplot(2,5,eee+1)\n",
    "    \n",
    "    if eee==0: \n",
    "        mask1 = ~np.isnan(SRR[0:e_pos[0]])\n",
    "        mask2 = ~np.isnan(GRR[0:e_pos[0]])\n",
    "                             \n",
    "        pl.hist(SRR[0:e_pos[0]][mask1],bins=50,\n",
    "                   label='elevation: '+str(gr_ele[0])[0:3],\n",
    "                   color='black', normed=1, facecolor=\"None\", alpha=1, edgecolor='black')\n",
    "        pl.hist(GRR[0:e_pos[0]][mask2],bins=50,\n",
    "                   color='blue', normed=1, facecolor=\"None\", alpha=1, edgecolor='blue')\n",
    "\n",
    "        pl.title('r: ' + corcor(SRR[0:e_pos[0]],GRR[0:e_pos[0]])[0:4]+\n",
    "                '\\n bias: ' +  str(stats(SRR[0:e_pos[0]],GRR[0:e_pos[0]])[0])[0:4]+\n",
    "                 '\\n rmse: ' + str(stats(SRR[0:e_pos[0]],GRR[0:e_pos[0]])[1])[0:4]\n",
    "                )\n",
    "        pl.legend(loc='upper left')\n",
    "\n",
    "        \n",
    "    else:\n",
    "        mask3 = ~np.isnan(SRR[e_pos[eee-1]:e_pos[eee]])\n",
    "        mask4 = ~np.isnan(GRR[e_pos[eee-1]:e_pos[eee]])\n",
    "        \n",
    "        pl.hist(SRR[e_pos[eee-1]:e_pos[eee]][mask3],bins=50,\n",
    "                   label='elevation: '+str(gr_ele[eee])[0:3],\n",
    "                   color='black', normed=1, facecolor=\"None\", alpha=1, edgecolor='black')\n",
    "        pl.hist(GRR[e_pos[eee-1]:e_pos[eee]][mask4],bins=50,\n",
    "                   color='blue', normed=1, facecolor=\"None\", alpha=1, edgecolor='blue')\n",
    "    \n",
    "        pl.title('r: ' + corcor(SRR[e_pos[eee-1]:e_pos[eee]],GRR[e_pos[eee-1]:e_pos[eee]])[0:4]+\n",
    "                '\\n bias: ' +  str(stats(SRR[e_pos[eee-1]:e_pos[eee]],GRR[e_pos[eee-1]:e_pos[eee]])[0])[0:4] +\n",
    "                 '\\n rmse: ' + str(stats(SRR[e_pos[eee-1]:e_pos[eee]],GRR[e_pos[eee-1]:e_pos[eee]])[1])[0:4]\n",
    "                )\n",
    "        pl.legend(loc='upper left',markerscale=0)\n",
    "\n",
    "        \n",
    "    #pl.xlabel('SR in dbz')\n",
    "    #pl.ylabel('GR in dbz')\n",
    "    pl.grid()\n",
    "    pl.xlim(10,45)\n",
    "    pl.ylim(0,0.4)\n",
    "    #pl.plot([0,45],[0,45], color='black', ls='--')\n",
    "\n",
    "pl.tight_layout()\n",
    "pl.savefig('/automount/ftp/velibor/validation/new_corr_elev_hight_num_BB_stats2_hist.pdf')\n",
    "pl.savefig('/automount/ftp/velibor/validation/new_corr_elev_hight_num_BB_stats2_hist.png')\n",
    "pl.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def show_swath():\n",
    "    \"\"\"\n",
    "    Simple Darstellung des DPR Swaths\n",
    "    \n",
    "    w.i.p MS fehlt\n",
    "    \"\"\"\n",
    "    import matplotlib.pyplot as pl\n",
    "    import numpy as np\n",
    "    swath_ku = np.arange(1,50,1) \n",
    "    pos_ku =  np.ones(swath_ku.shape)\n",
    "\n",
    "    swath_ka = np.arange(1+12,26+12,1) \n",
    "    pos_ka =  np.ones(swath_ka.shape)+0.1\n",
    "\n",
    "    print('Ka: ',len(swath_ka))\n",
    "    print('Ku: ',len(swath_ku))\n",
    "    print('Ku - Ka (DFR): ',len(swath_ku)-len(swath_ka))\n",
    "\n",
    "    pl.figure(figsize=(20,5))\n",
    "    pl.scatter(swath_ku,pos_ku, c='blue', s=200, label='Ku')\n",
    "    pl.scatter(swath_ka,pos_ka, c='red', s=200, label='Ka')\n",
    "\n",
    "    pl.scatter(swath_ku[12:26+11],pos_ku[12:26+11], c='black', s=200, label='DFR')\n",
    "    pl.scatter(swath_ku[12:26+11],pos_ku[12:26+11], c='red', s=20, marker='+')\n",
    "\n",
    "    pl.grid()\n",
    "    pl.legend()\n",
    "    pl.ylim(0.5,1.5)\n",
    "    pl.show()\n",
    "\n",
    "show_swath()"
   ]
  },
  {
   "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": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def stats(a, b):\n",
    "    \"\"\"beschriften und Dict!\"\"\"\n",
    "    mask = ~np.isnan(a) & ~np.isnan(b)\n",
    "    H = len(a[mask])\n",
    "    bias = np.nansum(a - b)/H\n",
    "    rmse = np.sqrt(np.nansum(((a -b )**2.0)/H))\n",
    "    \n",
    "    return bias, rmse\n",
    "    \n",
    "    \n",
    "    \n",
    "\n",
    "ZZ = Z.copy()\n",
    "SRR = SR.copy()\n",
    "GRR = GR.copy()\n",
    "GN = GRnum.copy()\n",
    "\n",
    "_num = 150\n",
    "_h1, _h2 = 2000, 3000\n",
    "\n",
    "\n",
    "# BB entfernen\n",
    "#ZZ[(Z>=_h1) & (Z<=_h2)] = np.nan\n",
    "#SRR[(Z>=_h1) & (Z<=_h2)] = np.nan\n",
    "#GRR[(Z>=_h1) & (Z<=_h2)] = np.nan\n",
    "#GN[(Z>=_h1) & (Z<=_h2)] = np.nan\n",
    "\n",
    "# NUM entfernen\n",
    "#ZZ[(GRnum < _num) ] = np.nan\n",
    "#SRR[(GRnum < _num) ] = np.nan\n",
    "#GRR[(GRnum < _num)] = np.nan\n",
    "#GN[(GRnum < _num) ] = np.nan\n",
    "\n",
    "cmap = pl.cm.Set1\n",
    "# extract all colors from the .jet map\n",
    "cmaplist = [cmap(i) for i in range(cmap.N)]\n",
    "# force the first color entry to be grey\n",
    "cmaplist[0] = (.5,.5,.5,1.0)\n",
    "# create the new map\n",
    "cmap = cmap.from_list('Custom cmap', cmaplist, cmap.N)\n",
    "\n",
    "# define the bins and normalize\n",
    "#bounds = np.linspace(0,5000,10)\n",
    "bounds = np.arange(0,7000,1000)\n",
    "norm = mpl.colors.BoundaryNorm(bounds, cmap.N)\n",
    "\n",
    "aa=1\n",
    "# Elevation unterteilen mit Cumulativer Summe\n",
    "e_pos = gr_nr_elev.cumsum().astype(int)\n",
    "\n",
    "pl.figure(figsize=(26,10))\n",
    "\n",
    "for eee in range(len(gr_nr_elev)):\n",
    "    pl.subplot(2,5,eee+1)\n",
    "    \n",
    "    if eee==0: \n",
    "        sc = pl.scatter(SRR[0:e_pos[0]],\n",
    "                   GRR[0:e_pos[0]],\n",
    "                   c=ZZ[0:e_pos[0]],\n",
    "                   alpha=aa,label=r'$\\Phi$: '+str(gr_ele[0])[0:4],\n",
    "                   color='black',\n",
    "                   s=40,#GN[0:e_pos[0]]/10,\n",
    "                   vmin=0.,\n",
    "                   vmax=7000,\n",
    "                   marker='o',\n",
    "                   cmap=cmap, norm=norm)\n",
    "        sc.set_facecolor(\"none\")\n",
    "        pl.title('r: ' + corcor(SRR[0:e_pos[0]],GRR[0:e_pos[0]])[0:4]+\n",
    "                '\\n bias: ' +  str(stats(SRR[0:e_pos[0]],GRR[0:e_pos[0]])[0])[0:4]#+\n",
    "                 #'\\n rmse: ' + str(stats(SRR[0:e_pos[0]],GRR[0:e_pos[0]])[1])[0:4]\n",
    "                ,fontsize=25)\n",
    "        pl.legend(loc='lower right',markerscale=0, fontsize=20, borderpad=0.1,frameon=False)\n",
    "        #cbar = pl.colorbar(shrink=0.7, cmap=cmap, norm=norm, ticks=bounds, boundaries=bounds)\n",
    "        #cbar.ax.set_ylabel('Hight in m', fontsize=10)\n",
    "        #cbar.ax.tick_params(labelsize=20)\n",
    "        \n",
    "    else:\n",
    "        sc = pl.scatter(SRR[e_pos[eee-1]:e_pos[eee]],\n",
    "                   GRR[e_pos[eee-1]:e_pos[eee]],\n",
    "                   c=ZZ[e_pos[eee-1]:e_pos[eee]],\n",
    "                   vmin = 0.,\n",
    "                   vmax = 7000.,\n",
    "                   alpha=aa,\n",
    "                   label=r'$\\Phi$: ' + str(gr_ele[eee])[0:4],\n",
    "                   color='black',\n",
    "                   s=40,#GN[e_pos[eee-1]:e_pos[eee]]/10,\n",
    "                   marker='o', \n",
    "                   cmap=cmap, norm=norm)\n",
    "        sc.set_facecolor(\"none\")\n",
    "        pl.title('r: ' + corcor(SRR[e_pos[eee-1]:e_pos[eee]],GRR[e_pos[eee-1]:e_pos[eee]])[0:4]+\n",
    "                '\\n bias: ' +  str(stats(SRR[e_pos[eee-1]:e_pos[eee]],GRR[e_pos[eee-1]:e_pos[eee]])[0])[0:4] #+\n",
    "                 #'\\n rmse: ' + str(stats(SRR[e_pos[eee-1]:e_pos[eee]],GRR[e_pos[eee-1]:e_pos[eee]])[1])[0:4]\n",
    "                ,fontsize=25)\n",
    "        pl.legend(loc='lower right',markerscale=0, fontsize=20, borderpad=0.1,frameon=False)\n",
    "        #cbar = pl.colorbar(shrink=0.7, cmap=cmap, norm=norm, ticks=bounds, boundaries=bounds)\n",
    "        #cbar.ax.set_ylabel('Hight in m', fontsize=20)\n",
    "        #cbar.ax.tick_params(labelsize=20)\n",
    "        \n",
    "    cbar = pl.colorbar(shrink=0.7, cmap=cmap, norm=norm, ticks=bounds, boundaries=bounds)\n",
    "    cbar.ax.set_ylabel('Hight in m', fontsize=20)\n",
    "    cbar.ax.tick_params(labelsize=20)    \n",
    "    pl.xlabel('SR in dbz', fontsize=20)\n",
    "    pl.ylabel('GR in dbz', fontsize=20)\n",
    "    pl.grid()\n",
    "    pl.xlim(15,45)\n",
    "    pl.ylim(15,45)\n",
    "    pl.plot([0,45],[0,45], color='black', ls='--')\n",
    "    pl.xticks(fontsize=20); pl.yticks(fontsize=20)\n",
    "\n",
    "pl.tight_layout()\n",
    "pl.savefig('/automount/ftp/velibor/validation/irs.pdf')\n",
    "pl.savefig('/automount/ftp/velibor/validation/irs.png')\n",
    "pl.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def stats(a, b):\n",
    "    \"\"\"beschriften und Dict!\"\"\"\n",
    "    mask = ~np.isnan(a) & ~np.isnan(b)\n",
    "    H = len(a[mask])\n",
    "    bias = np.nansum(a - b)/H\n",
    "    rmse = np.sqrt(np.nansum(((a -b )**2.0)/H))\n",
    "    \n",
    "    return bias, rmse\n",
    "    \n",
    "    \n",
    "    \n",
    "\n",
    "ZZ = Z.copy()/1000.\n",
    "SRR = SR.copy()\n",
    "GRR = GR.copy()\n",
    "GN = GRnum.copy()\n",
    "\n",
    "_num = 150\n",
    "_h1, _h2 = 2000, 3000\n",
    "\n",
    "\n",
    "# BB entfernen\n",
    "#ZZ[(Z>=_h1) & (Z<=_h2)] = np.nan\n",
    "#SRR[(Z>=_h1) & (Z<=_h2)] = np.nan\n",
    "#GRR[(Z>=_h1) & (Z<=_h2)] = np.nan\n",
    "#GN[(Z>=_h1) & (Z<=_h2)] = np.nan\n",
    "\n",
    "# NUM entfernen\n",
    "#ZZ[(GRnum < _num) ] = np.nan\n",
    "#SRR[(GRnum < _num) ] = np.nan\n",
    "#GRR[(GRnum < _num)] = np.nan\n",
    "#GN[(GRnum < _num) ] = np.nan\n",
    "\n",
    "cmap = pl.cm.jet\n",
    "# extract all colors from the .jet map\n",
    "cmaplist = [cmap(i) for i in range(cmap.N)]\n",
    "# force the first color entry to be grey\n",
    "cmaplist[0] = (.5,.5,.5,1.0)\n",
    "# create the new map\n",
    "cmap = cmap.from_list('Custom cmap', cmaplist, cmap.N)\n",
    "\n",
    "# define the bins and normalize\n",
    "#bounds = np.linspace(0,5000,10)\n",
    "bounds = np.arange(0,7000,500)/1000.\n",
    "norm = mpl.colors.BoundaryNorm(bounds, cmap.N)\n",
    "\n",
    "aa=1\n",
    "# Elevation unterteilen mit Cumulativer Summe\n",
    "e_pos = gr_nr_elev.cumsum().astype(int)\n",
    "\n",
    "pl.figure(figsize=(26,10))\n",
    "\n",
    "for eee in range(len(gr_nr_elev)):\n",
    "    pl.subplot(2,5,eee+1)\n",
    "    \n",
    "    if eee==0: \n",
    "        sc = pl.scatter(SRR,\n",
    "                   GRR,\n",
    "                   c=ZZ,\n",
    "                   alpha=aa,label=r'$\\Phi$: all',\n",
    "                   color='black',\n",
    "                   s=40,#GN[0:e_pos[0]]/10,\n",
    "                   vmin=0.,\n",
    "                   vmax=7000,\n",
    "                   marker='o',\n",
    "                   cmap=cmap, norm=norm)\n",
    "        sc.set_facecolor(\"none\")\n",
    "        pl.title('r: ' + corcor(SRR,GRR)[0:4]+\n",
    "                '\\n bias: ' +  str(stats(SRR,GRR)[0])[0:4]#+\n",
    "                 #'\\n rmse: ' + str(stats(SRR[0:e_pos[0]],GRR[0:e_pos[0]])[1])[0:4]\n",
    "                ,fontsize=25)\n",
    "        pl.legend(loc='lower right',markerscale=0, fontsize=20, borderpad=0.1,frameon=False)\n",
    "        #cbar = pl.colorbar(shrink=0.7, cmap=cmap, norm=norm, ticks=bounds, boundaries=bounds)\n",
    "        #cbar.ax.set_ylabel('Hight in m', fontsize=10)\n",
    "        #cbar.ax.tick_params(labelsize=20)\n",
    "        \n",
    "    else:\n",
    "        sc = pl.scatter(SRR[e_pos[eee-1]:e_pos[eee]],\n",
    "                   GRR[e_pos[eee-1]:e_pos[eee]],\n",
    "                   c=ZZ[e_pos[eee-1]:e_pos[eee]],\n",
    "                   vmin = 0.,\n",
    "                   vmax = 7000.,\n",
    "                   alpha=aa,\n",
    "                   label=r'$\\Phi$: ' + str(gr_ele[eee])[0:4],\n",
    "                   color='black',\n",
    "                   s=40,#GN[e_pos[eee-1]:e_pos[eee]]/10,\n",
    "                   marker='o', \n",
    "                   cmap=cmap, norm=norm)\n",
    "        sc.set_facecolor(\"none\")\n",
    "        pl.title('r: ' + corcor(SRR[e_pos[eee-1]:e_pos[eee]],GRR[e_pos[eee-1]:e_pos[eee]])[0:4]+\n",
    "                '\\n bias: ' +  str(stats(SRR[e_pos[eee-1]:e_pos[eee]],GRR[e_pos[eee-1]:e_pos[eee]])[0])[0:4] #+\n",
    "                 #'\\n rmse: ' + str(stats(SRR[e_pos[eee-1]:e_pos[eee]],GRR[e_pos[eee-1]:e_pos[eee]])[1])[0:4]\n",
    "                ,fontsize=25)\n",
    "        pl.legend(loc='lower right',markerscale=0, fontsize=20, borderpad=0.1,frameon=False)\n",
    "        #cbar = pl.colorbar(shrink=0.7, cmap=cmap, norm=norm, ticks=bounds, boundaries=bounds)\n",
    "        #cbar.ax.set_ylabel('Hight in m', fontsize=20)\n",
    "        #cbar.ax.tick_params(labelsize=20)\n",
    "        \n",
    "    cbar = pl.colorbar(shrink=0.7, cmap=cmap, norm=norm, ticks=bounds, boundaries=bounds, orientation='vertical')\n",
    "    cbar.ax.set_ylabel('Hight in km', fontsize=20)\n",
    "    cbar.ax.tick_params(labelsize=20)    \n",
    "    pl.xlabel('SR in dbz', fontsize=20)\n",
    "    pl.ylabel('GR in dbz', fontsize=20)\n",
    "    pl.grid()\n",
    "    pl.xlim(15,45)\n",
    "    pl.ylim(15,45)\n",
    "    pl.plot([0,45],[0,45], color='black', ls='--')\n",
    "    pl.xticks(fontsize=20); pl.yticks(fontsize=20)\n",
    "\n",
    "pl.tight_layout()\n",
    "pl.savefig('/automount/ftp/velibor/validation/irs22.pdf')\n",
    "pl.savefig('/automount/ftp/velibor/validation/irs22.png')\n",
    "pl.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "np.exp(-1000000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "pl.figure(figsize=(8,8))\n",
    "\n",
    "c = pl.scatter(SRR,\n",
    "                   GRR,\n",
    "                   c=ZZ,\n",
    "                   alpha=aa,label=r'$\\Phi$: all',\n",
    "                   color='black',\n",
    "                   s=40,#GN[0:e_pos[0]]/10,\n",
    "                   vmin=0.,\n",
    "                   vmax=7000,\n",
    "                   marker='o',\n",
    "                   cmap=cmap, norm=norm)\n",
    "sc.set_facecolor(\"none\")\n",
    "pl.title('r: ' + corcor(SRR,GRR)[0:4]+\n",
    "        '\\n bias: ' +  str(stats(SRR,GRR)[0])[0:4]#+\n",
    "         #'\\n rmse: ' + str(stats(SRR[0:e_pos[0]],GRR[0:e_pos[0]])[1])[0:4]\n",
    "        ,fontsize=25)\n",
    "pl.legend(loc='lower right',markerscale=0, fontsize=20, borderpad=0.1,frameon=False)\n",
    "cbar = pl.colorbar(shrink=0.7, cmap=cmap, norm=norm, ticks=bounds, boundaries=bounds, orientation='vertical')\n",
    "cbar.ax.set_ylabel('Hight in km', fontsize=20)\n",
    "cbar.ax.tick_params(labelsize=20)    \n",
    "pl.xlabel('SR in dbz', fontsize=20)\n",
    "pl.ylabel('GR in dbz', fontsize=20)\n",
    "pl.grid()\n",
    "pl.xlim(15,45)\n",
    "pl.ylim(15,45)\n",
    "pl.plot([0,45],[0,45], color='black', ls='--')\n",
    "pl.xticks(fontsize=20); pl.yticks(fontsize=20)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "Z"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "pl.figure(figsize=(8,8))\n",
    "srr, grr = SRR.copy(), GRR.copy()\n",
    "srr[Z>1000]= np.nan\n",
    "grr[Z>1000]= np.nan\n",
    "\n",
    "cmap = pl.cm.Set1\n",
    "# extract all colors from the .jet map\n",
    "cmaplist = [cmap(i) for i in range(cmap.N)]\n",
    "# force the first color entry to be grey\n",
    "cmaplist[0] = (.5,.5,.5,1.0)\n",
    "# create the new map\n",
    "cmap = cmap.from_list('Custom cmap', cmaplist, cmap.N)\n",
    "\n",
    "# define the bins and normalize\n",
    "#bounds = np.linspace(0,5000,10)\n",
    "bounds = np.arange(0,7000,1000)/1000.\n",
    "norm = mpl.colors.BoundaryNorm(bounds, cmap.N)\n",
    "\n",
    "aa=1\n",
    "# Elevation unterteilen mit Cumulativer Summe\n",
    "e_pos = gr_nr_elev.cumsum().astype(int)\n",
    "\n",
    "\n",
    "c = pl.scatter(srr,\n",
    "                   grr,\n",
    "                   c=ZZ,\n",
    "                   alpha=aa,label=r'$\\Phi$: all',\n",
    "                   color='black',\n",
    "                   s=40,#GN[0:e_pos[0]]/10,\n",
    "                   vmin=0.,\n",
    "                   vmax=7000,\n",
    "                   marker='o',\n",
    "                   cmap=cmap, norm=norm)\n",
    "sc.set_facecolor(\"none\")\n",
    "pl.title('r: ' + corcor(SRR,GRR)[0:4]+\n",
    "        '\\n bias: ' +  str(stats(SRR,GRR)[0])[0:4]#+\n",
    "         #'\\n rmse: ' + str(stats(SRR[0:e_pos[0]],GRR[0:e_pos[0]])[1])[0:4]\n",
    "        ,fontsize=25)\n",
    "pl.legend(loc='lower right',markerscale=0, fontsize=20, borderpad=0.1,frameon=False)\n",
    "cbar = pl.colorbar(shrink=0.7, cmap=cmap, norm=norm, ticks=bounds, boundaries=bounds, orientation='vertical')\n",
    "cbar.ax.set_ylabel('Hight in km', fontsize=20)\n",
    "cbar.ax.tick_params(labelsize=20)    \n",
    "pl.xlabel('SR in dbz', fontsize=20)\n",
    "pl.ylabel('GR in dbz', fontsize=20)\n",
    "pl.grid()\n",
    "pl.xlim(15,45)\n",
    "pl.ylim(15,45)\n",
    "pl.plot([0,45],[0,45], color='black', ls='--')\n",
    "pl.xticks(fontsize=20); pl.yticks(fontsize=20)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "pl.figure(figsize=(8,8))\n",
    "srr, grr = SRR.copy(), GRR.copy()\n",
    "srr[(Z<1000)]= np.nan\n",
    "grr[(Z<1000)]= np.nan\n",
    "\n",
    "srr[(Z>2000)]= np.nan\n",
    "grr[(Z>2000)]= np.nan\n",
    "\n",
    "cmap = pl.cm.Set1\n",
    "# extract all colors from the .jet map\n",
    "cmaplist = [cmap(i) for i in range(cmap.N)]\n",
    "# force the first color entry to be grey\n",
    "cmaplist[0] = (.5,.5,.5,1.0)\n",
    "# create the new map\n",
    "cmap = cmap.from_list('Custom cmap', cmaplist, cmap.N)\n",
    "\n",
    "# define the bins and normalize\n",
    "#bounds = np.linspace(0,5000,10)\n",
    "bounds = np.arange(0,7000,1000)/1000.\n",
    "norm = mpl.colors.BoundaryNorm(bounds, cmap.N)\n",
    "\n",
    "aa=1\n",
    "# Elevation unterteilen mit Cumulativer Summe\n",
    "e_pos = gr_nr_elev.cumsum().astype(int)\n",
    "\n",
    "\n",
    "c = pl.scatter(srr,\n",
    "                   grr,\n",
    "                   c=ZZ,\n",
    "                   alpha=aa,label=r'$\\Phi$: all',\n",
    "                   color='black',\n",
    "                   s=40,#GN[0:e_pos[0]]/10,\n",
    "                   vmin=0.,\n",
    "                   vmax=7000,\n",
    "                   marker='o',\n",
    "                   cmap=cmap, norm=norm)\n",
    "sc.set_facecolor(\"none\")\n",
    "pl.title('r: ' + corcor(SRR,GRR)[0:4]+\n",
    "        '\\n bias: ' +  str(stats(SRR,GRR)[0])[0:4]#+\n",
    "         #'\\n rmse: ' + str(stats(SRR[0:e_pos[0]],GRR[0:e_pos[0]])[1])[0:4]\n",
    "        ,fontsize=25)\n",
    "pl.legend(loc='lower right',markerscale=0, fontsize=20, borderpad=0.1,frameon=False)\n",
    "cbar = pl.colorbar(shrink=0.7, cmap=cmap, norm=norm, ticks=bounds, boundaries=bounds, orientation='vertical')\n",
    "cbar.ax.set_ylabel('Hight in km', fontsize=20)\n",
    "cbar.ax.tick_params(labelsize=20)    \n",
    "pl.xlabel('SR in dbz', fontsize=20)\n",
    "pl.ylabel('GR in dbz', fontsize=20)\n",
    "pl.grid()\n",
    "pl.xlim(15,45)\n",
    "pl.ylim(15,45)\n",
    "pl.plot([0,45],[0,45], color='black', ls='--')\n",
    "pl.xticks(fontsize=20); pl.yticks(fontsize=20)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "pl.figure(figsize=(8,8))\n",
    "srr, grr = SRR.copy(), GRR.copy()\n",
    "srr[(Z<2000)]= np.nan\n",
    "grr[(Z<2000)]= np.nan\n",
    "\n",
    "srr[(Z>3000)]= np.nan\n",
    "grr[(Z>3000)]= np.nan\n",
    "\n",
    "cmap = pl.cm.Set1\n",
    "# extract all colors from the .jet map\n",
    "cmaplist = [cmap(i) for i in range(cmap.N)]\n",
    "# force the first color entry to be grey\n",
    "cmaplist[0] = (.5,.5,.5,1.0)\n",
    "# create the new map\n",
    "cmap = cmap.from_list('Custom cmap', cmaplist, cmap.N)\n",
    "\n",
    "# define the bins and normalize\n",
    "#bounds = np.linspace(0,5000,10)\n",
    "bounds = np.arange(0,7000,1000)/1000.\n",
    "norm = mpl.colors.BoundaryNorm(bounds, cmap.N)\n",
    "\n",
    "aa=1\n",
    "# Elevation unterteilen mit Cumulativer Summe\n",
    "e_pos = gr_nr_elev.cumsum().astype(int)\n",
    "\n",
    "\n",
    "c = pl.scatter(srr,\n",
    "                   grr,\n",
    "                   c=ZZ,\n",
    "                   alpha=aa,label=r'$\\Phi$: all',\n",
    "                   color='black',\n",
    "                   s=40,#GN[0:e_pos[0]]/10,\n",
    "                   vmin=0.,\n",
    "                   vmax=7000,\n",
    "                   marker='o',\n",
    "                   cmap=cmap, norm=norm)\n",
    "sc.set_facecolor(\"none\")\n",
    "pl.title('r: ' + corcor(SRR,GRR)[0:4]+\n",
    "        '\\n bias: ' +  str(stats(SRR,GRR)[0])[0:4]#+\n",
    "         #'\\n rmse: ' + str(stats(SRR[0:e_pos[0]],GRR[0:e_pos[0]])[1])[0:4]\n",
    "        ,fontsize=25)\n",
    "pl.legend(loc='lower right',markerscale=0, fontsize=20, borderpad=0.1,frameon=False)\n",
    "cbar = pl.colorbar(shrink=0.7, cmap=cmap, norm=norm, ticks=bounds, boundaries=bounds, orientation='vertical')\n",
    "cbar.ax.set_ylabel('Hight in km', fontsize=20)\n",
    "cbar.ax.tick_params(labelsize=20)    \n",
    "pl.xlabel('SR in dbz', fontsize=20)\n",
    "pl.ylabel('GR in dbz', fontsize=20)\n",
    "pl.grid()\n",
    "pl.xlim(15,45)\n",
    "pl.ylim(15,45)\n",
    "pl.plot([0,45],[0,45], color='black', ls='--')\n",
    "pl.xticks(fontsize=20); pl.yticks(fontsize=20)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "pl.figure(figsize=(8,8))\n",
    "srr, grr = SRR.copy(), GRR.copy()\n",
    "srr[(Z<3000)]= np.nan\n",
    "grr[(Z<3000)]= np.nan\n",
    "\n",
    "\n",
    "cmap = pl.cm.Set1\n",
    "# extract all colors from the .jet map\n",
    "cmaplist = [cmap(i) for i in range(cmap.N)]\n",
    "# force the first color entry to be grey\n",
    "cmaplist[0] = (.5,.5,.5,1.0)\n",
    "# create the new map\n",
    "cmap = cmap.from_list('Custom cmap', cmaplist, cmap.N)\n",
    "\n",
    "# define the bins and normalize\n",
    "#bounds = np.linspace(0,5000,10)\n",
    "bounds = np.arange(0,7000,1000)/1000.\n",
    "norm = mpl.colors.BoundaryNorm(bounds, cmap.N)\n",
    "\n",
    "aa=1\n",
    "# Elevation unterteilen mit Cumulativer Summe\n",
    "e_pos = gr_nr_elev.cumsum().astype(int)\n",
    "\n",
    "\n",
    "c = pl.scatter(srr,\n",
    "                   grr,\n",
    "                   c=ZZ,\n",
    "                   alpha=aa,label=r'$\\Phi$: all',\n",
    "                   color='black',\n",
    "                   s=40,#GN[0:e_pos[0]]/10,\n",
    "                   vmin=0.,\n",
    "                   vmax=7000,\n",
    "                   marker='o',\n",
    "                   cmap=cmap, norm=norm)\n",
    "sc.set_facecolor(\"none\")\n",
    "pl.title('r: ' + corcor(SRR,GRR)[0:4]+\n",
    "        '\\n bias: ' +  str(stats(SRR,GRR)[0])[0:4]#+\n",
    "         #'\\n rmse: ' + str(stats(SRR[0:e_pos[0]],GRR[0:e_pos[0]])[1])[0:4]\n",
    "        ,fontsize=25)\n",
    "pl.legend(loc='lower right',markerscale=0, fontsize=20, borderpad=0.1,frameon=False)\n",
    "cbar = pl.colorbar(shrink=0.7, cmap=cmap, norm=norm, ticks=bounds, boundaries=bounds, orientation='vertical')\n",
    "cbar.ax.set_ylabel('Hight in km', fontsize=20)\n",
    "cbar.ax.tick_params(labelsize=20)    \n",
    "pl.xlabel('SR in dbz', fontsize=20)\n",
    "pl.ylabel('GR in dbz', fontsize=20)\n",
    "pl.grid()\n",
    "pl.xlim(15,45)\n",
    "pl.ylim(15,45)\n",
    "pl.plot([0,45],[0,45], color='black', ls='--')\n",
    "pl.xticks(fontsize=20); pl.yticks(fontsize=20)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "rho_gr.shape, alt.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "np.nanmean(SR-GR)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "pl.scatter(SR,GR)\n",
    "pl.xlim(0,50)\n",
    "pl.ylim(0,50)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "ref_sr.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "sr_data['refl'].shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "xyzp_sr.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "xyzp_sr[:,:,:,2].shape"
   ]
  },
  {
   "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": []
  },
  {
   "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
}