pcc_gpm_matching2-Copy1.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import warnings\n",
    "warnings.filterwarnings('ignore')\n",
    "import wradlib as wrl\n",
    "import matplotlib.pyplot as pl\n",
    "import matplotlib as mpl\n",
    "from matplotlib.collections import PatchCollection\n",
    "from matplotlib.colors import from_levels_and_colors\n",
    "from matplotlib.path import Path\n",
    "import matplotlib.patches as patches\n",
    "import matplotlib.cm as cm\n",
    "try:\n",
    "    get_ipython().magic(\"matplotlib inline\")\n",
    "except:\n",
    "    pl.ion()\n",
    "import numpy as np\n",
    "import datetime as dt\n",
    "from osgeo import osr\n",
    "warnings.filterwarnings('ignore')\n",
    "from satlib import corcor\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# define GPM data set\n",
    "# ----------------------\n",
    "gpm_file_path = '/automount/ags/velibor/gpmdata/dprV7/2A.GPM.DPR.V7-20170308.20141007-S015721-E032951.003445.V05A.HDF5'\n",
    "#gpm_file_path = '/automount/ags/velibor/gpmdata/dprV7/2A.GPM.DPR.V7-20170308.20141008-S084747-E102017.003465.V05A.HDF5'\n",
    "#gpm_file_path = '/automount/ags/velibor/gpmdata/dprV7/2A.GPM.DPR.V7-20170308.20160601-S171836-E185110.012837.V05A.HDF5'\n",
    "#gpm_file_path = '/automount/ags/velibor/gpmdata/dprV7/2A.GPM.DPR.V7-20170308.20150330-S224903-E002135.006166.V05A.HDF5'\n",
    "#gpm_file_path = '/automount/ags/velibor/gpmdata/dprV7/2A.GPM.DPR.V7-20170308.20151216-S014742-E032015.010213.V05A.HDF5'\n",
    "\n",
    "# define matching ground radar file\n",
    "# ----------------------\n",
    "gr_file_path = '/automount/radar-archiv/scans/2014/2014-10/2014-10-07/n_ppi_010deg/2014-10-07--02:37:44,00.mvol'\n",
    "#gr_file_path = '/automount/radar-archiv/scans/2014/2014-10/2014-10-08/n_ppi_010deg/2014-10-08--09:45:00,00.mvol'\n",
    "#gr_file_path = '/automount/radar/scans/2016/2016-06/2016-06-01/n_ppi_010deg/2016-06-01--17:59:50,00.mvol'\n",
    "#gr_file_path = '/automount/radar/scans/2015/2015-03/2015-03-30/n_ppi_010deg/2015-03-30--23:30:03,00.mvol'\n",
    "#gr_file_path = '/automount/radar/scans/2015/2015-12/2015-12-16/n_ppi_010deg/2015-12-16--02:45:01,00.mvol'\n",
    "\n",
    "gpm_file = wrl.util.get_wradlib_data_file(gpm_file_path)\n",
    "gr2gpm_file = wrl.util.get_wradlib_data_file(gr_file_path)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Space-born precipitation radar parameters\n",
    "sr_pars = {\"gpm\": {\n",
    "    \"zt\": 407000.,  # orbital height of GPM                 APPROXIMATION!\n",
    "    \"dr\": 125.,      # gate spacing of GPM\n",
    "    \"gr_file\": gr2gpm_file,\n",
    "}}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Set parameters for this procedure\n",
    "bw_sr = 0.71                  # SR beam width\n",
    "platf = \"gpm\"                 # SR platform/product: one out of [\"gpm\", \"trmm\"]\n",
    "zt = sr_pars[platf][\"zt\"]     # SR orbit height (meters)\n",
    "dr_sr = sr_pars[platf][\"dr\"]  # SR gate length (meters)\n",
    "gr_file = sr_pars[platf][\"gr_file\"]\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def read_gr_boxpol(filename, loaddata=True):\n",
    "    \"\"\"\n",
    "    Read ground radar boxpol scans\n",
    "    \n",
    "    in arbeit....\n",
    "    \n",
    "    \"\"\"\n",
    "    sdate = []\n",
    "    refl = []\n",
    "    \n",
    "    gr_data = wrl.io.read_generic_netcdf(gr_file)\n",
    "    grdata, grattrs = wrl.io.read_GAMIC_hdf5(gr_file)\n",
    "    print gr_file\n",
    "    \n",
    "    dat_ = gr_data['what']['date']\n",
    "    #print dat_\n",
    "    #dat = dat_[0:10]\n",
    "    #tim = dat_[11:-1]\n",
    "    #print (dat, tim)\n",
    "    # '2014-10-07T02:42:02Z'\n",
    "    #date = dt.datetime.strptime(dat + tim, \"%Y-%d-%m%H:%M:%S\")\n",
    "    #source = gr_data['what']['source']\n",
    "\n",
    "    lon = gr_data['where']['lon']\n",
    "    lat = gr_data['where']['lat']\n",
    "    alt = gr_data['where']['height']\n",
    "\n",
    "    import glob\n",
    "    c_liste = glob.glob(gr_file.replace('n_ppi_010deg','*ppi*'))\n",
    "    print c_liste\n",
    "    ntilt = len(c_liste)\n",
    "    \n",
    "    ref_data = [[]]*ntilt\n",
    "    rho_data = [[]]*ntilt\n",
    "    zdr_data = [[]]*ntilt\n",
    "\n",
    "\n",
    "    print ('ntilit:',ntilt)    \n",
    "    ngate = np.zeros(ntilt, dtype=np.int16)\n",
    "    nbeam = np.zeros(ntilt)\n",
    "    elang = np.zeros(ntilt)\n",
    "    r0 = np.zeros(ntilt)\n",
    "    dr = np.zeros(ntilt)\n",
    "    a0 = np.zeros(ntilt)\n",
    "    \n",
    "    for jj in range(ntilt):\n",
    "        print c_liste[jj]\n",
    "        \n",
    "        gr_data = wrl.io.read_generic_netcdf(c_liste[jj])\n",
    "        grdata, grattrs = wrl.io.read_GAMIC_hdf5(c_liste[jj])\n",
    "\n",
    "        #for i in range(0, ntilt):\n",
    "        dset = gr_data['scan0']\n",
    "        # Azimuth Start\n",
    "        a0[jj] = dset['how']['azi_start'] \n",
    "        # Elevations Winkel\n",
    "        elang[jj] = dset['how']['elevation']\n",
    "        # Anzahl bins\n",
    "        ngate[jj] = dset['how']['bin_count']\n",
    "        # Start und Step der range\n",
    "        r0[jj] = dset['how']['range_start']\n",
    "        dr[jj] = dset['how']['range_step']\n",
    "        # Anzahl Rays\n",
    "        nbeam[jj] = dset['how']['ray_count']\n",
    "\n",
    "        #print ('azi start: ', a0)    \n",
    "        #print ('elevation: ', elang)    \n",
    "        #print ('bin count: ', ngate)    \n",
    "        #print ('range start: ', r0)    \n",
    "        #print ('range step: ', dr)    \n",
    "        #print ('ray count: ', nbeam)    \n",
    "\n",
    "\n",
    "        gr_dict = {}\n",
    "        gr_dict.update({'date': dat_, 'lon': lon, 'lat': lat,\n",
    "                        'alt': alt, 'ngate': ngate, 'nbeam': nbeam, 'ntilt': ntilt,\n",
    "                        'r0': r0, 'dr': dr, 'a0': a0, 'elang': elang})\n",
    "        if not loaddata:\n",
    "            return gr_dict\n",
    "\n",
    "\n",
    "\n",
    "        #for i in range(0, ntilt):\n",
    "            #dset = gr_data['scan0']\n",
    "            #data = dset['variables']['moment_10']['data']\n",
    "        #dh = grdata['SCAN0']['ZH']['data']\n",
    "        #dv = grdata['SCAN0']['ZV']['data']\n",
    "        #data = (dh + dv)/2.\n",
    "        data = grdata['SCAN0']['ZH']['data']\n",
    "        \n",
    "        rho = grdata['SCAN0']['RHOHV']['data']\n",
    "        \n",
    "        zdr = grdata['SCAN0']['ZDR']['data']\n",
    "        \n",
    "\n",
    "        #refl.append(data)\n",
    "        ref_data[jj]=data\n",
    "        rho_data[jj]=rho\n",
    "        zdr_data[jj]=zdr\n",
    "        #ref_data[ref_data<=0]=np.nan\n",
    "\n",
    "        #refl1 = np.array(refl)\n",
    "    gr_dict.update({'refl': ref_data})\n",
    "    gr_dict.update({'rho': rho_data})\n",
    "    gr_dict.update({'zdr': zdr_data})\n",
    "    \n",
    "    return gr_dict, ntilt\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# read matching GR data\n",
    "gr_data, ntilt = read_gr_boxpol(gr_file)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Beginn der QnD shitty For schleife"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "SR, GR =  np.array([]), np.array([])\n",
    "X, Y, Z = np.array([]), np.array([]), np.array([])\n",
    "\n",
    "GR_zdr = np.array([])\n",
    "\n",
    "gr_nr_elev, gr_ele = np.array([]), np.array([])\n",
    "\n",
    "\n",
    "for ijk in range(ntilt):\n",
    "#for ijk in range(1,6,2):\n",
    "\n",
    "    ######### GR data and att\n",
    "    ee = ijk\n",
    "    # number of rays in gr sweep\n",
    "    nray_gr = gr_data['nbeam'].astype(\"i4\")[ee]\n",
    "    # number of gates in gr beam\n",
    "    ngate_gr = gr_data['ngate'].astype(\"i4\")[ee]\n",
    "    # number of sweeps\n",
    "    nelev = gr_data['ntilt']\n",
    "    # elevation of sweep (degree)\n",
    "    elev_gr = gr_data['elang'][ee]\n",
    "    # gate length (meters)\n",
    "    dr_gr = gr_data['dr'][ee]\n",
    "    # sweep datetime stamp\n",
    "    #date_gr = gr_data['sdate'][ee]\n",
    "    # range of first gate\n",
    "    r0_gr = gr_data['r0'][ee]\n",
    "    # azimuth angle of first beam\n",
    "    a0_gr = gr_data['a0'][ee]\n",
    "    # Longitude of GR\n",
    "    lon0_gr = gr_data['lon']\n",
    "    # Latitude of GR\n",
    "    lat0_gr = gr_data['lat']\n",
    "    # Altitude of GR (meters)\n",
    "    alt0_gr = gr_data['alt']\n",
    "    # Beam width of GR (degree)\n",
    "    bw_gr = 1.\n",
    "    # reflectivity array of sweep\n",
    "    ref_gr = gr_data['refl'][ee]\n",
    "    # rhohv array of sweep\n",
    "    rho_gr = gr_data['rho'][ee]\n",
    "    # zdr array of sweep\n",
    "    zdr_gr = gr_data['zdr'][ee]\n",
    "    \n",
    "    if elev_gr<0.01:\n",
    "        print (\"Elevation zu niedrig\")\n",
    "        pass\n",
    "    else:\n",
    "        ############################################################################## Korrektur der GR Daten\n",
    "        print (\"----------------Elevationswinkel: \", elev_gr)\n",
    "        print (\"________Offset______\")\n",
    "        #Hardcoded....passt in etwa\n",
    "        #ref_gr = ref_gr + 2\n",
    "\n",
    "\n",
    "        print (\"________ATTCORR______\")\n",
    "        \n",
    "        \"\"\"pia_harrison = wrl.atten.correctAttenuationHB(\n",
    "            ref_gr,\n",
    "            coefficients = dict(a=4.57e-5, b=0.731, gate_length=1.0),\n",
    "            mode=\"warn\",\n",
    "            thrs=59.)\n",
    "        pia_harrison[pia_harrison > 4.8] = 4.8\"\"\"\n",
    "\n",
    "        print (\"________ATTCORR2______\")\n",
    "        \"\"\"ref_gr = ref_gr + pia_harrison\"\"\"\n",
    "\n",
    "\n",
    "        print (\"________CLUTTER______\")\n",
    "        rho_th  = 0.85\n",
    "        ref_gr[rho_gr<= rho_th] = np.nan\n",
    "\n",
    "        print (\"________TH______\")\n",
    "        TH = 15\n",
    "        ref_gr[ref_gr<=TH]=np.nan\n",
    "        \n",
    "        # wichtige ZDR Werte nehmen\n",
    "        zdr_gr[zdr_gr > 5.]=np.nan\n",
    "        zdr_gr[zdr_gr < -5.]=np.nan\n",
    "\n",
    "        print ('----ref_gr: ',np.unique(ref_gr))\n",
    "\n",
    "        coord = wrl.georef.sweep_centroids(nray_gr, dr_gr, ngate_gr, elev_gr)\n",
    "        coords = wrl.georef.spherical_to_proj(coord[..., 0],\n",
    "                                              np.degrees(coord[..., 1]),\n",
    "                                              coord[..., 2],\n",
    "                                              (lon0_gr, lat0_gr, alt0_gr))\n",
    "        lon = coords[..., 0]\n",
    "        lat = coords[..., 1]\n",
    "        alt = coords[..., 2]\n",
    "        bbox = wrl.zonalstats.get_bbox(lon, lat)\n",
    "        print(\"Radar bounding box:\\n\\t%.2f\\n%.2f           %.2f\\n\\t%.2f\" %\n",
    "              (bbox['top'], bbox['left'], bbox['right'], bbox['bottom']))\n",
    "        sr_data = wrl.io.read_gpm(gpm_file, bbox)\n",
    "\n",
    "        # Longitudes of SR scans\n",
    "        sr_lon = sr_data['lon']\n",
    "        # Latitudes of SR scans\n",
    "        sr_lat = sr_data['lat']\n",
    "        # Precip flag\n",
    "        pflag = sr_data['pflag']\n",
    "        # Number of scans on SR data\n",
    "        nscan_sr= sr_data['nscan']\n",
    "        # Number of rays in one SR scan\n",
    "        nray_sr = sr_data['nray']\n",
    "        # Number of gates in one SR ray\n",
    "        ngate_sr = sr_data['nbin']\n",
    "\n",
    "        # Calculate equivalent earth radius\n",
    "        wgs84 = wrl.georef.get_default_projection()\n",
    "        re1 = wrl.georef.get_earth_radius(lat0_gr, wgs84)\n",
    "        #print(\"Earth radius 1:\", re1)\n",
    "        a = wgs84.GetSemiMajor()\n",
    "        b = wgs84.GetSemiMinor()\n",
    "        #print(\"SemiMajor, SemiMinor:\", a, b)\n",
    "\n",
    "        # Set up aeqd-projection gr-centered\n",
    "        rad = wrl.georef.proj4_to_osr(('+proj=aeqd +lon_0={lon:f} ' +\n",
    "                                           '+lat_0={lat:f} +a={a:f} ' +\n",
    "                                           '+b={b:f}').format(lon=lon0_gr,\n",
    "                                                              lat=lat0_gr,\n",
    "                                                              a=a, b=b))\n",
    "        re2 = wrl.georef.get_earth_radius(lat0_gr, rad)\n",
    "        #print(\"Earth radius 2:\", re2)\n",
    "\n",
    "        # create gr range and azimuth arrays\n",
    "        rmax_gr = r0_gr + ngate_gr * dr_gr\n",
    "        r_gr = np.arange(0, ngate_gr) * dr_gr + dr_gr/2.\n",
    "        az_gr = np.arange(0, nray_gr) - a0_gr\n",
    "        #print(\"Range/Azi-Shape:\", r_gr.shape, az_gr.shape)\n",
    "\n",
    "        # create gr polar grid and calculate aeqd-xyz coordinates\n",
    "        gr_polargrid = np.meshgrid(r_gr, az_gr)\n",
    "        gr_xyz, rad = wrl.georef.spherical_to_xyz(gr_polargrid[0], gr_polargrid[1], elev_gr, (lon0_gr, lat0_gr, alt0_gr ))\n",
    "        #print(\"XYZ-Grid-Shape:\", gr_xyz.shape)\n",
    "\n",
    "        # create gr poygon array in aeqd-xyz-coordinates\n",
    "        gr_poly, rad1 = wrl.georef.spherical_to_polyvert(r_gr, az_gr, elev_gr, (lon0_gr, lat0_gr, alt0_gr))\n",
    "        #print(gr_poly.shape, 360 * 600)\n",
    "        gr_poly.shape = (nray_gr, ngate_gr, 5, 3)\n",
    "\n",
    "        # get radar domain (outer ring)\n",
    "        gr_domain = gr_xyz[:,-1,0:2]\n",
    "        gr_domain = np.vstack((gr_domain, gr_domain[0]))\n",
    "        #print(\"Domain-Shape:\", gr_domain.shape)\n",
    "\n",
    "        sr_x, sr_y = wrl.georef.reproject(sr_lon, sr_lat,\n",
    "                                          projection_source=wgs84,\n",
    "                                          projection_target=rad)\n",
    "        sr_xy = np.dstack((sr_x, sr_y))\n",
    "        #print(\"SR-GRID-Shapes:\", sr_x.shape, sr_y.shape, sr_xy.shape)\n",
    "\n",
    "        precip_mask = (pflag >= 2) & wrl.zonalstats.get_clip_mask(sr_xy, gr_domain, rad)\n",
    "\n",
    "        #print(\"NRAY\", nray_sr)\n",
    "        #print(\"NBIN\", ngate_sr)\n",
    "\n",
    "        # use localZenith Angle\n",
    "        alpha = sr_data['zenith']\n",
    "        beta = abs(-17.04 + np.arange(nray_sr) * bw_sr)\n",
    "\n",
    "        # Correct for parallax, get 3D-XYZ-Array\n",
    "        #   xyzp_sr: Parallax corrected xyz coordinates\n",
    "        #   r_sr_inv: range array from ground to SR platform\n",
    "        #   zp: SR bin altitudes\n",
    "        xyp_sr, r_sr_inv, z_sr = wrl.georef.correct_parallax(sr_xy, ngate_sr, dr_sr, alpha)\n",
    "        print(xyp_sr.shape, r_sr_inv.shape, z_sr.shape)\n",
    "        xyzp_sr = np.concatenate((xyp_sr, z_sr[..., np.newaxis]),\n",
    "                           axis=-1)\n",
    "        r_sr, az_sr, elev_sr = wrl.georef.xyz_to_spherical(xyzp_sr, alt0_gr, proj=rad)\n",
    "        #TODO: hardcoded 1.0\n",
    "        mask = (elev_sr > (1.0 - bw_gr/2.)) & (elev_sr < (1.0 + bw_gr/2.))\n",
    "\n",
    "        rs = wrl.georef.dist_from_orbit(zt, alpha, beta, r_sr_inv, re1)\n",
    "\n",
    "        # Small anngle approximation\n",
    "        vol_sr2  = np.pi * dr_sr * rs**2 * np.radians(bw_sr / 2.)**2\n",
    "\n",
    "        # Or using wradlib's native function\n",
    "        vol_sr = wrl.qual.pulse_volume(rs, dr_sr, bw_sr)\n",
    "\n",
    "        # GR pulse volumes\n",
    "        #   along one beam\n",
    "        vol_gr = wrl.qual.pulse_volume(r_gr, dr_gr, bw_gr)\n",
    "        #   with shape (nray_gr, ngate_gr)\n",
    "        vol_gr = np.repeat(vol_gr, nray_gr).reshape((nray_gr, ngate_gr), order=\"F\")\n",
    "\n",
    "        Rs = 0.5 * (1 +  np.cos(np.radians(alpha)))[:,:,np.newaxis] * rs * np.tan(np.radians(bw_sr/2.))\n",
    "        Ds = dr_sr / np.cos(np.radians(alpha))\n",
    "        Ds = np.broadcast_to(Ds[..., np.newaxis], Rs.shape)\n",
    "\n",
    "        ratio, ibb = wrl.qual.get_bb_ratio(sr_data['zbb'], sr_data['bbwidth'], sr_data['quality'], z_sr)\n",
    "\n",
    "        zbb = sr_data['zbb'].copy()\n",
    "        zbb[~ibb] = np.nan\n",
    "\n",
    "        ref_sr = sr_data['refl'].filled(np.nan)\n",
    "\n",
    "\n",
    "        # First assumption: no valid SR bins (all False)\n",
    "        valid = np.asarray(elev_sr, dtype=np.bool)==False\n",
    "        print(valid.shape, precip_mask.shape)\n",
    "\n",
    "        # SR is inside GR range and is precipitating\n",
    "        iscan = precip_mask.nonzero()[0]\n",
    "        iray = precip_mask.nonzero()[1]\n",
    "        valid[iscan,iray] = True\n",
    "\n",
    "        # SR bins intersect with GR sweep\n",
    "        valid = valid & (elev_sr >= (elev_gr-bw_gr/2.)) & (elev_sr <= (elev_gr+bw_gr/2.))\n",
    "\n",
    "        # Number of matching SR bins per profile\n",
    "        nvalids = np.sum(valid, axis=2)\n",
    "\n",
    "        # scan and ray indices for profiles with at least one valid bin\n",
    "        vscan, vray = np.where(nvalids>0)\n",
    "\n",
    "        # number of profiles with at least one valid bin\n",
    "        nprof = len(vscan)\n",
    "\n",
    "        # average coordinates\n",
    "        xyz_v1 = xyzp_sr.copy()\n",
    "        print(xyz_v1.shape)\n",
    "        xyz_v1[~valid] = np.nan\n",
    "        xyz_c1 = xyzp_sr.copy()\n",
    "        xyz_c1[~valid] = 0\n",
    "        c = np.count_nonzero(xyz_c1[..., 0], axis=2)\n",
    "        ntotsr = c[vscan, vray]\n",
    "        xyz_m1 = np.nanmean(xyz_v1,axis=2)\n",
    "        xyz = xyz_m1[vscan, vray]\n",
    "        print(xyz.shape, c.shape)\n",
    "\n",
    "        # approximate Rs\n",
    "        rs_v1 = Rs.copy()\n",
    "        rs_v1[~valid] = np.nan\n",
    "        rs_m1 = np.nanmax(rs_v1, axis=2)\n",
    "        rs_prof = rs_m1[vscan, vray]\n",
    "        ds = rs_prof\n",
    "\n",
    "        # approximate Ds\n",
    "        ds_v1 = Ds.copy()\n",
    "        ds_v1[~valid] = np.nan\n",
    "        ds_m1 = np.nansum(ds_v1, axis=2)\n",
    "        ds_prof = ds_m1[vscan, vray]\n",
    "        dz = ds_prof\n",
    "\n",
    "        # approximate Vs\n",
    "        vs_v1 = vol_sr.copy()\n",
    "        vs_v1[~valid] = np.nan\n",
    "        vs_m1 = np.nansum(vs_v1, axis=2)\n",
    "        vs_prof = vs_m1[vscan, vray]\n",
    "        volsr1 = vs_prof\n",
    "\n",
    "        r_sr, az_sr, el_rs = wrl.georef.xyz_to_spherical(xyz, alt0_gr, proj=rad)\n",
    "\n",
    "        ref_sr_1 = wrl.trafo.idecibel(ref_sr)\n",
    "        ref_sr_1[~valid] = np.nan\n",
    "        refsr1a = np.nanmean(ref_sr_1, axis=2)[vscan,vray]\n",
    "        refsr1a = wrl.trafo.decibel(refsr1a)\n",
    "\n",
    "\n",
    "        zds = wrl.zonalstats.DataSource(xyz[:, 0:2].reshape(-1, 2), rad)\n",
    "\n",
    "\n",
    "        tmp_trg_lyr = zds.ds.GetLayer()\n",
    "        trg_poly = []\n",
    "        for i, feat in enumerate(tmp_trg_lyr):\n",
    "            geom = feat.GetGeometryRef().Buffer(rs_prof[i])\n",
    "            poly = wrl.georef.ogr_to_numpy(geom)\n",
    "            trg_poly.append(poly)    \n",
    "\n",
    "        zdp = wrl.zonalstats.ZonalDataPoly(gr_poly[..., 0:2].reshape(-1, 5, 2), trg_poly, srs=rad)\n",
    "        zdp.dump_vector('m3d_zonal_poly_{0}'.format(platf))\n",
    "\n",
    "        obj3 = wrl.zonalstats.GridCellsToPoly('m3d_zonal_poly_{0}'.format(platf))\n",
    "\n",
    "        print(obj3.ix.shape)\n",
    "        volgr1 = np.array([np.sum(vol_gr.ravel()[obj3.ix[i]])\n",
    "                           for i in np.arange(len(obj3.ix))[~obj3.check_empty()]])\n",
    "        print(volgr1.shape)\n",
    "\n",
    "        ref_gr_i = wrl.trafo.idecibel(ref_gr.ravel())\n",
    "        #ref_gr_ku_i = wrl.trafo.idecibel(ref_gr_ku.ravel())\n",
    "        refgr1a = np.array([np.nanmean(ref_gr_i[obj3.ix[i]])\n",
    "                     for i in np.arange(len(obj3.ix))[~obj3.check_empty()]])\n",
    "        refgr1a = wrl.trafo.decibel(refgr1a)\n",
    "        \n",
    "        \n",
    "        #zdr\n",
    "        zdr_gr_i = wrl.trafo.idecibel(zdr_gr.ravel())\n",
    "        #ref_gr_ku_i = wrl.trafo.idecibel(ref_gr_ku.ravel())\n",
    "        zdrgr1a = np.array([np.nanmedian(zdr_gr_i[obj3.ix[i]])\n",
    "                     for i in np.arange(len(obj3.ix))[~obj3.check_empty()]])\n",
    "        zdrgr1a = wrl.trafo.decibel(zdrgr1a)\n",
    "        \n",
    "\n",
    "        print ('Durchgang: ', ijk)\n",
    "\n",
    "        pl.figure(figsize=(18,10))\n",
    "        pl.subplot(1,3,1)\n",
    "        pl.scatter(xyz[...,0], xyz[...,1], c=refsr1a, label=str(elev_gr),vmin=0, vmax=40)\n",
    "        pl.colorbar()\n",
    "        pl.subplot(1,3,2)\n",
    "        pl.scatter(xyz[...,0], xyz[...,1], c=refgr1a, label=str(elev_gr),vmin=0, vmax=40)\n",
    "        pl.colorbar()\n",
    "        pl.subplot(1,3,3)\n",
    "        pl.scatter(xyz[...,0], xyz[...,1], c=zdrgr1a, label=str(elev_gr))\n",
    "        pl.colorbar()\n",
    "        \n",
    "        #pl.colorbar() \n",
    "        #pl.savefig('/automount/ftp/velibor/validation/match_ele'+str(elev_gr)+'.png')\n",
    "        #pl.close()\n",
    "        pl.show()\n",
    "        \n",
    "        X = np.append(X,xyz[...,0])\n",
    "        Y = np.append(Y,xyz[...,1])\n",
    "        Z = np.append(Z,xyz[...,2])\n",
    "\n",
    "        SR = np.append(SR, refsr1a)\n",
    "        GR = np.append(GR, refgr1a)\n",
    "        GR_zdr = np.append(GR_zdr, zdrgr1a)\n",
    "\n",
    "        gr_nr_elev = np.append(gr_nr_elev, refgr1a.shape)\n",
    "        gr_ele = np.append(gr_ele, elev_gr)\n",
    "\n",
    "\n",
    "\n",
    "    \n",
    "\n",
    "    \n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from satlib import corcor\n",
    "from pcc import get_miub_cmap\n",
    "TH = 15\n",
    "SR[SR<=TH]=np.nan\n",
    "GR[GR<=TH]=np.nan\n",
    "\n",
    "\n",
    "fig = pl.figure(figsize=(14,14))\n",
    "ax = fig.add_subplot(221, aspect=\"equal\")\n",
    "pl.scatter(GR, SR, marker=\"+\", c=\"black\")\n",
    "pl.plot([10,50],[10,50], linestyle=\"solid\", color=\"black\")\n",
    "#pl.xlim(10,50)\n",
    "#pl.ylim(10,50)\n",
    "pl.grid()\n",
    "pl.xlabel(\"GR reflectivity (dBZ)\")\n",
    "pl.ylabel(\"SR reflectivity (dBZ)\")\n",
    "pl.title(corcor(GR, SR))\n",
    "\n",
    "ax = fig.add_subplot(222)\n",
    "m0 = ~np.isnan(GR)& ~np.isnan(SR)\n",
    "\n",
    "pl.hist(GR[m0], bins=20, edgecolor=\"blue\", facecolor=\"grey\", label=\"GR\")\n",
    "pl.hist(SR[m0], bins=20, edgecolor=\"red\", facecolor=\"None\", label=\"SR\", lw=3)\n",
    "pl.xlabel(\"Reflectivity (dBZ)\")\n",
    "pl.legend()\n",
    "pl.grid()\n",
    "pl.title('Matched SR and GR Ref')\n",
    "pl.xlim(10,50)\n",
    "\n",
    "ax = fig.add_subplot(223)\n",
    "pl.hist2d(GR[m0], SR[m0],bins=20)\n",
    "\n",
    "pl.xlabel(\"GR reflectivity (dBZ)\")\n",
    "pl.ylabel(\"SR reflectivity (dBZ)\")\n",
    "pl.legend()\n",
    "pl.grid()\n",
    "pl.colorbar()\n",
    "#pl.xlim(10,50)\n",
    "\n",
    "\n",
    "\n",
    "pl.tight_layout()\n",
    "pl.savefig('/automount/ftp/velibor/validation/matching1.png')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "fig = pl.figure(figsize=(18,6))\n",
    "ax1 = fig.add_subplot(131, projection='3d')\n",
    "pppl = ax1.scatter(X, Y, Z, c=SR.ravel(), vmin=0, vmax=40, cmap=get_miub_cmap())\n",
    "fig.colorbar(pppl, shrink=0.5, aspect=5)\n",
    "ax1.set_title('Ref SR')\n",
    "\n",
    "ax2 = fig.add_subplot(132, projection='3d')\n",
    "pppl = ax2.scatter(X, Y, Z, c=GR.ravel(), vmin=0, vmax=40, cmap=get_miub_cmap())\n",
    "fig.colorbar(pppl, shrink=0.5, aspect=5)\n",
    "ax2.set_title('Ref GR')\n",
    "\n",
    "ax3 = fig.add_subplot(133, projection='3d')\n",
    "pppl = ax3.scatter(X, Y, Z, c=SR.ravel()-GR.ravel(),\n",
    "                   cmap='coolwarm',vmin=-10,vmax=10)\n",
    "fig.colorbar(pppl, shrink=0.5, aspect=5)\n",
    "ax3.set_title('diffRef SR-GR')\n",
    "\n",
    "from satlib import corcor\n",
    "print ('Corellation: ',corcor( SR, GR))\n",
    "pl.savefig('/automount/ftp/velibor/validation/matching2.png')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "fig = pl.figure(figsize=(12,8))\n",
    "ax = fig.add_subplot(131, aspect=\"auto\")\n",
    "pl.scatter(Y,Z, c=SR,  vmin=0, vmax=40, edgecolor=\"None\", cmap=get_miub_cmap())\n",
    "pl.title(\"SR reflectivity\")\n",
    "pl.xlim(np.nanmin(Y),np.nanmax(Y))\n",
    "pl.ylim(np.nanmin(Z),8000)\n",
    "pl.grid()\n",
    "pl.colorbar(shrink=0.5)\n",
    "\n",
    "ax = fig.add_subplot(132, aspect=\"auto\")\n",
    "pl.scatter(Y,Z, c=GR,  vmin=0, vmax=40, edgecolor=\"None\", cmap=get_miub_cmap())\n",
    "pl.title(\"GR reflectivity\")\n",
    "pl.xlim(np.nanmin(Y),np.nanmax(Y))\n",
    "pl.ylim(np.nanmin(Z),8000)\n",
    "pl.grid()\n",
    "pl.colorbar(shrink=0.5)\n",
    "#fig.suptitle(\"uncorrected GR vs uncorrected SR\")\n",
    "pl.tight_layout()\n",
    "\n",
    "ax = fig.add_subplot(133, aspect=\"auto\")\n",
    "pl.scatter(Y,Z, c=SR-GR, cmap='jet',vmin=-10,vmax=10 ,edgecolor=\"None\")\n",
    "pl.title(\"SR - GR reflectivity\")\n",
    "pl.xlim(np.nanmin(Y),np.nanmax(Y))\n",
    "pl.ylim(np.nanmin(Z),8000)\n",
    "pl.grid()\n",
    "pl.colorbar(shrink=0.5)\n",
    "#fig.suptitle(\"uncorrected SR - GR\")\n",
    "pl.tight_layout()\n",
    "pl.savefig('/automount/ftp/velibor/validation/matching3.png')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "fig = pl.figure(figsize=(12,8))\n",
    "ax = fig.add_subplot(131, aspect=\"auto\")\n",
    "pl.scatter(X,Z, c=SR, cmap=get_miub_cmap(), vmin=0, vmax=40, edgecolor=\"None\")\n",
    "pl.title(\"SR reflectivity\")\n",
    "pl.xlim(np.nanmin(X),np.nanmax(X))\n",
    "pl.ylim(np.nanmin(Z),8000)\n",
    "pl.grid()\n",
    "pl.colorbar(shrink=0.5)\n",
    "\n",
    "ax = fig.add_subplot(132, aspect=\"auto\")\n",
    "pl.scatter(X,Z, c=GR, cmap=get_miub_cmap(), vmin=0, vmax=40, edgecolor=\"None\")\n",
    "pl.title(\"GR reflectivity\")\n",
    "pl.xlim(np.nanmin(X),np.nanmax(X))\n",
    "pl.ylim(np.nanmin(Z),8000)\n",
    "pl.grid()\n",
    "pl.colorbar(shrink=0.5)\n",
    "#fig.suptitle(\"uncorrected GR vs uncorrected SR\")\n",
    "pl.tight_layout()\n",
    "\n",
    "ax = fig.add_subplot(133, aspect=\"auto\")\n",
    "pl.scatter(X,Z, c=SR-GR, cmap='jet',vmin=-10,vmax=10 ,edgecolor=\"None\")\n",
    "pl.title(\"SR - GR reflectivity\")\n",
    "pl.xlim(np.nanmin(X),np.nanmax(X))\n",
    "pl.ylim(np.nanmin(Z),8000)\n",
    "pl.grid()\n",
    "pl.colorbar(shrink=0.5)\n",
    "#fig.suptitle(\"uncorrected SR - GR\")\n",
    "pl.tight_layout()\n",
    "pl.savefig('/automount/ftp/velibor/validation/matching4.png')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "fig = pl.figure(figsize=(12,8))\n",
    "ax = fig.add_subplot(131, aspect=\"equal\")\n",
    "pl.scatter(X,Y, c=SR, cmap=get_miub_cmap(), vmin=0, vmax=40, edgecolor=\"None\",s=50)\n",
    "pl.title(\"SR reflectivity\")\n",
    "pl.xlim(np.nanmin(X),np.nanmax(X))\n",
    "pl.ylim(np.nanmin(Y),np.nanmax(Y))\n",
    "pl.grid()\n",
    "pl.colorbar(shrink=0.5)\n",
    "\n",
    "ax = fig.add_subplot(132, aspect=\"equal\")\n",
    "pl.scatter(X,Y, c=GR, cmap=get_miub_cmap(), vmin=0, vmax=40, edgecolor=\"None\",s=50)\n",
    "pl.title(\"GR reflectivity\")\n",
    "pl.xlim(np.nanmin(X),np.nanmax(X))\n",
    "pl.ylim(np.nanmin(Y),np.nanmax(Y))\n",
    "pl.grid()\n",
    "pl.colorbar(shrink=0.5)\n",
    "#fig.suptitle(\"uncorrected GR vs uncorrected SR\")\n",
    "pl.tight_layout()\n",
    "\n",
    "ax = fig.add_subplot(133, aspect=\"equal\")\n",
    "pl.scatter(X,Y, c=SR-GR, cmap='jet',vmin=-5,vmax=5 ,edgecolor=\"None\",s=50)\n",
    "pl.title(\"SR - GR reflectivity\")\n",
    "pl.xlim(np.nanmin(X),np.nanmax(X))\n",
    "pl.ylim(np.nanmin(Y),np.nanmax(Y))\n",
    "pl.grid()\n",
    "pl.colorbar(shrink=0.5)\n",
    "#fig.suptitle(\"uncorrected SR - GR\")\n",
    "pl.tight_layout()\n",
    "pl.savefig('/automount/ftp/velibor/validation/matching5.png')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def plot_cfad(x1, x2, ZP):\n",
    "    \"\"\"\n",
    "    Plot 4 CFAD\n",
    "    \"\"\"\n",
    "    import matplotlib.pyplot as plt\n",
    "    m1 = ~np.isnan(x1) & ~np.isnan(x2)\n",
    "    \n",
    "    x3 = np.array([]); x4 = np.array([]); x5 = np.array([])\n",
    "\n",
    "    hh = np.array([])\n",
    "    for jj in range(0,20000,100):\n",
    "        hh = np.append(hh,jj)\n",
    "        x3 = np.append(x3,np.nanmean(x1[np.where((x2>jj)&(x2<jj+100))]))\n",
    "        x4 = np.append(x4,np.nanmax(x1[np.where((x2>jj)&(x2<jj+100))]))\n",
    "    \n",
    "    plt.figure(figsize=(12,12)) \n",
    "    plt.hexbin(x1[m1], x2[m1],bins=10, cmap='inferno')\n",
    "    plt.colorbar()\n",
    "    plt.xlabel('ZH in DBz')\n",
    "    plt.ylabel('Height in m')\n",
    "    plt.title('CFAD: '+ ZP)\n",
    "    plt.plot(x3,hh, linewidth=2, label='mean')\n",
    "    plt.plot(x4,hh, linewidth=2, label='median')\n",
    "    plt.ylim(0,10000)\n",
    "    plt.xlim(0,60)\n",
    "    plt.legend(loc='upper right', fontsize=10)\n",
    "    plt.grid(color='white')\n",
    "    #plt.close()\n",
    "    #plt.ylim(0,)\n",
    "    #plt.show()\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plot_cfad(SR,Z,'SR')\n",
    "pl.savefig('/automount/ftp/velibor/validation/cfadSR.png')\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plot_cfad(GR,Z,'GR')\n",
    "pl.savefig('/automount/ftp/velibor/validation/CFADGR.png')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "hhh=0\n",
    "\n",
    "h = np.arange(0,15000,100)\n",
    "cor = np.array([])\n",
    "for ccc in h:\n",
    "    cor = np.append(cor,corcor(SR[Z>ccc],GR[Z>ccc]))\n",
    "\n",
    "                    \n",
    "#       \n",
    "pl.figure(figsize=(14,14))\n",
    "pl.subplot(2,2,1)\n",
    "pl.scatter(SR[Z>hhh],GR[Z>hhh],c=Z[Z>hhh],alpha=0.5)\n",
    "pl.xlabel('SR in dbz')\n",
    "pl.ylabel('GR in dbz')\n",
    "pl.title(corcor(SR[Z>hhh],GR[Z>hhh]))\n",
    "cb = pl.colorbar()\n",
    "cb.set_label('h in m')\n",
    "pl.grid()\n",
    "            \n",
    "pl.subplot(2,2,2)\n",
    "pl.plot(h,cor)\n",
    "pl.grid()\n",
    "pl.ylim(0,1)\n",
    "pl.xlabel('h in m')\n",
    "pl.ylabel('correlation')\n",
    "\n",
    "pl.savefig('/automount/ftp/velibor/validation/vertical_corr_Z.png')\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "\"\"\"aa=0.9\n",
    "\n",
    "pl.figure(figsize=(16,10))\n",
    "pl.subplot(2,5,1)\n",
    "pl.scatter(SR[0:629],GR[0:629],alpha=aa,label='ele: 0.10',color='black')\n",
    "pl.legend(loc='upper left')\n",
    "pl.xlabel('SR in dbz')\n",
    "pl.ylabel('GR in dbz')\n",
    "pl.title(corcor(SR[0:629],GR[0:629]))\n",
    "pl.grid()\n",
    "pl.xlim(0,45)\n",
    "pl.ylim(0,45)\n",
    "pl.plot([0,45],[0,45], color='black', ls='--')\n",
    "\n",
    "pl.subplot(2,5,5)\n",
    "pl.scatter(SR[629:1254],GR[629:1254],alpha=aa,label='ele: 4.50',color='blue')\n",
    "pl.legend(loc='upper left')\n",
    "pl.xlabel('SR in dbz')\n",
    "pl.ylabel('GR in dbz')\n",
    "pl.title(corcor(SR[629:1254],GR[629:1254]))\n",
    "pl.grid()\n",
    "pl.xlim(0,45)\n",
    "pl.ylim(0,45)\n",
    "pl.plot([0,45],[0,45], color='black', ls='--')\n",
    "\n",
    "pl.subplot(2,5,6)\n",
    "pl.scatter(SR[1254:2379],GR[1254:2379],alpha=aa,label='ele: 8.20',color='green')\n",
    "pl.legend(loc='upper left')\n",
    "pl.xlabel('SR in dbz')\n",
    "pl.ylabel('GR in dbz')\n",
    "pl.title(corcor(SR[1254:2379],GR[1254:2379]))\n",
    "pl.grid()\n",
    "pl.xlim(0,45)\n",
    "pl.ylim(0,45)\n",
    "pl.plot([0,45],[0,45], color='black', ls='--')\n",
    "\n",
    "pl.subplot(2,5,7)\n",
    "pl.scatter(SR[2379:2676],GR[2379:2676],alpha=aa,label='ele: 11.0',color='tan')\n",
    "pl.legend(loc='upper left')\n",
    "pl.xlabel('SR in dbz')\n",
    "pl.ylabel('GR in dbz')\n",
    "pl.title(corcor(SR[2379:2676],GR[2379:2676]))\n",
    "pl.grid()\n",
    "pl.xlim(0,45)\n",
    "pl.ylim(0,45)\n",
    "pl.plot([0,45],[0,45], color='black', ls='--')\n",
    "\n",
    "pl.subplot(2,5,8)\n",
    "pl.scatter(SR[2676:2861],GR[2676:2861],alpha=aa,label='ele: 14.0',color='red')\n",
    "pl.legend(loc='upper left')\n",
    "pl.xlabel('SR in dbz')\n",
    "pl.ylabel('GR in dbz')\n",
    "pl.title(corcor(SR[2676:2861],GR[2676:2861]))\n",
    "pl.grid()\n",
    "pl.xlim(0,45)\n",
    "pl.ylim(0,45)\n",
    "pl.plot([0,45],[0,45], color='black', ls='--')\n",
    "\n",
    "pl.subplot(2,5,9)\n",
    "pl.scatter(SR[2861:2944],GR[2861:2944],alpha=aa,label='ele: 18.0',color='gray')\n",
    "pl.legend(loc='upper left')\n",
    "pl.xlabel('SR in dbz')\n",
    "pl.ylabel('GR in dbz')\n",
    "pl.title(corcor(SR[2861:2944],GR[2861:2944]))\n",
    "pl.grid()\n",
    "pl.xlim(0,45)\n",
    "pl.ylim(0,45)\n",
    "pl.plot([0,45],[0,45], color='black', ls='--')\n",
    "\n",
    "pl.subplot(2,5,10)\n",
    "pl.scatter(SR[2944:2950],GR[2944:2950],alpha=aa,label='ele: 28.0',color='orange')\n",
    "pl.legend(loc='upper left')\n",
    "pl.xlabel('SR in dbz')\n",
    "pl.ylabel('GR in dbz')\n",
    "pl.title(corcor(SR[2944:2950],GR[2944:2950]))\n",
    "pl.grid()\n",
    "pl.xlim(0,45)\n",
    "pl.ylim(0,45)\n",
    "pl.plot([0,45],[0,45], color='black', ls='--')\n",
    "\n",
    "pl.subplot(2,5,2)\n",
    "pl.scatter(SR[2950:3252],GR[2950:3252],alpha=aa,label='ele: 1.5',color='brown')\n",
    "pl.legend(loc='upper left')\n",
    "pl.xlabel('SR in dbz')\n",
    "pl.ylabel('GR in dbz')\n",
    "pl.title(corcor(SR[2950:3252],GR[2950:3252]))\n",
    "pl.grid()\n",
    "pl.xlim(0,45)\n",
    "pl.ylim(0,45)\n",
    "pl.plot([0,45],[0,45], color='black', ls='--')\n",
    "\n",
    "pl.subplot(2,5,3)\n",
    "pl.scatter(SR[3252:3554],GR[3252:3554],alpha=aa,label='ele: 2.5',color='cyan')\n",
    "pl.legend(loc='upper left')\n",
    "pl.xlabel('SR in dbz')\n",
    "pl.ylabel('GR in dbz')\n",
    "pl.title(corcor(SR[3252:3554],GR[3252:3554]))\n",
    "pl.xlim(0,45)\n",
    "pl.ylim(0,45)\n",
    "pl.plot([0,45],[0,45], color='black', ls='--')\n",
    "pl.grid()\n",
    "\n",
    "pl.subplot(2,5,4)\n",
    "pl.scatter(SR[3554:3854],GR[3554:3854],alpha=aa,label='ele: 3.4',color='magenta')\n",
    "pl.legend(loc='upper left')\n",
    "pl.xlabel('SR in dbz')\n",
    "pl.ylabel('GR in dbz')\n",
    "pl.title(corcor(SR[3554:3854],GR[3554:3854]))\n",
    "pl.xlim(0,45)\n",
    "pl.ylim(0,45)\n",
    "pl.plot([0,45],[0,45], color='black', ls='--')\n",
    "\n",
    "pl.grid()\n",
    "pl.savefig('/automount/ftp/velibor/validation/corr_elev.png')\"\"\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "\"\"\"aa=0.9\n",
    "from pcc import get_miub_cmap\n",
    "pl.figure(figsize=(30,10))\n",
    "pl.subplot(2,5,1)\n",
    "pl.scatter(X[0:629],Y[0:629],c=GR[0:629], cmap=get_miub_cmap(), vmin=0, vmax=40)\n",
    "pl.grid()\n",
    "pl.colorbar()\n",
    "pl.xlim(-150000,150000)\n",
    "pl.ylim(-150000,150000)\n",
    "\n",
    "\n",
    "pl.subplot(2,5,5)\n",
    "pl.scatter(X[629:1254],Y[629:1254],c=GR[629:1254], cmap=get_miub_cmap(), vmin=0, vmax=40)\n",
    "pl.grid()\n",
    "pl.colorbar()\n",
    "pl.xlim(-150000,150000)\n",
    "pl.ylim(-150000,150000)\n",
    "\n",
    "pl.subplot(2,5,6)\n",
    "pl.scatter(X[1254:2379],Y[1254:2379],c=GR[1254:2379], cmap=get_miub_cmap(), vmin=0, vmax=40)\n",
    "pl.grid()\n",
    "pl.colorbar()\n",
    "pl.xlim(-150000,150000)\n",
    "pl.ylim(-150000,150000)\n",
    "\n",
    "pl.subplot(2,5,7)\n",
    "pl.scatter(X[2379:2676],Y[2379:2676],c=GR[2379:2676], cmap=get_miub_cmap(), vmin=0, vmax=40)\n",
    "pl.grid()\n",
    "pl.colorbar()\n",
    "pl.xlim(-150000,150000)\n",
    "pl.ylim(-150000,150000)\n",
    "\n",
    "pl.subplot(2,5,8)\n",
    "pl.scatter(X[2676:2861],Y[2676:2861],c=GR[2676:2861], cmap=get_miub_cmap(), vmin=0, vmax=40)\n",
    "pl.grid()\n",
    "pl.colorbar()\n",
    "pl.xlim(-150000,150000)\n",
    "pl.ylim(-150000,150000)\n",
    "\n",
    "pl.subplot(2,5,9)\n",
    "pl.scatter(X[2861:2944],Y[2861:2944],c=GR[2861:2944], cmap=get_miub_cmap(), vmin=0, vmax=40)\n",
    "pl.grid()\n",
    "pl.colorbar()\n",
    "pl.xlim(-150000,150000)\n",
    "pl.ylim(-150000,150000)\n",
    "\n",
    "pl.subplot(2,5,10)\n",
    "pl.scatter(X[2944:2950],Y[2944:2950],c=GR[2944:2950], cmap=get_miub_cmap(), vmin=0, vmax=40)\n",
    "pl.grid()\n",
    "pl.colorbar()\n",
    "pl.xlim(-150000,150000)\n",
    "pl.ylim(-150000,150000)\n",
    "\n",
    "pl.subplot(2,5,2)\n",
    "pl.scatter(X[2950:3252],Y[2950:3252],c=GR[2950:3252], cmap=get_miub_cmap(), vmin=0, vmax=40)\n",
    "pl.grid()\n",
    "pl.colorbar()\n",
    "pl.xlim(-150000,150000)\n",
    "pl.ylim(-150000,150000)\n",
    "\n",
    "\n",
    "pl.subplot(2,5,3)\n",
    "pl.scatter(X[3252:3554],Y[3252:3554],c=GR[3252:3554], cmap=get_miub_cmap(), vmin=0, vmax=40)\n",
    "pl.grid()\n",
    "pl.colorbar()\n",
    "pl.xlim(-150000,150000)\n",
    "pl.ylim(-150000,150000)\n",
    "\n",
    "\n",
    "pl.subplot(2,5,4)\n",
    "pl.scatter(X[3554:3854],Y[3554:3854],c=GR[3554:3854], cmap=get_miub_cmap(), vmin=0, vmax=40)\n",
    "pl.colorbar()\n",
    "pl.xlim(-150000,150000)\n",
    "pl.ylim(-150000,150000)\n",
    "pl.grid()\n",
    "pl.tight_layout()\n",
    "pl.savefig('/automount/ftp/velibor/validation/corr_elev_ppi.png')\"\"\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "mm = ~np.isnan(GR)& ~np.isnan(SR)\n",
    "pl.figure(figsize=(14,10))\n",
    "\n",
    "pl.hist(GR[mm]-SR[mm],bins=100, facecolor=\"grey\", alpha=0.6)\n",
    "pl.axvline(np.mean(GR[mm]-SR[mm]), color='blue', label='Mean diff')\n",
    "pl.axvline(np.median(GR[mm]-SR[mm]), color='green', label='Median diff')\n",
    "pl.axvline(0, color='black',ls='--',lw=3)\n",
    "\n",
    "xmin, xmax = np.nanmin(GR[mm]-SR[mm]), np.nanmax(GR[mm]-SR[mm])\n",
    "limits = np.nanmax([abs(xmin),abs(xmax)])\n",
    "pl.xlim(-limits,limits)\n",
    "\n",
    "pl.title('Diff GR-SR, bias: '+str( np.sum(GR[mm]-SR[mm])/len(GR[mm]-SR[mm]))+' dbz')\n",
    "pl.grid()\n",
    "pl.legend()\n",
    "pl.xlabel('SR overestimate----------------------------Diff in dbz----------------------------GR overestimate')\n",
    "pl.ylabel('number')\n",
    "pl.savefig('/automount/ftp/velibor/validation/matching_bias.png')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "\n",
    "m1 = ~np.isnan(GR)\n",
    "m2 = ~np.isnan(SR)\n",
    "\n",
    "pl.figure(figsize=(14,5))\n",
    "pl.hist(SR[mm],bins=50, facecolor=\"None\", alpha=1, edgecolor='black', histtype='step', label='SR')\n",
    "pl.xlim(left=0)\n",
    "\n",
    "\n",
    "pl.hist(GR[mm],bins=50, facecolor=\"None\", alpha=1, edgecolor='blue', histtype='step', label='GR')\n",
    "pl.axvline(0, color='black',ls='--',lw=3)\n",
    "pl.grid()\n",
    "pl.legend(fontsize=20)\n",
    "#pl.ylim(top=90)\n",
    "pl.ylabel('number', fontsize=20)\n",
    "pl.xlabel(r'$Z_H$ in dbz', fontsize=20)\n",
    "pl.xticks(fontsize=20); pl.yticks(fontsize=20)\n",
    "pl.savefig('/automount/ftp/velibor/validation/matching_hist.png')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "aa=0.9\n",
    "# Elevation unterteilen mit Cumulativer Summe\n",
    "e_pos = gr_nr_elev.cumsum().astype(int)\n",
    "\n",
    "pl.figure(figsize=(16,10))\n",
    "\n",
    "for eee in range(len(gr_nr_elev)):\n",
    "    pl.subplot(2,5,eee+1)\n",
    "    \n",
    "    if eee==0: \n",
    "        pl.scatter(SR[0:e_pos[0]],GR[0:e_pos[0]],alpha=aa,label='elevation: '+str(gr_ele[0])[0:5],color='black',s=10)\n",
    "        pl.title('r: ' + corcor(SR[0:e_pos[0]],GR[0:e_pos[0]])[0:4])\n",
    "        pl.legend(loc='upper left')\n",
    "    else:\n",
    "        pl.scatter(SR[e_pos[eee-1]:e_pos[eee]],GR[e_pos[eee-1]:e_pos[eee]],alpha=aa,\n",
    "                   label='elevation: ' + str(gr_ele[eee])[0:5],color='black',s=10)\n",
    "        pl.title('r: ' + corcor(SR[e_pos[eee-1]:e_pos[eee]],GR[e_pos[eee-1]:e_pos[eee]])[0:6])\n",
    "        pl.legend(loc='upper left')\n",
    "    pl.xlabel('SR in dbz')\n",
    "    pl.ylabel('GR in dbz')\n",
    "    pl.grid()\n",
    "    pl.xlim(0,45)\n",
    "    pl.ylim(0,45)\n",
    "    pl.plot([0,45],[0,45], color='black', ls='--')\n",
    "\n",
    "pl.tight_layout()\n",
    "pl.savefig('/automount/ftp/velibor/validation/corr_elev.png')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Elevation unterteilen mit Cumulativer Summe\n",
    "e_pos = gr_nr_elev.cumsum().astype(int)\n",
    "from pcc import get_miub_cmap\n",
    "\n",
    "pl.figure(figsize=(30,10))\n",
    "\n",
    "for eee in range(len(gr_nr_elev)):\n",
    "    pl.subplot(2,5,eee+1)\n",
    "    \n",
    "    if eee==0: \n",
    "        pl.scatter(X[0:e_pos[0]],Y[0:e_pos[0]],c=GR[0:e_pos[0]],label=str(gr_ele[0]),\n",
    "                   cmap=get_miub_cmap(), vmin=0, vmax=40)\n",
    "        pl.legend(loc='upper left')\n",
    "    else:\n",
    "        pl.scatter(X[e_pos[eee-1]:e_pos[eee]],Y[e_pos[eee-1]:e_pos[eee]],c=GR[e_pos[eee-1]:e_pos[eee]],\n",
    "                   label=str(gr_ele[eee]), cmap=get_miub_cmap(), vmin=0, vmax=40)\n",
    "        pl.legend(loc='upper left')\n",
    "    pl.xlabel('SR in dbz')\n",
    "    pl.ylabel('GR in dbz')\n",
    "    pl.grid()\n",
    "    pl.xlim(-150000,150000)\n",
    "    pl.ylim(-150000,150000)\n",
    "    pl.colorbar()\n",
    "\n",
    "\n",
    "pl.tight_layout()\n",
    "pl.savefig('/automount/ftp/velibor/validation/corr_ppiGR.png')\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Elevation unterteilen mit Cumulativer Summe\n",
    "e_pos = gr_nr_elev.cumsum().astype(int)\n",
    "from pcc import get_miub_cmap\n",
    "\n",
    "pl.figure(figsize=(30,10))\n",
    "\n",
    "for eee in range(len(gr_nr_elev)):\n",
    "    pl.subplot(2,5,eee+1)\n",
    "    \n",
    "    if eee==0: \n",
    "        pl.scatter(X[0:e_pos[0]],Y[0:e_pos[0]],c=SR[0:e_pos[0]],label=str(gr_ele[0]),\n",
    "                   cmap=get_miub_cmap(), vmin=0, vmax=40)\n",
    "        pl.legend(loc='upper left')\n",
    "    else:\n",
    "        pl.scatter(X[e_pos[eee-1]:e_pos[eee]],Y[e_pos[eee-1]:e_pos[eee]],c=SR[e_pos[eee-1]:e_pos[eee]],\n",
    "                   label=str(gr_ele[eee]), cmap=get_miub_cmap(), vmin=0, vmax=40)\n",
    "        pl.legend(loc='upper left')\n",
    "    pl.xlabel('SR in dbz')\n",
    "    pl.ylabel('GR in dbz')\n",
    "    pl.grid()\n",
    "    pl.xlim(-150000,150000)\n",
    "    pl.ylim(-150000,150000)\n",
    "    pl.colorbar()\n",
    "\n",
    "\n",
    "pl.tight_layout()\n",
    "pl.savefig('/automount/ftp/velibor/validation/corr_ppiSR.png')\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Elevation unterteilen mit Cumulativer Summe\n",
    "e_pos = gr_nr_elev.cumsum().astype(int)\n",
    "from pcc import get_miub_cmap\n",
    "\n",
    "pl.figure(figsize=(30,10))\n",
    "\n",
    "for eee in range(len(gr_nr_elev)):\n",
    "    pl.subplot(2,5,eee+1)\n",
    "    \n",
    "    if eee==0: \n",
    "        pl.scatter(X[0:e_pos[0]],Z[0:e_pos[0]],c=SR[0:e_pos[0]],label=str(gr_ele[0]),\n",
    "                   cmap=get_miub_cmap(), vmin=0, vmax=40)\n",
    "        pl.legend(loc='upper left')\n",
    "    else:\n",
    "        pl.scatter(X[e_pos[eee-1]:e_pos[eee]],Z[e_pos[eee-1]:e_pos[eee]],c=SR[e_pos[eee-1]:e_pos[eee]],\n",
    "                   label=str(gr_ele[eee]), cmap=get_miub_cmap(), vmin=0, vmax=40)\n",
    "        pl.legend(loc='upper left')\n",
    "    pl.xlabel('SR in dbz')\n",
    "    pl.ylabel('GR in dbz')\n",
    "    pl.grid()\n",
    "    pl.xlim(-150000,150000)\n",
    "    pl.ylim(0,7000)\n",
    "    pl.colorbar()\n",
    "\n",
    "\n",
    "pl.tight_layout()\n",
    "pl.savefig('/automount/ftp/velibor/validation/corr_XppiSR.png')\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Elevation unterteilen mit Cumulativer Summe\n",
    "e_pos = gr_nr_elev.cumsum().astype(int)\n",
    "from pcc import get_miub_cmap\n",
    "\n",
    "pl.figure(figsize=(30,10))\n",
    "\n",
    "for eee in range(len(gr_nr_elev)):\n",
    "    pl.subplot(2,5,eee+1)\n",
    "    \n",
    "    if eee==0: \n",
    "        pl.scatter(X[0:e_pos[0]],Z[0:e_pos[0]],c=GR[0:e_pos[0]],label=str(gr_ele[0]),\n",
    "                   cmap=get_miub_cmap(), vmin=0, vmax=40)\n",
    "        pl.legend(loc='upper left')\n",
    "    else:\n",
    "        pl.scatter(X[e_pos[eee-1]:e_pos[eee]],Z[e_pos[eee-1]:e_pos[eee]],c=GR[e_pos[eee-1]:e_pos[eee]],\n",
    "                   label=str(gr_ele[eee]), cmap=get_miub_cmap(), vmin=0, vmax=40)\n",
    "        pl.legend(loc='upper left')\n",
    "    pl.xlabel('SR in dbz')\n",
    "    pl.ylabel('GR in dbz')\n",
    "    pl.grid()\n",
    "    pl.xlim(-150000,150000)\n",
    "    pl.ylim(0,7000)\n",
    "    pl.colorbar()\n",
    "\n",
    "\n",
    "pl.tight_layout()\n",
    "pl.savefig('/automount/ftp/velibor/validation/corr_XppiGR.png')\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "np.nanmin(ref_sr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Elevation unterteilen mit Cumulativer Summe\n",
    "e_pos = gr_nr_elev.cumsum().astype(int)\n",
    "from pcc import get_miub_cmap\n",
    "\n",
    "pl.figure(figsize=(30,10))\n",
    "\n",
    "for eee in range(len(gr_nr_elev)):\n",
    "    pl.subplot(2,5,eee+1)\n",
    "    \n",
    "    if eee==0: \n",
    "        pl.scatter(X[0:e_pos[0]],Y[0:e_pos[0]],c=SR[0:e_pos[0]]-GR[0:e_pos[0]],label=str(gr_ele[0]),\n",
    "                   cmap='jet', vmin=-10, vmax=10)\n",
    "        pl.legend(loc='upper left')\n",
    "    else:\n",
    "        pl.scatter(X[e_pos[eee-1]:e_pos[eee]],Y[e_pos[eee-1]:e_pos[eee]],c=SR[e_pos[eee-1]:e_pos[eee]]-GR[e_pos[eee-1]:e_pos[eee]],\n",
    "                   label=str(gr_ele[eee]), cmap='jet', vmin=-10, vmax=10)\n",
    "        pl.legend(loc='upper left')\n",
    "    pl.xlabel('X')\n",
    "    pl.ylabel('Y')\n",
    "    pl.title('SR-GR')\n",
    "    pl.grid()\n",
    "    pl.xlim(-150000,150000)\n",
    "    pl.ylim(-150000,150000)\n",
    "    pl.colorbar()\n",
    "\n",
    "\n",
    "pl.tight_layout()\n",
    "pl.savefig('/automount/ftp/velibor/validation/corr_ppiSR.png')\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Elevation unterteilen mit Cumulativer Summe\n",
    "e_pos = gr_nr_elev.cumsum().astype(int)\n",
    "from pcc import get_miub_cmap\n",
    "\n",
    "pl.figure(figsize=(30,10))\n",
    "\n",
    "for eee in range(len(gr_nr_elev)):\n",
    "    pl.subplot(2,5,eee+1)\n",
    "    \n",
    "    if eee==0: \n",
    "        pl.scatter(X[0:e_pos[0]],Z[0:e_pos[0]],c=SR[0:e_pos[0]]-GR[0:e_pos[0]],label=str(gr_ele[0]),\n",
    "                   cmap='jet', vmin=-10, vmax=10)\n",
    "        pl.legend(loc='upper left')\n",
    "    else:\n",
    "        pl.scatter(X[e_pos[eee-1]:e_pos[eee]],Z[e_pos[eee-1]:e_pos[eee]],c=SR[e_pos[eee-1]:e_pos[eee]]-GR[e_pos[eee-1]:e_pos[eee]],\n",
    "                   label=str(gr_ele[eee]), cmap='jet', vmin=-10, vmax=10)\n",
    "        pl.legend(loc='upper left')\n",
    "    pl.xlabel('X')\n",
    "    pl.ylabel('Z')\n",
    "    pl.grid()\n",
    "    pl.xlim(-150000,150000)\n",
    "    pl.ylim(0,7000)\n",
    "    pl.colorbar()\n",
    "\n",
    "\n",
    "pl.tight_layout()\n",
    "pl.savefig('/automount/ftp/velibor/validation/corr_XppiSR.png')\n",
    "\n"
   ]
  },
  {
   "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": [
    "SRi, GRi = SR.copy(), GR.copy()\n",
    "\n",
    "SRi[SRi>36]=np.nan\n",
    "SRi[SRi<24]=np.nan\n",
    "GRi[GRi>36]=np.nan\n",
    "GRi[GRi<24]=np.nan\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "fig = pl.figure(figsize=(14,14))\n",
    "ax = fig.add_subplot(221, aspect=\"equal\")\n",
    "pl.scatter(GRi, SRi, marker=\"+\", c=\"black\")\n",
    "pl.plot([10,50],[10,50], linestyle=\"solid\", color=\"black\")\n",
    "#pl.xlim(10,50)\n",
    "#pl.ylim(10,50)\n",
    "pl.grid()\n",
    "pl.xlabel(\"GR reflectivity (dBZ)\")\n",
    "pl.ylabel(\"SR reflectivity (dBZ)\")\n",
    "pl.title(corcor(GRi, SRi))\n",
    "\n",
    "ax = fig.add_subplot(222)\n",
    "m0i = ~np.isnan(GRi)& ~np.isnan(SRi)\n",
    "\n",
    "pl.hist(GRi[m0i], bins=20, edgecolor=\"blue\", facecolor=\"None\", label=\"GR\", histtype='step')\n",
    "pl.hist(SRi[m0i], bins=20, edgecolor=\"red\", facecolor=\"None\", label=\"SR\", lw=3, histtype='step')\n",
    "pl.xlabel(\"Reflectivity (dBZ)\")\n",
    "pl.legend()\n",
    "pl.grid()\n",
    "pl.title('Matched SR and GR Ref')\n",
    "\n",
    "ax = fig.add_subplot(223)\n",
    "pl.hist2d(GRi[m0i], SRi[m0i],bins=20)\n",
    "\n",
    "pl.xlabel(\"GR reflectivity (dBZ)\")\n",
    "pl.ylabel(\"SR reflectivity (dBZ)\")\n",
    "pl.legend()\n",
    "pl.grid()\n",
    "pl.colorbar()\n",
    "#pl.xlim(10,50)\n",
    "\n",
    "\n",
    "\n",
    "ax = fig.add_subplot(224)\n",
    "\n",
    "mmi = ~np.isnan(GRi)& ~np.isnan(SRi)\n",
    "\n",
    "pl.hist(SRi[mmi]-GRi[mmi],bins=50, facecolor=\"grey\", alpha=0.6)\n",
    "pl.axvline(np.mean(SRi[mmi]-GRi[mmi]), color='blue', label='Mean diff')\n",
    "pl.axvline(np.median(SRi[mmi]-GRi[mmi]), color='green', label='Median diff')\n",
    "pl.axvline(0, color='black',ls='--',lw=3)\n",
    "\n",
    "xmin, xmax = np.nanmin(SRi[mmi]-GRi[mmi]), np.nanmax(SRi[mmi]-GRi[mmi])\n",
    "limits = np.nanmax([abs(xmin),abs(xmax)])\n",
    "pl.xlim(-limits,limits)\n",
    "\n",
    "pl.title('Diff SR-GR, bias: '+str( np.sum(SRi[mmi]-GRi[mmi])/len(SRi[mmi]-GRi[mmi]))+' dbz')\n",
    "pl.grid()\n",
    "pl.legend()\n",
    "pl.xlabel('GR overestimate----------------------------Diff in dbz----------------------------SR overestimate')\n",
    "pl.ylabel('number')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "print ('von',mmi.shape, 'Pixel matschen sich ', sum(mmi==True))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "print sum(~np.isnan(GR)), sum(~np.isnan(SR))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "150./600.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "\n",
    "pl.plot(sr_data['zbb']-sr_data['bbwidth'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def subplot_cfad(x1, x2, ZP):\n",
    "    \"\"\"\n",
    "    Plot 4 CFAD\n",
    "    \"\"\"\n",
    "    import matplotlib.pyplot as plt\n",
    "    m1 = ~np.isnan(x1) & ~np.isnan(x2)\n",
    "    \n",
    "    x3 = np.array([]); x4 = np.array([]); x5 = np.array([])\n",
    "\n",
    "    hh = np.array([])\n",
    "    for jj in range(0,20000,100):\n",
    "        hh = np.append(hh,jj)\n",
    "        x3 = np.append(x3,np.nanmean(x1[np.where((x2>jj)&(x2<jj+100))]))\n",
    "        x4 = np.append(x4,np.nanmedian(x1[np.where((x2>jj)&(x2<jj+100))]))\n",
    "    \n",
    "    plt.hist2d(x1[m1], x2[m1],bins=30, cmap='inferno')\n",
    "    #plt.hexbin(x1[m1], x2[m1],bins=10, cmap='inferno')\n",
    "    plt.colorbar()\n",
    "    plt.xlabel('ZH in DBz')\n",
    "    plt.ylabel('Height in m')\n",
    "    plt.title('CFAD: '+ ZP)\n",
    "    plt.plot(x3,hh, linewidth=2, label='mean')\n",
    "    plt.plot(x4,hh, linewidth=2, label='median')\n",
    "\n",
    "    plt.legend(loc='upper right', fontsize=10)\n",
    "    plt.grid(color='white')\n",
    "    #plt.close()\n",
    "    #plt.ylim(0,)\n",
    "    #plt.show()\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "diff_SRGR = SR - GR\n",
    "pl.figure(figsize=(12,6)) \n",
    "pl.subplot(1,3,1)\n",
    "subplot_cfad(SR,Z,'SR')\n",
    "#pl.xlim(15,45)\n",
    "\n",
    "pl.subplot(1,3,2)\n",
    "subplot_cfad(GR,Z,'GR')\n",
    "#pl.xlim(15,45)\n",
    "\n",
    "pl.subplot(1,3,3)\n",
    "subplot_cfad(diff_SRGR,Z,'SR - GR')\n",
    "pl.axvline(0, color='white')\n",
    "#pl.xlim(-15,20)\n",
    "pl.tight_layout()\n",
    "pl.savefig('/automount/ftp/velibor/validation/cfadSR-GR.pdf')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "pl.figure(figsize=(12,6)) \n",
    "pl.subplot(1,3,1)\n",
    "subplot_cfad(sr_data['refl'].reshape(28* 49* 176),xyzp_sr[:,:,:,2].reshape(28* 49* 176),'SR_real')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "print xyz_v1[:,:,:,2].shape\n",
    "print sr_data['refl'].shape\n",
    "print xyzp_sr[:,:,:,2].shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "A = np.nanmean(sr_data['refl'], axis=(0,1))\n",
    "h = np.arange(0,176,1)*250.\n",
    "h = h-3000.\n",
    "print A.shape, h.shape\n",
    "\n",
    "pl.plot(A,h)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "A\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "# Distanz Berechnen\n",
    "x1, x2 = 0, X\n",
    "y1, y2 = 0, Y\n",
    "\n",
    "dist = np.hypot(x2 - x1, y2 - y1)\n",
    "\n",
    "#plt.pcolormesh(dpr_lon, dpr_lat, Z_boxpol_ipoli)\n",
    "pl.scatter(X, Y, c=dist)\n",
    "pl.colorbar()\n",
    "\n",
    "\n",
    "pl.figure(figsize=(18,16))\n",
    "sc = pl.scatter(SR,GR,c=Z, marker='o', s=3**(dist/10000), cmap='spectral')\n",
    "sc.set_facecolor(\"none\")\n",
    "pl.plot([15,50],[15,50])\n",
    "pl.colorbar()\n",
    "pl.grid()\n",
    "pl.xlim(15,50); pl.ylim(15,50)\n",
    "pl.xlabel('SR in dbz', fontsize=18); pl.ylabel('GR in dbz', fontsize=18)\n",
    "#pl.savefig('/automount/ftp/velibor/validation/007elev.pdf')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "print abs(np.arange(0,15,1)-15)\n",
    "print np.arange(0,15,1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "# Distanz Berechnen\n",
    "x1, x2 = 0, X\n",
    "y1, y2 = 0, Y\n",
    "\n",
    "dist = np.hypot(x2 - x1, y2 - y1)\n",
    "\n",
    "#plt.pcolormesh(dpr_lon, dpr_lat, Z_boxpol_ipoli)\n",
    "pl.scatter(X, Y, c=dist)\n",
    "pl.colorbar()\n",
    "\n",
    "\n",
    "pl.figure(figsize=(18,16))\n",
    "sc = pl.scatter(SR,GR,c=GR_zdr, marker='o', s=dist/1000, cmap='jet')\n",
    "sc.set_facecolor(\"none\")\n",
    "pl.plot([15,50],[15,50])\n",
    "pl.colorbar()\n",
    "pl.grid()\n",
    "pl.xlim(15,50); pl.ylim(15,50)\n",
    "pl.xlabel('SR in dbz', fontsize=18); pl.ylabel('GR in dbz', fontsize=18)\n",
    "pl.savefig('/automount/ftp/velibor/validation/007elev.pdf')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "m = ~np.isnan(GR_zdr)\n",
    "pl.hist(GR_zdr[m])\n",
    "pl.yscale('log')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "1e5 * 1e4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "(1./5e9) * ((-8.*5e4) - (4./np.cos(np.deg2rad(30))*1e5 - 10.*5e4 + 4.*1e5/np.cos(np.deg2rad(30))))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "(1./5e9) * ((-9.*5e4) - (4./np.cos(np.deg2rad(30))*1e5 - 9.*5e4 + 4.*1e5/np.cos(np.deg2rad(30))))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "0.000094 *3600"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "#bounds2 = np.array([-3,0,1,1.5,2.5,3,3.5,4,4.5])\n",
    "\n",
    "print (bounds2)\n",
    "\n",
    "# Distanz Berechnen\n",
    "x1, x2 = 0, X\n",
    "y1, y2 = 0, Y\n",
    "\n",
    "dist = np.hypot(x2 - x1, y2 - y1)\n",
    "\n",
    "v1, v2 =0, 4000 \n",
    "bounds2 = np.linspace(np.nanmin(Z), np.nanmax(Z), 20)\n",
    "bounds2 = np.arange(v1,v2, 1000)\n",
    "\n",
    "#plt.pcolormesh(dpr_lon, dpr_lat, Z_boxpol_ipoli)\n",
    "pl.scatter(X, Y, c=dist)\n",
    "pl.colorbar()\n",
    "#pl.colorbar(shrink=0.8,extend='both',ticks=bounds2, boundaries=bounds2)\n",
    "\n",
    "pl.figure(figsize=(18,16))\n",
    "sc = pl.scatter(SR,GR,c=Z, marker='o', s=dist/1000, cmap='cool',vmin = 0, vmax = 4000)\n",
    "sc.set_facecolor(\"none\")\n",
    "pl.plot([15,50],[15,50])\n",
    "pl.colorbar(shrink=0.8,extend='both')\n",
    "pl.grid()\n",
    "pl.xlim(15,50); pl.ylim(15,50)\n",
    "pl.xlabel('SR in dbz', fontsize=18); pl.ylabel('GR in dbz', fontsize=18)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "#bounds2 = np.array([-3,0,1,1.5,2.5,3,3.5,4,4.5])\n",
    "\n",
    "print (bounds2)\n",
    "\n",
    "# Distanz Berechnen\n",
    "x1, x2 = 0, X\n",
    "y1, y2 = 0, Y\n",
    "\n",
    "dist = np.hypot(x2 - x1, y2 - y1)\n",
    "\n",
    "v1, v2 =0, 4000 \n",
    "bounds2 = np.linspace(np.nanmin(Z), np.nanmax(Z), 20)\n",
    "bounds2 = np.arange(v1,v2, 1000)\n",
    "\n",
    "#plt.pcolormesh(dpr_lon, dpr_lat, Z_boxpol_ipoli)\n",
    "pl.scatter(X, Y, c=dist)\n",
    "pl.colorbar()\n",
    "#pl.colorbar(shrink=0.8,extend='both',ticks=bounds2, boundaries=bounds2)\n",
    "\n",
    "pl.figure(figsize=(18,16))\n",
    "sc = pl.scatter(SR,GR,c=GR_zdr, marker='o', s=Z/100, cmap='jet', vmin=np.nanmin(GR_zdr), vmax=3)\n",
    "sc.set_facecolor(\"none\")\n",
    "pl.plot([15,50],[15,50])\n",
    "pl.colorbar(shrink=0.8,extend='both')\n",
    "pl.grid()\n",
    "pl.xlim(15,50); pl.ylim(15,50)\n",
    "pl.xlabel('SR in dbz', fontsize=18); pl.ylabel('GR in dbz', fontsize=18)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "\n",
    "pl.figure(figsize=(18,16))\n",
    "sc = pl.scatter(SR,GR,c=GR_zdr, marker='o', s=80, cmap='jet')\n",
    "pl.colorbar()\n",
    "\n",
    "#pl.colorbar(shrink=0.8,extend='both',ticks=bounds2, boundaries=bounds2)\n",
    "\n",
    "pl.grid()\n",
    "\n",
    "pl.xlabel('SR'); pl.ylabel('GR')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "aa=0.9\n",
    "# Elevation unterteilen mit Cumulativer Summe\n",
    "e_pos = gr_nr_elev.cumsum().astype(int)\n",
    "\n",
    "pl.figure(figsize=(20,10))\n",
    "\n",
    "for eee in range(len(gr_nr_elev)):\n",
    "    pl.subplot(2,5,eee+1)\n",
    "    \n",
    "    if eee==0: \n",
    "        pl.scatter(SR[0:e_pos[0]],GR[0:e_pos[0]],c=Z[0:e_pos[0]],alpha=aa,\n",
    "                   label='elevation: '+str(gr_ele[0])[0:5],color='black',s=10, vmin=0, vmax=4000)\n",
    "        pl.title('r: ' + corcor(SR[0:e_pos[0]],GR[0:e_pos[0]])[0:4])\n",
    "        pl.legend(loc='upper left')\n",
    "        pl.colorbar()\n",
    "    else:\n",
    "        pl.scatter(SR[e_pos[eee-1]:e_pos[eee]],GR[e_pos[eee-1]:e_pos[eee]],c=Z[e_pos[eee-1]:e_pos[eee]],alpha=aa,\n",
    "                   label='elevation: ' + str(gr_ele[eee])[0:5],color='black',s=10, vmin=0, vmax=4000)\n",
    "        pl.title('r: ' + corcor(SR[e_pos[eee-1]:e_pos[eee]],GR[e_pos[eee-1]:e_pos[eee]])[0:6])\n",
    "        pl.legend(loc='upper left')\n",
    "        pl.colorbar()\n",
    "    pl.xlabel('SR in dbz')\n",
    "    pl.ylabel('GR in dbz')\n",
    "    pl.grid()\n",
    "    pl.xlim(0,45)\n",
    "    pl.ylim(0,45)\n",
    "    pl.plot([0,45],[0,45], color='black', ls='--')\n",
    "\n",
    "pl.tight_layout()\n",
    "pl.savefig('/automount/ftp/velibor/validation/ZDRcorr_elev.pdf')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "np.unique(GR_zdr)"
   ]
  },
  {
   "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": [
    "fig = pl.figure(figsize=(12,8))\n",
    "ax = fig.add_subplot(131, aspect=\"auto\")\n",
    "pl.scatter(Y,Z, c=SR,  vmin=0, vmax=40, edgecolor=\"None\", cmap=get_miub_cmap())\n",
    "pl.title(\"SR reflectivity\")\n",
    "pl.xlim(np.nanmin(Y),np.nanmax(Y))\n",
    "pl.ylim(np.nanmin(Z),8000)\n",
    "pl.grid()\n",
    "pl.colorbar(shrink=0.5)\n",
    "\n",
    "ax = fig.add_subplot(132, aspect=\"auto\")\n",
    "pl.scatter(Y,Z, c=GR,  vmin=0, vmax=40, edgecolor=\"None\", cmap=get_miub_cmap())\n",
    "pl.title(\"GR reflectivity\")\n",
    "pl.xlim(np.nanmin(Y),np.nanmax(Y))\n",
    "pl.ylim(np.nanmin(Z),8000)\n",
    "pl.grid()\n",
    "pl.colorbar(shrink=0.5)\n",
    "#fig.suptitle(\"uncorrected GR vs uncorrected SR\")\n",
    "pl.tight_layout()\n",
    "\n",
    "ax = fig.add_subplot(133, aspect=\"auto\")\n",
    "pl.scatter(Y,Z, c=GR_zdr, cmap='jet',vmin=-1,vmax=3 ,edgecolor=\"None\")\n",
    "pl.title(\"ZDR reflectivity\")\n",
    "pl.xlim(np.nanmin(Y),np.nanmax(Y))\n",
    "pl.ylim(np.nanmin(Z),8000)\n",
    "pl.grid()\n",
    "pl.colorbar(shrink=0.5)\n",
    "#fig.suptitle(\"uncorrected SR - GR\")\n",
    "pl.tight_layout()\n",
    "pl.savefig('/automount/ftp/velibor/validation/matching3.png')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "zbb.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "28*49"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "aa=0.9\n",
    "# Elevation unterteilen mit Cumulativer Summe\n",
    "e_pos = gr_nr_elev.cumsum().astype(int)\n",
    "\n",
    "pl.figure(figsize=(16,10))\n",
    "\n",
    "for eee in range(len(gr_nr_elev)):\n",
    "    pl.subplot(2,5,eee+1)\n",
    "    \n",
    "    if eee==0: \n",
    "        pl.scatter(SR[0:e_pos[0]],GR[0:e_pos[0]],c=GR_zdr[0:e_pos[0]],alpha=aa,label='elevation: '+str(gr_ele[0])[0:5],color='black',s=10, vmin=-0.2, vmax=2)\n",
    "        pl.title('r: ' + corcor(SR[0:e_pos[0]],GR[0:e_pos[0]])[0:4])\n",
    "        pl.legend(loc='upper left')\n",
    "        pl.colorbar()\n",
    "    else:\n",
    "        pl.scatter(SR[e_pos[eee-1]:e_pos[eee]],GR[e_pos[eee-1]:e_pos[eee]],c=GR_zdr[e_pos[eee-1]:e_pos[eee]],alpha=aa,\n",
    "                   label='elevation: ' + str(gr_ele[eee])[0:5],color='black',s=10, vmin=-0.2, vmax=2)\n",
    "        pl.title('r: ' + corcor(SR[e_pos[eee-1]:e_pos[eee]],GR[e_pos[eee-1]:e_pos[eee]])[0:6])\n",
    "        pl.legend(loc='upper left')\n",
    "        pl.colorbar()\n",
    "    pl.xlabel('SR in dbz')\n",
    "    pl.ylabel('GR in dbz')\n",
    "    pl.grid()\n",
    "    pl.xlim(0,45)\n",
    "    pl.ylim(0,45)\n",
    "    pl.plot([0,45],[0,45], color='black', ls='--')\n",
    "\n",
    "pl.tight_layout()\n",
    "##pl.savefig('/automount/ftp/velibor/validation/corr_elev.png')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "xyz_v1[:,:,:,2].shape\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "('r:', 49.627307582066095)\n",
      "('h:', 6.093467170257374)\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "print ('r:',np.cos(np.deg2rad(7.))*50000/1000.)\n",
    "print ('h:' ,np.sin(np.deg2rad(7.))*50000/1000.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "500./0.02"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "16.0"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "4000./250"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "27.04"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(26**2)/25."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "75000"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "75*1000"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "16.0"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "4000./250.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.0"
      ]
     },
     "execution_count": 38,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.cos(np.deg2rad(0.))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.0"
      ]
     },
     "execution_count": 42,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.cos(np.deg2rad(180))**2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "8"
      ]
     },
     "execution_count": 37,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data = np.array([[1., 2., 3., 4., np.nan],[1,2,3,4,np.nan]])\n",
    "np.count_nonzero(~np.isnan(data))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "6.123233995736766e-17"
      ]
     },
     "execution_count": 39,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.cos()"
   ]
  },
  {
   "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
}