_3_PAPER_stat.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import h5py\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as pl\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "import wradlib\n",
    "from osgeo import osr\n",
    "import datetime as dt\n",
    "import glob\n",
    "from osgeo import gdal\n",
    "from scipy.spatial import cKDTree\n",
    "import wradlib as wrl\n",
    "import matplotlib as mpl\n",
    "\n",
    "#% matplotlib inline \n",
    "\n",
    "import datetime as dt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def corcorcor(A,B):\n",
    "    from scipy import stats, linspace\n",
    "\n",
    "        \n",
    "    try:\n",
    "        mask = ~np.isnan(A) & ~np.isnan(B)\n",
    "        diffi = A[mask]-B[mask]\n",
    "        bias = np.nansum(diffi)/len(diffi)\n",
    "        slope, intercept, r_value, p_value, std_err = stats.linregress(A[mask], B[mask])\n",
    "        line = slope * A +intercept\n",
    "        corr, bias = np.round(r_value,2),  np.round(bias,2)\n",
    "    except:\n",
    "        corr, bias = np.nan, np.nan\n",
    "    \n",
    "    return corr, bias"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def skill_score(estimate, reference, th):\n",
    "    \"\"\"\n",
    "    estimate\n",
    "    \n",
    "    reference\n",
    "    \n",
    "    threshold\n",
    "    \"\"\"\n",
    "    # SkillScore berechnung\n",
    "    #contigency tables Lit: Tang et al. 2015\n",
    "    from scipy import stats, linspace\n",
    "    \n",
    "    val1 = estimate.copy()\n",
    "    val2 = reference.copy()\n",
    "\n",
    "    H_pos = np.array((val1 > th) & (val2 > th))\n",
    "    M_pos = np.array((val1 <= th) & (val2 > th))\n",
    "    F_pos = np.array((val1 > th) & (val2 <= th))\n",
    "    C_pos = np.array((val1 <= th) & (val2 <= th))\n",
    "\n",
    "    H = float(len(H_pos[H_pos==True]))\n",
    "    M = float(len(M_pos[M_pos==True]))\n",
    "    F = float(len(F_pos[F_pos==True]))\n",
    "    C = float(len(C_pos[C_pos==True]))\n",
    "\n",
    "    #maske = ~np.isnan(val1) & ~np.isnan(val2)\n",
    "    #N = float(len(maske[maske==True]))\n",
    "    N = np.nansum([H, M, F, C])\n",
    "    \n",
    "    # Rates\n",
    "    try:\n",
    "        E = 1./N * (((H + M)*(H + F)) + ((C + M) * (C + F)))\n",
    "    except:\n",
    "        E = np.nan\n",
    "\n",
    "    try:\n",
    "        POD = H/(H + M)\n",
    "    except:\n",
    "        POD = np.nan   \n",
    "\n",
    "    try:    \n",
    "        FAR = F/(H + F)\n",
    "    except:\n",
    "        FAR = np.nan\n",
    "\n",
    "    try:\n",
    "        BID = (H +F)/(H + M)\n",
    "    except:\n",
    "        BID = np.nan\n",
    "\n",
    "    try:    \n",
    "        HSS = (H + C - E) / (N - E)\n",
    "    except:\n",
    "        HSS = np.nan\n",
    "    try:\n",
    "        HR = (H+C)/N\n",
    "    except:\n",
    "        HR = np.nan\n",
    "    #print (H, M, F, C, N, E)\n",
    "    #print (HSS, FAR, BID, POD)\n",
    "\n",
    "\n",
    "    #\n",
    "    val1[val1<th]=np.nan\n",
    "    val2[val2<th]=np.nan\n",
    "    maske2 = ~np.isnan(val1) & ~np.isnan(val2)\n",
    "    diffi = val1[maske2]-val2[maske2]\n",
    "    summi = np.nansum(val2[maske2])\n",
    "    #bias\n",
    "    bias = np.nansum(diffi)/len(diffi)\n",
    "    rmse = np.sqrt(np.nansum(((diffi)**2.0)/len(diffi)))\n",
    "\n",
    "    slope, intercept, r_value, p_value, std_err = stats.linregress(val1[maske2], val2[maske2])\n",
    "    \n",
    "    PE = 100. * np.nansum(diffi)/summi\n",
    "    \n",
    "    #print(r_value, bias, rmse)\n",
    "    result = {'H': H,\n",
    "              'M': M,\n",
    "              'F': F,\n",
    "              'C': C,\n",
    "              'N': N,\n",
    "              'POD': POD,\n",
    "              'FAR': FAR,\n",
    "              'BID': BID,\n",
    "              'HSS': HSS,\n",
    "              'bias': bias,\n",
    "              'RMSE': rmse,\n",
    "              'corr':r_value,\n",
    "             'HR':HR,\n",
    "             'PE':PE}\n",
    "\n",
    "    #for key, value in result.items():\n",
    "    #    print(key + ':', value)\n",
    "\n",
    "    return result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def plotplot(ax, A, B, tit):\n",
    "    ff=20\n",
    "    from matplotlib.colors import LogNorm\n",
    "\n",
    "\n",
    "\n",
    "    #pl.hist2d(ry1[(rhmin1>bbh1+bbw1/2)&(dprh1>bbh1+bbw1/2)],\n",
    "    #           dpr1[(rhmin1>bbh1+bbw1/2)&(dprh1>bbh1+bbw1/2)],\n",
    "    #          bins=100, vmin=1, vmax=100, norm=LogNorm())\n",
    "\n",
    "    #fig, ax = plt.subplots(figsize=(10,10))\n",
    "    #ax = pl.subplot()\n",
    "    xbins = 10**np.arange(-3, 3, 0.1)\n",
    "    ybins = 10**np.arange(-3, 3, 0.1)\n",
    "\n",
    "    counts, _, _ = np.histogram2d(A,\n",
    "                                  B,\n",
    "                                  bins=(xbins, ybins))\n",
    "\n",
    "    a = np.ma.masked_where(counts.T == 0., counts.T)\n",
    "    pl.pcolormesh(xbins, ybins, a, cmap='jet', vmin=0)\n",
    "    print str(corcorcor(A,B))\n",
    "    \n",
    "    cb = plt.colorbar(shrink=0.5)\n",
    "    cb.set_label('number of samples', fontsize=ff)\n",
    "    cb.ax.tick_params(labelsize=20)\n",
    "\n",
    "    ax.set_xscale('log')\n",
    "    ax.set_yscale('log')\n",
    "\n",
    "    pl.title(tit,fontsize=ff)\n",
    "    pl.ylabel('DPR Rainrate in mm/h',fontsize=ff)\n",
    "    pl.xlabel('Radolan Rainrate in mm/h',fontsize=ff)\n",
    "    \n",
    "    pl.xticks(fontsize=ff)\n",
    "    pl.yticks(fontsize=ff)\n",
    "    pl.plot([0,1000], [0,1000], color='black', label='r: '+ str(corcorcor(A,B)[0])+'\\n bias: '+ str(corcorcor(A,B)[1]))\n",
    "    pl.legend(loc='lower right', markerscale=0)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 0 ns, sys: 13.1 ms, total: 13.1 ms\n",
      "Wall time: 27.5 ms\n"
     ]
    }
   ],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Es gibt 39 Overpasses\n"
     ]
    }
   ],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'%%time\\npfad = glob.glob(\"/automount/ags/velibor/gpmdata/dumpdataV7/dpr_ry_NS/dprrado*\")\\n\\nprint (\\'Es gibt \\' + str(len(pfad)) + \\' Overpasses\\')\\n\\npara_dict = {\\'sr_x\\': [],\\n             \\'sr_y\\': [],\\n             \\'gr_h_max\\': [],\\n             \\'gr_h_mean\\': [],\\n             \\'gr_h_min\\': [],\\n             \\'gr_h_std\\': [],\\n             ##!memory!##\\'gr_pp\\': [],\\n             \\'gr_pp_ipoli\\': [],\\n             ##### \\'gr_time\\': np.array([]),\\n             ##!memory!##\\'gr_xy\\': [],\\n             \\'sr_bbh\\': [],\\n             \\'sr_bbw\\': [],\\n             \\'sr_cfb\\': [],\\n             \\'sr_ele\\': [],\\n             \\'sr_lst\\': [],\\n             \\'sr_pav\\': [],\\n             \\'sr_pes\\': [],\\n             \\'sr_phase\\': [],\\n             \\'sr_pns\\': [],\\n             \\'sr_sth\\': [],\\n             ##### \\'sr_time\\': np.array([]),\\n             \\'sr_times\\': [],\\n             \\'sr_type\\': []\\n}   \\n\\na=0\\n\\nfor dict_name in para_dict.keys():\\n    # DR als liste\\n    print(dict_name)\\n    \\n    if \\'gr\\' in dict_name:\\n        gr_name = \\'GR\\'\\n    else:\\n        gr_name = \\'SR\\'\\n    \\n    for p in pfad:\\n\\n        DR = h5py.File(p, \\'r\\')\\n        # for \\xc3\\xbcber liste\\n        a+=1\\n        # Listen append x.append(x)\\n        para_dict[dict_name].append(DR[gr_name][dict_name][:])\\n        \\n        # ! Colse hdf5 file ! \\n        DR.close()\\n        del(DR)\\n        \\n        \\n    ac = np.concatenate(para_dict[dict_name],axis=0)\\n    print (ac.shape)\\n\\n    np.save(\\'/automount/ags/velibor/gpmdata/dumpdataV7/dpr_ry_NS_compress/\\'+dict_name+\\'.npy\\',ac) \\n    para_dict.pop(dict_name,None)'"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "%%time\n",
    "pfad = glob.glob(\"/automount/ags/velibor/gpmdata/dumpdataV7/dpr_ry_NS/dprrado*\")[1:40]\n",
    "print ('Es gibt ' + str(len(pfad)) + ' Overpasses')\n",
    "\n",
    "para_dict = {'sr_x': [],\n",
    "             'sr_y': [],\n",
    "             'gr_h_max': [],\n",
    "             'gr_h_mean': [],\n",
    "             'gr_h_min': [],\n",
    "             'gr_h_std': [],\n",
    "             ##!memory!##'gr_pp': [],\n",
    "             'gr_pp_ipoli': [],\n",
    "             ##### 'gr_time': np.array([]),\n",
    "             ##!memory!##'gr_xy': [],\n",
    "             'sr_bbh': [],\n",
    "             'sr_bbw': [],\n",
    "             'sr_cfb': [],\n",
    "             'sr_ele': [],\n",
    "             'sr_lst': [],\n",
    "             'sr_pav': [],\n",
    "             'sr_pes': [],\n",
    "             'sr_phase': [],\n",
    "             'sr_pns': [],\n",
    "             'sr_sth': [],\n",
    "             ##### 'sr_time': np.array([]),\n",
    "             'sr_times': [],\n",
    "             'sr_type': []\n",
    "}   \n",
    "\n",
    "\n",
    "time_delta = np.array([])\n",
    "\n",
    "a=0\n",
    "\n",
    "for dict_name in para_dict.keys():\n",
    "    # DR als liste\n",
    "    print(dict_name)\n",
    "    \n",
    "    if 'gr' in dict_name:\n",
    "        gr_name = 'GR'\n",
    "    else:\n",
    "        gr_name = 'SR'\n",
    "    \n",
    "    for p in pfad:\n",
    "\n",
    "        DR = h5py.File(p, 'r')\n",
    "        # for über liste\n",
    "        a+=1\n",
    "        # Listen append x.append(x)\n",
    "        para_dict[dict_name].append(DR[gr_name][dict_name][:])\n",
    "        \n",
    "        if dict_name == 'sr_times':\n",
    "            \n",
    "        \n",
    "        # ! Colse hdf5 file ! \n",
    "        DR.close()\n",
    "        del(DR)\n",
    "        \n",
    "        \n",
    "    ac = np.concatenate(para_dict[dict_name],axis=0)\n",
    "    print (ac.shape)\n",
    "\n",
    "    #np.save('/automount/ags/velibor/gpmdata/dumpdataV7/dpr_ry_NS_compress/'+dict_name+'.npy',ac) \n",
    "    para_dict.pop(dict_name,None)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'%%time\\npfad = glob.glob(\"/automount/ags/velibor/gpmdata/dumpdataV7/dpr_ry_NS/dprrado*\")\\n\\nprint (\\'Es gibt \\' + str(len(pfad)) + \\' Overpasses\\')\\n\\npara_dict = {\\'sr_x\\': [],\\n             \\'sr_y\\': [],\\n             \\'gr_h_max\\': [],\\n             \\'gr_h_mean\\': [],\\n             \\'gr_h_min\\': [],\\n             \\'gr_h_std\\': [],\\n             ##!memory!##\\'gr_pp\\': [],\\n             \\'gr_pp_ipoli\\': [],\\n             ##### \\'gr_time\\': np.array([]),\\n             ##!memory!##\\'gr_xy\\': [],\\n             \\'sr_bbh\\': [],\\n             \\'sr_bbw\\': [],\\n             \\'sr_cfb\\': [],\\n             \\'sr_ele\\': [],\\n             \\'sr_lst\\': [],\\n             \\'sr_pav\\': [],\\n             \\'sr_pes\\': [],\\n             \\'sr_phase\\': [],\\n             \\'sr_pns\\': [],\\n             \\'sr_sth\\': [],\\n             ##### \\'sr_time\\': np.array([]),\\n             \\'sr_times\\': [],\\n             \\'sr_type\\': []\\n}   \\n\\na=0\\n\\nfor dict_name in para_dict.keys():\\n    # DR als liste\\n    print(dict_name)\\n    \\n    if \\'gr\\' in dict_name:\\n        gr_name = \\'GR\\'\\n    else:\\n        gr_name = \\'SR\\'\\n    \\n    for p in pfad:\\n\\n        DR = h5py.File(p, \\'r\\')\\n        # for \\xc3\\xbcber liste\\n        a+=1\\n        # Listen append x.append(x)\\n        para_dict[dict_name].append(DR[gr_name][dict_name][:])\\n        \\n        # ! Colse hdf5 file ! \\n        DR.close()\\n        del(DR)\\n        \\n        \\n    ac = np.concatenate(para_dict[dict_name],axis=0)\\n    print (ac.shape)\\n\\n    np.save(\\'/automount/ags/velibor/gpmdata/dumpdataV7/dpr_ry_NS_compress/\\'+dict_name+\\'.npy\\',ac) \\n    para_dict.pop(dict_name,None)'"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\"\"\"%%time\n",
    "pfad = glob.glob(\"/automount/ags/velibor/gpmdata/dumpdataV7/dpr_ry_NS/dprrado*\")\n",
    "\n",
    "print ('Es gibt ' + str(len(pfad)) + ' Overpasses')\n",
    "\n",
    "para_dict = {'sr_x': [],\n",
    "             'sr_y': [],\n",
    "             'gr_h_max': [],\n",
    "             'gr_h_mean': [],\n",
    "             'gr_h_min': [],\n",
    "             'gr_h_std': [],\n",
    "             ##!memory!##'gr_pp': [],\n",
    "             'gr_pp_ipoli': [],\n",
    "             ##### 'gr_time': np.array([]),\n",
    "             ##!memory!##'gr_xy': [],\n",
    "             'sr_bbh': [],\n",
    "             'sr_bbw': [],\n",
    "             'sr_cfb': [],\n",
    "             'sr_ele': [],\n",
    "             'sr_lst': [],\n",
    "             'sr_pav': [],\n",
    "             'sr_pes': [],\n",
    "             'sr_phase': [],\n",
    "             'sr_pns': [],\n",
    "             'sr_sth': [],\n",
    "             ##### 'sr_time': np.array([]),\n",
    "             'sr_times': [],\n",
    "             'sr_type': []\n",
    "}   \n",
    "\n",
    "a=0\n",
    "\n",
    "for dict_name in para_dict.keys():\n",
    "    # DR als liste\n",
    "    print(dict_name)\n",
    "    \n",
    "    if 'gr' in dict_name:\n",
    "        gr_name = 'GR'\n",
    "    else:\n",
    "        gr_name = 'SR'\n",
    "    \n",
    "    for p in pfad:\n",
    "\n",
    "        DR = h5py.File(p, 'r')\n",
    "        # for über liste\n",
    "        a+=1\n",
    "        # Listen append x.append(x)\n",
    "        para_dict[dict_name].append(DR[gr_name][dict_name][:])\n",
    "        \n",
    "        # ! Colse hdf5 file ! \n",
    "        DR.close()\n",
    "        del(DR)\n",
    "        \n",
    "        \n",
    "    ac = np.concatenate(para_dict[dict_name],axis=0)\n",
    "    print (ac.shape)\n",
    "\n",
    "    np.save('/automount/ags/velibor/gpmdata/dumpdataV7/dpr_ry_NS_compress/'+dict_name+'.npy',ac) \n",
    "    para_dict.pop(dict_name,None)\"\"\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "xl = np.load('/automount/ags/velibor/gpmdata/dumpdataV7/dpr_ry_NS_compress/sr_x.npy')\n",
    "yl = np.load('/automount/ags/velibor/gpmdata/dumpdataV7/dpr_ry_NS_compress/sr_y.npy')\n",
    "\n",
    "sr_pns = np.load('/automount/ags/velibor/gpmdata/dumpdataV7/dpr_ry_NS_compress/sr_pns.npy')\n",
    "sr_pav = np.load('/automount/ags/velibor/gpmdata/dumpdataV7/dpr_ry_NS_compress/sr_pav.npy')\n",
    "gr_pns = np.load('/automount/ags/velibor/gpmdata/dumpdataV7/dpr_ry_NS_compress/gr_pp_ipoli.npy')\n",
    "\n",
    "gr_hmin = np.load('/automount/ags/velibor/gpmdata/dumpdataV7/dpr_ry_NS_compress/gr_h_min.npy')\n",
    "gr_hmax = np.load('/automount/ags/velibor/gpmdata/dumpdataV7/dpr_ry_NS_compress/gr_h_max.npy')\n",
    "sr_ele = np.load('/automount/ags/velibor/gpmdata/dumpdataV7/dpr_ry_NS_compress/sr_ele.npy')\n",
    "\n",
    "sr_typ = np.load('/automount/ags/velibor/gpmdata/dumpdataV7/dpr_ry_NS_compress/sr_type.npy')\n",
    "sr_pha = np.load('/automount/ags/velibor/gpmdata/dumpdataV7/dpr_ry_NS_compress/sr_phase.npy')\n",
    "sr_cfb = np.load('/automount/ags/velibor/gpmdata/dumpdataV7/dpr_ry_NS_compress/sr_cfb.npy')\n",
    "\n",
    "sr_bbh = np.load('/automount/ags/velibor/gpmdata/dumpdataV7/dpr_ry_NS_compress/sr_bbh.npy')\n",
    "sr_bbw = np.load('/automount/ags/velibor/gpmdata/dumpdataV7/dpr_ry_NS_compress/sr_bbw.npy')\n",
    "sr_sth = np.load('/automount/ags/velibor/gpmdata/dumpdataV7/dpr_ry_NS_compress/sr_sth.npy')\n",
    "\n",
    "sr_lst = np.load('/automount/ags/velibor/gpmdata/dumpdataV7/dpr_ry_NS_compress/sr_lst.npy')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "skill_score(gr_pns, sr_pav, 0.5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "th_set = np.arange(0,1,0.01)\n",
    "#h_set = np.array([0,1,0.1])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "res1 = res2 = np.array([])\n",
    "\n",
    "for i in th_set:\n",
    "    #print i\n",
    "    res1 = np.append(res1, skill_score(sr_pns[sr_pha==2], gr_pns[sr_pha==2], th=i))\n",
    "    #print skill_score(sr_pns, gr_pns, th=i)['corr']\n",
    "    res2 = np.append(res2, skill_score(sr_pav[sr_pha==2], gr_pns[sr_pha==2], th=i))\n",
    "    #print skill_score(sr_pav, gr_pns, th=i)['corr']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "res2[1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "pl.figure(figsize=(15,10))\n",
    "for i in range(len(res1)):\n",
    "    #print th_set[i], res1[i][wort] \n",
    "    pl.subplot(3,4,1)\n",
    "    wort = 'HR'\n",
    "    #pl.yscale('log')\n",
    "    pl.scatter(th_set[i], res1[i][wort], color='green')\n",
    "    #pl.scatter(th_set[i], res2[i][wort], color = 'blue')\n",
    "    pl.axhline(0)\n",
    "    pl.axvline(0.2, color='red')\n",
    "    pl.axvline(0.5, color='black')\n",
    "    pl.ylabel(wort + ' - sr-gr')\n",
    "    pl.xlabel('Threshold in mm/h')\n",
    "    pl.title('pns: green, pav:blue')\n",
    "    pl.grid()\n",
    "    \n",
    "    pl.subplot(3,4,2)\n",
    "    wort = 'corr'\n",
    "    pl.scatter(th_set[i], res1[i][wort], color='green')\n",
    "    #pl.scatter(th_set[i], res2[i][wort], color = 'blue')\n",
    "    pl.axhline(0)\n",
    "    pl.axvline(0.2, color='red')\n",
    "    pl.axvline(0.5, color='black')\n",
    "    pl.ylabel(wort + ' - sr-gr')\n",
    "    pl.xlabel('Threshold in mm/h')\n",
    "    pl.grid()\n",
    "    \n",
    "    pl.subplot(3,4,3)\n",
    "    wort = 'bias'\n",
    "    pl.scatter(th_set[i], res1[i][wort], color='green')\n",
    "    #pl.scatter(th_set[i], res2[i][wort], color = 'blue')\n",
    "    pl.axhline(0)\n",
    "    pl.axvline(0.2, color='red')\n",
    "    pl.axvline(0.5, color='black')\n",
    "    pl.ylabel(wort + ' - sr-gr')\n",
    "    pl.xlabel('Threshold in mm/h')\n",
    "    pl.grid()\n",
    "    \n",
    "    pl.subplot(3,4,4)\n",
    "    wort = 'RMSE'\n",
    "    pl.scatter(th_set[i], res1[i][wort], color='green')\n",
    "    #pl.scatter(th_set[i], res2[i][wort], color = 'blue')\n",
    "    pl.axhline(0)\n",
    "    pl.axvline(0.2, color='red')\n",
    "    pl.axvline(0.5, color='black')\n",
    "    pl.ylabel(wort + ' - sr-gr')\n",
    "    pl.xlabel('Threshold in mm/h')\n",
    "    pl.grid()\n",
    "    \n",
    "    pl.subplot(3,4,5)\n",
    "    wort = 'POD'\n",
    "    pl.scatter(th_set[i], res1[i][wort], color='green')\n",
    "    #pl.scatter(th_set[i], res2[i][wort], color = 'blue')\n",
    "    pl.axhline(0)\n",
    "    pl.axvline(0.2, color='red')\n",
    "    pl.axvline(0.5, color='black')\n",
    "    pl.ylabel(wort + ' - sr-gr')\n",
    "    pl.xlabel('Threshold in mm/h')\n",
    "    pl.grid()\n",
    "    \n",
    "    pl.subplot(3,4,6)\n",
    "    wort = 'BID'\n",
    "    pl.scatter(th_set[i], res1[i][wort], color='green')\n",
    "    #pl.scatter(th_set[i], res2[i][wort], color = 'blue')\n",
    "    pl.axhline(0)\n",
    "    pl.axvline(0.2, color='red')\n",
    "    pl.axvline(0.5, color='black')\n",
    "    pl.ylabel(wort + ' - sr-gr')\n",
    "    pl.xlabel('Threshold in mm/h')\n",
    "    pl.grid()\n",
    "    \n",
    "    pl.subplot(3,4,7)\n",
    "    wort = 'FAR'\n",
    "    pl.scatter(th_set[i], res1[i][wort], color='green')\n",
    "    #pl.scatter(th_set[i], res2[i][wort], color = 'blue')\n",
    "    pl.axhline(0)\n",
    "    pl.axvline(0.2, color='red')\n",
    "    pl.axvline(0.5, color='black')\n",
    "    pl.ylabel(wort + ' - sr-gr')\n",
    "    pl.xlabel('Threshold in mm/h')\n",
    "    pl.grid()\n",
    "    \n",
    "    pl.subplot(3,4,8)\n",
    "    wort = 'HSS'\n",
    "    pl.scatter(th_set[i], res1[i][wort], color='green')\n",
    "    #pl.scatter(th_set[i], res2[i][wort], color = 'blue')\n",
    "    pl.axhline(0)\n",
    "    pl.axvline(0.2, color='red')\n",
    "    pl.axvline(0.5, color='black')\n",
    "    pl.ylabel(wort + ' - sr-gr')\n",
    "    pl.xlabel('Threshold in mm/h')\n",
    "    pl.grid()\n",
    "    \n",
    "    pl.subplot(3,4,9)\n",
    "    wort = 'H'\n",
    "    pl.scatter(th_set[i], res1[i][wort], color='green')\n",
    "    #pl.scatter(th_set[i], res2[i][wort], color = 'blue')\n",
    "    pl.axhline(0)\n",
    "    pl.axvline(0.2, color='red')\n",
    "    pl.axvline(0.5, color='black')\n",
    "    pl.ylabel(wort + ' - sr-gr')\n",
    "    pl.xlabel('Threshold in mm/h')\n",
    "    pl.grid()\n",
    "    \n",
    "    pl.subplot(3,4,10)\n",
    "    wort = 'F'\n",
    "    pl.scatter(th_set[i], res1[i][wort], color='green')\n",
    "    #pl.scatter(th_set[i], res2[i][wort], color = 'blue')\n",
    "    pl.axhline(0)\n",
    "    pl.axvline(0.2, color='red')\n",
    "    pl.axvline(0.5, color='black')\n",
    "    pl.ylabel(wort + ' - sr-gr')\n",
    "    pl.xlabel('Threshold in mm/h')\n",
    "    pl.grid()\n",
    "    \n",
    "    pl.subplot(3,4,11)\n",
    "    wort = 'M'\n",
    "    pl.scatter(th_set[i], res1[i][wort], color='green')\n",
    "    #pl.scatter(th_set[i], res2[i][wort], color = 'blue')\n",
    "    pl.axhline(0)\n",
    "    pl.axvline(0.2, color='red')\n",
    "    pl.axvline(0.5, color='black')\n",
    "    pl.ylabel(wort + ' - sr-gr')\n",
    "    pl.xlabel('Threshold in mm/h')\n",
    "    pl.grid()\n",
    "    \n",
    "    pl.subplot(3,4,12)\n",
    "    wort = 'PE'\n",
    "    pl.scatter(th_set[i], res1[i][wort], color='green')\n",
    "    #pl.scatter(th_set[i], res2[i][wort], color = 'blue')\n",
    "    pl.axhline(0)\n",
    "    pl.axvline(0.2, color='red')\n",
    "    pl.axvline(0.5, color='black')\n",
    "    pl.ylabel(wort + ' - sr-gr')\n",
    "    pl.xlabel('Threshold in mm/h')\n",
    "    pl.grid()\n",
    "\n",
    "pl.tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0 2]\n"
     ]
    }
   ],
   "source": [
    "#Umwandeln von typ\n",
    "\"\"\"\n",
    "a: Main rain type. (a=1,2,3),\n",
    "b: 0,\n",
    "c: 0,\n",
    "d: V rain type,\n",
    "e: H rain type,\n",
    "f: BB in outter sath detected,\n",
    "g: Shallow rain,\n",
    "h: Small size cell\n",
    "\n",
    "Major rain type = typePrecip/10000000\n",
    "    = 1    stratiform\n",
    "    = 2    convective\n",
    "    = 3    other\n",
    "\n",
    "b = 0: BB not detected,\n",
    "    1: BB detected.\n",
    "\n",
    "g = 0: No shallow rain,\n",
    "    1: Shallow isolated,\n",
    "    3: Shallow non-isolated.\n",
    "\n",
    "h = 0: No small size cell,\n",
    "    1: Single cell,\n",
    "    2: Small size cell consisting of two adjacent pixels.\n",
    "\"\"\"\n",
    "typ = sr_typ.copy()\n",
    "\n",
    "typ = typ//10000000\n",
    "\n",
    "# Umwandeln von Phase.. Achtung es ist schon Umgewandelt\n",
    "\"\"\"\n",
    "phaseNearSurface < 100 Temperature(C)=phaseNearSurface-100\n",
    "phaseNearSurface > 200 Temperature(C)=phaseNearSurface-200\n",
    "phaseNearSurface = 100 Top of the bright band\n",
    "phaseNearSurface = 200 Bottom of the bright band\n",
    "phaseNearSurface = 125 is used for the range bins between\n",
    "                    the top and peak of bright band\n",
    "phaseNearSurface = 175 is used for the range bins between\n",
    "                    the peak and bottom of bright band\n",
    "                    \n",
    "Integer values of phaseNearSurface/100 =\n",
    "0 - solid\n",
    "1 - mixed phase\n",
    "2 - liquid\n",
    "255 - Missing\n",
    "\"\"\"\n",
    "pha = sr_pha.copy()\n",
    "pha = pha/100\n",
    "\n",
    "print(np.unique(pha))\n",
    "# Umwandeln von clutterfree bottom\n",
    "#cfb = sr_cfb.copy()\n",
    "#cfb[cfb<0]=np.nan\n",
    "#cfb = 125.*(176. - cfb)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "pl.figure(figsize=(18,14))\n",
    "xbins = 10**np.arange(-3, 3, 0.1)\n",
    "\n",
    "pl.subplot(4,3,1)\n",
    "pl.xscale('log')\n",
    "pl.hist(sr_pns[sr_pns>0], bins=xbins)\n",
    "pl.xlabel('Rainrate (mm/h)', fontsize=20)\n",
    "pl.title('PNS', fontsize=20)\n",
    "pl.xticks(fontsize=20)\n",
    "pl.yticks(fontsize=20)\n",
    "\n",
    "pl.subplot(4,3,2)\n",
    "pl.xscale('log')\n",
    "pl.hist(sr_pav[sr_pav>0], bins=xbins)\n",
    "pl.xlabel('Rainrate (mm/h)', fontsize=20)\n",
    "pl.title('PAV', fontsize=20)\n",
    "pl.xticks(fontsize=20)\n",
    "pl.yticks(fontsize=20)\n",
    "\n",
    "pl.subplot(4,3,3)\n",
    "pl.xscale('log')\n",
    "pl.hist(gr_pns[gr_pns>0], bins=xbins)\n",
    "pl.xlabel('Rainrate (mm/h)', fontsize=20)\n",
    "pl.title('GRP', fontsize=20)\n",
    "pl.xticks(fontsize=20)\n",
    "pl.yticks(fontsize=20)\n",
    "\n",
    "pl.subplot(4,3,4)\n",
    "pl.hist(sr_bbh[sr_bbh>0]/1000, bins=50)\n",
    "pl.xlabel('Height (km)', fontsize=20)\n",
    "pl.title('BBH', fontsize=20)\n",
    "pl.xticks(fontsize=20)\n",
    "pl.yticks(fontsize=20)\n",
    "\n",
    "pl.subplot(4,3,5)\n",
    "pl.hist(sr_bbw[sr_bbw>0]/1000, bins=50)\n",
    "pl.xlabel('Width (km)', fontsize=20)\n",
    "pl.title('BBW', fontsize=20)\n",
    "pl.xticks(fontsize=20)\n",
    "pl.yticks(fontsize=20)\n",
    "\n",
    "pl.subplot(4,3,6)\n",
    "pl.hist(sr_sth[sr_sth>=0]/1000, bins=50)\n",
    "pl.xlabel('Height (km)', fontsize=20)\n",
    "pl.title('STH', fontsize=20)\n",
    "pl.xticks(fontsize=20)\n",
    "pl.yticks(fontsize=20)\n",
    "\n",
    "pl.subplot(4,3,7)\n",
    "pl.hist(sr_lst[(sr_lst>=0) & (sr_lst!=-9999)], bins=30)\n",
    "pl.xlim(0,400)\n",
    "pl.xlabel('Ocean ----  Land  ---- Coast ---- InlandWater')\n",
    "pl.title('LST', fontsize=20)\n",
    "pl.xticks(fontsize=20)\n",
    "pl.yticks(fontsize=20)\n",
    "\n",
    "pl.subplot(4,3,8)\n",
    "pl.hist(typ[typ>0])\n",
    "pl.xlabel('Strati ---- Convectiv ---- Other')\n",
    "pl.title('TYP', fontsize=20)\n",
    "pl.xticks(fontsize=20)\n",
    "pl.yticks(fontsize=20)\n",
    "\n",
    "pl.subplot(4,3,9)\n",
    "pl.yscale('log')\n",
    "pl.hist(sr_pha[sr_pha>=0], bins=5)\n",
    "pl.xlabel('Solid ---- Mixed  ---- Liquid')\n",
    "pl.title('PHA', fontsize=20)\n",
    "pl.xticks(fontsize=20)\n",
    "pl.yticks(fontsize=20)\n",
    "\n",
    "pl.subplot(4,3,10)\n",
    "srh = sr_cfb - sr_ele\n",
    "pl.hist(srh[srh>=0]/1000, bins=30)\n",
    "pl.xlim(0,6)\n",
    "pl.xlabel('Height (km)', fontsize=20)\n",
    "pl.title('CFB', fontsize=20)\n",
    "pl.xticks(fontsize=20)\n",
    "pl.yticks(fontsize=20)\n",
    "\n",
    "pl.subplot(4,3,12)\n",
    "pl.hist(gr_hmin[gr_hmin>=0]/1000., bins=30)\n",
    "pl.xlim(0,6)\n",
    "pl.xlabel('Height (km)', fontsize=20)\n",
    "pl.title('GRH', fontsize=20)\n",
    "pl.xticks(fontsize=20)\n",
    "pl.yticks(fontsize=20)\n",
    "\n",
    "pl.subplot(4,3,11)\n",
    "pl.hist(sr_ele[sr_ele>0]/1000, bins=30)\n",
    "pl.xlim(0,6)\n",
    "pl.xlabel('Height (km)', fontsize=20)\n",
    "pl.title('ELE', fontsize=20)\n",
    "pl.xticks(fontsize=20)\n",
    "pl.yticks(fontsize=20)\n",
    "\n",
    "pl.tight_layout()\n",
    "pl.show()"
   ]
  },
  {
   "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": false
   },
   "outputs": [],
   "source": [
    "# new"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Hoehen unterschied zwischen Sr un GR der 3 Jahre\n",
    "-------------------------------------------------------------"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "srh = sr_cfb-sr_ele\n",
    "grh = gr_hmin \n",
    "ff = 20\n",
    "def get_miub_cmap():\n",
    "    import matplotlib.colors as col\n",
    "    startcolor = 'white'  # a dark olive\n",
    "    color1 = '#8ec7ff'#'cyan'    # a bright yellow\n",
    "    color2 = 'dodgerblue'\n",
    "    color3 = 'lime'\n",
    "    color4 = 'yellow'\n",
    "    color5 = 'darkorange'\n",
    "    color6 = 'red'\n",
    "    color7 = 'purple'\n",
    "    #color6 = 'grey'\n",
    "    endcolor = 'darkmagenta'    # medium dark red\n",
    "    colors = [startcolor, color1, color2, color3, color4, color5, color6, endcolor]\n",
    "    return col.LinearSegmentedColormap.from_list('miub1',colors)\n",
    "\n",
    "\n",
    "mh = ~np.isnan(grh) & ~np.isnan(srh)\n",
    "\n",
    "pl.hist2d(srh[mh], grh[mh], bins=200, cmap=get_miub_cmap() )\n",
    "pl.xlabel(\"SR measuring height in m\", fontsize=ff)\n",
    "pl.ylabel(\"GR measuring height in m\", fontsize=ff)\n",
    "cb = pl.colorbar(shrink=0.5)\n",
    "cb.set_label('number of samples', fontsize=ff)\n",
    "cb.ax.tick_params(labelsize=20)\n",
    "pl.xlim(0,4000)\n",
    "pl.ylim(0,4000)\n",
    "pl.xticks(fontsize=ff)\n",
    "pl.yticks(fontsize=ff)\n",
    "\n",
    "pl.plot([0,6000],[0,6000], linestyle=\"solid\", color=\"black\")\n",
    "pl.tight_layout()\n",
    "pl.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "pl.hist(grh[mh], bins=200, histtype='step', normed=False,label='GR', color='blue')\n",
    "pl.hist(srh[mh], bins=200, histtype='step', normed=False, label='SR', color='green')\n",
    "pl.legend()\n",
    "pl.xlabel(\"Measuring height in m\", fontsize=ff)\n",
    "pl.grid()\n",
    "pl.xticks(fontsize=ff)\n",
    "pl.yticks(fontsize=ff)\n",
    "pl.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "srh = sr_cfb-sr_ele\n",
    "grh = gr_hmin \n",
    "\n",
    "hdiff = grh - srh\n",
    "\n",
    "ff = 20\n",
    "def get_miub_cmap():\n",
    "    import matplotlib.colors as col\n",
    "    startcolor = 'white'  # a dark olive\n",
    "    color1 = '#8ec7ff'#'cyan'    # a bright yellow\n",
    "    color2 = 'dodgerblue'\n",
    "    color3 = 'lime'\n",
    "    color4 = 'yellow'\n",
    "    color5 = 'darkorange'\n",
    "    color6 = 'red'\n",
    "    color7 = 'purple'\n",
    "    #color6 = 'grey'\n",
    "    endcolor = 'darkmagenta'    # medium dark red\n",
    "    colors = [startcolor, color1, color2, color3, color4, color5, color6, endcolor]\n",
    "    return col.LinearSegmentedColormap.from_list('miub1',colors)\n",
    "\n",
    "\n",
    "mh = ~np.isnan(hdiff) & ~np.isnan(yl)\n",
    "\n",
    "pl.figure(figsize=(12,12))\n",
    "pl.hist2d(hdiff[mh], yl[mh], bins=200, cmap=get_miub_cmap() )\n",
    "pl.xlabel(\"GR - SR measuring height in m\", fontsize=ff)\n",
    "pl.ylabel(\"y in km\", fontsize=ff)\n",
    "cb = pl.colorbar(shrink=0.5)\n",
    "cb.set_label('number of samples', fontsize=ff)\n",
    "cb.ax.tick_params(labelsize=20)\n",
    "pl.xlim(-4000,4000)\n",
    "\n",
    "pl.xticks(fontsize=ff)\n",
    "pl.yticks(fontsize=ff)\n",
    "pl.axvline(0, linestyle=\"solid\", color=\"black\")\n",
    "pl.axvline(-1000, linestyle=\"-.\", color=\"black\")\n",
    "pl.axvline(1000, linestyle=\"-.\", color=\"black\")\n",
    "\n",
    "pl.tight_layout()\n",
    "pl.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "srh = sr_cfb-sr_ele\n",
    "grh = gr_hmin \n",
    "\n",
    "hdiff = grh - srh\n",
    "\n",
    "ff = 20\n",
    "def get_miub_cmap():\n",
    "    import matplotlib.colors as col\n",
    "    startcolor = 'white'  # a dark olive\n",
    "    color1 = '#8ec7ff'#'cyan'    # a bright yellow\n",
    "    color2 = 'dodgerblue'\n",
    "    color3 = 'lime'\n",
    "    color4 = 'yellow'\n",
    "    color5 = 'darkorange'\n",
    "    color6 = 'red'\n",
    "    color7 = 'purple'\n",
    "    #color6 = 'grey'\n",
    "    endcolor = 'darkmagenta'    # medium dark red\n",
    "    colors = [startcolor, color1, color2, color3, color4, color5, color6, endcolor]\n",
    "    return col.LinearSegmentedColormap.from_list('miub1',colors)\n",
    "\n",
    "\n",
    "mh = ~np.isnan(hdiff) & ~np.isnan(xl)\n",
    "\n",
    "pl.figure(figsize=(12,12))\n",
    "pl.hist2d(xl[mh], hdiff[mh], bins=200, cmap=get_miub_cmap() )\n",
    "pl.ylabel(\"GR - SR measuring height in m\", fontsize=ff)\n",
    "pl.xlabel(\"x in km\", fontsize=ff)\n",
    "cb = pl.colorbar(shrink=0.5)\n",
    "cb.set_label('number of samples', fontsize=ff)\n",
    "cb.ax.tick_params(labelsize=20)\n",
    "#pl.xlim(-4000,4000)\n",
    "\n",
    "pl.xticks(fontsize=ff)\n",
    "pl.yticks(fontsize=ff)\n",
    "pl.axhline(0, linestyle=\"solid\", color=\"black\")\n",
    "pl.axhline(-1000, linestyle=\"-.\", color=\"black\")\n",
    "pl.axhline(1000, linestyle=\"-.\", color=\"black\")\n",
    "\n",
    "pl.tight_layout()\n",
    "pl.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Höhenabhängigkeitsbeweis?\n",
    "-------------------------------------------"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/velibor/anaconda2/envs/wradlib/lib/python2.7/site-packages/numpy/core/_methods.py:29: RuntimeWarning: invalid value encountered in reduce\n",
      "  return umr_minimum(a, axis, None, out, keepdims)\n",
      "/home/velibor/anaconda2/envs/wradlib/lib/python2.7/site-packages/numpy/core/_methods.py:26: RuntimeWarning: invalid value encountered in reduce\n",
      "  return umr_maximum(a, axis, None, out, keepdims)\n",
      "/home/velibor/anaconda2/envs/wradlib/lib/python2.7/site-packages/numpy/lib/function_base.py:991: RuntimeWarning: invalid value encountered in greater_equal\n",
      "  not_smaller_than_edge = (sample[:, i] >= edges[i][-1])\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(0.62, 0.3)\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/velibor/anaconda2/envs/wradlib/lib/python2.7/site-packages/ipykernel/__main__.py:16: RuntimeWarning: invalid value encountered in greater\n",
      "/home/velibor/anaconda2/envs/wradlib/lib/python2.7/site-packages/ipykernel/__main__.py:17: RuntimeWarning: invalid value encountered in less_equal\n",
      "/home/velibor/anaconda2/envs/wradlib/lib/python2.7/site-packages/ipykernel/__main__.py:17: RuntimeWarning: invalid value encountered in greater\n",
      "/home/velibor/anaconda2/envs/wradlib/lib/python2.7/site-packages/ipykernel/__main__.py:18: RuntimeWarning: invalid value encountered in greater\n",
      "/home/velibor/anaconda2/envs/wradlib/lib/python2.7/site-packages/ipykernel/__main__.py:18: RuntimeWarning: invalid value encountered in less_equal\n",
      "/home/velibor/anaconda2/envs/wradlib/lib/python2.7/site-packages/ipykernel/__main__.py:19: RuntimeWarning: invalid value encountered in less_equal\n",
      "/home/velibor/anaconda2/envs/wradlib/lib/python2.7/site-packages/ipykernel/__main__.py:64: RuntimeWarning: invalid value encountered in less\n",
      "/home/velibor/anaconda2/envs/wradlib/lib/python2.7/site-packages/ipykernel/__main__.py:65: RuntimeWarning: invalid value encountered in less\n",
      "/home/velibor/anaconda2/envs/wradlib/lib/python2.7/site-packages/ipykernel/__main__.py:35: RuntimeWarning: invalid value encountered in less\n",
      "/home/velibor/anaconda2/envs/wradlib/lib/python2.7/site-packages/ipykernel/__main__.py:36: RuntimeWarning: invalid value encountered in less\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(0.71, 0.44)\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/velibor/anaconda2/envs/wradlib/lib/python2.7/site-packages/ipykernel/__main__.py:38: RuntimeWarning: invalid value encountered in less\n"
     ]
    }
   ],
   "source": [
    "ff=20\n",
    "\n",
    "gr, sr = gr_pns.copy(), sr_pns.copy()\n",
    "hmin, cfb, ele = gr_hmin.copy(), sr_cfb.copy(), sr_ele.copy()\n",
    "\n",
    "gr[gr==0.]=np.nan\n",
    "sr[sr==0.]=np.nan\n",
    "\n",
    "#stratiform 1\n",
    "gr[typ!=1.]=np.nan\n",
    "sr[typ!=1.]=np.nan\n",
    "\n",
    "#liquid 2\n",
    "gr[pha!=2.]=np.nan\n",
    "sr[pha!=2.]=np.nan\n",
    "\n",
    "#pos = (gr_hmin - (sr_cfb - sr_ele)) < 500.\n",
    "pl.figure(figsize=(15,5))\n",
    "ax1 = pl.subplot(1,2,1, aspect='equal')\n",
    "#pl.xscale('log')\n",
    "#pl.yscale('log')\n",
    "#pl.xlim(1e-3,1000)\n",
    "#pl.ylim(1e-3,1000)\n",
    "plotplot(ax1, gr, sr, 'All Measurements')\n",
    "#pl.plot([0,1000], [0,1000], color='black')\n",
    "skill_score(gr, sr, 0.5)\n",
    "\n",
    "hdiff = 100.\n",
    "\n",
    "ax2 = pl.subplot(1,2,2, aspect='equal')\n",
    "#pl.xscale('log')\n",
    "#pl.yscale('log')\n",
    "#pl.xlim(1e-3,1000)\n",
    "#pl.ylim(1e-3,1000)\n",
    "plotplot(ax2, gr[abs((hmin - (cfb - ele))) < hdiff],\n",
    "         sr[abs((hmin - (cfb - ele))) < hdiff], 'Only Measurements with diff \\n in height: '+str(hdiff)+' m')\n",
    "#pl.plot([0,1000], [0,1000], color='black')\n",
    "skill_score(gr[abs((hmin - (cfb - ele))) < hdiff], sr[abs((hmin - (cfb - ele))) < hdiff], 0.5)\n",
    "\n",
    "pl.tight_layout()\n",
    "pl.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "hdiff = 100.\n",
    "\n",
    "for hdiff in np.arange(100,2000,100):\n",
    "    print hdiff\n",
    "    SS = skill_score(gr[abs((hmin - (cfb - ele))) < hdiff], sr[abs((hmin - (cfb - ele))) < hdiff], 0.5)\n",
    "    plt.scatter(hdiff, SS['corr'], color='blue')\n",
    "    plt.scatter(hdiff, SS['POD'], color='red')\n",
    "    plt.scatter(hdiff, SS['HR'], color='green')\n",
    "    plt.scatter(hdiff, SS['HSS'], color='black')\n",
    "    plt.scatter(hdiff, SS['FAR'], color='grey')\n",
    "    \n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "hdiff = 100.\n",
    "\n",
    "for hdiff in np.arange(100,2000,100):\n",
    "    print hdiff\n",
    "    SS = skill_score(gr[(abs((hmin - (cfb - ele))) < hdiff) & (abs((hmin - (cfb - ele))) > hdiff-100.)],\n",
    "                     sr[(abs((hmin - (cfb - ele))) < hdiff) & (abs((hmin - (cfb - ele))) > hdiff-100.)],\n",
    "                     0.5)\n",
    "    plt.scatter(hdiff, SS['N'], color='blue')\n",
    "    plt.scatter(hdiff, SS['H'], color='red')\n",
    "    plt.scatter(hdiff, SS['F'], color='green')\n",
    "    plt.scatter(hdiff, SS['M'], color='black')\n",
    "    plt.scatter(hdiff, SS['C'], color='grey')\n",
    "    \n",
    "\n",
    "plt.xlabel('hdiff in m')    \n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "125.0\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/velibor/anaconda2/envs/wradlib/lib/python2.7/site-packages/ipykernel/__main__.py:11: RuntimeWarning: invalid value encountered in less\n",
      "/home/velibor/anaconda2/envs/wradlib/lib/python2.7/site-packages/ipykernel/__main__.py:11: RuntimeWarning: invalid value encountered in greater\n",
      "/home/velibor/anaconda2/envs/wradlib/lib/python2.7/site-packages/ipykernel/__main__.py:12: RuntimeWarning: invalid value encountered in less\n",
      "/home/velibor/anaconda2/envs/wradlib/lib/python2.7/site-packages/ipykernel/__main__.py:12: RuntimeWarning: invalid value encountered in greater\n",
      "/home/velibor/anaconda2/envs/wradlib/lib/python2.7/site-packages/ipykernel/__main__.py:16: RuntimeWarning: invalid value encountered in greater\n",
      "/home/velibor/anaconda2/envs/wradlib/lib/python2.7/site-packages/ipykernel/__main__.py:17: RuntimeWarning: invalid value encountered in less_equal\n",
      "/home/velibor/anaconda2/envs/wradlib/lib/python2.7/site-packages/ipykernel/__main__.py:17: RuntimeWarning: invalid value encountered in greater\n",
      "/home/velibor/anaconda2/envs/wradlib/lib/python2.7/site-packages/ipykernel/__main__.py:18: RuntimeWarning: invalid value encountered in greater\n",
      "/home/velibor/anaconda2/envs/wradlib/lib/python2.7/site-packages/ipykernel/__main__.py:18: RuntimeWarning: invalid value encountered in less_equal\n",
      "/home/velibor/anaconda2/envs/wradlib/lib/python2.7/site-packages/ipykernel/__main__.py:19: RuntimeWarning: invalid value encountered in less_equal\n",
      "/home/velibor/anaconda2/envs/wradlib/lib/python2.7/site-packages/ipykernel/__main__.py:64: RuntimeWarning: invalid value encountered in less\n",
      "/home/velibor/anaconda2/envs/wradlib/lib/python2.7/site-packages/ipykernel/__main__.py:65: RuntimeWarning: invalid value encountered in less\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "250.0\n",
      "375.0\n",
      "500.0\n",
      "625.0\n",
      "750.0\n",
      "875.0\n",
      "1000.0\n",
      "1125.0\n",
      "1250.0\n",
      "1375.0\n",
      "1500.0\n",
      "1625.0\n",
      "1750.0\n",
      "1875.0\n",
      "2000.0\n",
      "2125.0\n",
      "2250.0\n",
      "2375.0\n",
      "2500.0\n",
      "2625.0\n",
      "2750.0\n",
      "2875.0\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "\"    plt.subplot(3,1,2)\\n    plt.scatter(hdiff, SS['N'],  marker='o', s=40)\\n    plt.scatter(hdiff, SS['H'],  marker='^', s=40)\\n    plt.scatter(hdiff, SS['F'],  marker='h', s=40)\\n    plt.scatter(hdiff, SS['M'],  marker='*', s=40)\\n    plt.scatter(hdiff, SS['C'],  marker='+', s=40)\\n    \\n    plt.legend()\\n    plt.grid()    \\n    plt.xlabel('hdiff in m')\\n    plt.ylabel('o:N, ^:H, h:F, *:M, +:C ') \\n    \\n    plt.subplot(3,1,3)\\n    plt.scatter(hdiff, SS['POD'],  marker='o', s=40)\\n    plt.scatter(hdiff, SS['FAR'],  marker='^', s=40)\\n    plt.scatter(hdiff, SS['BID'],  marker='h', s=40)\\n    plt.scatter(hdiff, SS['HSS'],  marker='*', s=40)\\n    plt.scatter(hdiff, SS['PE']/100,  marker='+', s=40)\\n    plt.scatter(hdiff, SS['HR'],  marker='d', s=40)\\n    plt.grid()    \\n    plt.xlabel('hdiff in m')\\n    plt.ylabel('o:POD, ^:FAR, h:BID, *:HSS, +:PR , d:HR') \\n    \\n    plt.subplot(3,1,1)\\n    plt.scatter(hdiff, SS['bias'], c=SS['corr'], vmin=0.3, vmax=0.7, marker='o', s=50)\\n    plt.scatter(hdiff, SS['RMSE'], c=SS['corr'], vmin=0.3, vmax=0.7, marker='^', s=50)\\n    plt.axhline(0)\\n    plt.grid()    \\n    plt.xlabel('hdiff in m')\\n    plt.ylabel('bias o ; RMSE ^ ') \\n\\nplt.colorbar()\\nplt.tight_layout()\\nplt.show()\\n\""
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "step = 125.\n",
    "corr, bias, RMSE = np.array([]), np.array([]), np.array([])\n",
    "N, H, F, M, C = np.array([]), np.array([]), np.array([]), np.array([]), np.array([])\n",
    "POD, HSS, HR, BID, PE, FAR = np.array([]), np.array([]), np.array([]), np.array([]), np.array([]), np.array([])\n",
    "\n",
    "\n",
    "_mean, _std = np.array([]), np.array([])\n",
    "\n",
    "for hdiff in np.arange(step,3000,step):\n",
    "    print hdiff\n",
    "    SS = skill_score(gr[(abs((hmin - (cfb - ele))) < hdiff) & (abs((hmin - (cfb - ele))) > hdiff-step)],\n",
    "                     sr[(abs((hmin - (cfb - ele))) < hdiff) & (abs((hmin - (cfb - ele))) > hdiff-step)],\n",
    "                     0.5) \n",
    "    \n",
    "    corr = np.append(corr, SS['corr'])\n",
    "    bias = np.append(bias, SS['bias'])\n",
    "    RMSE = np.append(RMSE, SS['RMSE'])\n",
    "    POD= np.append(POD, SS['POD'])\n",
    "    HSS = np.append(HSS, SS['HSS'])\n",
    "    HR = np.append(HR, SS['HR'])\n",
    "    BID = np.append(BID, SS['BID'])\n",
    "    PE = np.append(PE, SS['PE'])\n",
    "    FAR = np.append(FAR, SS['FAR'])\n",
    "    N = np.append(N, SS['N'])\n",
    "    H = np.append(H, SS['H'])\n",
    "    F = np.append(F, SS['F'])\n",
    "    M = np.append(M, SS['M'])\n",
    "    C = np.append(C, SS['C'])\n",
    "\n",
    "    \n",
    "\"\"\"    plt.subplot(3,1,2)\n",
    "    plt.scatter(hdiff, SS['N'],  marker='o', s=40)\n",
    "    plt.scatter(hdiff, SS['H'],  marker='^', s=40)\n",
    "    plt.scatter(hdiff, SS['F'],  marker='h', s=40)\n",
    "    plt.scatter(hdiff, SS['M'],  marker='*', s=40)\n",
    "    plt.scatter(hdiff, SS['C'],  marker='+', s=40)\n",
    "    \n",
    "    plt.legend()\n",
    "    plt.grid()    \n",
    "    plt.xlabel('hdiff in m')\n",
    "    plt.ylabel('o:N, ^:H, h:F, *:M, +:C ') \n",
    "    \n",
    "    plt.subplot(3,1,3)\n",
    "    plt.scatter(hdiff, SS['POD'],  marker='o', s=40)\n",
    "    plt.scatter(hdiff, SS['FAR'],  marker='^', s=40)\n",
    "    plt.scatter(hdiff, SS['BID'],  marker='h', s=40)\n",
    "    plt.scatter(hdiff, SS['HSS'],  marker='*', s=40)\n",
    "    plt.scatter(hdiff, SS['PE']/100,  marker='+', s=40)\n",
    "    plt.scatter(hdiff, SS['HR'],  marker='d', s=40)\n",
    "    plt.grid()    \n",
    "    plt.xlabel('hdiff in m')\n",
    "    plt.ylabel('o:POD, ^:FAR, h:BID, *:HSS, +:PR , d:HR') \n",
    "    \n",
    "    plt.subplot(3,1,1)\n",
    "    plt.scatter(hdiff, SS['bias'], c=SS['corr'], vmin=0.3, vmax=0.7, marker='o', s=50)\n",
    "    plt.scatter(hdiff, SS['RMSE'], c=SS['corr'], vmin=0.3, vmax=0.7, marker='^', s=50)\n",
    "    plt.axhline(0)\n",
    "    plt.grid()    \n",
    "    plt.xlabel('hdiff in m')\n",
    "    plt.ylabel('bias o ; RMSE ^ ') \n",
    "\n",
    "plt.colorbar()\n",
    "plt.tight_layout()\n",
    "plt.show()\n",
    "\"\"\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "293906\n",
      "293393.0\n"
     ]
    }
   ],
   "source": [
    "print np.count_nonzero(~np.isnan(sr- gr))\n",
    "print np.nansum(N)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "hdiff = np.arange(step,3000,step)\n",
    "ff=20\n",
    "plt.figure(figsize=(20,14))\n",
    "plt.subplot(3,1,2,aspect='auto')\n",
    "plt.plot(hdiff, N,  marker='o', label='N')\n",
    "plt.plot(hdiff, H,  marker='^', label='H')\n",
    "plt.plot(hdiff, F,  marker='h', label='F')\n",
    "plt.plot(hdiff, M,  marker='*', label='M')\n",
    "plt.plot(hdiff, C,  marker='+', label='C')\n",
    "plt.xlim(0, 4000)\n",
    "plt.legend(loc='upper right',fontsize=ff)\n",
    "plt.grid()    \n",
    "plt.xticks(fontsize=ff)\n",
    "plt.yticks(fontsize=ff)\n",
    "\n",
    "# -------------------------------\n",
    "plt.subplot(3,1,3,aspect='auto')\n",
    "plt.plot(hdiff, corr,  marker='s', label='Corr', lw=3)\n",
    "plt.plot(hdiff, POD,  marker='o', label='POD')\n",
    "plt.plot(hdiff, FAR,  marker='^', label='FAR')\n",
    "plt.plot(hdiff, BID,  marker='h', label='BID')\n",
    "plt.plot(hdiff, HSS,  marker='*', label='HSS')\n",
    "plt.plot(hdiff, PE/100,  marker='D', label='PE')\n",
    "plt.plot(hdiff, HR,  marker='d', label='HR')\n",
    "plt.axhline(0, color='grey', lw=4, alpha = 0.5)\n",
    "\n",
    "plt.legend(loc='upper right',fontsize=ff)\n",
    "plt.grid()    \n",
    "plt.xlabel('hdiff in m',fontsize=ff)\n",
    "plt.xlim(0, 4000)\n",
    "plt.xticks(fontsize=ff)\n",
    "plt.yticks(fontsize=ff)\n",
    "\n",
    "# ------------------------------\n",
    "plt.subplot(3,1,1,aspect='auto')\n",
    "plt.plot(hdiff, bias, marker='o', color='blue', label='bias')\n",
    "plt.plot(hdiff, RMSE, marker='^', color='black', label='RMSE')\n",
    "plt.axhline(0, color='grey', lw=4)\n",
    "plt.grid()    \n",
    "plt.legend(loc='upper right',fontsize=ff)\n",
    "plt.xticks(fontsize=ff)\n",
    "plt.yticks(fontsize=ff)\n",
    "plt.title('DPR NS vs. RY -  Liquid Stratiform 2014-2018, TH=0.5 mm/h, N_all:' + str(np.nansum(N)), fontsize=ff)\n",
    "\n",
    "\n",
    "plt.xlim(0, 4000)\n",
    "\n",
    "#plt.colorbar()\n",
    "#plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%%time\n",
    "\n",
    "#pl.figure(figsize=(10,10))\n",
    "#pl.scatter(para1, para2, c=para3, edgecolors='none', s=10)\n",
    "#pl.xlim(-1,60)\n",
    "#pl.ylim(-1,60)\n",
    "#pl.colorbar()\n",
    "#pl.plot([0,60],[0,60], linestyle=\"solid\", color=\"black\")\n",
    "#pl.grid()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "dpr = sr_pns.copy() # para2 = GPM \n",
    "dpr2 = sr_pav.copy() # para2 = GPM \n",
    "ry = gr_pns.copy() # para1 = RY\n",
    "xx = xl.copy()\n",
    "yy = yl.copy()\n",
    "#pha = np.array(sr_pha.copy(),dtype=float)\n",
    "typ1 = typ.copy()\n",
    "#dpr_h = ele_offest - sr_ele\n",
    "#dpr_h = (178. - sr_cfb.copy())*125\n",
    "#ry_h = gr_hmin.copy()\n",
    "\n",
    "ff=20\n",
    "#yy[(ry == 0) & (dpr > 0)] = np.nan\n",
    "#xx[(ry == 0) & (dpr > 0)] = np.nan\n",
    "#\"\"\"\n",
    "th = 0.\n",
    "\n",
    "#ry[ry <= th] = np.nan\n",
    "#dpr[dpr <= th] = np.nan\n",
    "#dpr2[dpr2 <= th] = np.nan\n",
    "\n",
    "#typ[(ry <= th) & (dpr <= th) & (dpr2 <= th) ] = np.nan\n",
    "#pha[(ry <= th) & (dpr <= th) & (dpr2 <= th) ] = np.nan\n",
    "\n",
    "#mask = ~np.isnan(ry) & ~np.isnan(dpr) & ~np.isnan(dpr2) & ~np.isnan(typ) & ~np.isnan(pha)\n",
    "\n",
    "\n",
    "#\"\"\"\n",
    "#dpr_h[(ry <= 0) & (dpr <= 0) & (dpr2 <= 0)] = np.nan\n",
    "#ry_h[(ry <= 0) & (dpr <= 0) & (dpr2 <= 0)] = np.nan\n",
    "\n",
    "#ry[abs(ry_h-dpr_h)<500] = np.nan\n",
    "#dpr[abs(ry_h-dpr_h)<500] = np.nan\n",
    "\n",
    "#sth[(ry <= 0) & (dpr <= 0) & (dpr2 <= 0)] = np.nan\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "np.unique(typ)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "pl.hist(typ[typ!=-1], bins=3,histtype='step', align='mid')\n",
    "pl.xlim(0,4)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%%time\n",
    "ff=20\n",
    "pl.figure(figsize=(10,6))\n",
    "pl.xscale('log')\n",
    "pl.yscale('log')\n",
    "pl.xlim(1e-3,1000)\n",
    "pl.ylim(1e-3,1000)\n",
    "#pl.scatter(ry[mask],dpr[mask], c=sth[mask], edgecolors='none')\n",
    "pl.scatter(ry[mask],dpr2[mask], c=typ[mask], edgecolors='none', cmap='jet')\n",
    "print(corcorcor(ry[mask],dpr2[mask]))\n",
    "\n",
    "cb = pl.colorbar()\n",
    "cb.set_label('GR - SR height in m', fontsize=ff)\n",
    "cb.ax.tick_params(labelsize=20)\n",
    "\n",
    "pl.grid()\n",
    "\n",
    "pl.ylabel('DPR Rainrate in mm/h',fontsize=ff)\n",
    "pl.xlabel('Radolan Rainrate in mm/h',fontsize=ff)\n",
    "pl.xticks(fontsize=ff)\n",
    "pl.yticks(fontsize=ff)\n",
    "pl.plot([0,1000], [0,1000], color='black')\n",
    "#pl.savefig(\"/automount/ftp/velibor/test1.png\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "typ1= typ.astype(int)/10000000\n",
    "print('___Corr______Err______BIAS____')\n",
    "\n",
    "print(corcorcor(ry[pha==0.],dpr2[pha==0.]))\n",
    "print(corcorcor(ry[pha==2.],dpr2[pha==2.]))\n",
    "\n",
    "print('___________________')\n",
    "print(corcorcor(ry[typ1==1.],dpr2[typ1==1.]))\n",
    "print(corcorcor(ry[typ1==2.],dpr2[typ1==2.]))\n",
    "print(corcorcor(ry[typ1==3.],dpr2[typ1==3.]))\n",
    "\n",
    "print('___Min__________Max______')\n",
    "\n",
    "print(np.nanmin(ry[ry>0]), np.nanmax(ry))\n",
    "print(np.nanmin(dpr[dpr>0]), np.nanmax(dpr))\n",
    "print(np.nanmin(dpr2[dpr2>0]), np.nanmax(dpr2))\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "np.unique(sr_pha) # Keine mixed phase ? "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%%time\n",
    "\n",
    "ff=20\n",
    "pl.figure(figsize=(15,15))\n",
    "pl.subplot(3,2,1, aspect='equal')\n",
    "pl.xscale('log')\n",
    "pl.yscale('log')\n",
    "pl.xlim(1e-3,1000)\n",
    "pl.ylim(1e-3,1000)\n",
    "#pl.scatter(ry[mask],dpr[mask], c=sth[mask], edgecolors='none')\n",
    "pl.scatter(ry[(pha==0.) & (typ1==1.)],dpr[(pha==0.) & (typ1==1.)],  edgecolors='none', color='green', \n",
    "          label=str(corcorcor(ry[(pha==0.) & (typ1==1.)],dpr[(pha==0.) & (typ1==1.)])))\n",
    "pl.plot([0,1000], [0,1000], color='black')\n",
    "pl.ylabel('DPR RR in mm/h \\n Stratiform',fontsize=ff)\n",
    "pl.title('Solid',fontsize=ff, color='green')\n",
    "pl.xticks(fontsize=ff)\n",
    "pl.yticks(fontsize=ff)\n",
    "pl.legend()\n",
    "pl.plot([0,1000], [0,1000], color='black')\n",
    "\n",
    "pl.subplot(3,2,2, aspect='equal')\n",
    "pl.xscale('log')\n",
    "pl.yscale('log')\n",
    "pl.xlim(1e-3,1000)\n",
    "pl.ylim(1e-3,1000)\n",
    "#pl.scatter(ry[mask],dpr[mask], c=sth[mask], edgecolors='none')\n",
    "pl.scatter(ry[(pha==2.) & (typ1==1.)],dpr[(pha==2.) & (typ1==1.)],  edgecolors='none',\n",
    "          label=str(corcorcor(ry[(pha==2.) & (typ1==1.)],dpr[(pha==2.) & (typ1==1.)])))\n",
    "pl.plot([0,1000], [0,1000], color='black')\n",
    "pl.xticks(fontsize=ff)\n",
    "pl.yticks(fontsize=ff)\n",
    "pl.title('Liquid',fontsize=ff, color='blue')\n",
    "pl.legend()\n",
    "pl.plot([0,1000], [0,1000], color='black')\n",
    "\n",
    "pl.subplot(3,2,3, aspect='equal')\n",
    "pl.xscale('log')\n",
    "pl.yscale('log')\n",
    "pl.xlim(1e-3,1000)\n",
    "pl.ylim(1e-3,1000)\n",
    "#pl.scatter(ry[mask],dpr[mask], c=sth[mask], edgecolors='none')\n",
    "pl.scatter(ry[(pha==0.) & (typ1==2.)],dpr[(pha==0.) & (typ1==2.)],  edgecolors='none', color='green',\n",
    "          label=str(corcorcor(ry[(pha==0.) & (typ1==2.)],dpr[(pha==0.) & (typ1==2.)])))\n",
    "pl.plot([0,1000], [0,1000], color='black')\n",
    "pl.ylabel('DPR RR in mm/h \\n Convective',fontsize=ff)\n",
    "pl.xticks(fontsize=ff)\n",
    "pl.yticks(fontsize=ff)\n",
    "pl.legend()\n",
    "pl.plot([0,1000], [0,1000], color='black')\n",
    "\n",
    "pl.subplot(3,2,4, aspect='equal')\n",
    "pl.xscale('log')\n",
    "pl.yscale('log')\n",
    "pl.xlim(1e-3,1000)\n",
    "pl.ylim(1e-3,1000)\n",
    "#pl.scatter(ry[mask],dpr[mask], c=sth[mask], edgecolors='none')\n",
    "pl.scatter(ry[(pha==2.) & (typ1==2.)],dpr[(pha==2.) & (typ1==2.)],  edgecolors='none',\n",
    "          label=str(corcorcor(ry[(pha==2.) & (typ1==2.)],dpr[(pha==2.) & (typ1==2.)])))\n",
    "pl.plot([0,1000], [0,1000], color='black')\n",
    "pl.xticks(fontsize=ff)\n",
    "pl.yticks(fontsize=ff)\n",
    "pl.legend()\n",
    "pl.plot([0,1000], [0,1000], color='black')\n",
    "\n",
    "pl.subplot(3,2,5, aspect='equal')\n",
    "pl.xscale('log')\n",
    "pl.yscale('log')\n",
    "pl.xlim(1e-3,1000)\n",
    "pl.ylim(1e-3,1000)\n",
    "#pl.scatter(ry[mask],dpr[mask], c=sth[mask], edgecolors='none')\n",
    "pl.scatter(ry[(pha==0.) & (typ1==3.)],dpr[(pha==0.) & (typ1==3.)],  edgecolors='none', color='green',\n",
    "          label=str(corcorcor(ry[(pha==0.) & (typ1==3.)],dpr[(pha==0.) & (typ1==3.)])))\n",
    "pl.plot([0,1000], [0,1000], color='black')\n",
    "pl.xlabel('Radolan RR in mm/h',fontsize=ff)\n",
    "pl.ylabel('DPR RR in mm/h \\n Other',fontsize=ff)\n",
    "pl.xticks(fontsize=ff)\n",
    "pl.yticks(fontsize=ff)\n",
    "pl.legend()\n",
    "pl.plot([0,1000], [0,1000], color='black')\n",
    "\n",
    "pl.subplot(3,2,6, aspect='equal')\n",
    "pl.xscale('log')\n",
    "pl.yscale('log')\n",
    "pl.xlim(1e-3,1000)\n",
    "pl.ylim(1e-3,1000)\n",
    "#pl.scatter(ry[mask],dpr[mask], c=sth[mask], edgecolors='none')\n",
    "pl.scatter(ry[(pha==2.) & (typ1==3.)],dpr[(pha==2.) & (typ1==3.)],  edgecolors='none', \n",
    "          label=str(corcorcor(ry[(pha==2.) & (typ1==3.)],dpr[(pha==2.) & (typ1==3.)])))\n",
    "pl.plot([0,1000], [0,1000], color='black')\n",
    "pl.xlabel('Radolan RR in mm/h',fontsize=ff)\n",
    "pl.xticks(fontsize=ff)\n",
    "pl.yticks(fontsize=ff)\n",
    "pl.legend()\n",
    "pl.plot([0,1000], [0,1000], color='black')\n",
    "\n",
    "\n",
    "pl.tight_layout()\n",
    "#pl.savefig(\"/automount/ftp/velibor/test1.png\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%%time\n",
    "\n",
    "srh = sr_cfb-sr_ele\n",
    "grh = gr_hmin \n",
    "\n",
    "hdiff = abs(grh - srh)\n",
    "\n",
    "ff=20\n",
    "pl.figure(figsize=(20,20))\n",
    "ax1 = pl.subplot(3,2,1, aspect='equal')\n",
    "#pl.xscale('log')\n",
    "#pl.yscale('log')\n",
    "#pl.xlim(1e-3,1000)\n",
    "#pl.ylim(1e-3,1000)\n",
    "plotplot(ax1, ry[(pha==0.) & (typ1==1.) & (hdiff<=100)],dpr[(pha==0.) & (typ1==1.)& (hdiff<=100)], 'Solid, Stratiform')\n",
    "#pl.plot([0,1000], [0,1000], color='black')\n",
    "#pl.ylabel('DPR RR in mm/h \\n Stratiform',fontsize=ff)\n",
    "#pl.title('Solid',fontsize=ff, color='green')\n",
    "#pl.xticks(fontsize=ff)\n",
    "#pl.yticks(fontsize=ff)\n",
    "#pl.legend()\n",
    "#pl.plot([0,1000], [0,1000], color='black')\n",
    "\n",
    "ax2 = pl.subplot(3,2,2, aspect='equal')\n",
    "#pl.xscale('log')\n",
    "#pl.yscale('log')\n",
    "#pl.xlim(1e-3,1000)\n",
    "#pl.ylim(1e-3,1000)\n",
    "#pl.scatter(ry[mask],dpr[mask], c=sth[mask], edgecolors='none')\n",
    "plotplot(ax2,ry[(pha==2.) & (typ1==1.)& (hdiff<=100)],dpr[(pha==2.) & (typ1==1.)& (hdiff<=100)],'Liquid, Stratiform')\n",
    "#pl.plot([0,1000], [0,1000], color='black')\n",
    "#pl.xticks(fontsize=ff)\n",
    "#pl.yticks(fontsize=ff)\n",
    "#pl.title('Liquid',fontsize=ff, color='blue')\n",
    "#pl.legend()\n",
    "#pl.plot([0,1000], [0,1000], color='black')\n",
    "\n",
    "ax3 = pl.subplot(3,2,3, aspect='equal')\n",
    "#pl.xscale('log')\n",
    "#pl.yscale('log')\n",
    "#pl.xlim(1e-3,1000)\n",
    "#pl.ylim(1e-3,1000)\n",
    "#pl.scatter(ry[mask],dpr[mask], c=sth[mask], edgecolors='none')\n",
    "plotplot(ax3, ry[(pha==0.) & (typ1==2.)& (hdiff<=100)],dpr[(pha==0.) & (typ1==2.)& (hdiff<=100)], 'Solid, Convective')\n",
    "#pl.plot([0,1000], [0,1000], color='black')\n",
    "#pl.ylabel('DPR RR in mm/h \\n Convective',fontsize=ff)\n",
    "#pl.xticks(fontsize=ff)\n",
    "#pl.yticks(fontsize=ff)\n",
    "#pl.legend()\n",
    "#pl.plot([0,1000], [0,1000], color='black')\n",
    "\n",
    "ax4 = pl.subplot(3,2,4, aspect='equal')\n",
    "#pl.xscale('log')\n",
    "#pl.yscale('log')\n",
    "#pl.xlim(1e-3,1000)\n",
    "#pl.ylim(1e-3,1000)\n",
    "#pl.scatter(ry[mask],dpr[mask], c=sth[mask], edgecolors='none')\n",
    "plotplot(ax4,ry[(pha==2.) & (typ1==2.)& (hdiff<=100)],dpr[(pha==2.) & (typ1==2.)& (hdiff<=100)], 'Liquid Convective')\n",
    "#pl.plot([0,1000], [0,1000], color='black')\n",
    "#pl.xticks(fontsize=ff)\n",
    "#pl.yticks(fontsize=ff)\n",
    "#pl.legend()\n",
    "#pl.plot([0,1000], [0,1000], color='black')\n",
    "\n",
    "ax5 = pl.subplot(3,2,5, aspect='equal')\n",
    "\n",
    "#pl.scatter(ry[mask],dpr[mask], c=sth[mask], edgecolors='none')\n",
    "plotplot(ax5, ry[(pha==0.) & (typ1==3.)& (hdiff<=100)],dpr[(pha==0.) & (typ1==3.)& (hdiff<=100)], 'Solid, Other')\n",
    "\n",
    "\n",
    "ax6 = pl.subplot(3,2,6, aspect='equal')\n",
    "\n",
    "#pl.scatter(ry[mask],dpr[mask], c=sth[mask], edgecolors='none')\n",
    "plotplot(ax6,ry[(pha==2.) & (typ1==3.)& (hdiff<=100)],dpr[(pha==2.) & (typ1==3.)& (hdiff<=100)],  'Liquid, Other')\n",
    "\n",
    "\n",
    "\n",
    "pl.tight_layout()\n",
    "pl.savefig(\"/automount/ftp/velibor/with_hdiff.png\")\n",
    "pl.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "pl.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "pl.figure(figsize=(20,20))\n",
    "ax1 = pl.subplot(1,1,1, aspect='equal')\n",
    "plotplot(ax1,ry[(pha==0.) & (typ1==1.)],dpr[(pha==0.) & (typ1==1.)], \n",
    "          str(corcorcor(ry[(pha==0.) & (typ1==1.)],dpr[(pha==0.) & (typ1==1.)])))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Test GR und SR Messhoehen in abhängigkeit zum BB\n",
    "--------------------------------------------"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "dpr1 = sr_pns.copy() # para2 = GPM \n",
    "ry1 = gr_pns.copy() # para1 = RY\n",
    "\n",
    "rhmax = gr_hmax.copy()\n",
    "rhmin1 = gr_hmin.copy()\n",
    "dprh1 = sr_cfb.copy()\n",
    "\n",
    "bbh1 = sr_bbh.copy()\n",
    "bbw1 = sr_bbw.copy()\n",
    "\n",
    "bbh1[bbh1<=0]=np.nan\n",
    "#bbh[bbh==0]=np.nan\n",
    "dprh1[dprh1 < 0] = np.nan"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "dprh1 = 125.*(176. - dprh1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "dprh1[typ!=1]= np.nan # 1 =strati ; 2 =conv "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "dprh1[sr_pha!=2]= np.nan # 0=solid, 2=liquid"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "np.unique(typ), np.unique(sr_pha) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# for hist2d only values without nan\n",
    "bbh = bbh1[~np.isnan(bbh1) & ~np.isnan(bbw1) & ~np.isnan(rhmin1) & ~np.isnan(dprh1) & ~np.isnan(ry1) & ~np.isnan(dpr1)]\n",
    "bbw = bbw1[~np.isnan(bbh1) & ~np.isnan(bbw1) & ~np.isnan(rhmin1) & ~np.isnan(dprh1) & ~np.isnan(ry1) & ~np.isnan(dpr1)]\n",
    "dprh = dprh1[~np.isnan(bbh1) & ~np.isnan(bbw1) & ~np.isnan(rhmin1) & ~np.isnan(dprh1) & ~np.isnan(ry1) & ~np.isnan(dpr1)]\n",
    "rhmin = rhmin1[~np.isnan(bbh1) & ~np.isnan(bbw1) & ~np.isnan(rhmin1) & ~np.isnan(dprh1) & ~np.isnan(ry1) & ~np.isnan(dpr1)]\n",
    "ry = ry1[~np.isnan(bbh1) & ~np.isnan(bbw1) & ~np.isnan(rhmin1) & ~np.isnan(dprh1) & ~np.isnan(ry1) & ~np.isnan(dpr1)]\n",
    "dpr = dpr1[~np.isnan(bbh1) & ~np.isnan(bbw1) & ~np.isnan(rhmin1) & ~np.isnan(dprh1) & ~np.isnan(ry1) & ~np.isnan(dpr1)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Wie oft misst RADOLAN über unter oder im BB\n",
    "-------------------------"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "pl.figure(figsize=(20,20))\n",
    "ax1 = pl.subplot(3,3,1, aspect='equal')\n",
    "plotplot(ax1,ry[(dprh<bbh+bbw/2)&(dprh>bbh-bbw/2)&(rhmin<bbh+bbw/2)&(rhmin>bbh-bbw/2)],\n",
    "         dpr[(dprh<bbh+bbw/2)&(dprh>bbh-bbw/2)&(rhmin<bbh+bbw/2)&(rhmin>bbh-bbw/2)],\n",
    "         'SR and GR in BB')\n",
    "\n",
    "ax2 = pl.subplot(3,3,2, aspect='equal')\n",
    "plotplot(ax2, ry[(dprh<bbh+bbw/2)&(dprh>bbh-bbw/2)&(rhmin>bbh+bbw/2)],\n",
    "           dpr[(dprh<bbh+bbw/2)&(dprh>bbh-bbw/2)&(rhmin>bbh+bbw/2)], 'GR over & SR in BB')\n",
    "\n",
    "ax3 = pl.subplot(3,3,8, aspect='equal')\n",
    "plotplot(ax3, ry[(rhmin<bbh+bbw/2)&(rhmin>bbh-bbw/2)&(dprh>bbh+bbw/2)],\n",
    "           dpr[(rhmin<bbh+bbw/2)&(rhmin>bbh-bbw/2)&(dprh>bbh+bbw/2)], 'GR in & SR over BB')\n",
    "\n",
    "ax4 = pl.subplot(3,3,4, aspect='equal')\n",
    "plotplot(ax4,ry[(dprh<bbh+bbw/2)&(dprh>bbh-bbw/2)&(rhmin<bbh-bbw/2)],\n",
    "           dpr[(dprh<bbh+bbw/2)&(dprh>bbh-bbw/2)&(rhmin<bbh-bbw/2)], 'GR under & SR in BB')\n",
    "\n",
    "ax5 = pl.subplot(3,3,5, aspect='equal')\n",
    "plotplot(ax5, ry[(rhmin<bbh+bbw/2)&(rhmin>bbh-bbw/2)&(dprh<bbh-bbw/2)],\n",
    "           dpr[(rhmin<bbh+bbw/2)&(rhmin>bbh-bbw/2)&(dprh<bbh-bbw/2)], 'GR in & SR under BB')\n",
    "\n",
    "ax6 = pl.subplot(3,3,6, aspect='equal')\n",
    "plotplot(ax6, ry[(rhmin>bbh+bbw/2)&(dprh<bbh-bbw/2)],\n",
    "           dpr[(rhmin>bbh+bbw/2)&(dprh<bbh-bbw/2)], 'GR over & SR under BB')\n",
    "\n",
    "ax7 = pl.subplot(3,3,7, aspect='equal')\n",
    "plotplot(ax7, ry[(rhmin<bbh-bbw/2)&(dprh>bbh+bbw/2)],\n",
    "               dpr[(rhmin<bbh-bbw/2)&(dprh>bbh+bbw/2)], 'GR under & SR over BB' )\n",
    "\n",
    "ax8 = pl.subplot(3,3,3, aspect='equal')\n",
    "plotplot(ax8, ry[(rhmin<bbh-bbw/2)&(dprh<bbh-bbw/2)],\n",
    "           dpr[(rhmin<bbh-bbw/2)&(dprh<bbh-bbw/2)],'GR under & SR under BB')\n",
    "\n",
    "ax9 = pl.subplot(3,3,9, aspect='equal')\n",
    "plotplot(ax9, ry[(rhmin>bbh+bbw/2)&(dprh>bbh+bbw/2)],\n",
    "           dpr[(rhmin>bbh+bbw/2)&(dprh>bbh+bbw/2)], 'GR over & SR over BB')\n",
    "\n",
    "pl.tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#pl.yscale('log')\n",
    "ff=20\n",
    "pl.figure(figsize=(8,6))\n",
    "pl.hist(rhmin[(rhmin<bbh+bbw/2)&(rhmin>bbh-bbw/2)], bins=30, histtype='step', label='in BB')\n",
    "pl.hist(rhmin[(rhmin>bbh+bbw/2)], bins=30, histtype='step', label='over BB')\n",
    "pl.hist(rhmin[(rhmin<bbh-bbw/2)], bins=30, histtype='step', label='under BB')\n",
    "pl.legend(loc='upper right', fontsize=20)\n",
    "pl.xlabel('RADOLAN measured height in m', fontsize=20)\n",
    "pl.ylabel('number', fontsize=20)\n",
    "pl.xticks(fontsize=20)\n",
    "pl.yticks(fontsize=20)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#pl.yscale('log')\n",
    "ff=20\n",
    "pl.figure(figsize=(8,6))\n",
    "pl.hist(dprh[(dprh<bbh+bbw/2)&(dprh>bbh-bbw/2)], bins=21, histtype='step', label='in BB')\n",
    "pl.hist(dprh[(dprh>bbh+bbw/2)], bins=21, histtype='step', label='over BB')\n",
    "pl.hist(dprh[(dprh<bbh-bbw/2)], bins=21, histtype='step', label='under BB')\n",
    "pl.legend(loc='upper right', fontsize=20)\n",
    "pl.xlabel('DPR measured height in m', fontsize=20)\n",
    "pl.ylabel('number', fontsize=20)\n",
    "pl.xticks(fontsize=20)\n",
    "pl.yticks(fontsize=20)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "ff=20"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def plotplot(ax, A, B, tit):\n",
    "    ff=20\n",
    "    from matplotlib.colors import LogNorm\n",
    "\n",
    "\n",
    "\n",
    "    #pl.hist2d(ry1[(rhmin1>bbh1+bbw1/2)&(dprh1>bbh1+bbw1/2)],\n",
    "    #           dpr1[(rhmin1>bbh1+bbw1/2)&(dprh1>bbh1+bbw1/2)],\n",
    "    #          bins=100, vmin=1, vmax=100, norm=LogNorm())\n",
    "\n",
    "    #fig, ax = plt.subplots(figsize=(10,10))\n",
    "    #ax = pl.subplot()\n",
    "    xbins = 10**np.arange(-3, 3, 0.1)\n",
    "    ybins = 10**np.arange(-3, 3, 0.1)\n",
    "\n",
    "    counts, _, _ = np.histogram2d(A,\n",
    "                                  B,\n",
    "                                  bins=(xbins, ybins))\n",
    "\n",
    "    a = np.ma.masked_where(counts.T == 0., counts.T)\n",
    "    pl.pcolormesh(xbins, ybins, a, cmap='jet', vmin=0)\n",
    "    print str(corcorcor(A,B))\n",
    "    \n",
    "    cb = plt.colorbar(shrink=0.5)\n",
    "    cb.set_label('number of samples', fontsize=ff)\n",
    "    cb.ax.tick_params(labelsize=20)\n",
    "\n",
    "    ax.set_xscale('log')\n",
    "    ax.set_yscale('log')\n",
    "\n",
    "    pl.title(tit,fontsize=ff)\n",
    "    pl.ylabel('DPR Rainrate in mm/h',fontsize=ff)\n",
    "    pl.xlabel('Radolan Rainrate in mm/h',fontsize=ff)\n",
    "    \n",
    "    pl.xticks(fontsize=ff)\n",
    "    pl.yticks(fontsize=ff)\n",
    "    pl.plot([0,1000], [0,1000], color='black', label='r: '+ str(corcorcor(A,B)[0])+'\\n bias: '+ str(corcorcor(A,B)[1]))\n",
    "    pl.legend(loc='lower right', markerscale=0)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "np.unique(dpr1[(rhmin1>bbh1+bbw1/2)&(dprh1>bbh1+bbw1/2)])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "len(np.array([])) +  len(np.array([0]))"
   ]
  },
  {
   "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": false
   },
   "outputs": [],
   "source": [
    "pl.xscale('log')\n",
    "pl.yscale('log')\n",
    "\n",
    "pl.xlim(th,1000)\n",
    "pl.ylim(th,1000)\n",
    "pl.scatter(ry[typ1==1.],dpr[typ1==1.])\n",
    "pl.plot([0,1000], [0,1000], color='black')\n",
    "corcorcor(ry[typ1==1.],dpr[typ1==1.])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "pl.xscale('log')\n",
    "pl.yscale('log')\n",
    "pl.xlim(th,1000)\n",
    "pl.ylim(th,1000)\n",
    "pl.scatter(ry[typ1==2.],dpr[typ1==2.])\n",
    "pl.plot([0,1000], [0,1000], color='black')\n",
    "corcorcor(ry[typ1==2.],dpr[typ1==2.])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "pl.xscale('log')\n",
    "pl.yscale('log')\n",
    "pl.xlim(th,1000)\n",
    "pl.ylim(th,1000)\n",
    "pl.scatter(ry[typ1==3.],dpr[typ1==3.])\n",
    "pl.plot([0,1000], [0,1000], color='black')\n",
    "corcorcor(ry[typ1==3.],dpr[typ1==3.])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "pl.xscale('log')\n",
    "pl.yscale('log')\n",
    "pl.xlim(th,1000)\n",
    "pl.ylim(th,1000)\n",
    "pl.scatter(ry[pha==2.],dpr[pha==2.])\n",
    "pl.plot([0,1000], [0,1000], color='black')\n",
    "corcorcor(ry[pha==2.],dpr[pha==2.])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "pl.xscale('log')\n",
    "pl.yscale('log')\n",
    "pl.xlim(th,1000)\n",
    "pl.ylim(th,1000)\n",
    "pl.scatter(ry[pha==0.],dpr[pha==0.])\n",
    "pl.plot([0,1000], [0,1000], color='black')\n",
    "corcorcor(ry[pha==0.],dpr[pha==0.])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "np.unique(typ1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "mask = ~np.isnan(A) & ~np.isnan(B)\n",
    "diffi = A[mask]-B[mask]\n",
    "bias = np.nansum(diffi)/len(diffi)\n",
    "slope, intercept, r_value, p_value, std_err = stats.linregress(A[mask], B[mask])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%%time\n",
    "ff=20\n",
    "pl.figure(figsize=(10,6))\n",
    "pl.xscale('log')\n",
    "pl.yscale('log')\n",
    "pl.xlim(1e-3,1000)\n",
    "pl.ylim(1e-3,1000)\n",
    "#pl.scatter(ry[mask],dpr[mask], c=sth[mask], edgecolors='none')\n",
    "pl.scatter(ry[mask],dpr[mask], c=ry_h[mask]-dpr_h[mask], edgecolors='none', cmap='jet', vmin=-2000, vmax=2000)\n",
    "print(corcorcor(ry[mask],dpr[mask]))\n",
    "cb = pl.colorbar()\n",
    "cb.set_label('storm height in m', fontsize=ff)\n",
    "cb.ax.tick_params(labelsize=20)\n",
    "\n",
    "pl.grid()\n",
    "\n",
    "pl.ylabel('DPR Rainrate in mm/h',fontsize=ff)\n",
    "pl.xlabel('Radolan Rainrate in mm/h',fontsize=ff)\n",
    "pl.xticks(fontsize=ff)\n",
    "pl.yticks(fontsize=ff)\n",
    "pl.plot([0,1000], [0,1000], color='black')\n",
    "#pl.savefig(\"/automount/ftp/velibor/test2.png\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "pl.figure(figsize=(10,6))\n",
    "pl.hist(ry[mask],bins=np.logspace(np.log10(0.1), np.log10(100),1250),\n",
    "        normed=True, cumulative=True, histtype='step', label='RY')\n",
    "pl.hist(dpr[mask],bins=np.logspace(np.log10(0.1), np.log10(100),1250),\n",
    "        normed=True, cumulative=True, histtype='step', label='PNS')\n",
    "pl.hist(dpr2[mask],bins=np.logspace(np.log10(0.1), np.log10(100),1250),\n",
    "        normed=True, cumulative=True, histtype='step', label='PAV')\n",
    "pl.legend(fontsize=20, loc='lower right')\n",
    "pl.xscale('log', nonposy='clip')\n",
    "pl.grid()\n",
    "pl.xlabel('RR (mm/h)', fontsize=20)\n",
    "pl.ylabel('cdf', fontsize=20)\n",
    "pl.xticks(fontsize=20);pl.yticks(fontsize=20)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "pl.figure(figsize=(10,6))\n",
    "pl.hist(ry[mask],bins=np.logspace(np.log10(0.001), np.log10(100),250),\n",
    "        normed=True, histtype='step', label='RY')\n",
    "pl.hist(dpr[mask],bins=np.logspace(np.log10(0.001), np.log10(100),250),\n",
    "        normed=True, histtype='step', label='PNS')\n",
    "pl.hist(dpr2[mask],bins=np.logspace(np.log10(0.001), np.log10(100),250),\n",
    "        normed=True, histtype='step', label='PAV')\n",
    "pl.legend(fontsize=20, loc='lower right')\n",
    "pl.xscale('log', nonposy='clip')\n",
    "pl.grid()\n",
    "pl.xlabel('RR (mm/h)', fontsize=20)\n",
    "pl.ylabel('pdf', fontsize=20)\n",
    "pl.xticks(fontsize=20);pl.yticks(fontsize=20)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "A, B = ry.copy(), dpr.copy()\n",
    "C = sth.copy()\n",
    "\n",
    "from pcc import get_my_cmap2\n",
    "\n",
    "\n",
    "\n",
    "fig = pl.figure()\n",
    "ax = pl.gca()\n",
    "ax.scatter(A[mask],B[mask], c=C[mask], cmap=get_my_cmap2())\n",
    "cb = pl.colorbar()\n",
    "cb.set_label('sth', fontsize=ff)\n",
    "#pl.plot([0,1000], [0,1000], lw=2, color='black')\n",
    "ax.set_xscale('log')\n",
    "ax.set_yscale('log')\n",
    "pl.ylabel('DPR Rainrate in mm/h',fontsize=ff)\n",
    "pl.xlabel('Radolan Rainrate in mm/h',fontsize=ff)\n",
    "pl.xticks(fontsize=ff)\n",
    "pl.yticks(fontsize=ff)\n",
    "pl.xlim(0,1000)\n",
    "pl.ylim(0,1000)\n",
    "#plt.title('Corr: ' + str(round(corr,3)) + r'$\\pm$'+  str(round(eror,3))+\n",
    "#            '\\nBias: '+ str(round(bias,3)), fontsize=20)\n",
    "\n",
    "pl.grid()\n",
    "pl.show()\n",
    "print corcorcor(A,B)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "A, B = ry.copy(), dpr.copy()\n",
    "from pcc import get_my_cmap2\n",
    "\n",
    "xbins = 10**np.arange(-1, 2.08, 0.08)\n",
    "ybins = 10**np.arange(-1, 2.08, 0.08)\n",
    "\n",
    "counts, _, _ = np.histogram2d(A[mask],B[mask], bins=(xbins, ybins))\n",
    "\n",
    "fig, ax = pl.subplots(figsize=(10,10))\n",
    "pl.pcolormesh(xbins, ybins, counts.T, cmap=get_my_cmap2(),vmin=0.1)\n",
    "cb = pl.colorbar()\n",
    "cb.set_label('number of samples', fontsize=ff)\n",
    "pl.plot([0,1000], [0,1000], lw=2, color='black')\n",
    "ax.set_xscale('log')\n",
    "ax.set_yscale('log')\n",
    "pl.ylabel('DPR Rainrate in mm/h',fontsize=ff)\n",
    "pl.xlabel('Radolan Rainrate in mm/h',fontsize=ff)\n",
    "pl.xticks(fontsize=ff)\n",
    "pl.yticks(fontsize=ff)\n",
    "\n",
    "#plt.title('Corr: ' + str(round(corr,3)) + r'$\\pm$'+  str(round(eror,3))+\n",
    "#            '\\nBias: '+ str(round(bias,3)), fontsize=20)\n",
    "\n",
    "pl.grid()\n",
    "pl.show()\n",
    "print corcorcor(A,B)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "A, B = ry.copy(), dpr2.copy()\n",
    "from pcc import get_my_cmap2\n",
    "\n",
    "xbins = 10**np.arange(-1, 2.08, 0.08)\n",
    "ybins = 10**np.arange(-1, 2.08, 0.08)\n",
    "\n",
    "counts, _, _ = np.histogram2d(A[mask],B[mask], bins=(xbins, ybins))\n",
    "\n",
    "fig, ax = pl.subplots(figsize=(10,10))\n",
    "pl.pcolormesh(xbins, ybins, counts.T, cmap=get_my_cmap2(),vmin=0.1)\n",
    "cb = pl.colorbar()\n",
    "cb.set_label('number of samples', fontsize=ff)\n",
    "pl.plot([0,1000], [0,1000], lw=2, color='black')\n",
    "ax.set_xscale('log')\n",
    "ax.set_yscale('log')\n",
    "pl.ylabel('DPR Rainrate in mm/h',fontsize=ff)\n",
    "pl.xlabel('Radolan Rainrate in mm/h',fontsize=ff)\n",
    "pl.xticks(fontsize=ff)\n",
    "pl.yticks(fontsize=ff)\n",
    "\n",
    "#plt.title('Corr: ' + str(round(corr,3)) + r'$\\pm$'+  str(round(eror,3))+\n",
    "#            '\\nBias: '+ str(round(bias,3)), fontsize=20)\n",
    "\n",
    "pl.grid()\n",
    "pl.show()\n",
    "print corcorcor(A,B)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "### Overpass statistics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def get_radar_locations():\n",
    "\n",
    "    radars = {}\n",
    "    radar = {}\n",
    "    radar['name'] = 'ASR Dresden'\n",
    "    radar['wmo'] = 10487\n",
    "    radar['lon'] = 13.76347\n",
    "    radar['lat'] = 51.12404\n",
    "    radar['alt'] = 261\n",
    "    radars['ASD'] = radar\n",
    "    \n",
    "    radar = {}\n",
    "    #### Kordinaten un hoehe checken! Neu\n",
    "    # Hoehe üNN 2m\n",
    "    # Turm 55 m\n",
    "    #radar['name'] = 'Borkum'\n",
    "    #radar['wmo'] = 10103\n",
    "    #radar['lon'] = 6.748292\n",
    "    #radar['lat'] = 53.564011\n",
    "    #radar['alt'] = 57.\n",
    "    #radars['ASB'] = radar\n",
    "    \n",
    "    #Fake Borkum/Emden\n",
    "    radar['name'] = 'Emden'\n",
    "    radar['wmo'] = 10103\n",
    "    radar['lon'] = 7.02377\n",
    "    radar['lat'] = 53.33872\n",
    "    radar['alt'] = 58.\n",
    "    radars['ASB'] = radar\n",
    "\n",
    "    radar = {}\n",
    "    radar['name'] = 'Boostedt'\n",
    "    radar['wmo'] = 10132\n",
    "    radar['lon'] = 10.04687\n",
    "    radar['lat'] = 54.00438\n",
    "    radar['alt'] = 124.56\n",
    "    radars['BOO'] = radar\n",
    "\n",
    "    radar = {}\n",
    "    radar['name'] = 'Dresden'\n",
    "    radar['wmo'] = 10488\n",
    "    radar['lon'] = 13.76865\n",
    "    radar['lat'] = 51.12465\n",
    "    radar['alt'] = 263.36\n",
    "    radars['DRS'] = radar\n",
    "\n",
    "    radar = {}\n",
    "    radar['name'] = 'Eisberg'\n",
    "    radar['wmo'] = 10780\n",
    "    radar['lon'] = 12.40278\n",
    "    radar['lat'] = 49.54066\n",
    "    radar['alt'] = 798.79\n",
    "    radars['EIS'] = radar\n",
    "\n",
    "    radar = {}\n",
    "    radar['name'] = 'Emden'\n",
    "    radar['wmo'] = 10204\n",
    "    radar['lon'] = 7.02377\n",
    "    radar['lat'] = 53.33872\n",
    "    radar['alt'] = 58\n",
    "    radars['EMD'] = radar\n",
    "\n",
    "    radar = {}\n",
    "    radar['name'] = 'Essen'\n",
    "    radar['wmo'] = 10410\n",
    "    radar['lon'] = 6.96712\n",
    "    radar['lat'] = 51.40563\n",
    "    radar['alt'] = 185.10\n",
    "    radars['ESS'] = radar\n",
    "\n",
    "    radar = {}\n",
    "    radar['name'] = 'Feldberg'\n",
    "    radar['wmo'] = 10908\n",
    "    radar['lon'] = 8.00361\n",
    "    radar['lat'] = 47.87361\n",
    "    radar['alt'] = 1516.10\n",
    "    radars['FBG'] = radar\n",
    "\n",
    "    radar = {}\n",
    "    radar['name'] = 'Flechtdorf'\n",
    "    radar['wmo'] = 10440\n",
    "    radar['lon'] = 8.802\n",
    "    radar['lat'] = 51.3112\n",
    "    radar['alt'] = 627.88\n",
    "    radars['FLD'] = radar\n",
    "\n",
    "    radar = {}\n",
    "    radar['name'] = 'Hannover'\n",
    "    radar['wmo'] = 10339\n",
    "    radar['lon'] = 9.69452\n",
    "    radar['lat'] = 52.46008\n",
    "    radar['alt'] = 97.66\n",
    "    radars['HNR'] = radar\n",
    "\n",
    "    radar = {}\n",
    "    radar['name'] = 'Neuhaus'\n",
    "    radar['wmo'] = 10557\n",
    "    radar['lon'] = 11.13504\n",
    "    radar['lat'] = 50.50012\n",
    "    radar['alt'] = 878.04\n",
    "    radars['NEU'] = radar\n",
    "\n",
    "    radar = {}\n",
    "    radar['name'] = 'Neuheilenbach'\n",
    "    radar['wmo'] = 10605\n",
    "    radar['lon'] = 6.54853\n",
    "    radar['lat'] = 50.10965\n",
    "    radar['alt'] = 585.84\n",
    "    radars['NHB'] = radar\n",
    "\n",
    "    radar = {}\n",
    "    radar['name'] = 'Offenthal'\n",
    "    radar['wmo'] = 10629\n",
    "    radar['lon'] = 8.71293\n",
    "    radar['lat'] = 49.9847\n",
    "    radar['alt'] = 245.80\n",
    "    radars['OFT'] = radar\n",
    "\n",
    "    radar = {}\n",
    "    radar['name'] = 'Proetzel'\n",
    "    radar['wmo'] = 10392\n",
    "    radar['lon'] = 13.85821\n",
    "    radar['lat'] = 52.64867\n",
    "    radar['alt'] = 193.92\n",
    "    radars['PRO'] = radar\n",
    "\n",
    "    radar = {}\n",
    "    radar['name'] = 'Memmingen'\n",
    "    radar['wmo'] = 10950\n",
    "    radar['lon'] = 10.21924\n",
    "    radar['lat'] = 48.04214\n",
    "    radar['alt'] = 724.40\n",
    "    radars['MEM'] = radar\n",
    "\n",
    "    radar = {}\n",
    "    radar['name'] = 'Rostock'\n",
    "    radar['wmo'] = 10169\n",
    "    radar['lon'] = 12.05808\n",
    "    radar['lat'] = 54.17566\n",
    "    radar['alt'] = 37\n",
    "    radars['ROS'] = radar\n",
    "\n",
    "    radar = {}\n",
    "    radar['name'] = 'Isen'\n",
    "    radar['wmo'] = 10873\n",
    "    radar['lon'] = 12.10177\n",
    "    radar['lat'] = 48.1747\n",
    "    radar['alt'] = 677.77\n",
    "    radars['ISN'] = radar\n",
    "\n",
    "    radar = {}\n",
    "    radar['name'] = 'Tuerkheim'\n",
    "    radar['wmo'] = 10832\n",
    "    radar['lon'] = 9.78278\n",
    "    radar['lat'] = 48.58528\n",
    "    radar['alt'] = 767.62\n",
    "    radars['TUR'] = radar\n",
    "\n",
    "    radar = {}\n",
    "    radar['name'] = 'Ummendorf'\n",
    "    radar['wmo'] = 10356\n",
    "    radar['lon'] = 11.17609\n",
    "    radar['lat'] = 52.16009\n",
    "    radar['alt'] = 183\n",
    "    radars['UMM'] = radar\n",
    "\n",
    "    \n",
    "    return radars\n",
    "\n",
    "def plot_radar(radar, ax, proj):\n",
    "\n",
    "    site = (radar['lon'], radar['lat'], radar['alt'] )\n",
    "    \n",
    "    r = np.arange(1, 151) * 1000\n",
    "\n",
    "    # azimuth array 1 degree spacing\n",
    "    az = np.linspace(0, 360, 361)[0:-1]\n",
    "\n",
    "    # build polygons for maxrange rangering\n",
    "    polygons = wrl.georef.spherical_to_polyvert(r, az, 0,\n",
    "                                                site,\n",
    "                                                proj=proj)\n",
    "    polygons = polygons[..., 0:2]\n",
    "    polygons.shape = (len(az), len(r), 5, 2)\n",
    "    polygons = polygons[:, -1, :, :]\n",
    "\n",
    "    x_loc, y_loc = wrl.georef.reproject(site[0], site[1],\n",
    "                                        projection_source=proj_wgs,\n",
    "                                        projection_target=proj)\n",
    "\n",
    "\n",
    "    # create PolyCollections and add to respective axes\n",
    "    polycoll = mpl.collections.PolyCollection(polygons, closed=True,\n",
    "                                              edgecolors='k',\n",
    "                                              facecolors='k',\n",
    "                                             linewidths=0.1)\n",
    "    ax.add_collection(polycoll, autolim=True)\n",
    "\n",
    "    # plot radar location and information text\n",
    "    #ax.scatter(x_loc, y_loc, s=40, color='red', edgecolors='black')\n",
    "    #ax.text(x_loc, y_loc, radar['name'], color='k')\n",
    "\n",
    "def plot_all_cband2(ax):\n",
    "    for i in get_radar_locations().keys():\n",
    "\n",
    "        if i =='EMD':\n",
    "            pass\n",
    "        else:\n",
    "            plot_radar(get_radar_locations()[i], ax , proj_stereo)\n",
    "            \n",
    "# create radolan projection object\n",
    "proj_stereo = wrl.georef.create_osr(\"dwd-radolan\")\n",
    "\n",
    "# create wgs84 projection object\n",
    "proj_wgs = osr.SpatialReference()\n",
    "proj_wgs.ImportFromEPSG(4326)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "np.unique(sr_cfb)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%%time\n",
    "\n",
    "cc=0.5\n",
    "ff=20\n",
    "\n",
    "def get_miub_cmap():\n",
    "    import matplotlib.colors as col\n",
    "    startcolor = 'white'  # a dark olive\n",
    "    color1 = '#8ec7ff'#'cyan'    # a bright yellow\n",
    "    color2 = 'dodgerblue'\n",
    "    color3 = 'lime'\n",
    "    color4 = 'yellow'\n",
    "    color5 = 'darkorange'\n",
    "    color6 = 'red'\n",
    "    color7 = 'purple'\n",
    "    #color6 = 'grey'\n",
    "    endcolor = 'darkmagenta'    # medium dark red\n",
    "    colors = [startcolor, color1, color2, color3, color4, color5, color6, endcolor]\n",
    "    return col.LinearSegmentedColormap.from_list('miub1',colors)\n",
    "\n",
    "dpr = sr_pns.copy() # para2 = GPM \n",
    "ry = gr_pns.copy() # para1 = RY\n",
    "xx = xl.copy()\n",
    "yy = yl.copy()\n",
    "\n",
    "#yy[(ry == 0) & (dpr > 0)] = np.nan\n",
    "#xx[(ry == 0) & (dpr > 0)] = np.nan\n",
    "\n",
    "#yy[(ry > 0) & (dpr == 0)] = np.nan\n",
    "#xx[(ry > 0) & (dpr == 0)] = np.nan\n",
    "\n",
    "#yy[(ry > 0) & (dpr > 0)] = np.nan\n",
    "#xx[(ry > 0) & (dpr > 0)] = np.nan\n",
    "\n",
    "#hd =gr_hmin - (sr_cfb - sr_ele) \n",
    "#yy[hd<=1000] = np.nan\n",
    "#xx[hd<=1000] = np.nan\n",
    "\n",
    "#yy[(sr_cfb-sr_ele)<=2000] = np.nan\n",
    "#xx[(sr_cfb-sr_ele)<=2000] = np.nan\n",
    "\n",
    "#yy[sr_bbh < gr_hmin] = np.nan\n",
    "#xx[sr_bbh < gr_hmin] = np.nan\n",
    "\n",
    "yy[sr_sth < 5000] = np.nan\n",
    "xx[sr_sth < 5000] = np.nan\n",
    "\n",
    "#yy[sr_cfb<=2000] = np.nan\n",
    "#xx[sr_cfb<=2000] = np.nan\n",
    "\n",
    "#yy[ry<=0] = np.nan\n",
    "#xx[ry<=0] = np.nan\n",
    "#yy[dpr<=0] = np.nan\n",
    "#xx[dpr<=0] = np.nan\n",
    "#yy[sr_bbh<=2] = np.nan\n",
    "#xx[sr_bbh<=2] = np.nan\n",
    "#yy[typ!=1] = np.nan\n",
    "#xx[typ!=1] = np.nan\n",
    "#yy[(ry <= 0) & (dpr <= 0)] = np.nan\n",
    "#xx[(ry <= 0) & (dpr <= 0)] = np.nan\n",
    "\n",
    "mask = ~np.isnan(yy) & ~np.isnan(xx)\n",
    "\n",
    "mask = ~np.isnan(yy) & ~np.isnan(xx)\n",
    "fig = pl.figure(figsize=(12,12))\n",
    "#plt.subplots_adjust(left=0.1, bottom=0.1, right=0.8, top=0.8, wspace=0.2, hspace=0.2)\n",
    "\n",
    "ax1 = fig.add_subplot(111, aspect='equal')\n",
    "pl.hist2d(xx[mask],yy[mask], bins=100, cmap=get_miub_cmap())\n",
    "cb = pl.colorbar(shrink=cc)\n",
    "cb.set_label(\"#\", fontsize=ff)\n",
    "cb.ax.tick_params(labelsize=ff)\n",
    "pl.xlabel('x (km)',fontsize=ff)\n",
    "pl.ylabel('y (km)',fontsize=ff)\n",
    "#pl.title('Number of Overpasses 2014-2017 \\n DPR rain:yes - RADOLAN rain:yes',fontsize=ff)\n",
    "\n",
    "pl.title('Number of Overpasses 2014-2018 \\n CFB SR >= 1km',fontsize=ff)\n",
    "pl.xticks(fontsize=ff)\n",
    "pl.yticks(fontsize=ff)\n",
    "plot_all_cband2(ax1)\n",
    "pl.xlim(-420, 377);\n",
    "pl.ylim(-4660, -3759)\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": [
    "%%time\n",
    "ff=20; cc=0.5\n",
    "dpr = sr_pns.copy() # para2 = GPM \n",
    "ry = gr_pns.copy() # para1 = RY\n",
    "xx = xl.copy()\n",
    "yy = yl.copy()\n",
    "\n",
    "#yy[(ry == 0) & (dpr > 0)] = np.nan\n",
    "#xx[(ry == 0) & (dpr > 0)] = np.nan\n",
    "\n",
    "#yy[(ry <= 0) & (dpr <= 0)] = np.nan\n",
    "#xx[(ry <= 0) & (dpr <= 0)] = np.nan\n",
    "\n",
    "mask = ~np.isnan(yy) & ~np.isnan(xx)\n",
    "pl.figure(figsize=(10,10))\n",
    "pl.hist2d(xx[mask],yy[mask], bins=200)\n",
    "cb = pl.colorbar(shrink=cc)\n",
    "cb.set_label(\"#\", fontsize=ff)\n",
    "cb.ax.tick_params(labelsize=ff)\n",
    "pl.xlabel('x (km)',fontsize=ff)\n",
    "pl.ylabel('y (km)',fontsize=ff)\n",
    "pl.title('Number of Overpasses 2014-2017',fontsize=ff)\n",
    "pl.xticks(fontsize=ff)\n",
    "pl.yticks(fontsize=ff)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "# Richtige einteilung wäre \n",
    "----------------------------------------------\n",
    "##### Aufgrund er Regenraten berechnung des GPM\n",
    "\n",
    "$$R = \\epsilon ^r \\cdot p D_m ^q$$\n",
    "\n",
    "$\\epsilon -> k-Ze$\n",
    "\n",
    "\n",
    "$r, q, p -> raintyp$\n",
    "\n",
    "\n",
    "$Dm -> phase$\n",
    "\n",
    "Tabelleneinteilung \n",
    "raintyp ::: -  stratifrom S convective C\n",
    "\n",
    "\n",
    "phase ::: -  Liquid L und Solid S\n",
    "\n",
    "\n",
    "rainintensity ::: -  light l moderate m heavy h\n",
    "\n",
    "> l ::: 0.5 mm/h <= rr < 1 mm/h \n",
    "\n",
    "> m ::: 1 mm/h <= rr < 10 mm/h\n",
    "\n",
    "> h ::: rr > 10 mm/h\n",
    "\n",
    "    = 1    stratiform\n",
    "    = 2    convective\n",
    "    = 3    other\n",
    "\n",
    "\n",
    "0 - solid\n",
    "1 - mixed phase\n",
    "2 - liquid\n",
    "255 - Missing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "srp, grp = sr_pns.copy(), gr_pns.copy()\n",
    "\n",
    "SR_SIL = srp[(typ==1.) & (pha==0) & (grp>=0.5) & (grp<1.0)]\n",
    "SR_SIM = srp[(typ==1.) & (pha==0) & (grp>=1.0) & (grp<10.)]\n",
    "SR_SIH = srp[(typ==1.) & (pha==0) & (grp>10.)]\n",
    "\n",
    "SR_SLL = srp[(typ==1.) & (pha==2) & (grp>=0.5) & (grp<1.0)]\n",
    "SR_SLM = srp[(typ==1.) & (pha==2) & (grp>=1.0) & (grp<10.)]\n",
    "SR_SLH = srp[(typ==1.) & (pha==2) & (grp>10.)]\n",
    "\n",
    "SR_CIL = srp[(typ==2.) & (pha==0) & (grp>=0.5) & (grp<1.0)]\n",
    "SR_CIM = srp[(typ==2.) & (pha==0) & (grp>=1.0) & (grp<10.)]\n",
    "SR_CIH = srp[(typ==2.) & (pha==0) & (grp>10.)]\n",
    "\n",
    "SR_CLL = srp[(typ==2.) & (pha==2) & (grp>=0.5) & (grp<1.0)]\n",
    "SR_CLM = srp[(typ==2.) & (pha==2) & (grp>=1.0) & (grp<10.)]\n",
    "SR_CLH = srp[(typ==2.) & (pha==2) & (grp>10.)]\n",
    "\n",
    "\n",
    "GR_SIL = grp[(typ==1.) & (pha==0) & (grp>=0.5) & (grp<1.0)]\n",
    "GR_SIM = grp[(typ==1.) & (pha==0) & (grp>=1.0) & (grp<10.)]\n",
    "GR_SIH = grp[(typ==1.) & (pha==0) & (grp>10.)]\n",
    "\n",
    "GR_SLL = grp[(typ==1.) & (pha==2) & (grp>=0.5) & (grp<1.0)]\n",
    "GR_SLM = grp[(typ==1.) & (pha==2) & (grp>=1.0) & (grp<10.)]\n",
    "GR_SLH = grp[(typ==1.) & (pha==2) & (grp>10.)]\n",
    "\n",
    "GR_CIL = grp[(typ==2.) & (pha==0) & (grp>=0.5) & (grp<1.0)]\n",
    "GR_CIM = grp[(typ==2.) & (pha==0) & (grp>=1.0) & (grp<10.)]\n",
    "GR_CIH = grp[(typ==2.) & (pha==0) & (grp>10.)]\n",
    "\n",
    "GR_CLL = grp[(typ==2.) & (pha==2) & (grp>=0.5) & (grp<1.0)]\n",
    "GR_CLM = grp[(typ==2.) & (pha==2) & (grp>=1.0) & (grp<10.)]\n",
    "GR_CLH = grp[(typ==2.) & (pha==2) & (grp>10.)]\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "srh = sr_cfb-sr_ele\n",
    "grh = gr_hmin \n",
    "\n",
    "hdiff = abs(grh - srh)\n",
    "\n",
    "srp, grp = sr_pns.copy(), gr_pns.copy()\n",
    "\n",
    "SR_SIL = srp[(typ==1.) & (pha==0) & (grp<1.0)]\n",
    "SR_SIM = srp[(typ==1.) & (pha==0) & (grp>=1.0) & (grp<10.)]\n",
    "SR_SIH = srp[(typ==1.) & (pha==0) & (grp>10.)]\n",
    "\n",
    "SR_SLL = srp[(typ==1.) & (pha==2) &  (grp<1.0)]\n",
    "SR_SLM = srp[(typ==1.) & (pha==2) & (grp>=1.0) & (grp<10.)]\n",
    "SR_SLH = srp[(typ==1.) & (pha==2) & (grp>10.)]\n",
    "\n",
    "SR_CIL = srp[(typ==2.) & (pha==0) & (grp<1.0)]\n",
    "SR_CIM = srp[(typ==2.) & (pha==0) & (grp>=1.0) & (grp<10.)]\n",
    "SR_CIH = srp[(typ==2.) & (pha==0) & (grp>10.)]\n",
    "\n",
    "SR_CLL = srp[(typ==2.) & (pha==2) & (grp<1.0)]\n",
    "SR_CLM = srp[(typ==2.) & (pha==2) & (grp>=1.0) & (grp<10.)]\n",
    "SR_CLH = srp[(typ==2.) & (pha==2) & (grp>10.)]\n",
    "\n",
    "\n",
    "GR_SIL = grp[(typ==1.) & (pha==0) & (grp<1.0)]\n",
    "GR_SIM = grp[(typ==1.) & (pha==0) & (grp>=1.0) & (grp<10.)]\n",
    "GR_SIH = grp[(typ==1.) & (pha==0) & (grp>10.)]\n",
    "\n",
    "GR_SLL = grp[(typ==1.) & (pha==2) & (grp<1.0)]\n",
    "GR_SLM = grp[(typ==1.) & (pha==2) & (grp>=1.0) & (grp<10.)]\n",
    "GR_SLH = grp[(typ==1.) & (pha==2) & (grp>10.)]\n",
    "\n",
    "GR_CIL = grp[(typ==2.) & (pha==0) & (grp<1.0)]\n",
    "GR_CIM = grp[(typ==2.) & (pha==0) & (grp>=1.0) & (grp<10.)]\n",
    "GR_CIH = grp[(typ==2.) & (pha==0) & (grp>10.)]\n",
    "\n",
    "GR_CLL = grp[(typ==2.) & (pha==2) & (grp<1.0)]\n",
    "GR_CLM = grp[(typ==2.) & (pha==2) & (grp>=1.0) & (grp<10.)]\n",
    "GR_CLH = grp[(typ==2.) & (pha==2) & (grp>10.)]\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "srh = sr_cfb-sr_ele\n",
    "grh = gr_hmin \n",
    "\n",
    "hdiff = abs(grh - srh)\n",
    "\n",
    "srp, grp = sr_pns.copy(), gr_pns.copy()\n",
    "\n",
    "SR_SIL = srp[(typ==1.) & (pha==0) & (grp<1.0)& (hdiff<=100)]\n",
    "SR_SIM = srp[(typ==1.) & (pha==0) & (grp>=1.0) & (grp<10.)& (hdiff<=100)]\n",
    "SR_SIH = srp[(typ==1.) & (pha==0) & (grp>10.)& (hdiff<=100)]\n",
    "\n",
    "SR_SLL = srp[(typ==1.) & (pha==2) &  (grp<1.0)& (hdiff<=100)]\n",
    "SR_SLM = srp[(typ==1.) & (pha==2) & (grp>=1.0) & (grp<10.)& (hdiff<=100)]\n",
    "SR_SLH = srp[(typ==1.) & (pha==2) & (grp>10.)& (hdiff<=100)]\n",
    "\n",
    "SR_CIL = srp[(typ==2.) & (pha==0) & (grp<1.0)& (hdiff<=100)]\n",
    "SR_CIM = srp[(typ==2.) & (pha==0) & (grp>=1.0) & (grp<10.)& (hdiff<=100)]\n",
    "SR_CIH = srp[(typ==2.) & (pha==0) & (grp>10.)& (hdiff<=100)]\n",
    "\n",
    "SR_CLL = srp[(typ==2.) & (pha==2) & (grp<1.0)& (hdiff<=100)]\n",
    "SR_CLM = srp[(typ==2.) & (pha==2) & (grp>=1.0) & (grp<10.)& (hdiff<=100)]\n",
    "SR_CLH = srp[(typ==2.) & (pha==2) & (grp>10.)& (hdiff<=100)]\n",
    "\n",
    "\n",
    "GR_SIL = grp[(typ==1.) & (pha==0) & (grp<1.0)& (hdiff<=100)]\n",
    "GR_SIM = grp[(typ==1.) & (pha==0) & (grp>=1.0) & (grp<10.)& (hdiff<=100)]\n",
    "GR_SIH = grp[(typ==1.) & (pha==0) & (grp>10.)& (hdiff<=100)]\n",
    "\n",
    "GR_SLL = grp[(typ==1.) & (pha==2) & (grp<1.0)& (hdiff<=100)]\n",
    "GR_SLM = grp[(typ==1.) & (pha==2) & (grp>=1.0) & (grp<10.)& (hdiff<=100)]\n",
    "GR_SLH = grp[(typ==1.) & (pha==2) & (grp>10.)& (hdiff<=100)]\n",
    "\n",
    "GR_CIL = grp[(typ==2.) & (pha==0) & (grp<1.0)& (hdiff<=100)]\n",
    "GR_CIM = grp[(typ==2.) & (pha==0) & (grp>=1.0) & (grp<10.)& (hdiff<=100)]\n",
    "GR_CIH = grp[(typ==2.) & (pha==0) & (grp>10.)& (hdiff<=100)]\n",
    "\n",
    "GR_CLL = grp[(typ==2.) & (pha==2) & (grp<1.0)& (hdiff<=100)]\n",
    "GR_CLM = grp[(typ==2.) & (pha==2) & (grp>=1.0) & (grp<10.)& (hdiff<=100)]\n",
    "GR_CLH = grp[(typ==2.) & (pha==2) & (grp>10.)& (hdiff<=100)]\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "print (SR_SLL.shape ,SR_SLM.shape , SR_SLH.shape )\n",
    "print (GR_CIL.shape ,GR_CIM.shape , GR_CIH.shape )\n",
    "print (GR_CLL.shape ,GR_CLM.shape , GR_CLH.shape )\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "sr_cfb, gr_hmin"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "ff=20\n",
    "pl.figure(figsize=(25,25))\n",
    "ax1 = pl.subplot(4,3,1, aspect='equal')\n",
    "plotplot(ax1,GR_SLL ,SR_SLL, 'SLL')\n",
    "\n",
    "\n",
    "ax1 = pl.subplot(4,3,2, aspect='equal')\n",
    "plotplot(ax1,GR_SLM ,SR_SLM, 'SLM')\n",
    "\n",
    "ax1 = pl.subplot(4,3,3, aspect='equal')\n",
    "plotplot(ax1,GR_SLH ,SR_SLH, 'SLH')\n",
    "\n",
    "ax1 = pl.subplot(4,3,4, aspect='equal')\n",
    "plotplot(ax1,GR_SIL ,SR_SIL, 'SIL')\n",
    "\n",
    "\n",
    "ax1 = pl.subplot(4,3,5, aspect='equal')\n",
    "plotplot(ax1,GR_SIM ,SR_SIM, 'SIM')\n",
    "\n",
    "ax1 = pl.subplot(4,3,6, aspect='equal')\n",
    "plotplot(ax1,GR_SIH ,SR_SIH, 'SIH')\n",
    "\n",
    "ax1 = pl.subplot(4,3,7, aspect='equal')\n",
    "plotplot(ax1,GR_CLL ,SR_CLL, 'CLL')\n",
    "\n",
    "\n",
    "ax1 = pl.subplot(4,3,8, aspect='equal')\n",
    "plotplot(ax1,GR_CLM ,SR_CLM, 'CLM')\n",
    "\n",
    "ax1 = pl.subplot(4,3,9, aspect='equal')\n",
    "plotplot(ax1,GR_CLH ,SR_CLH, 'CLH')\n",
    "\n",
    "ax1 = pl.subplot(4,3,10, aspect='equal')\n",
    "plotplot(ax1,GR_CIL ,SR_CIL, 'CIL')\n",
    "\n",
    "\n",
    "ax1 = pl.subplot(4,3,11, aspect='equal')\n",
    "plotplot(ax1,GR_CIM ,SR_CIM, 'CIM')\n",
    "\n",
    "ax1 = pl.subplot(4,3,12, aspect='equal')\n",
    "plotplot(ax1,GR_CIH ,SR_CIH, 'CIH')\n",
    "pl.tight_layout()\n",
    "pl.savefig(\"/automount/ftp/velibor/test2.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": [
    "pd.DataFrame.from_dict(skill_score(GR_SLM, SR_SLM,0.1), orient='index')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "GR_SLM.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "# Times\n",
    "----------------------"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import pandas as pd\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "sr_time = np.load('/automount/ags/velibor/gpmdata/dumpdataV7/dpr_ry_NS_compress/sr_times.npy')\n",
    "sr_pns = np.load('/automount/ags/velibor/gpmdata/dumpdataV7/dpr_ry_NS_compress/sr_pns.npy')\n",
    "gr_pns = np.load('/automount/ags/velibor/gpmdata/dumpdataV7/dpr_ry_NS_compress/gr_pp_ipoli.npy')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "sr_time.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "df = pd.DataFrame(dict(sr_pns=sr_pns,\n",
    "                       gr_pns=gr_pns,\n",
    "                       sr_gr = sr_pns-gr_pns,time=sr_time))\n",
    "\n",
    "formater = '%Y-%m-%d %H:%M:%S'\n",
    "df['Datetime'] = pd.to_datetime(df['time'], format=formater)\n",
    "df = df.set_index(pd.DatetimeIndex(df['Datetime']))\n",
    "df = df.drop(['Datetime','time'], axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#a = df['sr_gr'].groupby([df.index, pd.TimeGrouper(freq='M')])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "df['sr_gr'].plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plt.show()"
   ]
  },
  {
   "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": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "df['time'] = pd.to_datetime(df['time'])  \n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "df = df.set_index('time')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plt.subplot(3,1,1)\n",
    "df['gr_pns'][0:10000].plot()\n",
    "plt.subplot(3,1,2)\n",
    "\n",
    "df['sr_pns'][0:10000].plot()\n",
    "\n",
    "plt.subplot(3,1,3)\n",
    "\n",
    "df['sr_gr'][0:10000].plot()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "df['sr_gr'].groupby(df.index).mean()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "df['sr_gr'].groupby([df.index, pd.TimeGrouper(freq='Y')])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "df.index"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plt.semilogy(sr_time, sr_pns - gr_pns)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "type(time)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import datetime\n",
    "time = np.array(sr_time,dtype='datetime64')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "time\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plt.plot(time, sr_pns - gr_pns)\n",
    "#pl.ylim(1e-6,1e6)\n",
    "plt.grid()\n",
    "plt.axhline(0,linewidth=4, color='r' )\n",
    "plt.show()"
   ]
  },
  {
   "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
}