_1_Paper_data_preparation_RY.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import h5py\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as pl\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",
    "\n",
    "% matplotlib inline \n",
    "\n",
    "#from pcc import get_time_of_gpm\n",
    "import datetime as dt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Used Functions\n",
    "-----------------"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def dpr_swath_contour(sr_grid):\n",
    "    \"\"\"\n",
    "\n",
    "    Function for dpr swath grid contour\n",
    "    \n",
    "    ## sr_grid[:, 0]        # - complete left column\n",
    "    ## sr_grid[-1,1:-1]     # - upper edge without left and right corner\n",
    "    ## sr_grid[:, -1][::-1] # - complete right column (backward)\n",
    "    ## sr_grid[0,0:-1]      # - lower edge (without right corner)\n",
    "    \n",
    "    Input:\n",
    "    ------\n",
    "    sr_grid ::: Space-borne grid x,y\n",
    "    \n",
    "    Output:\n",
    "    -------\n",
    "    trg ::: Contour of the space-borne grid\n",
    "    \n",
    "    \"\"\"\n",
    "    \n",
    "    trg = np.r_[sr_grid[:, 0], sr_grid[-1, 1:-1], sr_grid[:, -1][::-1], sr_grid[0, 0:-1][::-1]]\n",
    "    \n",
    "    return trg\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def idx_sr_in_gr_area(_bingrid, _gpm_xy):\n",
    "    \"\"\"\n",
    "    Funktion:\n",
    "    ---------\n",
    "    Search for SR Footprint (Index) that are located in the GR scan area \n",
    "    \n",
    "    Input: \n",
    "    ------\n",
    "    _bingrid ::: Binary-grid of the GR scan area\n",
    "        \n",
    "    _gpm_xy ::: GPM footptint coordinates\n",
    "    \n",
    "    Output:\n",
    "    -------\n",
    "    \n",
    "    _gpm_xy_outer_idx ::: index of all footprints not included in radolan scan area\n",
    "    _gpm_xy_inner_idx ::: idex of all footprints included in radolan area but\n",
    "                           not scanned by ground radars \n",
    "    \n",
    "    \"\"\"\n",
    "    _gpm_xy = _gpm_xy.copy()\n",
    "    # Bestimmen von Eckpunkten bei RADOLAN RY GRID!\n",
    "    # HARDCODED: but ok for regula RADOLAN Produkts\n",
    "    xmin, xmax = -523.4621669218559, 375.5378330781441\n",
    "    ymin, ymax = -4658.644724265572, -3759.644724265572\n",
    "    \n",
    "    # Radolangitter um 1 gridpoit zu allen seiten erweitern\n",
    "    _xex, _yex = np.arange(xmin-1, xmax+2, 1), np.arange(ymin-1, ymax+2, 1)\n",
    "    _xxx, _yyy = np.meshgrid(_xex,_yex)\n",
    "    \n",
    "    # Bin grid bestimmen\n",
    "    _bingrid = _bingrid.copy()\n",
    "    _bingrid[_bingrid >= 0]=1\n",
    "    _bingrid[_bingrid < 0]=0\n",
    "    \n",
    "    # Gittererweiterungen zusammenstellen\n",
    "    ## rn_tb ::: bin gird top and bottom\n",
    "    ## rn_lr ::: bin gird left and right\n",
    "    rn_tb = np.zeros(900)\n",
    "    rn_lr = np.zeros(902)\n",
    "    \n",
    "    # Griderweiterungen einsetzen\n",
    "    # erst oben und unten\n",
    "    rn1 = np.c_[rn_tb, _bingrid, rn_tb]\n",
    "    # dann rechts\n",
    "    rn2 = np.vstack((rn1, rn_lr))\n",
    "    # dann rechts\n",
    "    rn3 = np.vstack((rn_lr,rn2))\n",
    "    \n",
    "    from skimage import measure\n",
    "    # Contouren des Bin-grids erstellen\n",
    "    contours = measure.find_contours(rn3, 0, positive_orientation='high',\n",
    "                                     fully_connected='high' )\n",
    "    ### print(len(contours), ' detected Polygons')\n",
    "    #print(\"---------------------\")\n",
    "    #for j in range(len(contours)):\n",
    "    #    print(len(contours[j]))\n",
    "    #print(\"---------------------\")\n",
    "\n",
    "    # bestiimung der Listen für Polygon Points\n",
    "    _xx = []\n",
    "    _yy = []\n",
    "\n",
    "    # suche nach dem größten Polygon... dieser sollte der RADOLAN UMRISS sein\n",
    "    for i in range(len(contours)):\n",
    "        # Polygone mit 3 oder weniger points ergeben keine nutzbare Fläche\n",
    "        if len(_xxx[contours[i][:,0].astype(int),contours[i][:,1].astype(int)]) <= 3:\n",
    "            print('not relevant polygon removed')\n",
    "        else:   \n",
    "            # PolygonPoints bestimmen (RADOLAN RAND POINTS)\n",
    "            _xx.append(_xxx[contours[i][:,0].astype(int),contours[i][:,1].astype(int)])\n",
    "            _yy.append(_yyy[contours[i][:,0].astype(int),contours[i][:,1].astype(int)])\n",
    "            \n",
    "    # sortiern nach der länge, das letzt element ist das größte und sollte somit der Radolarand sein\n",
    "    _xx.sort(key=len)\n",
    "    _yy.sort(key=len)\n",
    "\n",
    "    # Outter polygon erstellen\n",
    "    _xy = np.vstack((_xx[-1].ravel(), _yy[-1].ravel())).transpose()    \n",
    "    \n",
    "    # Suche nach den GPM footprints im Scanngebiet von RADOLAN mit ZONALSTATS\n",
    "    zdpoly = wradlib.zonalstats.ZonalDataPoly(_gpm_xy, [_xy])\n",
    "    _gpm_xy_outer_idx = zdpoly.get_source_index(0)\n",
    "    \n",
    "    # Array for inner Polygons\n",
    "    _gpm_xy_inner_idx = np.array([])\n",
    "    \n",
    "    # Entferne Inner Polygons\n",
    "    # Auch nur wenn es weitere Polygone gibt \n",
    "    if len(_xx)>1:\n",
    "\n",
    "        for inner_poly_index in range(len(_xx)-1):\n",
    "\n",
    "            ### print ('Polygon size: ', len(_xx[inner_poly_index].ravel()))\n",
    "\n",
    "            _xy_inner = np.vstack((_xx[inner_poly_index].ravel(),\n",
    "                                  _yy[inner_poly_index].ravel())).transpose()\n",
    "\n",
    "            _zdpoly = wradlib.zonalstats.ZonalDataPoly(_gpm_xy, [_xy_inner])\n",
    "\n",
    "            _inner_idx = _zdpoly.get_source_index(0)\n",
    "\n",
    "            _gpm_xy_inner_idx = np.append(_gpm_xy_inner_idx, _inner_idx)\n",
    "\n",
    "\n",
    "            ### if _gpm_xy_inner_idx.size==0:\n",
    "                ### print ('inner polygons do not match with SR grid')\n",
    "\n",
    "            ### else:\n",
    "                ### print('match idx with inner polygon: ', _gpm_xy_inner_idx.shape)\n",
    "                    \n",
    "\n",
    "    \n",
    "    _idx_r = ~np.isin(_gpm_xy_outer_idx, _gpm_xy_inner_idx)\n",
    "    \n",
    "    return _gpm_xy_outer_idx[_idx_r]\n",
    "        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def dpr_antenna_weighting(r):\n",
    "    \"\"\"\n",
    "    Funktion: explained in \n",
    "    Watters, D., A. Battaglia, K. Mroz, and F. Tridon, 0: \n",
    "    Validation of the GPM Version-5 Surface Rainfall \n",
    "    Products over Great Britain and Ireland. \n",
    "    J. Hydrometeor., 0, https://doi.org/10.1175/JHM-D-18-0051.1\n",
    "    \n",
    "    and\n",
    "    \n",
    "    Mroz, K., A. Battaglia, T. J. Lang, D. J. Cecil, \n",
    "    S. Tanelli, and F. Tridon, 2017:  Hail-detection\n",
    "    algorithm for the GPM Core Observatory satellite sensors.\n",
    "    Journal of Applied Meteorology and Climatology,56 (7), 1939–1957\n",
    "    \n",
    "    Input:\n",
    "    ------\n",
    "    r ::: Distance to the center point of the footprint\n",
    "    \n",
    "    Output:\n",
    "    -------\n",
    "    omega ::: Weighting for the point with the distance\n",
    "    \n",
    "    \"\"\"\n",
    "    \n",
    "    omega = np.exp(-(r/2.5)**2. * np.log10(4.))\n",
    "    \n",
    "    return omega"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def ipoli_radi_w(xy_cut, gpm_xy,rwdata_cut, dpr_footprint, k=25):\n",
    "    \"\"\"\n",
    "    ...\n",
    "    \"\"\"\n",
    "    # cKDTree radolan\n",
    "    tree = cKDTree(xy_cut, balanced_tree=False)\n",
    "\n",
    "    # cKDTree gpm dpr\n",
    "    tree_gpm = cKDTree(gpm_xy, balanced_tree=True)\n",
    "\n",
    "    dists, ix = tree.query(gpm_xy, k=k, n_jobs=-1) # k maximal possible ry pixel in dpr footprint\n",
    "    ix2 = tree.query_ball_point(gpm_xy, dpr_footprint)\n",
    "\n",
    "    ry_pns_w = []\n",
    "\n",
    "    for i in range(ix.shape[0]):\n",
    "        #i is all points in one dpr footprint\n",
    "        #distancen for all i\n",
    "        index = np.isin(ix[i,:],ix2[i] )\n",
    "        res1 = np.nansum( dpr_antenna_weighting(dists[i,:][index]) * rwdata_cut.ravel()[ix[i,:][index]])\n",
    "        res2 = np.nansum(dpr_antenna_weighting(dists[i,:][index]))\n",
    "\n",
    "        ry_pns_w.append(res1/res2)\n",
    "\n",
    "\n",
    "    return np.array(ry_pns_w)   "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def ipoli_radi(xy_cut, gpm_xy,rwdata_cut, dpr_footprint, k=25, calc='mean'):\n",
    "    \"\"\"\n",
    "    calc: mean, max, min\n",
    "    \"\"\"\n",
    "    # cKDTree radolan\n",
    "    tree = cKDTree(xy_cut, balanced_tree=False)\n",
    "\n",
    "    # cKDTree gpm dpr\n",
    "    tree_gpm = cKDTree(gpm_xy, balanced_tree=True)\n",
    "\n",
    "    dists, ix = tree.query(gpm_xy, k=k, n_jobs=-1) # k maximal possible ry pixel in dpr footprint\n",
    "    ix2 = tree.query_ball_point(gpm_xy, dpr_footprint)\n",
    "\n",
    "    ry_par = []\n",
    "    \n",
    "    if calc=='mean':\n",
    "\n",
    "        for i in range(ix.shape[0]):\n",
    "            #i is all points in one dpr footprint\n",
    "            #distancen for all i\n",
    "            index = np.isin(ix[i,:],ix2[i] )\n",
    "            res1 = np.nanmean(rwdata_cut.ravel()[ix[i,:][index]])\n",
    "            ry_par.append(res1)\n",
    "            \n",
    "    elif calc=='max':\n",
    "\n",
    "            for i in range(ix.shape[0]):\n",
    "                #i is all points in one dpr footprint\n",
    "                #distancen for all i\n",
    "                index = np.isin(ix[i,:],ix2[i] )\n",
    "                res1 = np.nanmax(rwdata_cut.ravel()[ix[i,:][index]])\n",
    "                ry_par.append(res1)\n",
    "              \n",
    "    elif calc=='min':\n",
    "\n",
    "            for i in range(ix.shape[0]):\n",
    "                #i is all points in one dpr footprint\n",
    "                #distancen for all i\n",
    "                index = np.isin(ix[i,:],ix2[i] )\n",
    "                res1 = np.nanmin(rwdata_cut.ravel()[ix[i,:][index]])\n",
    "                ry_par.append(res1)\n",
    "    \n",
    "    elif calc=='std':\n",
    "\n",
    "            for i in range(ix.shape[0]):\n",
    "                #i is all points in one dpr footprint\n",
    "                #distancen for all i\n",
    "                index = np.isin(ix[i,:],ix2[i] )\n",
    "                res1 = np.nanstd(rwdata_cut.ravel()[ix[i,:][index]])\n",
    "                ry_par.append(res1)            \n",
    "                \n",
    "    else:\n",
    "        print ('Wrong calc parameter declaration!')\n",
    "\n",
    "\n",
    "    return np.array(ry_par)   "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def gpm_scan_time(filename, scanswath='NS'):\n",
    "    from netCDF4 import Dataset\n",
    "    pr_data = Dataset(filename, mode=\"r\")\n",
    "    year = pr_data[scanswath]['ScanTime'].variables['Year']\n",
    "    month = pr_data[scanswath]['ScanTime'].variables['Month']\n",
    "    dayofmonth = pr_data[scanswath]['ScanTime'].variables['DayOfMonth']\n",
    "    # dayofyear = pr_data['NS']['ScanTime'].variables['DayOfYear'][mask]\n",
    "    hour = pr_data[scanswath]['ScanTime'].variables['Hour']\n",
    "    minute = pr_data[scanswath]['ScanTime'].variables['Minute']\n",
    "    second = pr_data[scanswath]['ScanTime'].variables['Second']\n",
    "    # secondofday = pr_data['NS']['ScanTime'].variables['SecondOfDay'][mask]\n",
    "    millisecond = pr_data[scanswath]['ScanTime'].variables['MilliSecond']\n",
    "    date_array = zip(year, month, dayofmonth,\n",
    "                     hour, minute, second,\n",
    "                     millisecond)\n",
    "    pr_time = np.array(\n",
    "        [dt.datetime(d[0], d[1], d[2], d[3], d[4], d[5], d[6]) for d in date_array])\n",
    "    \n",
    "    #ttt = gpm_scan_time(gpm_file)\n",
    "    s_n = {'NS':49, 'HS':24, 'MS':25}\n",
    "    \n",
    "    time_array = np.array(s_n[scanswath] * [pr_time]).T\n",
    "    \n",
    "    return time_array"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def get_time_of_gpm2(gpm_time):\n",
    "    \"\"\"einfache mittlere zeit\"\"\"\n",
    "    iii = int(len(np.array(gpm_time['Year']))/2)\n",
    "    gpm_year = np.array(gpm_time['Year'])[iii]\n",
    "    gpm_month = np.array(gpm_time['Month'])[iii]\n",
    "    gpm_day = np.array(gpm_time['DayOfMonth'])[iii]\n",
    "    gpm_hour = np.array(gpm_time['Hour'])[iii]\n",
    "    gpm_min = np.array(gpm_time['Minute'])[iii]\n",
    "    gpm_sek = np.array(gpm_time['Second'])[iii]\n",
    "    gpm_dt = dt.datetime(gpm_year,gpm_month, gpm_day, gpm_hour, gpm_min, gpm_sek).strftime(\"%Y.%m.%d -- %H:%M:%S\")\n",
    "    return gpm_dt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def write2hdf(name, x, y, sr_pns, sr_pes, sr_pav, sr_bbh, sr_bbw, sr_phase, sr_type, sr_sth,\n",
    "             sr_cfb, sr_lst, sr_ele, sr_tim, gr_pp_ipoli, gr_xy, gr_pp, gr_h_max, gr_h_min, gr_h_mean, gr_h_std,\n",
    "             gr_beam_width, \n",
    "             sr_time, gr_time, gpm_lbn, gpm_ibn, gpm_mbn, sr_bbb, sr_bbt, sr_lza, sr_flh, sr_hip, sr_shr, sr_anv,\n",
    "              sr_flb, gr_dem_mean, gr_dem_std, gr_slope_std, gr_slope_mean):\n",
    "    \"\"\"\n",
    "    \n",
    "    Write data in HDF5 file\n",
    "    \n",
    "    \"\"\"\n",
    "    hdat = h5py.File(str(name) + '.hdf5', 'w')\n",
    "    \n",
    "    ## Create Groups SR and GR estimates\n",
    "    grp_sr = hdat.create_group(\"SR\")\n",
    "    grp_sr.attrs['description']=\" space-borne estimates (for more information see gpm ATBD)\"\n",
    "    grp_gr = hdat.create_group(\"GR\")\n",
    "    grp_gr.attrs['description']=\" ground-based estimates (for more information see RADOLAN Docs)\"\n",
    "    \n",
    "    ## SR\n",
    "    g_x=grp_sr.create_dataset('sr_x', data=x, compression=\"gzip\")\n",
    "    g_x.attrs['description'] = 'space-borne x-coordinates in km' \n",
    "    \n",
    "    g_y = grp_sr.create_dataset('sr_y', data=y, compression=\"gzip\")\n",
    "    g_y.attrs['description'] = 'space-borne y-coordinates in km'\n",
    "    \n",
    "    g_sr_pns = grp_sr.create_dataset('sr_pns', data=sr_pns, compression=\"gzip\")\n",
    "    g_sr_pns.attrs['description'] = 'precip near surface (mm/h)'\n",
    "    \n",
    "    g_sr_pes = grp_sr.create_dataset('sr_pes', data=sr_pes, compression=\"gzip\")\n",
    "    g_sr_pes.attrs['description'] = 'adjusted precip estimated observ. height (mm/h)'\n",
    "    \n",
    "    g_sr_pav = grp_sr.create_dataset('sr_pav', data=sr_pav, compression=\"gzip\")\n",
    "    g_sr_pav.attrs['description'] = 'precip mean between 2 and 4 km (mm/h)'\n",
    "    \n",
    "    g_sr_bbh = grp_sr.create_dataset('sr_bbh', data=sr_bbh, compression=\"gzip\")\n",
    "    g_sr_bbh.attrs['description'] = 'bright band height (m)'\n",
    "    \n",
    "    g_sr_bbw = grp_sr.create_dataset('sr_bbw', data=sr_bbw, compression=\"gzip\")\n",
    "    g_sr_bbw.attrs['description'] = 'bright band width (m)'\n",
    "    \n",
    "    g_sr_pha = grp_sr.create_dataset('sr_phase', data=sr_phase, compression=\"gzip\")\n",
    "    g_sr_pha.attrs['description'] = 'hydrometeor phase near surface //100 ->, 0=solid, 1=mixed, 2=liquid'\n",
    "    \n",
    "    g_sr_typ = grp_sr.create_dataset('sr_type', data=sr_type, compression=\"gzip\")\n",
    "    g_sr_typ.attrs['description'] = 'precip type, 1=strat, 2=conv, 3=other'\n",
    "    \n",
    "    g_sr_sth = grp_sr.create_dataset('sr_sth', data=sr_sth, compression=\"gzip\")\n",
    "    g_sr_sth.attrs['description'] = 'storm top hight (m)'\n",
    "    \n",
    "    g_sr_cfb = grp_sr.create_dataset('sr_cfb', data=sr_cfb, compression=\"gzip\")\n",
    "    g_sr_cfb.attrs['description'] = 'range bin number for clutter free bottom'\n",
    "    \n",
    "    g_sr_lst = grp_sr.create_dataset('sr_lst', data=sr_lst, compression=\"gzip\")\n",
    "    g_sr_lst.attrs['description'] = 'land surface type, O-99=Ocean, 100-199=Land,200-299=Coast, 300-399=InlandWater'\n",
    "    \n",
    "    g_sr_ele = grp_sr.create_dataset('sr_ele', data=sr_ele, compression=\"gzip\")\n",
    "    g_sr_ele.attrs['description'] = 'elevation of measurement point, Copy of DEMHmean (m)'\n",
    "    \n",
    "    g_sr_times = grp_sr.create_dataset('sr_times', data=sr_tim, compression=\"gzip\")\n",
    "    g_sr_times.attrs['description'] = 'scan times'\n",
    "    \n",
    "    g_sr_aph = grp_sr.create_dataset('sr_lbn', data=gpm_lbn, compression=\"gzip\")\n",
    "    g_sr_aph.attrs['description'] = 'number of liquid phase bins'\n",
    "    \n",
    "    g_sr_aph = grp_sr.create_dataset('sr_ibn', data=gpm_ibn, compression=\"gzip\")\n",
    "    g_sr_aph.attrs['description'] = 'number of iced phase bins'\n",
    "    \n",
    "    g_sr_aph = grp_sr.create_dataset('sr_mbn', data=gpm_mbn, compression=\"gzip\")\n",
    "    g_sr_aph.attrs['description'] = 'number of mixed phase bins'\n",
    "    \n",
    "    # new\n",
    "    g_sr_bbb = grp_sr.create_dataset('sr_bbb', data=sr_bbb, compression=\"gzip\")\n",
    "    g_sr_bbb.attrs['description'] = 'bright band bottom in m above NN'\n",
    "    \n",
    "    g_sr_bbt = grp_sr.create_dataset('sr_bbt', data=sr_bbt, compression=\"gzip\")\n",
    "    g_sr_bbt.attrs['description'] = 'bright band top in m above NN'\n",
    "    \n",
    "    g_sr_lza = grp_sr.create_dataset('sr_lza', data=sr_lza, compression=\"gzip\")\n",
    "    g_sr_lza.attrs['description'] = 'Locan Zenith Angle'\n",
    "    \n",
    "    g_sr_flh = grp_sr.create_dataset('sr_flh', data=sr_flh, compression=\"gzip\")\n",
    "    g_sr_flh.attrs['description'] = 'freezing level height in m above NN'\n",
    "    \n",
    "    g_sr_hip = grp_sr.create_dataset('sr_hip', data=sr_hip, compression=\"gzip\")\n",
    "    g_sr_hip.attrs['description'] = 'Heavy Ice Precip Flag'\n",
    "    \n",
    "    g_sr_shr = grp_sr.create_dataset('sr_shr', data=sr_shr, compression=\"gzip\")\n",
    "    g_sr_shr.attrs['description'] = 'Flag for Shallow rain'\n",
    "    \n",
    "    g_sr_anv = grp_sr.create_dataset('sr_anv', data=sr_anv, compression=\"gzip\")\n",
    "    g_sr_anv.attrs['description'] = 'Flag for Anvil'\n",
    "    \n",
    "    g_sr_flb = grp_sr.create_dataset('sr_flb', data=sr_flb, compression=\"gzip\")\n",
    "    g_sr_flb.attrs['description'] = 'Flag for BB'\n",
    "    \n",
    "    ## GR\n",
    "    g_gr_ipp = grp_gr.create_dataset('gr_pp_ipoli', data=gr_pp_ipoli, compression=\"gzip\")\n",
    "    g_gr_ipp.attrs['description'] = 'ground precip ipol on sr grid'\n",
    "    g_gr_xyc = grp_gr.create_dataset('gr_xy', data=gr_xy, compression=\"gzip\")\n",
    "    g_gr_xyc.attrs['description'] = 'ground based xy coordinates (km)'\n",
    "    g_gr_npp = grp_gr.create_dataset('gr_pp', data=gr_pp, compression=\"gzip\")\n",
    "    g_gr_npp.attrs['description'] = 'ground precip on gr grid'\n",
    "    g_gr_hmax = grp_gr.create_dataset('gr_h_max', data=gr_h_max, compression=\"gzip\")\n",
    "    g_gr_hmax.attrs['description'] = 'max measured height over ground on sr grid (m)'\n",
    "    g_gr_hmin = grp_gr.create_dataset('gr_h_min', data=gr_h_min, compression=\"gzip\")\n",
    "    g_gr_hmin.attrs['description'] = 'min measured height over ground on sr grid (m)'\n",
    "    g_gr_hmea = grp_gr.create_dataset('gr_h_mean', data=gr_h_mean, compression=\"gzip\")\n",
    "    g_gr_hmea.attrs['description'] = 'mean measured height over ground on sr grid (m)'\n",
    "    g_gr_hstd = grp_gr.create_dataset('gr_h_std', data=gr_h_std, compression=\"gzip\")\n",
    "    g_gr_hstd.attrs['description'] = 'std measured height over ground on sr grid (m)'\n",
    "    \n",
    "    g_gr_cbw = grp_gr.create_dataset('gr_c_beamwidth', data=gr_beam_width, compression=\"gzip\")\n",
    "    g_gr_cbw.attrs['description'] = 'gr beam width measured on sr grid (m)'\n",
    "    \n",
    "    ## DEM altitude and slope\n",
    "    g_gr_dems = grp_gr.create_dataset('gr_dem_std', data=gr_dem_std, compression=\"gzip\")\n",
    "    g_gr_dems.attrs['description'] = 'Std DEM in m'\n",
    "    \n",
    "    g_gr_dem = grp_gr.create_dataset('gr_dem_mean', data=gr_dem_mean, compression=\"gzip\")\n",
    "    g_gr_dem.attrs['description'] = 'Mean DEM in m'\n",
    "    \n",
    "    g_gr_slopedems = grp_gr.create_dataset('gr_slope_std', data=gr_slope_std, compression=\"gzip\")\n",
    "    g_gr_slopedems.attrs['description'] = 'Std DEM slope in deg'\n",
    "    \n",
    "    g_gr_slopedem = grp_gr.create_dataset('gr_slope_mean', data=gr_slope_mean, compression=\"gzip\")\n",
    "    g_gr_slopedem.attrs['description'] = 'Mean DEM slope in deg'\n",
    "    \n",
    "    # Times\n",
    "    g_sr_otime = grp_sr.create_dataset('sr_time', data=sr_time)\n",
    "    g_sr_otime.attrs['description'] = 'SR Overpass time'\n",
    "    g_gr_otime = grp_gr.create_dataset('gr_time', data=gr_time)\n",
    "    g_gr_otime.attrs['description'] = 'GR Meassured time'\n",
    "\n",
    "    hdat.close()\n",
    "    \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def calc_height(bin_range_number, elipsoid_offset, localZenithAngle, scan_swath='NS'):\n",
    "    # Function for height calculation from bins.....\n",
    "    # GPM DPR ATBD Level 2 2017 Awaka et al S. 21\n",
    "\n",
    "    if scan_swath in ('NS', 'MS'):\n",
    "        rangeBinSize = 125.\n",
    "        binEllipsoid = 176.\n",
    "    else:\n",
    "        rangeBinSize = 250.\n",
    "        binEllipsoid = 88.\n",
    "\n",
    "    calc_1 = ((binEllipsoid - bin_range_number) * rangeBinSize + elipsoid_offset)\n",
    "\n",
    "    calc_h = calc_1  * np.cos(np.deg2rad(localZenithAngle))\n",
    "\n",
    "    return calc_h\n",
    "\n",
    "\n",
    "\n",
    "def calc_bin(height, elipsoid_offset, localZenithAngle, scan_swath='NS'):\n",
    "\n",
    "    # Function for height calculation from bins.....\n",
    "    # GPM DPR ATBD Level 2 2017 Awaka et al S. 21\n",
    "\n",
    "    if scan_swath in ('NS', 'MS'):\n",
    "        rangeBinSize = 125.\n",
    "        binEllipsoid = 176.\n",
    "    else:\n",
    "        rangeBinSize = 250.\n",
    "        binEllipsoid = 88.\n",
    "\n",
    "    res0 = (height/np.cos(np.deg2rad(localZenithAngle))) - elipsoid_offset \n",
    "    res1 = res0 / rangeBinSize \n",
    "    res_bin = binEllipsoid - res1\n",
    "\n",
    "    return res_bin"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Import other parameters\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%%time\n",
    "\n",
    "# Import radolan grid src\n",
    "radolan_xy = wradlib.georef.get_radolan_grid(900,900)\n",
    "radolan_xy = radolan_xy.reshape(-1, radolan_xy.shape[-1])\n",
    "zd = wradlib.zonalstats.DataSource(radolan_xy, name='src')\n",
    "\n",
    "\n",
    "# Determine projection\n",
    "proj_stereo = wradlib.georef.create_osr(\"dwd-radolan\")\n",
    "proj_wgs = osr.SpatialReference()\n",
    "proj_wgs.ImportFromEPSG(4326)\n",
    "\n",
    "# Import RADOLAN Heights\n",
    "rhmax = np.load(\"/automount/ags/velibor/data/radolan_dx/RY_H_max.npy\")\n",
    "rhmin = np.load(\"/automount/ags/velibor/data/radolan_dx/RY_H_min.npy\")\n",
    "rhmean = np.load(\"/automount/ags/velibor/data/radolan_dx/RY_H_mean.npy\")\n",
    "#rhstd = np.load(\"/automount/ags/velibor/data/radolan_dx/RADOLAN_H_xystd_900x900.npy\")[2]\n",
    "\n",
    "# Radolan hoehen ausgewaehlt mit nearest\n",
    "#rhmean = np.load('/automount/ags/velibor/data/radolan_dx/RADOLAN_H_XYNEAREST_with_DEM_900x900.npy')\n",
    "#rhmean = np.load(\"/automount/ags/velibor/data/radolan_dx/RADOLAN_H_xymean_900x900_dem1x1.npy\")[2]\n",
    "\n",
    "#radolan_height_std_cut = rhstd.ravel()[ry_idx2]\n",
    "radolan_beam_width = np.load('/automount/ags/velibor/data/radolan_dx/RY_C_min.npy')\n",
    "\n",
    "dem_slope = np.load('/automount/ags/velibor/data/radolan_dx/slope_array.npy')\n",
    "\n",
    "\n",
    "## DEM\n",
    "radolan_tif_path = '/automount/db01/python/data/DEM/radolan/radolan_900x900.tif'\n",
    "ds = wradlib.io.open_raster(radolan_tif_path)\n",
    "ds = wradlib.georef.reproject_raster_dataset(ds, spacing=1)\n",
    "cor_dem = wradlib.georef.raster.read_gdal_coordinates(ds,mode='centers')\n",
    "val_dem = wradlib.georef.raster.read_gdal_values(ds)\n",
    "DEM = np.ma.masked_invalid(val_dem[::-1, :])\n",
    "\n",
    "# Masken von rh mean und dem entfernen\n",
    "DEM = np.ma.filled(DEM.astype(float), np.nan)\n",
    "rhmean = np.ma.filled(rhmean.astype(float), np.nan)\n",
    "rhmin = np.ma.filled(rhmin.astype(float), np.nan)\n",
    "rhmax = np.ma.filled(rhmax.astype(float), np.nan)\n",
    "radolan_beam_width = np.ma.filled(radolan_beam_width.astype(float), np.nan)\n",
    "\n",
    "DEM[DEM==-32768.] = np.nan \n",
    "\n",
    "# Other Parameters\n",
    "# dpr_footprint\n",
    "# dpr swath NS MS HS\n",
    "# Scan swath\n",
    "sc = 'NS'\n",
    "dpr_para = {'NS':176, 'HS':88, 'MS':176}\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "fig = pl.figure(figsize=(10,10))\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(221, aspect='equal')\n",
    "pl.pcolormesh(np.ma.masked_invalid(rhmin))\n",
    "pl.colorbar(shrink=0.5)\n",
    "\n",
    "ax1 = fig.add_subplot(222, aspect='equal')\n",
    "pl.pcolormesh(np.ma.masked_invalid(radolan_beam_width))\n",
    "pl.colorbar(shrink=0.5)\n",
    "\n",
    "ax1 = fig.add_subplot(223, aspect='equal')\n",
    "pl.pcolormesh(dem_slope)\n",
    "pl.colorbar(shrink=0.5)\n",
    "\n",
    "ax1 = fig.add_subplot(224, aspect='equal')\n",
    "pl.pcolormesh(np.ma.masked_invalid(DEM))\n",
    "pl.colorbar(shrink=0.5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#np.save('/automount/ags/velibor/data/radolan_dx/0ryh.npy', rhmin)\n",
    "#np.save('/automount/ags/velibor/data/radolan_dx/0ryc.npy', radolan_beam_width)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#%%time\n",
    "\n",
    "gpm_dir =sorted(glob.glob('/automount/ags/velibor/gpmdata/dprV7_small/*radolan.GPM*.HDF5'))\n",
    "#gpm_dir = gpm_dir[200:1003]\n",
    "#gpm_dir = gpm_dir[::-1]\n",
    "\n",
    "for gpm_file in gpm_dir:\n",
    "    # read gpm file with importatnt parameters\n",
    "    gpmdpr = h5py.File(gpm_file, 'r')\n",
    "\n",
    "    # gpm scantimes\n",
    "    gpm_scan_times = gpm_scan_time(gpm_file, scanswath=sc)\n",
    "\n",
    "    # GPM navigation and time\n",
    "    gpm_lat = np.array(gpmdpr[sc]['Latitude'])\n",
    "    gpm_lon = np.array(gpmdpr[sc]['Longitude'])\n",
    "    gpm_time = gpmdpr[sc]['ScanTime']\n",
    "\n",
    "\n",
    "    try:\n",
    "        #gpm_zeit = get_time_of_gpm(gpm_lon, gpm_lat, gpm_time)\n",
    "        gpm_zeit = get_time_of_gpm2(gpm_time)\n",
    "\n",
    "    except ValueError:\n",
    "        pass\n",
    "        print ('____________ValueError____________')\n",
    "    else:\n",
    "        #gpm_zeit = get_time_of_gpm(gpm_lon, gpm_lat, gpm_time)\n",
    "        print gpm_zeit\n",
    "        ht, mt = gpm_zeit[14:16], str(int(round(float(gpm_zeit[17:19])/5.0)*5.0))\n",
    "        year, ye, m, d = gpm_zeit[0:4], gpm_zeit[2:4], gpm_zeit[5:7], gpm_zeit[8:10]\n",
    "\n",
    "        if mt == '0':\n",
    "            mt = '00'\n",
    "        if mt == '5':\n",
    "            mt = '05'\n",
    "        if mt == '60':\n",
    "            mt = '55'\n",
    "            ht = str(int(ht)+1)\n",
    "            if ht == '24':\n",
    "                d = str(int(d)+1)\n",
    "\n",
    "        print gpm_zeit\n",
    "\n",
    "        r_pro = 'ry'\n",
    "        #try:\n",
    "\n",
    "\n",
    "        r_pfad = ('/automount/radar/dwd/'+ r_pro +'/'+str(year)+'/'+str(year)+'-'+\n",
    "                str(m)+'/'+ str(year)+'-'+str(m)+'-'+str(d)+'/raa01-'+r_pro+'_10000-'+\n",
    "                str(ye)+str(m)+ str(d)+str(ht)+str(mt)+'-dwd---bi*')\n",
    "        \n",
    "        print r_pfad\n",
    "\n",
    "\n",
    "        ## TRY: gibt es eine RADOLAN Datei zum Overpass!?\n",
    "        ## read radolan ry file data and attributes\n",
    "        try:\n",
    "            rwdata, rwattrs = wradlib.io.read_radolan_composite(glob.glob(r_pfad)[0])\n",
    "\n",
    "            radolan_grid_xy = wradlib.georef.get_radolan_grid(900,900)\n",
    "            x = radolan_grid_xy[:,:,0]\n",
    "            y = radolan_grid_xy[:,:,1]\n",
    "            rwdata = np.ma.masked_equal(rwdata, -9999) *12 # Um auf mm/h zu kommen, beim Einlesen sind es mm/5min\n",
    "\n",
    "            radolan_scan_time = rwattrs['datetime']\n",
    "        except:\n",
    "            print('__________________________________________RY File Missing!')\n",
    "            \n",
    "        #if len(rwattrs['radarlocations'])<17:\n",
    "            #print (\"------ Not all 17 Radars included in Composite ------\")\n",
    "            #pass    \n",
    "            \n",
    "        else:    \n",
    "            # copy data for binary grid\n",
    "            bingrid = rwdata.copy()\n",
    "\n",
    "            # Precipitation data in RADOLAN region\n",
    "            # Posible precip products:\n",
    "            \n",
    "            # precipRateNearSurface\n",
    "            dpr_pns = np.array(gpmdpr[sc]['SLV']['precipRateNearSurface'])\n",
    "            dpr_pns[dpr_pns ==-9999.9]= np.nan\n",
    "            # precipRateESurface\n",
    "            #dpr_pes = np.array(gpmdpr[sc]['SLV']['precipRateESurface'])\n",
    "            #dpr_pes[dpr_pes ==-9999.9]= np.nan\n",
    "            # precipRateAve24\n",
    "            dpr_pav = np.array(gpmdpr[sc]['SLV']['precipRateAve24'])\n",
    "            dpr_pav[dpr_pav ==-9999.9]= np.nan\n",
    "            # precip rate 3d\n",
    "            dpr_p3d = np.array(gpmdpr[sc]['SLV']['precipRate'])\n",
    "            dpr_p3d[dpr_p3d ==-9999.9]= np.nan\n",
    "\n",
    "            # Brightband height\n",
    "            dpr_bbh = np.array(gpmdpr[sc]['CSF']['heightBB'], dtype=float)\n",
    "            dpr_bbh[dpr_bbh ==-9999.9] = np.nan\n",
    "\n",
    "            # Bright-band width\n",
    "            dpr_bbw = np.array(gpmdpr[sc]['CSF']['widthBB'], dtype=float)\n",
    "            dpr_bbw[dpr_bbw ==-9999.9] = np.nan\n",
    "            \n",
    "            # Bright-band Top/Bot\n",
    "            dpr_bot = np.array(gpmdpr[sc]['CSF']['binBBBottom'], dtype=float)\n",
    "            dpr_bot[dpr_bot ==-9999] = np.nan\n",
    "            dpr_bot[dpr_bot ==-1111] = np.nan\n",
    "            dpr_bot[dpr_bot ==0] = np.nan\n",
    "            \n",
    "            dpr_top = np.array(gpmdpr[sc]['CSF']['binBBTop'], dtype=float)\n",
    "            dpr_top[dpr_top ==-9999] = np.nan\n",
    "            dpr_top[dpr_top ==-1111] = np.nan\n",
    "            dpr_top[dpr_top ==0] = np.nan\n",
    "            \n",
    "            # FreezingLevelHeight\n",
    "            dpr_flh = np.array(gpmdpr[sc]['VER']['heightZeroDeg'], dtype=float)\n",
    "            dpr_flh[dpr_flh ==-9999.9] = np.nan\n",
    "\n",
    "            # PhaseNearSurface\n",
    "            dpr_pha = np.array(gpmdpr[sc]['SLV']['phaseNearSurface'], dtype=float)\n",
    "            #dpr_pha = dpr_pha/100\n",
    "            dpr_pha[dpr_pha==255]=np.nan\n",
    "\n",
    "            #typePrecip\n",
    "            dpr_typ = np.array(gpmdpr[sc]['CSF']['typePrecip'], dtype=float)\n",
    "\n",
    "            #StromTopHeight\n",
    "            dpr_sth = np.array(gpmdpr[sc]['PRE']['heightStormTop'], dtype=float)\n",
    "            dpr_sth[dpr_sth ==-9999.9] = np.nan\n",
    "\n",
    "            #Elevation\n",
    "            dpr_ele = np.array(gpmdpr[sc]['PRE']['elevation'], dtype=float)\n",
    "            dpr_ele[dpr_ele == -9999] = np.nan\n",
    "\n",
    "            #Clutterfreebottom\n",
    "            dpr_cfb_bin = np.array(gpmdpr[sc]['PRE']['binClutterFreeBottom'], dtype=float)\n",
    "            dpr_cfb_bin[dpr_cfb_bin==-9999]=np.nan\n",
    "            \n",
    "            dpr_eof = np.array(gpmdpr[sc]['PRE']['ellipsoidBinOffset'], dtype=float)\n",
    "            dpr_eof[dpr_eof==-9999]=np.nan\n",
    "            \n",
    "            dpr_lza = np.array(gpmdpr[sc]['PRE']['localZenithAngle'], dtype=float)\n",
    "            dpr_lza[dpr_lza==-9999]=np.nan\n",
    "            \n",
    "            ## Berechnung der clutterfree bottom height\n",
    "            dpr_cfb = calc_height(dpr_cfb_bin, dpr_eof, dpr_lza, scan_swath=sc)\n",
    "            dpr_bbt = calc_height(dpr_top, dpr_eof, dpr_lza, scan_swath=sc)\n",
    "            dpr_bbb = calc_height(dpr_bot, dpr_eof, dpr_lza, scan_swath=sc)\n",
    "            \n",
    "            # Landsurfacetyp\n",
    "            dpr_lst = np.array(gpmdpr[sc]['PRE']['landSurfaceType'], dtype=float)\n",
    "            dpr_lst[dpr_lst==-9999]=np.nan\n",
    "            \n",
    "            ## Flags\n",
    "            # HeavyIcePrecipitatio\n",
    "            dpr_hip = np.array(gpmdpr[sc]['CSF']['flagHeavyIcePrecip'], dtype=float)\n",
    "            dpr_hip[dpr_hip==-99]=np.nan\n",
    "            # flagShallowRain\n",
    "            dpr_shr = np.array(gpmdpr[sc]['CSF']['flagShallowRain'], dtype=float)\n",
    "            dpr_shr[dpr_shr==-9999]=np.nan\n",
    "            \n",
    "            # flagBB\n",
    "            dpr_flb = np.array(gpmdpr[sc]['CSF']['flagBB'], dtype=float)\n",
    "            \n",
    "            \n",
    "            # Anvil\n",
    "            dpr_anv = np.array(gpmdpr[sc]['CSF']['flagAnvil'], dtype=float)\n",
    "            dpr_anv[dpr_anv==-99]=np.nan\n",
    "            \n",
    "\n",
    "            #3D Phase\n",
    "            dpr_3dpha = np.array(gpmdpr[sc]['DSD']['phase'], dtype=float)\n",
    "            dpr_3dpha[dpr_3dpha==255]=np.nan\n",
    "            dpr_3dpha = dpr_3dpha//100\n",
    "            \n",
    "            #dpr_3dpha[dpr_3dpha == -9999] = np.nan\n",
    "            \n",
    "            # PiaFinal\n",
    "            #dpr_pia = np.array(gpmdpr[sc]['SLV']['piaFinal'], dtype=float)\n",
    "            #dpr_pia[dpr_pia == -9999] = np.nan\n",
    "            \n",
    "            # flagHeavyIcePrecip FLAG!\n",
    "            # dpr_ice = np.array(gpmdpr[sc]['PRE']['flagHeavyIcePrecip'], dtype=float)\n",
    "            # dpr_ice[dpr_ice == -99] = np.nan\n",
    "            \n",
    "\n",
    "            # GPM coordinate projection\n",
    "            gpm_x, gpm_y = wradlib.georef.reproject(gpm_lon, gpm_lat,\n",
    "                                                    projection_target=proj_stereo,\n",
    "                                                    projection_source=proj_wgs)\n",
    "\n",
    "\n",
    "            ## Remove all GPM footprint that are not in scan area of radolan\n",
    "            ## -------------------------------------------------------------\n",
    "\n",
    "            # GPM Koordinaten ravel\n",
    "            gpm_xy = np.vstack((gpm_x.ravel(), gpm_y.ravel())).transpose()\n",
    "\n",
    "            # index for GPM footprints in scan area\n",
    "            oi_idx = idx_sr_in_gr_area(bingrid, gpm_xy)\n",
    "            \n",
    "            # check if idx is empty\n",
    "            if len(oi_idx)==0:\n",
    "                print (\"____No Index of Overpass in Radolan region____\")\n",
    "                pass\n",
    "            \n",
    "            \n",
    "            else:\n",
    "                # entferne alle Indizes auserhalb des äußeren Polygons\n",
    "                gpm_xy = gpm_xy[oi_idx]\n",
    "                gpm_pns = dpr_pns.ravel()[oi_idx]\n",
    "                #gpm_pes = dpr_pes.ravel()[oi_idx]\n",
    "                gpm_pav = dpr_pav.ravel()[oi_idx]\n",
    "                gpm_bbh = dpr_bbh.ravel()[oi_idx]\n",
    "                gpm_bbw = dpr_bbw.ravel()[oi_idx]\n",
    "                gpm_pha = dpr_pha.ravel()[oi_idx]\n",
    "                gpm_typ = dpr_typ.ravel()[oi_idx]\n",
    "                gpm_sth = dpr_sth.ravel()[oi_idx]\n",
    "                gpm_ele = dpr_ele.ravel()[oi_idx]\n",
    "                gpm_cfb = dpr_cfb.ravel()[oi_idx]\n",
    "                gpm_lst = dpr_lst.ravel()[oi_idx]\n",
    "                gpm_tim = gpm_scan_times.ravel()[oi_idx]\n",
    "                # Parameters only for new pns product\n",
    "                dpr_p3d = np.reshape(dpr_p3d,[dpr_pns.shape[0]*dpr_pns.shape[1],dpr_para[sc]])\n",
    "                gpm_p3d = dpr_p3d[oi_idx,:]\n",
    "                # 3d phase for adjusting\n",
    "                dpr_3dpha = np.reshape(dpr_3dpha,[dpr_3dpha.shape[0]*dpr_3dpha.shape[1],dpr_para[sc]])\n",
    "                gpm_ph3d = dpr_3dpha[oi_idx,:]\n",
    "                \n",
    "                gpm_cfb = dpr_cfb.ravel()[oi_idx]\n",
    "                gpm_eof = dpr_eof.ravel()[oi_idx]\n",
    "                gpm_cfb_bin = dpr_cfb_bin.ravel()[oi_idx]\n",
    "                \n",
    "                gpm_lza = dpr_lza.ravel()[oi_idx]\n",
    "                gpm_bbb = dpr_bbb.ravel()[oi_idx]\n",
    "                gpm_bbt = dpr_bbt.ravel()[oi_idx]\n",
    "                gpm_flh = dpr_flh.ravel()[oi_idx]\n",
    "                gpm_hip = dpr_hip.ravel()[oi_idx]\n",
    "                gpm_shr = dpr_shr.ravel()[oi_idx]\n",
    "                gpm_anv = dpr_anv.ravel()[oi_idx]\n",
    "                gpm_flb = dpr_flb.ravel()[oi_idx]\n",
    "\n",
    "\n",
    "                ## Remove all RADOLAN Points not included in GPM DPR swath\n",
    "                ## -------------------------------------------------------\n",
    "\n",
    "                # Contour from ORIGINAL DPR Overpass\n",
    "                dpr_contour_x = dpr_swath_contour(gpm_x)\n",
    "                dpr_contour_y = dpr_swath_contour(gpm_y)\n",
    "                print(dpr_contour_x.shape)\n",
    "\n",
    "                # create dpr polygon of dpr xy contours \n",
    "                dpr_xy_poly = np.vstack((dpr_contour_x.ravel(), dpr_contour_y.ravel())).transpose()\n",
    "\n",
    "                #load regular radolan grid and overlie with dpr contours\n",
    "                zd_poly_radolan = wradlib.zonalstats.ZonalDataPoint(zd, [dpr_xy_poly], buf=2.5)\n",
    "\n",
    "                # get radolan index in dpr swath\n",
    "                ry_idx = zd_poly_radolan.get_source_index(0)\n",
    "\n",
    "                # get radolan index in dpr swath AND in binary grid\n",
    "                rwtest = ~np.ma.masked_less(rwdata, 0).mask\n",
    "\n",
    "                radolan_idx = np.flatnonzero(rwtest.ravel())\n",
    "\n",
    "                ry_idx2 = np.intersect1d(ry_idx, radolan_idx, assume_unique=True)\n",
    "\n",
    "                # Radolan xy array\n",
    "                xy_radolan = np.vstack((x.ravel(), y.ravel())).transpose()\n",
    "\n",
    "                # Extract all affected RADOLAN grid points\n",
    "                rwdata_cut = rwdata.ravel()[ry_idx2]\n",
    "                x_cut = xy_radolan[...,0].ravel()[ry_idx2].copy()\n",
    "                y_cut = xy_radolan[...,1].ravel()[ry_idx2].copy()\n",
    "\n",
    "\n",
    "                # extrect all Height points (need to have the same shape as radolan grid)\n",
    "                radolan_height_min_cut = rhmin.ravel()[ry_idx2]\n",
    "                radolan_height_max_cut = rhmax.ravel()[ry_idx2]\n",
    "                #radolan_height_std_cut = rhstd.ravel()[ry_idx2]\n",
    "                radolan_height_mean_cut = rhmean.ravel()[ry_idx2]\n",
    "\n",
    "                radolan_beam_width_cut = radolan_beam_width.ravel()[ry_idx2]\n",
    "                \n",
    "                DEM_cut = DEM.ravel()[ry_idx2]\n",
    "                dem_slope_cut = dem_slope.ravel()[ry_idx2]\n",
    "\n",
    "\n",
    "                # Compare all cutted radolan x and y \n",
    "                xy_cut = np.vstack((x_cut.ravel(), y_cut.ravel())).transpose()\n",
    "\n",
    "                #Dpr footprint radius \n",
    "                dpr_footprint = 2.5 #radius in km\n",
    "\n",
    "                # Interpolation with antena weighted function\n",
    "                ## -------------------------------\n",
    "                ry_pns_w = ipoli_radi_w(xy_cut, gpm_xy, rwdata_cut, dpr_footprint, k=25)\n",
    "\n",
    "                # Interpolation of other parameters on dpr grid\n",
    "                radolan_height_min_cut_ipoli = ipoli_radi(xy_cut, gpm_xy, radolan_height_min_cut,\n",
    "                                                          dpr_footprint, k=25, calc='mean')\n",
    "                radolan_height_max_cut_ipoli = ipoli_radi(xy_cut, gpm_xy, radolan_height_max_cut,\n",
    "                                                          dpr_footprint, k=25, calc='mean')\n",
    "                radolan_height_mean_cut_ipoli = ipoli_radi(xy_cut, gpm_xy, radolan_height_mean_cut,\n",
    "                                                           dpr_footprint, k=25, calc='mean')\n",
    "                radolan_height_std_cut_ipoli = ipoli_radi(xy_cut, gpm_xy, radolan_height_min_cut,\n",
    "                                                          dpr_footprint, k=25, calc='std')\n",
    "                ## DEM Altitude\n",
    "                DEM_cut_ipoli_mean = ipoli_radi(xy_cut, gpm_xy, DEM_cut,\n",
    "                                                          dpr_footprint, k=25, calc='mean')\n",
    "                DEM_cut_ipoli_std = ipoli_radi(xy_cut, gpm_xy, DEM_cut,\n",
    "                                                          dpr_footprint, k=25, calc='std')\n",
    "                ## DEM Slope\n",
    "                slop_cut_ipoli_mean = ipoli_radi(xy_cut, gpm_xy, dem_slope_cut,\n",
    "                                                          dpr_footprint, k=25, calc='mean')\n",
    "                \n",
    "                slop_cut_ipoli_std = ipoli_radi(xy_cut, gpm_xy, dem_slope_cut,\n",
    "                                                          dpr_footprint, k=25, calc='std')\n",
    "                \n",
    "                \n",
    "                radolan_beam_width_cut_ipoli = ipoli_radi(xy_cut, gpm_xy, radolan_beam_width_cut,\n",
    "                                                           dpr_footprint, k=25, calc='mean')\n",
    "\n",
    "                \n",
    "                ## Create adjusted pns and phase product from p3d\n",
    "                ## ----------------------------------------------\n",
    "\n",
    "                # delta h calculation\n",
    "                #delta_h = radolan_height_min_cut_ipoli - gpm_cfb\n",
    "                # Height over NN\n",
    "                delta_h = radolan_height_min_cut_ipoli - (gpm_cfb - gpm_ele)\n",
    "                \n",
    "                # new bin height  \n",
    "                new_bin = calc_bin(gpm_cfb + delta_h, gpm_eof, gpm_lza)\n",
    "                # round for bin\n",
    "                new_bin = np.ceil(new_bin)\n",
    "                # error calc with nan removal\n",
    "                new_bin[new_bin<0]=np.nan\n",
    "                \n",
    "                # Beachtung des local zenith angle\n",
    "                # gpm_bin_in_beam = beamwidth / sin(lza)\n",
    "                #rbwi = radolan_beam_width_cut_ipoli / np.cos(np.deg2rad(gpm_lza))\n",
    "                rbwi = radolan_beam_width_cut_ipoli\n",
    "                # column calculation ! Half column!\n",
    "                rcbin = np.round((rbwi/2.)/125.,0)\n",
    "                \n",
    "                # adjusted precip product \n",
    "                gpm_app = []\n",
    "                \n",
    "                # adjusted phase \n",
    "                # gpm_aph = []\n",
    "                \n",
    "                # number of bins with same phase (ice,liquid, mixed) in one column \n",
    "                gpm_lbn = []\n",
    "                gpm_ibn = []\n",
    "                gpm_mbn = []\n",
    "                    \n",
    "                    \n",
    "                for j in range(len(new_bin)):\n",
    "\n",
    "                    if np.isnan(new_bin[j]):\n",
    "                        # Check if bin is nan \n",
    "                        gpm_app.append(np.nan)\n",
    "                        gpm_lbn.append(np.nan)\n",
    "                        gpm_mbn.append(np.nan)\n",
    "                        gpm_ibn.append(np.nan)\n",
    "\n",
    "\n",
    "                    else:  \n",
    "                            if (new_bin[j]-1)>=dpr_para[sc]:\n",
    "                            # Check if bin index bigger than array\n",
    "                            # happen by round bin\n",
    "                                gpm_app.append(np.nan)\n",
    "                                gpm_lbn.append(np.nan)\n",
    "                                gpm_mbn.append(np.nan)\n",
    "                                gpm_ibn.append(np.nan)\n",
    "                            else:\n",
    "                                nbin = int(new_bin[j])-1\n",
    "                                cbin = rcbin[j]\n",
    "                                \n",
    "                                # Averaged Rainrate in gpm cpolumn (later Antenna)\n",
    "                                start_bin = nbin-int(cbin)\n",
    "                                end_bin = nbin+int(cbin)+1\n",
    "                                if start_bin<0:\n",
    "                                    start_bin = 0 \n",
    "                                if end_bin >=dpr_para[sc]:\n",
    "                                    end_bin = dpr_para[sc]-1\n",
    "                                    \n",
    "                                # Mean Precip rate in colum... later antena function C-band\n",
    "                                rr = gpm_p3d[j, start_bin : end_bin]\n",
    "                                rr = np.nanmean(rr)\n",
    "                                phase = gpm_ph3d[j, start_bin : end_bin]\n",
    "                                \n",
    "                                # Phase bin ratios... later also weighting!\n",
    "                                idx_l = phase == 2\n",
    "                                idx_i = phase == 0\n",
    "                                idx_m = phase == 1\n",
    "                                \n",
    "                                \n",
    "                                gpm_app.append(rr)\n",
    "                                #gpm_aph.append(gpm_ph3d[j,nbin])                                \n",
    "                                gpm_lbn.append(sum(idx_l))\n",
    "                                gpm_ibn.append(sum(idx_i))\n",
    "                                gpm_mbn.append(sum(idx_m))\n",
    "                                \n",
    "\n",
    "                gpm_app = np.array(gpm_app)\n",
    "                gpm_lbn = np.array(gpm_lbn)\n",
    "                gpm_ibn = np.array(gpm_ibn)\n",
    "                gpm_mbn = np.array(gpm_mbn)\n",
    "\n",
    "                ## Write and save to HDF5 file\n",
    "                ## ----------------------------\n",
    "                name = '/automount/ags/velibor/gpmdata/dumpdataV7/' \\\n",
    "                                   'dpr_ry_'+sc+'/dprrado_'+ sc + str(gpm_zeit)\n",
    "                \n",
    "                # array with datetime in array with datimestrings\n",
    "                gpm_tim_str = [dt.datetime.strftime(date, \"%Y-%m-%d %H:%M:%S\") for date in gpm_tim]    \n",
    "\n",
    "                write2hdf(name, gpm_xy[:,0], gpm_xy[:,1],\n",
    "                          gpm_pns, gpm_app, gpm_pav,\n",
    "                          gpm_bbh, gpm_bbw, gpm_pha,\n",
    "                          gpm_typ, gpm_sth, gpm_cfb,\n",
    "                          gpm_lst, gpm_ele, gpm_tim_str,\n",
    "                          ry_pns_w, xy_cut, rwdata_cut,\n",
    "                          radolan_height_max_cut_ipoli,\n",
    "                          radolan_height_min_cut_ipoli, \n",
    "                          radolan_height_mean_cut_ipoli,\n",
    "                          radolan_height_std_cut_ipoli,\n",
    "                          radolan_beam_width_cut_ipoli,\n",
    "                          gpm_zeit,\n",
    "                          dt.datetime.strftime(radolan_scan_time, \"%Y-%m-%d %H:%M:%S\"),\n",
    "                          gpm_lbn, gpm_ibn, gpm_mbn,\n",
    "                          gpm_bbb, gpm_bbt, gpm_lza, gpm_flh, gpm_hip, gpm_shr, gpm_anv, gpm_flb,\n",
    "                          DEM_cut_ipoli_mean, DEM_cut_ipoli_std, slop_cut_ipoli_std, slop_cut_ipoli_mean)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "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\n",
    "\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/ftp/velibor/data/'+dict_name+'.npy',ac) \n",
    "    para_dict.pop(dict_name,None)\"\"\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "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",
    "###### TODO calc also Timedifferenze!!!!!!!!!!! ##############\n",
    "\n",
    "import datetime as dt\n",
    "\n",
    "pfad = glob.glob(\"/automount/ags/velibor/gpmdata/dumpdataV7/dpr_ry_NS/dprrado*\")\n",
    "\n",
    "pfad = sorted(pfad)\n",
    "\n",
    "print ('Es gibt ' + str(len(pfad)) + ' Overpasses')\n",
    "\n",
    "para_dict = {#'sr_x': [],\n",
    "             #'sr_y': [],\n",
    "             #'gr_c_beamwidth':[],\n",
    "             #'sr_mbn': [],\n",
    "             #'sr_lbn': [],\n",
    "             #'sr_ibn': [],\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_bbb': [],\n",
    "             #'sr_bbt': [],\n",
    "             #'sr_lza': [],\n",
    "             #'sr_flh': [],\n",
    "             #'sr_hip': [],\n",
    "             #'sr_shr': [],\n",
    "             #'sr_anv': [],\n",
    "             #'sr_type': [],\n",
    "             #'sr_flb':[],\n",
    "             #'gr_dem_std':[],\n",
    "             #'gr_dem_mean':[],\n",
    "             #'gr_slope_std':[],\n",
    "             'gr_slope_mean':[]\n",
    "}   \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/ftp/velibor/data/'+dict_name+'.npy',ac) \n",
    "    para_dict.pop(dict_name,None)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "sorted(para_dict)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "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",
    "###### TODO calc also Timedifferenze!!!!!!!!!!! ##############\n",
    "\n",
    "import datetime as dt\n",
    "\n",
    "pfad = glob.glob(\"/automount/ags/velibor/gpmdata/dumpdataV7/dpr_ry_NS/dprrado*\")\n",
    "pfad = sorted(pfad)\n",
    "\n",
    "print ('Es gibt ' + str(len(pfad)) + ' Overpasses')\n",
    "\n",
    "para_dict = {'sr_times': []}   # 'gr_time': np.array([]),\n",
    "\n",
    "a=0\n",
    "\n",
    "for p in sorted(pfad):\n",
    "\n",
    "    DR = h5py.File(p, 'r')\n",
    "    # for über liste\n",
    "    a+=1\n",
    "    # Listen append x.append(x)\n",
    "    T1 = np.array(DR['SR']['sr_times'][:], dtype='datetime64')\n",
    "    T2 = np.array(DR['GR']['gr_time'].value, dtype='datetime64' )\n",
    "                \n",
    "    para_dict['sr_times'].append( T1 - T2 )\n",
    "\n",
    "    # ! Colse hdf5 file ! \n",
    "    DR.close()\n",
    "    del(DR)\n",
    "        \n",
    "        \n",
    "    ac = np.concatenate(para_dict['sr_times'],axis=0)\n",
    "    #print (ac.shape)\n",
    "\n",
    "    np.save('/automount/ftp/velibor/data/'+'timedelta'+'.npy',ac) \n",
    "    #para_dict.pop(dict_name,None)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "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",
    "###### TODO calc also Timedifferenze!!!!!!!!!!! ##############\n",
    "\n",
    "import datetime as dt\n",
    "\n",
    "pfad = glob.glob(\"/automount/ags/velibor/gpmdata/dumpdataV7/dpr_ry_NS/dprrado*\")\n",
    "\n",
    "pfad = sorted(pfad)\n",
    "\n",
    "print ('Es gibt ' + str(len(pfad)) + ' Overpasses')\n",
    "\n",
    "para_dict = {'sr_times': []}   # 'gr_time': np.array([]),\n",
    "\n",
    "a=0\n",
    "\n",
    "for p in sorted(pfad):\n",
    "\n",
    "    DR = h5py.File(p, 'r')\n",
    "    # for über liste\n",
    "    a+=1\n",
    "    # Listen append x.append(x)\n",
    "    T2 = DR['GR']['gr_time'].value\n",
    "                \n",
    "    para_dict['sr_times'].append(T2)\n",
    "\n",
    "    # ! Colse hdf5 file ! \n",
    "    DR.close()\n",
    "    del(DR)\n",
    "        \n",
    "        \n",
    "    #ac = np.concatenate(para_dict['sr_times'],axis=0)\n",
    "    #print (ac.shape)\n",
    "\n",
    "    np.save('/automount/ftp/velibor/data/'+'overpasstimes'+'.npy',ac) \n",
    "    #para_dict.pop(dict_name,None)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "default",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}