_0_PAPER_pcc_dx_dpr_RY.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 plt\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",
    "import numpy as np\n",
    "import glob"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def get_radar_locations():\n",
    "\n",
    "    radars = {}\n",
    "    radar = {}\n",
    "    radar['name'] = 'ASR Dresden'\n",
    "    radar['wmo'] = 10487\n",
    "    radar['lon'] = 13.76347\n",
    "    radar['lat'] = 51.12404\n",
    "    radar['alt'] = 261\n",
    "    radars['ASD'] = radar\n",
    "    \n",
    "    radar = {}\n",
    "    #### Kordinaten un hoehe checken! Neu\n",
    "    # Hoehe üNN 2m\n",
    "    # Turm 55 m\n",
    "    #radar['name'] = 'Borkum'\n",
    "    #radar['wmo'] = 10103\n",
    "    #radar['lon'] = 6.748292\n",
    "    #radar['lat'] = 53.564011\n",
    "    #radar['alt'] = 57.\n",
    "    #radars['ASB'] = radar\n",
    "    \n",
    "    #Fake Borkum/Emden\n",
    "    radar['name'] = 'Emden'\n",
    "    radar['wmo'] = 10103\n",
    "    radar['lon'] = 7.02377\n",
    "    radar['lat'] = 53.33872\n",
    "    radar['alt'] = 58.\n",
    "    radars['ASB'] = radar\n",
    "\n",
    "    radar = {}\n",
    "    radar['name'] = 'Boostedt'\n",
    "    radar['wmo'] = 10132\n",
    "    radar['lon'] = 10.04687\n",
    "    radar['lat'] = 54.00438\n",
    "    radar['alt'] = 124.56\n",
    "    radars['BOO'] = radar\n",
    "\n",
    "    radar = {}\n",
    "    radar['name'] = 'Dresden'\n",
    "    radar['wmo'] = 10488\n",
    "    radar['lon'] = 13.76865\n",
    "    radar['lat'] = 51.12465\n",
    "    radar['alt'] = 263.36\n",
    "    radars['DRS'] = radar\n",
    "\n",
    "    radar = {}\n",
    "    radar['name'] = 'Eisberg'\n",
    "    radar['wmo'] = 10780\n",
    "    radar['lon'] = 12.40278\n",
    "    radar['lat'] = 49.54066\n",
    "    radar['alt'] = 798.79\n",
    "    radars['EIS'] = radar\n",
    "\n",
    "    radar = {}\n",
    "    radar['name'] = 'Emden'\n",
    "    radar['wmo'] = 10204\n",
    "    radar['lon'] = 7.02377\n",
    "    radar['lat'] = 53.33872\n",
    "    radar['alt'] = 58\n",
    "    radars['EMD'] = radar\n",
    "\n",
    "    radar = {}\n",
    "    radar['name'] = 'Essen'\n",
    "    radar['wmo'] = 10410\n",
    "    radar['lon'] = 6.96712\n",
    "    radar['lat'] = 51.40563\n",
    "    radar['alt'] = 185.10\n",
    "    radars['ESS'] = radar\n",
    "\n",
    "    radar = {}\n",
    "    radar['name'] = 'Feldberg'\n",
    "    radar['wmo'] = 10908\n",
    "    radar['lon'] = 8.00361\n",
    "    radar['lat'] = 47.87361\n",
    "    radar['alt'] = 1516.10\n",
    "    radars['FBG'] = radar\n",
    "\n",
    "    radar = {}\n",
    "    radar['name'] = 'Flechtdorf'\n",
    "    radar['wmo'] = 10440\n",
    "    radar['lon'] = 8.802\n",
    "    radar['lat'] = 51.3112\n",
    "    radar['alt'] = 627.88\n",
    "    radars['FLD'] = radar\n",
    "\n",
    "    radar = {}\n",
    "    radar['name'] = 'Hannover'\n",
    "    radar['wmo'] = 10339\n",
    "    radar['lon'] = 9.69452\n",
    "    radar['lat'] = 52.46008\n",
    "    radar['alt'] = 97.66\n",
    "    radars['HNR'] = radar\n",
    "\n",
    "    radar = {}\n",
    "    radar['name'] = 'Neuhaus'\n",
    "    radar['wmo'] = 10557\n",
    "    radar['lon'] = 11.13504\n",
    "    radar['lat'] = 50.50012\n",
    "    radar['alt'] = 878.04\n",
    "    radars['NEU'] = radar\n",
    "\n",
    "    radar = {}\n",
    "    radar['name'] = 'Neuheilenbach'\n",
    "    radar['wmo'] = 10605\n",
    "    radar['lon'] = 6.54853\n",
    "    radar['lat'] = 50.10965\n",
    "    radar['alt'] = 585.84\n",
    "    radars['NHB'] = radar\n",
    "\n",
    "    radar = {}\n",
    "    radar['name'] = 'Offenthal'\n",
    "    radar['wmo'] = 10629\n",
    "    radar['lon'] = 8.71293\n",
    "    radar['lat'] = 49.9847\n",
    "    radar['alt'] = 245.80\n",
    "    radars['OFT'] = radar\n",
    "\n",
    "    radar = {}\n",
    "    radar['name'] = 'Proetzel'\n",
    "    radar['wmo'] = 10392\n",
    "    radar['lon'] = 13.85821\n",
    "    radar['lat'] = 52.64867\n",
    "    radar['alt'] = 193.92\n",
    "    radars['PRO'] = radar\n",
    "\n",
    "    radar = {}\n",
    "    radar['name'] = 'Memmingen'\n",
    "    radar['wmo'] = 10950\n",
    "    radar['lon'] = 10.21924\n",
    "    radar['lat'] = 48.04214\n",
    "    radar['alt'] = 724.40\n",
    "    radars['MEM'] = radar\n",
    "\n",
    "    radar = {}\n",
    "    radar['name'] = 'Rostock'\n",
    "    radar['wmo'] = 10169\n",
    "    radar['lon'] = 12.05808\n",
    "    radar['lat'] = 54.17566\n",
    "    radar['alt'] = 37\n",
    "    radars['ROS'] = radar\n",
    "\n",
    "    radar = {}\n",
    "    radar['name'] = 'Isen'\n",
    "    radar['wmo'] = 10873\n",
    "    radar['lon'] = 12.10177\n",
    "    radar['lat'] = 48.1747\n",
    "    radar['alt'] = 677.77\n",
    "    radars['ISN'] = radar\n",
    "\n",
    "    radar = {}\n",
    "    radar['name'] = 'Tuerkheim'\n",
    "    radar['wmo'] = 10832\n",
    "    radar['lon'] = 9.78278\n",
    "    radar['lat'] = 48.58528\n",
    "    radar['alt'] = 767.62\n",
    "    radars['TUR'] = radar\n",
    "\n",
    "    radar = {}\n",
    "    radar['name'] = 'Ummendorf'\n",
    "    radar['wmo'] = 10356\n",
    "    radar['lon'] = 11.17609\n",
    "    radar['lat'] = 52.16009\n",
    "    radar['alt'] = 183\n",
    "    radars['UMM'] = radar\n",
    "\n",
    "    \n",
    "    return radars\n"
   ]
  },
  {
   "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": [
    "#ZP = '201410070235'\n",
    "ZP = '201803072145'\n",
    "#ZP = '201409072140'\n",
    "\n",
    "def read_all_dx(ZP):\n",
    "    \n",
    "    radar_range = 150\n",
    "\n",
    "    YYYY, YY, MM, DD, hh, mm = ZP[0:4], ZP[2:4], ZP[4:6], ZP[6:8], ZP[8:10], ZP[10:12]\n",
    "    #path = '/automount/radar/dwd/dx/'+YYYY+'/'+YYYY+'-'+MM+'/'+YYYY+'-'+MM+'-'+DD+'/'+'*'+YY+MM+DD+hh+mm+'*'\n",
    "    #path = '/automount/radar/dwd/dx/'+YYYY+'/'+YYYY+'-'+MM+'/'+YYYY+'-'+MM+'-'+DD+'/'+'*'+YY+MM+DD+hh+mm+'*'\n",
    "    #path = '/automount/ags/velibor/data/radolan_dx/'+YYYY+'/'+YYYY+'-'+MM+'/'+YYYY+'-'+MM+'-'+DD+'/'+'*'+YY+MM+DD+hh+mm+'*'\n",
    "    path = '/automount/ags/velibor/data/radolan_dx/test_for_height/*'\n",
    "\n",
    "    #print (path)\n",
    "    c_liste = sorted(glob.glob(path))\n",
    "    #print c_liste\n",
    "    #print (c_liste)\n",
    "    nr_rad = len(c_liste)\n",
    "    #print (nr_rad)\n",
    "    # Daten Array mit Reflektivitaeten\n",
    "    all_scans = np.empty((nr_rad,360,128))\n",
    "    # Atribute\n",
    "    metadata = [[None] ] * nr_rad\n",
    "\n",
    "    ## Alle Positionen lon lat alt\n",
    "    all_lon = np.empty((nr_rad,360,radar_range))\n",
    "    all_lat = np.empty((nr_rad,360,radar_range))\n",
    "    all_alt = np.empty((nr_rad,360,radar_range))\n",
    "    ## Beambreite\n",
    "    all_kkk = np.empty((nr_rad,360,radar_range))\n",
    "    ## Radar\n",
    "    all_rad = np.empty((nr_rad,360,radar_range))\n",
    "    ## Entfernung zum Radar\n",
    "    all_dis = np.empty((nr_rad,360,radar_range))\n",
    "\n",
    "\n",
    "    for jj in range(nr_rad):\n",
    "        filename = wrl.util.get_wradlib_data_file(c_liste[jj])\n",
    "        all_scans[jj], metadata[jj] = wrl.io.read_dx(filename)\n",
    "\n",
    "    ## Radar Location\n",
    "    radloc = get_radar_locations()\n",
    "    \n",
    "    lon0 = np.empty(nr_rad)\n",
    "    lat0 = np.empty(nr_rad)\n",
    "    alt0 = np.empty(nr_rad)\n",
    "    radname = []\n",
    "\n",
    "    for rr in range(nr_rad):\n",
    "        \n",
    "        for ir in radloc.keys():\n",
    "\n",
    "            if str(radloc[ir]['wmo'])==metadata[rr]['radarid']:\n",
    "                print(str(radloc[ir]['wmo']),metadata[rr]['radarid'],radloc[ir]['lon'],radloc[ir]['lat'],radloc[ir]['name'])\n",
    "                lon0[rr] = radloc[ir]['lon']\n",
    "                lat0[rr] = radloc[ir]['lat']\n",
    "                alt0[rr] = radloc[ir]['alt']\n",
    "                radname.append(radloc[ir]['name'])\n",
    "            else:\n",
    "\n",
    "                pass\n",
    "\n",
    "    #print radname\n",
    "    proj_stereo = wrl.georef.create_osr(\"dwd-radolan\")\n",
    "    proj_wgs = osr.SpatialReference()\n",
    "    proj_wgs.ImportFromEPSG(4326)\n",
    "\n",
    "    for rrr in range(nr_rad):\n",
    "\n",
    "        radar_location = (lon0[rrr], lat0[rrr], alt0[rrr])\n",
    "        \n",
    "        azimuths = metadata[rrr]['azim']\n",
    "        \n",
    "        ranges = np.arange(0,radar_range)*1000 \n",
    "        \n",
    "        \n",
    "        elevations = metadata[rrr]['elev']\n",
    "        \n",
    "        if (lon0[rrr]==7.02377) & (lat0[rrr]==53.33872) & (alt0[rrr]==58.):\n",
    "            print (\"Elevationen bei Emden eingestellt!\")\n",
    "            # Hier werden einfach alle Elevation auf eins gesetzt und mit\n",
    "            # der 0.9 Elevation wie in Emden eingestellt\n",
    "            elevations = (elevations / elevations) * 0.9\n",
    "            \n",
    "        print('--->>> Elevations: ', np.unique(elevations))\n",
    "\n",
    "        polargrid = np.meshgrid(ranges, azimuths)\n",
    "        elevgrid = np.meshgrid(ranges, elevations)\n",
    "        \n",
    "        r = polargrid[0]\n",
    "\n",
    "        theta = metadata[rrr]['elev']\n",
    "\n",
    "        # Ray distance\n",
    "        rrrc = r.T * np.cos(np.deg2rad(theta))\n",
    "        rrrc = rrrc.T\n",
    "        \n",
    "        # Keulenbreite\n",
    "        alpha = 0.9 # C-band Offnungswinkel\n",
    "        kkk =rrrc * np.sin(np.deg2rad(alpha)) \n",
    "        \n",
    "        ## Polar to Proj\n",
    "        grid = wrl.georef.polar.spherical_to_proj(rrrc, polargrid[1],elevgrid[1], radar_location)\n",
    "\n",
    "        lon, lat, alt = grid[:,:,0], grid[:,:,1], grid[:,:,2]\n",
    "\n",
    "        lon, lat = wrl.georef.reproject(lon, lat, projection_target=proj_stereo ,projection_source=proj_wgs)        \n",
    "        \n",
    "        all_lon[rrr,:,:] = lon\n",
    "        all_lat[rrr,:,:] = lat\n",
    "        all_alt[rrr,:,:] = alt\n",
    "        all_kkk[rrr,:,:] = kkk\n",
    "        all_rad[rrr,:,:] = (kkk/kkk)*rrr\n",
    "        all_dis[rrr,:,:] = rrrc\n",
    "        \n",
    "        #print ('Max Alt:', np.nanmax(alt))\n",
    "        #print ('Shapes: ', lon.shape, lat.shape, alt.shape)\n",
    "\n",
    "\n",
    "    return all_lon, all_lat, all_alt,all_kkk, all_scans, all_rad, all_dis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#Einlesen aller DX Standorte\n",
    "#all_lon, all_lat, all_alt, all_scans = read_all_dx(ZP)\n",
    "all_lon, all_lat, all_alt,all_kkk, all_scans, all_rad, all_dis = read_all_dx(ZP)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plt.figure(figsize=(18,10))\n",
    "\n",
    "for plpl in range(all_lon.shape[0]):\n",
    "    #print plpl\n",
    "    plt.pcolormesh(all_lon[plpl,...], all_lat[plpl,...], all_dis[plpl,...]/1000., vmin=0. , vmax=150)\n",
    "\n",
    "plt.xlim(-600,400)\n",
    "plt.ylim(-4700, -3700)\n",
    "plt.colorbar()\n",
    "plt.grid()\n",
    "plt.title('Distance in km')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plt.figure(figsize=(18,10))\n",
    "\n",
    "for plpl in range(all_lon.shape[0]):\n",
    "    #print plpl\n",
    "    plt.pcolormesh(all_lon[plpl,...], all_lat[plpl,...], all_kkk[plpl,...], vmin=0. )\n",
    "\n",
    "plt.xlim(-600,400)\n",
    "plt.ylim(-4700, -3700)\n",
    "plt.colorbar()\n",
    "plt.grid()\n",
    "plt.title('Ray Width')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plt.figure(figsize=(18,10))\n",
    "\n",
    "for plpl in range(all_lon.shape[0]):\n",
    "    #print plpl\n",
    "    plt.pcolormesh(all_lon[plpl,...], all_lat[plpl,...], all_alt[plpl,...], vmin=0, vmax=4500)\n",
    "    print(np.nanmax(all_alt[plpl,...]), np.nanmax(all_kkk[plpl,...]), np.nanmax(all_dis[plpl,...]))\n",
    "\n",
    "plt.xlim(-600,400)\n",
    "plt.ylim(-4700, -3700)\n",
    "plt.colorbar()\n",
    "plt.grid()\n",
    "plt.title('Scan Height')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plt.figure(figsize=(18,10))\n",
    "\n",
    "color=('black', 'gray', 'magenta', 'cyan', 'blue', 'red', 'yellow', 'darkgreen', 'green',\n",
    "      'violet', 'coral', 'darkblue', 'teal', 'steelblue', 'lightgray', 'darkcyan', 'orange')\n",
    "\n",
    "for plpl in range(all_lon.shape[0]):\n",
    "    #print plpl\n",
    "    plt.scatter(all_alt[plpl,...], all_kkk[plpl,...], color=color[plpl], s=0.1, marker='o', alpha=1,\n",
    "               edgecolors=None)\n",
    "\n",
    "plt.xlim(0,7000)\n",
    "plt.ylim(0, 2500)\n",
    "#plt.colorbar()\n",
    "plt.xlabel('Beam Altitude in m (without DEM)', fontsize=20)\n",
    "plt.ylabel('Beam Width in m', fontsize=20)\n",
    "plt.xticks(fontsize=20)\n",
    "plt.yticks(fontsize=20)\n",
    "\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# RADOLAN GRID LADEN\n",
    "print(wrl.georef.get_radolan_grid().shape)\n",
    "x = wrl.georef.get_radolan_grid(900,900)[:,:,0]\n",
    "y = wrl.georef.get_radolan_grid(900,900)[:,:,1]\n",
    "print (x.shape, y.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## Radolan tif\n",
    "## ----\n",
    "\n",
    "radolan_tif_path = '/automount/db01/python/data/DEM/radolan/radolan_900x900.tif'\n",
    "ds = wrl.io.open_raster(radolan_tif_path)\n",
    "ds = wrl.georef.reproject_raster_dataset(ds, spacing=1)\n",
    "cor_dem = wrl.georef.raster.read_gdal_coordinates(ds,mode='centers')\n",
    "val_dem = wrl.georef.raster.read_gdal_values(ds)\n",
    "val_dem[val_dem== -32768.]=0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "## RADOLAN Maske\n",
    "\n",
    "rm = np.load(\"/automount/ags/velibor/data/radolan_dx/RADOLAN_H_xystd_900x900.npy\")[2]\n",
    "rm = rm/rm\n",
    "print(np.ma.masked_invalid(np.unique(rm)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from scipy.spatial import cKDTree\n",
    "\n",
    "def ipoli_footprint(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": "markdown",
   "metadata": {},
   "source": [
    "### Interpolate dx on RADOLAN gird\n",
    "-------------------------------"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def ipoli_radi(gr_grid, gr_heig, sr_grid,radius, gr_widt, gr_dist):\n",
    "    \"\"\"\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    gr_grid ::: grid of the ground radar\n",
    "    gr_high ::: hight of the ground radar\n",
    "    gr_widt ::: ground radar beam width\n",
    "    gr_dist ::: distance of ground radar\n",
    "    sr_grid ::: grid oft the Spaceborn Radar\n",
    "    radius  ::: radius of the dpr foot prints\n",
    "\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    gs_grid ::: Interpolated Groundradar Data on Spaceborne Grid\n",
    "    \n",
    "    Heights: Max(H), Min(H), Mean(H), Std(H)\n",
    "    Widths : Max(W), Min(W), Mean(W), Std(W)\n",
    "    Distance : Max(D), Min(D), Mean(D), Std(D)\n",
    "\n",
    "    \"\"\"\n",
    "    \n",
    "    import sys\n",
    "    \n",
    "    ## Beam Hights\n",
    "    gr_ipoli_high_min = np.zeros((sr_grid.shape[0]))\n",
    "    gr_ipoli_high_max = np.zeros((sr_grid.shape[0]))\n",
    "    gr_ipoli_high_mean = np.zeros((sr_grid.shape[0]))\n",
    "    gr_ipoli_high_std = np.zeros((sr_grid.shape[0]))\n",
    "    \n",
    "    ## Beam Widths\n",
    "    gr_ipoli_widt_min = np.zeros((sr_grid.shape[0]))\n",
    "    gr_ipoli_widt_max = np.zeros((sr_grid.shape[0]))\n",
    "    gr_ipoli_widt_mean = np.zeros((sr_grid.shape[0]))\n",
    "    gr_ipoli_widt_std = np.zeros((sr_grid.shape[0]))\n",
    "    \n",
    "    ## Beam Distance\n",
    "    gr_ipoli_dist_min = np.zeros((sr_grid.shape[0]))\n",
    "    gr_ipoli_dist_max = np.zeros((sr_grid.shape[0]))\n",
    "    gr_ipoli_dist_mean = np.zeros((sr_grid.shape[0]))\n",
    "    gr_ipoli_dist_std = np.zeros((sr_grid.shape[0]))\n",
    "    \n",
    "\n",
    "    for i in range(sr_grid.shape[0]):\n",
    "        \n",
    "\n",
    "        ## Determine X0 Y0\n",
    "        x0, y0 = sr_grid[i,0], sr_grid[i,1]  \n",
    "        \n",
    "        ## Distance to X0 Y0\n",
    "        rr = np.sqrt((gr_grid[:,0] - x0)**2 + (gr_grid[:,1] - y0)**2)\n",
    "\n",
    "        if gr_data[rr < radius].size!=0:\n",
    "            \n",
    "            # Beam Height Min Max Mean Std\n",
    "            gr_ipoli_high_max[i] = np.nanmax(gr_high[rr < radius])\n",
    "            gr_ipoli_high_min[i] = np.nanmin(gr_high[rr < radius])\n",
    "            gr_ipoli_high_mean[i] = np.nanmean(gr_high[rr < radius])\n",
    "            gr_ipoli_high_std[i] = np.nanstd(gr_high[rr < radius])\n",
    "            # Bemwidth Min Max Mean Std\n",
    "            gr_ipoli_widt_max[i] = np.nanmax(gr_widt[rr < radius])\n",
    "            gr_ipoli_widt_min[i] = np.nanmin(gr_widt[rr < radius])\n",
    "            gr_ipoli_widt_mean[i] = np.nanmean(gr_widt[rr < radius])\n",
    "            gr_ipoli_widt_std[i] = np.nanstd(gr_widt[rr < radius])\n",
    "            # Distance Min Max Mean Std\n",
    "            gr_ipoli_dist_max[i] = np.nanmax(gr_widt[rr < radius])\n",
    "            gr_ipoli_dist_min[i] = np.nanmin(gr_widt[rr < radius])\n",
    "            gr_ipoli_dist_mean[i] = np.nanmean(gr_widt[rr < radius])\n",
    "            gr_ipoli_dist_std[i] = np.nanstd(gr_widt[rr < radius])\n",
    "            \n",
    "            \n",
    "        \n",
    "        else:\n",
    "            gr_ipoli_high_max[i] = np.nan\n",
    "            gr_ipoli_high_min[i] = np.nan\n",
    "            gr_ipoli_high_mean[i] = np.nan\n",
    "            gr_ipoli_high_std[i] = np.nan\n",
    "            \n",
    "            gr_ipoli_widt_mean[i] = np.nan\n",
    "            gr_ipoli_widt_min[i] = np.nan\n",
    "            gr_ipoli_widt_max[i] = np.nan\n",
    "            gr_ipoli_widt_std[i] = np.nan\n",
    "            \n",
    "            gr_ipoli_dist_mean[i] = np.nan\n",
    "            gr_ipoli_dist_min[i] = np.nan\n",
    "            gr_ipoli_dist_max[i] = np.nan\n",
    "            gr_ipoli_dist_std[i] = np.nan\n",
    "            \n",
    "\n",
    "\n",
    "    RES = [gr_ipoli_high_max, gr_ipoli_high_min, gr_ipoli_high_mean, gr_ipoli_high_std,\n",
    "           gr_ipoli_widt_max, gr_ipoli_widt_min, gr_ipoli_widt_mean, gr_ipoli_widt_std, \n",
    "           gr_ipoli_dist_max, gr_ipoli_dist_min, gr_ipoli_dist_mean, gr_ipoli_dist_std]\n",
    "    \n",
    "    \n",
    "    return RES"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import sys\n",
    "def loadingBar(count,total,size):\n",
    "    percent = float(count)/float(total)*100\n",
    "    sys.stdout.write(\"\\r\" + str(int(count)).rjust(3,'0')+\"/\"+str(int(total)).rjust(3,'0') + ' [' + '='*int(percent/10)*size + ' '*(10-int(percent/10))*size + ']')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def ipoli_radi(gr_grid, gr_high, sr_grid,radius, gr_widt, gr_dist):\n",
    "    \"\"\"\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    gr_grid ::: grid of the ground radar\n",
    "    gr_high ::: hight of the ground radar\n",
    "    gr_widt ::: ground radar beam width\n",
    "    gr_dist ::: distance of ground radar\n",
    "    sr_grid ::: grid oft the Spaceborn Radar\n",
    "    radius  ::: radius of the dpr foot prints\n",
    "\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    gs_grid ::: Interpolated Groundradar Data on Spaceborne Grid\n",
    "    \n",
    "    Heights: Max(H), Min(H), Mean(H), Std(H)\n",
    "    Widths : Max(W), Min(W), Mean(W), Std(W)\n",
    "    Distance : Max(D), Min(D), Mean(D), Std(D)\n",
    "\n",
    "    \"\"\"\n",
    "    \n",
    "    import sys\n",
    "    \n",
    "    ## Beam Hights\n",
    "    gr_ipoli_high = np.zeros((sr_grid.shape[0]))\n",
    "\n",
    "    ## Beam Widths\n",
    "    gr_ipoli_widt = np.zeros((sr_grid.shape[0]))\n",
    "\n",
    "    ## Beam Distance\n",
    "    gr_ipoli_dist = np.zeros((sr_grid.shape[0]))\n",
    "\n",
    "    \n",
    "\n",
    "    for i in range(sr_grid.shape[0]):\n",
    "        loadingBar(i,sr_grid.shape[0],2)\n",
    "        \n",
    "\n",
    "        ## Determine X0 Y0\n",
    "        x0, y0 = sr_grid[i,0], sr_grid[i,1]  \n",
    "        \n",
    "        ## Distance to X0 Y0\n",
    "        rr = np.sqrt((gr_grid[:,0] - x0)**2 + (gr_grid[:,1] - y0)**2)\n",
    "\n",
    "        if gr_high[rr < radius].size!=0:\n",
    "            \n",
    "            H = gr_high[rr < radius]\n",
    "            W = gr_widt[rr < radius]\n",
    "            D = gr_dist[rr < radius]\n",
    "            \n",
    "            \n",
    "            # Beam Height \n",
    "            gr_ipoli_high[i] = H[np.nanargmin(D)]\n",
    "\n",
    "            # Bemwidth \n",
    "            gr_ipoli_widt[i] = W[np.nanargmin(D)]\n",
    "\n",
    "            # Distance \n",
    "            gr_ipoli_dist[i] = D[np.nanargmin(D)]\n",
    "        \n",
    "        else:\n",
    "            gr_ipoli_high[i] = np.nan\n",
    "            \n",
    "            gr_ipoli_widt[i] = np.nan\n",
    "            \n",
    "            gr_ipoli_dist[i] = np.nan\n",
    "\n",
    "            \n",
    "\n",
    "\n",
    "    RES = [gr_ipoli_high, \n",
    "           gr_ipoli_widt,\n",
    "           gr_ipoli_dist]\n",
    "    \n",
    "    \n",
    "    return RES"
   ]
  },
  {
   "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": [
    "%%time\n",
    "\n",
    "#dx\n",
    "xy_dx = np.vstack((all_lon.ravel(), all_lat.ravel())).transpose()\n",
    "\n",
    "#radolan\n",
    "grid_xy = np.vstack((x.ravel(), y.ravel())).transpose()\n",
    "\n",
    "# interpolation\n",
    "RES  = ipoli_radi(xy_dx, all_alt.ravel(), grid_xy, 1.5, all_kkk.ravel(), all_dis.ravel())\n",
    "\n",
    "\n",
    "h_on_radolan = RES[0].reshape(x.shape)\n",
    "\n",
    "c_on_radolan = RES[1].reshape(x.shape)\n",
    "\n",
    "d_on_radolan = RES[2].reshape(x.shape)\n",
    "\n",
    "\n",
    "\"\"\"h_on_radolan_max = RES[0].reshape(x.shape)\n",
    "h_on_radolan_min = RES[1].reshape(x.shape)\n",
    "h_on_radolan_mean = RES[2].reshape(x.shape)\n",
    "h_on_radolan_std = RES[3].reshape(x.shape)\n",
    "\n",
    "c_on_radolan_max = RES[4].reshape(x.shape)\n",
    "c_on_radolan_min = RES[5].reshape(x.shape)\n",
    "c_on_radolan_mean = RES[6].reshape(x.shape)\n",
    "c_on_radolan_std = RES[7].reshape(x.shape)\n",
    "\n",
    "d_on_radolan_max = RES[8].reshape(x.shape)\n",
    "d_on_radolan_min = RES[9].reshape(x.shape)\n",
    "d_on_radolan_mean = RES[10].reshape(x.shape)\n",
    "d_on_radolan_std = RES[11].reshape(x.shape)\"\"\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%%time\n",
    "\n",
    "#dx\n",
    "xy_dx = np.vstack((all_lon.ravel(), all_lat.ravel())).transpose()\n",
    "\n",
    "#radolan\n",
    "grid_xy = np.vstack((x.ravel(), y.ravel())).transpose()\n",
    "\n",
    "# interpolation\n",
    "RES  = ipoli_footprint(xy_dx, grid_xy, all_alt.ravel(), dpr_footprint=1.5, k=30)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "h_on_radolan = RES[0].reshape(x.shape)\n",
    "\n",
    "c_on_radolan = RES[1].reshape(x.shape)\n",
    "\n",
    "d_on_radolan = RES[2].reshape(x.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(12,12))\n",
    "#plt.subplots_adjust(left=0.1, bottom=0.1, right=0.8, top=0.8, wspace=0.2, hspace=0.2)\n",
    "\n",
    "ax1 = fig.add_subplot(111, aspect='equal')\n",
    "plt.pcolormesh(x,y,np.ma.masked_invalid(c_on_radolan))\n",
    "plt.colorbar(shrink=0.5)\n",
    "plot_all_cband2(ax1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(12,12))\n",
    "#plt.subplots_adjust(left=0.1, bottom=0.1, right=0.8, top=0.8, wspace=0.2, hspace=0.2)\n",
    "\n",
    "ax1 = fig.add_subplot(111, aspect='equal')\n",
    "plt.pcolormesh(x,y,np.ma.masked_invalid(h_on_radolan), vmax=4000)\n",
    "plt.colorbar(shrink=0.5)\n",
    "plot_all_cband2(ax1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(16,16))\n",
    "#plt.subplots_adjust(left=0.1, bottom=0.1, right=0.8, top=0.8, wspace=0.2, hspace=0.2)\n",
    "\n",
    "ax1 = fig.add_subplot(111, aspect='equal')\n",
    "plt.pcolormesh(x,y,np.ma.masked_invalid(d_on_radolan))\n",
    "plt.colorbar()\n",
    "plot_all_cband2(ax1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plt.figure(figsize=(20,20))\n",
    "plt.pcolormesh(np.ma.masked_invalid(c_on_radolan))\n",
    "plt.colorbar()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plt.pcolormesh(np.ma.masked_invalid(res))\n",
    "plt.colorbar()\n",
    "plt.figure(figsize=(20,20))\n",
    "plt.pcolormesh(res, vmin=0)\n",
    "plt.colorbar()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "\n",
    "ry_h_min = ((np.ma.masked_invalid(h_on_radolan_min))* np.ma.masked_invalid(rm))- np.ma.masked_invalid(val_dem[::-1, :])\n",
    "\n",
    "ry_h_mean = ((np.ma.masked_invalid(h_on_radolan_mean))* np.ma.masked_invalid(rm))- np.ma.masked_invalid(val_dem[::-1, :])\n",
    "\n",
    "ry_h_max = ((np.ma.masked_invalid(h_on_radolan_max))* np.ma.masked_invalid(rm))- np.ma.masked_invalid(val_dem[::-1, :])\n",
    "\n",
    "ry_c_min = np.ma.masked_invalid(c_on_radolan_min) * np.ma.masked_invalid(rm)\n",
    "\n",
    "ry_c_mean =np.ma.masked_invalid(c_on_radolan_mean) * np.ma.masked_invalid(rm)\n",
    "\n",
    "ry_c_max =np.ma.masked_invalid(c_on_radolan_max) * np.ma.masked_invalid(rm)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## Save\n",
    "ry_h_min.dump('/automount/ags/velibor/data/radolan_dx/RY_H_min.npy')\n",
    "ry_h_mean.dump('/automount/ags/velibor/data/radolan_dx/RY_H_mean.npy')\n",
    "ry_h_max.dump('/automount/ags/velibor/data/radolan_dx/RY_H_max.npy')\n",
    "\n",
    "ry_c_min.dump('/automount/ags/velibor/data/radolan_dx/RY_C_min.npy')\n",
    "ry_c_mean.dump('/automount/ags/velibor/data/radolan_dx/RY_C_mean.npy')\n",
    "ry_c_max.dump('/automount/ags/velibor/data/radolan_dx/RY_C_max.npy')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "cc = 0.8\n",
    "plt.figure(figsize=(16,10))\n",
    "plt.subplot(2,3,1)\n",
    "plt.pcolormesh(ry_h_min, vmax=6000)\n",
    "cb = plt.colorbar(shrink=cc)\n",
    "cb.set_label(\"Min Beam Height in m\", fontsize=20)\n",
    "cb.ax.tick_params(labelsize=20)\n",
    "\n",
    "plt.subplot(2,3,2)\n",
    "plt.pcolormesh(ry_h_mean, vmax=6000)\n",
    "cb = plt.colorbar(shrink=cc)\n",
    "cb.set_label(\"Mean Beam Height in m\", fontsize=20)\n",
    "cb.ax.tick_params(labelsize=20)\n",
    "\n",
    "plt.subplot(2,3,3)\n",
    "plt.pcolormesh(ry_h_max, vmax=6000)\n",
    "cb = plt.colorbar(shrink=cc)\n",
    "cb.set_label(\"Max Beam Height in m\", fontsize=20)\n",
    "cb.ax.tick_params(labelsize=20)\n",
    "\n",
    "plt.subplot(2,3,4)\n",
    "plt.pcolormesh(ry_c_min, vmax=2000)\n",
    "cb = plt.colorbar(shrink=cc)\n",
    "cb.set_label(\"Min Beam Width in m\", fontsize=20)\n",
    "cb.ax.tick_params(labelsize=20)\n",
    "\n",
    "\n",
    "plt.subplot(2,3,5)\n",
    "plt.pcolormesh(ry_c_mean, vmax=2000)\n",
    "cb = plt.colorbar(shrink=cc)\n",
    "cb.set_label(\"Mean Beam Width in m\", fontsize=20)\n",
    "cb.ax.tick_params(labelsize=20)\n",
    "\n",
    "\n",
    "plt.subplot(2,3,6)\n",
    "plt.pcolormesh(ry_c_max, vmax=2000)\n",
    "cb = plt.colorbar(shrink=cc)\n",
    "cb.set_label(\"Max Beam Width in m\", fontsize=20)\n",
    "cb.ax.tick_params(labelsize=20)\n",
    "\n",
    "plt.tight_layout()\n",
    "#plt.savefig('/automount/ftp/velibor/gpm/nearest_4deci_mehrfach_')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "ry_h_mean"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plt.figure(figsize=(16,8))\n",
    "plt.subplot(1,2,1)\n",
    "plt.pcolormesh((np.ma.masked_invalid(h_on_radolan_max)-np.ma.masked_invalid(h_on_radolan_min)) \n",
    "               * np.ma.masked_invalid(rm), vmax=4000)\n",
    "plt.colorbar()\n",
    "\n",
    "\n",
    "\n",
    "plt.subplot(1,2,2)\n",
    "plt.pcolormesh((np.ma.masked_invalid(c_on_radolan_max)-np.ma.masked_invalid(c_on_radolan_min)) \n",
    "               * np.ma.masked_invalid(rm))\n",
    "plt.colorbar()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plt.pcolormesh(ry_h_min)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "ry_h_min.dump('/automount/ags/velibor/data/radolan_dx/RY_H_min.npy')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "a = np.load('/automount/ags/velibor/data/radolan_dx/RY_C_min.npy')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "a = np.ma.filled(a.astype(float), np.nan)\n",
    "\n",
    "plt.pcolormesh(a)\n",
    "plt.colorbar()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "np.save('/automount/ftp/velibor/RY_C_mean.npy', a)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def ipoli_radi_dem(gr_grid, gr_data,sr_grid,radius):\n",
    "    \"\"\"\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    gr_grid ::: grid of the Ground Radar\n",
    "    gr_data ::: data of the Ground Radar\n",
    "    sr_grid ::: grid oft the Spaceborn Radar\n",
    "    radius  ::: radius of the dpr foot prints\n",
    "\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    gs_grid ::: Interpolated Groundradar Data on Spaceborne Grid\n",
    "\n",
    "    \"\"\"\n",
    "\n",
    "    gr_ipoli_data_min = np.zeros((sr_grid.shape[0]))\n",
    "    gr_ipoli_data_max = np.zeros((sr_grid.shape[0]))\n",
    "\n",
    "    for i in range(sr_grid.shape[0]):\n",
    "\n",
    "        x0, y0 = sr_grid[i,0], sr_grid[i,1]  ###########x y richtig?\n",
    "\n",
    "        rr = np.sqrt((gr_grid[:,0] - x0)**2 + (gr_grid[:,1] - y0)**2)\n",
    "        ## Todo: hier fehlt die Wichtung\n",
    "        #print (gr_data[rr < radius])\n",
    "        #gr_ipoli_data[i] = np.nanmean(gr_data[rr < radius])\n",
    "        if gr_data[rr < radius].size!=0:\n",
    "\n",
    "            gr_ipoli_data_max[i] = np.nanmax(gr_data[rr < radius])\n",
    "            gr_ipoli_data_min[i] = np.nanmin(gr_data[rr < radius])\n",
    "        else:\n",
    "            gr_ipoli_data_max[i] = np.nan\n",
    "            gr_ipoli_data_min[i] = np.nan\n",
    "            # mein Pull Verfahren\n",
    "            # pull_range = 1 #km\n",
    "            # gr_ipoli_data[i] = np.nanmax(gr_data[rr < radius + pull_range])\n",
    "\n",
    "    return gr_ipoli_data_max, gr_ipoli_data_min"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def ipoli_radi_fast(xy_cut, gpm_xy,rwdata_cut, dpr_footprint, k=25):\n",
    "    \"\"\"\n",
    "    ...\n",
    "    \"\"\"\n",
    "    from scipy.spatial import cKDTree\n",
    "\n",
    "    # cKDTree radolan\n",
    "    tree = cKDTree(xy_cut)\n",
    "\n",
    "    # cKDTree gpm dpr\n",
    "    tree_gpm = cKDTree(gpm_xy)\n",
    "\n",
    "    dists, ix = tree.query(gpm_xy, k=k) # 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.nanmean(rwdata_cut.ravel()[ix[i,:][index]])\n",
    "    \n",
    "\n",
    "        ry_pns_w.append(res1)\n",
    "\n",
    "\n",
    "    return np.array(ry_pns_w)   "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def ipoli_radi_nearest(gr_grid, gr_data,sr_grid, gr_beam_width):\n",
    "    \"\"\"\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    gr_grid ::: grid of the Ground Radar\n",
    "    gr_data ::: data of the Ground Radar\n",
    "    sr_grid ::: grid oft the Spaceborn Radar\n",
    "    gr_beam_width ::: additional gr data like beam width\n",
    "\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    gr_grid ::: Interpolated Groundradar Data on Spaceborne Grid\n",
    "    gr_beam_width ::: Interpolated additional ground data on sr grid\n",
    "\n",
    "    \"\"\"\n",
    "\n",
    "    gr_ipoli_data = np.zeros((sr_grid.shape[0]))\n",
    "    gr_ipoli_beam_width = np.zeros((sr_grid.shape[0]))\n",
    "\n",
    "\n",
    "    c = 0\n",
    "    e = 0\n",
    "    for i in range(sr_grid.shape[0]):\n",
    "\n",
    "        x0, y0 = sr_grid[i,0], sr_grid[i,1] \n",
    "\n",
    "        rr = np.sqrt((gr_grid[:,0] - x0)**2 + (gr_grid[:,1] - y0)**2)\n",
    "\n",
    "        # an welchen stellen befinden sich die nahesten Punkte\n",
    "        idx = np.where(rr==np.nanmin(rr))\n",
    "        \n",
    "        # bei mehrfacher belegung nan setzen, weil Hoehe nicht eindeutig\n",
    "        if len(idx[0]) != 1:\n",
    "            gr_ipoli_data[i] = np.nan\n",
    "            gr_ipoli_beam_width[i] = np.nan\n",
    "            c = c + 1\n",
    "        # bei eindeutiger belegung Hoehe eindeutig\n",
    "        else:    \n",
    "            gr_ipoli_data[i] = gr_data[idx]\n",
    "            gr_ipoli_beam_width[i] = gr_beam_width[idx]\n",
    "            e = e + 1\n",
    "        \n",
    "        \n",
    "    print('Mehrfachbelegung: ', c)\n",
    "    print('Einfachbelegung: ', e)\n",
    "\n",
    "    \n",
    "    return gr_ipoli_data, gr_ipoli_beam_width"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "%%time\n",
    "\n",
    "#dx koordinates\n",
    "xy_dx = np.vstack((all_lon.ravel(), all_lat.ravel())).transpose()\n",
    "#radolan grid\n",
    "grid_xy = np.vstack((x.ravel(), y.ravel())).transpose()\n",
    "\n",
    "h_on_radolan, c_on_radolan   = ipoli_radi_nearest(xy_dx, all_alt.ravel(), grid_xy, all_kkk.ravel())\n",
    "h_on_radolan = h_on_radolan.reshape(x.shape)\n",
    "c_on_radolan = c_on_radolan.reshape(x.shape)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plt.figure(figsize=(12,6))\n",
    "plt.subplot(1,2,1)\n",
    "plt.pcolormesh(np.ma.masked_invalid(h_on_radolan) * np.ma.masked_invalid(rm))\n",
    "plt.colorbar()\n",
    "\n",
    "plt.subplot(1,2,2)\n",
    "plt.pcolormesh(np.ma.masked_invalid(c_on_radolan) * np.ma.masked_invalid(rm))\n",
    "plt.colorbar()\n",
    "#plt.savefig('/automount/ftp/velibor/gpm/nearest_4deci_mehrfach_')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "np.ma.masked_invalid(rm)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plt.pcolormesh(np.ma.masked_invalid(c_on_radolan) * np.ma.masked_invalid(rm))\n",
    "plt.colorbar()\n",
    "print(h_on_radolan.shape)\n",
    "#plt.savefig(\"/automount/ftp/velibor/radplot.png\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "rm = np.load(\"/automount/ags/velibor/data/radolan_dx/RADOLAN_H_xystd_900x900.npy\")[2]\n",
    "rm = rm/rm\n",
    "print(np.ma.masked_invalid(np.unique(rm)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# RADOLAN H without DEM\n",
    "rh_without_dem = np.ma.masked_invalid(h_on_radolan) * np.ma.masked_invalid(rm)\n",
    "rh_without_dem.dump('/automount/ags/velibor/data/radolan_dx/RADOLAN_H_XYNEAREST_without_DEM_900x900.npy')\n",
    "\n",
    "\n",
    "a = np.load('/automount/ags/velibor/data/radolan_dx/RADOLAN_H_XYNEAREST_without_DEM_900x900.npy')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# RADOLAN C without DEM\n",
    "rc_without_dem = np.ma.masked_invalid(c_on_radolan) * np.ma.masked_invalid(rm)\n",
    "\n",
    "rcwd = np.ma.filled(rc_without_dem .astype(float), np.nan)\n",
    "np.save('/automount/ftp/velibor/data/dem/rhwd',rcwd)\n",
    "\n",
    "#rc_without_dem.dump('/automount/ags/velibor/data/radolan_dx/RADOLAN_H_XYNEAREST_columns_900x900.npy')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "np.load('/automount/ags/velibor/data/radolan_dx/RADOLAN_H_XYNEAREST_columns_900x900.npy')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plt.figure(figsize=(12,12))\n",
    "plt.pcolormesh(a- np.ma.masked_invalid(val_dem[::-1, :]), vmin=0, vmax=4000)\n",
    "plt.colorbar()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "rh_with_dem = a- np.ma.masked_invalid(val_dem[::-1, :])\n",
    "\n",
    "rh_with_dem.dump('/automount/ags/velibor/data/radolan_dx/RADOLAN_H_XYNEAREST_with_DEM_900x900.npy')\n",
    "b = np.load('/automount/ags/velibor/data/radolan_dx/RADOLAN_H_XYNEAREST_with_DEM_900x900.npy')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "rhwd = np.ma.filled(rh_with_dem .astype(float), np.nan)\n",
    "np.save('/automount/ftp/velibor/data/dem/rhwd',rhwd)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "np.unique(rhwd)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plt.figure(figsize=(12,12))\n",
    "plt.pcolormesh(rhwd, vmin=0, vmax=4000)\n",
    "plt.colorbar()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "np.unique(rh_without_dem)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "## Save\n",
    "## ----\n",
    "\n",
    "#radolan_tif_path = '/automount/db01/python/data/DEM/radolan/radolan_1100x900.tif'\n",
    "#ds = wrl.io.open_raster(radolan_tif_path)\n",
    "#ds = wrl.georef.reproject_raster_dataset(ds, spacing=1)\n",
    "#cor_dem = wrl.georef.raster.read_gdal_coordinates(ds,mode='centers')\n",
    "#val_dem = wrl.georef.raster.read_gdal_values(ds)\n",
    "\n",
    "#np.save(\"/automount/ags/velibor/data/radolan_dx/RADOLAN_H_xymax_1100x900.npy\",[x,y, h_on_radolan_max])\n",
    "#np.save(\"/automount/ags/velibor/data/radolan_dx/RADOLAN_H_xymin_1100x900.npy\",[x,y, h_on_radolan_min])\n",
    "#np.save(\"/automount/ags/velibor/data/radolan_dx/RADOLAN_H_xymean_1100x900.npy\",[x,y, h_on_radolan_mean])\n",
    "#np.save(\"/automount/ags/velibor/data/radolan_dx/RADOLAN_H_xystd_1100x900.npy\",[x,y, h_on_radolan_std])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## Radolan tif\n",
    "## ----\n",
    "\n",
    "radolan_tif_path = '/automount/db01/python/data/DEM/radolan/radolan_900x900.tif'\n",
    "ds = wrl.io.open_raster(radolan_tif_path)\n",
    "ds = wrl.georef.reproject_raster_dataset(ds, spacing=1)\n",
    "cor_dem = wrl.georef.raster.read_gdal_coordinates(ds,mode='centers')\n",
    "val_dem = wrl.georef.raster.read_gdal_values(ds)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "rstd = np.load(\"/automount/ags/velibor/data/radolan_dx/RADOLAN_H_xystd_900x900.npy\")[2]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plt.imshow(rstd)\n",
    "plt.colorbar()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "-np.ma.masked_invalid(val_dem[::-1, :])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#h_on_radolan_max, h_on_radolan_min, h_on_radolan_mean, h_on_radolan_std\n",
    "plt.imshow(h_on_radolan_std)\n",
    "plt.colorbar()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Save\n",
    "#np.save(\"/automount/ags/velibor/data/radolan_dx/RY_maxHxy_b2017.npy\",[x,y, h_on_radolan_max])\n",
    "#np.save(\"/automount/ags/velibor/data/radolan_dx/RY_minHxy_b2017.npy\",[x,y, h_on_radolan_min])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "rhmax = np.load(\"/automount/ags/velibor/data/radolan_dx/RY_maxHxy_b2017.npy\")\n",
    "rhmin = np.load(\"/automount/ags/velibor/data/radolan_dx/RY_minHxy_b2017.npy\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from pcc import plot_borders\n",
    "from pcc import plot_radar\n",
    "from pcc import plot_radar2\n",
    "from pcc import plot_dem"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def plot_radar(radar, ax, proj):\n",
    "\n",
    "    site = (radar['lon'], radar['lat'], radar['alt'] )\n",
    "    \n",
    "    r = np.arange(1, 151) * 1000\n",
    "\n",
    "    # azimuth array 1 degree spacing\n",
    "    az = np.linspace(0, 360, 361)[0:-1]\n",
    "\n",
    "    # build polygons for maxrange rangering\n",
    "    polygons = wrl.georef.spherical_to_polyvert(r, az, 0,\n",
    "                                                site,\n",
    "                                                proj=proj)\n",
    "    polygons = polygons[..., 0:2]\n",
    "    polygons.shape = (len(az), len(r), 5, 2)\n",
    "    polygons = polygons[:, -1, :, :]\n",
    "\n",
    "    x_loc, y_loc = wrl.georef.reproject(site[0], site[1],\n",
    "                                        projection_source=proj_wgs,\n",
    "                                        projection_target=proj)\n",
    "\n",
    "\n",
    "    # create PolyCollections and add to respective axes\n",
    "    polycoll = mpl.collections.PolyCollection(polygons, closed=True,\n",
    "                                              edgecolors='k',\n",
    "                                              facecolors='k')\n",
    "    ax.add_collection(polycoll, autolim=True)\n",
    "\n",
    "    # plot radar location and information text\n",
    "    #ax.scatter(x_loc, y_loc, s=40, color='red', edgecolors='black')\n",
    "    #ax.text(x_loc, y_loc, radar['name'], color='k')\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def plot_all_cband2(ax):\n",
    "    for i in get_radar_locations().keys():\n",
    "\n",
    "        if i =='EMD':\n",
    "            pass\n",
    "        else:\n",
    "            plot_radar(get_radar_locations()[i], ax , proj_stereo)\n",
    "            \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# create radolan projection object\n",
    "proj_stereo = wrl.georef.create_osr(\"dwd-radolan\")\n",
    "\n",
    "# create wgs84 projection object\n",
    "proj_wgs = osr.SpatialReference()\n",
    "proj_wgs.ImportFromEPSG(4326)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Read DEM Data for RADOLAN\n",
    "-------------------------------"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "rhmax = np.load(\"/automount/ags/velibor/data/radolan_dx/RY_maxHxy_b2017.npy\")\n",
    "rhmin = np.load(\"/automount/ags/velibor/data/radolan_dx/RY_minHxy_b2017.npy\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "radolan_tif_path = '/automount/db01/python/data/DEM/radolan/radolan_900x900.tif'\n",
    "ds = wrl.io.open_raster(radolan_tif_path)\n",
    "ds = wrl.georef.reproject_raster_dataset(ds, spacing=1)\n",
    "cor_dem = wrl.georef.raster.read_gdal_coordinates(ds,mode='centers')\n",
    "val_dem = wrl.georef.raster.read_gdal_values(ds)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "cor_dem = wrl.georef.raster.read_gdal_coordinates(ds,mode='centers')\n",
    "val_dem = wrl.georef.raster.read_gdal_values(ds)\n",
    "print(cor_dem.shape, val_dem.shape)\n",
    "#cor_dem = cor_dem[:,0:900, :]\n",
    "#val_dem = val_dem[:,0:900]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "val_dem0 = val_dem.copy()\n",
    "val_dem[val_dem<0]=0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plt.figure(figsize=(15,5))\n",
    "plt.subplot(1,3,1)\n",
    "plt.imshow(val_dem[::-1, ::-1],  cmap=cm.terrain,vmin=0, vmax=1500)\n",
    "plt.colorbar(shrink=cc)\n",
    "plt.subplot(1,3,2)\n",
    "plt.imshow(np.ma.masked_invalid(h_on_radolan_min),  cmap=cm.jet,vmin=0, vmax=4000)\n",
    "plt.colorbar(shrink=cc)\n",
    "plt.subplot(1,3,3)\n",
    "plt.imshow(h_on_radolan_min- np.ma.masked_invalid(val_dem[::-1, :]),  cmap=cm.jet,vmin=0, vmax=4000)\n",
    "plt.colorbar(shrink=cc)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "a = val_dem[::-1, :]\n",
    "\n",
    "np.save(\"/automount/ftp/velibor/radolan_dem.npy\", a)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "a.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "cc=0.5\n",
    "def get_miub_cmap():\n",
    "    import matplotlib.colors as col\n",
    "    #startcolor = 'blue'  # a dark olive\n",
    "    #color1 = 'blue'\n",
    "    color3 = 'lightgreen'\n",
    "    color4 = 'green'\n",
    "    color5 = 'yellow'\n",
    "    color6 = 'orange'\n",
    "    color7 = 'brown'\n",
    "    color8 = 'saddlebrown'\n",
    "    #color6 = 'grey'\n",
    "    endcolor = 'white'    # medium dark red\n",
    "    colors = [ color3, color4, color5, color6, color7, color8, endcolor]\n",
    "    return col.LinearSegmentedColormap.from_list('miub1',colors)\n",
    "\n",
    "\n",
    "fig = plt.figure(figsize=(12,12))\n",
    "\n",
    "ax4 = fig.add_subplot(111, aspect='equal')\n",
    "\n",
    "cmap3 = get_miub_cmap()\n",
    "# extract all colors from the .jet map\n",
    "cmaplist3 = [cmap3(i) for i in range(cmap3.N)]\n",
    "# force the first color entry to be grey\n",
    "cmaplist3[0] = (.5,.5,.5,1.0)\n",
    "# create the new map\n",
    "cmap3 = cmap3.from_list('Custom cmap', cmaplist3, cmap3.N)\n",
    "\n",
    "# define the bins and normalize\n",
    "#bounds = np.linspace(0,5000,10)\n",
    "bounds3 = np.array([0,50,500,1200,1600,2400,2800,3400])\n",
    "norm3 = mpl.colors.BoundaryNorm(bounds3, cmap3.N)\n",
    "\n",
    "\n",
    "filename = wrl.util.get_wradlib_data_file('geo/radolan_900x900_cr_500.tif')\n",
    "# pixel_spacing is in output units (lonlat)\n",
    "ds = wrl.io.open_raster(filename)\n",
    "# pixel_spacing is in output units (lonlat)\n",
    "from matplotlib.colors import LogNorm\n",
    "ds = wrl.georef.reproject_raster_dataset(ds, spacing=1)\n",
    "rastervalues, rastercoords, proj = wrl.georef.extract_raster_dataset(ds)\n",
    "# specify kwargs for plotting, using terrain colormap and LogNorm\n",
    "dem = ax4.pcolormesh(rastercoords[..., 0], rastercoords[..., 1],\n",
    "                    val_dem0+8,  cmap=cmap3,\n",
    "                    vmin=0, vmax=1500)\n",
    "\n",
    "#plot_all_cband2(ax4)\n",
    "a = plot_borders(ax4)\n",
    "\n",
    "cb = plt.colorbar(dem, shrink=cc)\n",
    "cb.set_label(\"Height [m]\", fontsize=20)\n",
    "cb.ax.tick_params(labelsize=20)\n",
    "\n",
    "plt.grid()\n",
    "plt.tick_params(axis='both', which='both', bottom='off', top='off',\n",
    "                labelbottom='off',right='off', left='off', labelleft='off')\n",
    "#plt.axis('off')\n",
    "plt.xlim(-450,350)\n",
    "plt.ylim(-4650, -3770)\n",
    "\n",
    "plt.tight_layout()\n",
    "# pad\n",
    "# W-pad\n",
    "# H_pad\n",
    "#plt.subplots_adjust(hspace=0.1,wspace = 0.1 )\n",
    "#plt.savefig(\"/automount/ftp/velibor/RY_H2_maxmindiff_RADOLANNETWORK.png\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def get_borders():\n",
    "    # Landesgrenzen Deutschlands\n",
    "    from osgeo import osr\n",
    "    import wradlib as wrl\n",
    "    import wradlib\n",
    "    import numpy as np\n",
    "    \n",
    "    gxy = []\n",
    "    wgs84 = osr.SpatialReference()\n",
    "    wgs84.ImportFromEPSG(4326)\n",
    "    india = osr.SpatialReference()\n",
    "    # asia south albers equal area conic\n",
    "    india.ImportFromEPSG(102028)\n",
    "\n",
    "    proj_gk = osr.SpatialReference()\n",
    "    proj_gk.ImportFromEPSG(31466)\n",
    "    proj_ll = osr.SpatialReference()\n",
    "    proj_ll.ImportFromEPSG(4326)\n",
    "    gk3 = wradlib.georef.epsg_to_osr(31467)\n",
    "    proj_stereo = wrl.georef.create_osr(\"dwd-radolan\")\n",
    "    proj_wgs = osr.SpatialReference()\n",
    "    proj_wgs.ImportFromEPSG(4326)\n",
    "\n",
    "    # country list\n",
    "    countries = ['Germany']#,'France','Denmark', 'Netherlands', 'Poland']\n",
    "    # open the input data source and get the layer\n",
    "    filename = wradlib.util.get_wradlib_data_file('/automount/db01/python/data/NED/10m/cultural/10m_cultural/10m_cultural/ne_10m_admin_0_countries.shp')\n",
    "    dataset, inLayer = wradlib.io.open_vector(filename)\n",
    "    # iterate over countries, filter accordingly, get coordinates and plot\n",
    "    for item in countries:\n",
    "        #print item\n",
    "        # SQL-like selection syntax\n",
    "        fattr = \"(name='\"+item+\"')\"\n",
    "        inLayer.SetAttributeFilter(fattr)\n",
    "        # get borders and names\n",
    "        borders, keys = wradlib.georef.get_vector_coordinates(inLayer, key='name')\n",
    "\n",
    "        for j in range(borders.shape[0]):\n",
    "            bu = np.array(borders[j].shape)\n",
    "            a = np.array(bu.shape)\n",
    "\n",
    "            if a==1:\n",
    "                for i in range(0,borders[j].shape[0],1):\n",
    "                    bordx, bordy = wrl.georef.reproject(borders[j][i][:,0], borders[j][i][:,1], projection_source=proj_wgs, projection_target=proj_stereo)\n",
    "                    bord_xy = np.vstack((bordx.ravel(), bordy.ravel())).transpose()\n",
    "\n",
    "                    gxy.append(bord_xy)\n",
    "                    \n",
    "            if a==2:    #richtig\n",
    "                bordx, bordy = wrl.georef.reproject(borders[j][:,0], borders[j][:,1], projection_source=proj_wgs, projection_target=proj_stereo)\n",
    "                bord_xy = np.vstack((bordx.ravel(), bordy.ravel())).transpose()\n",
    "                gxy.append(bord_xy)\n",
    "\n",
    "            #bord_xy = np.vstack((bordx.ravel(), bordy.ravel())).transpose()\n",
    "\n",
    "            #gxy.append(bord_xy)\n",
    "            \n",
    "    return gxy    \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "gbx = []\n",
    "gby = []\n",
    "\n",
    "borders = get_borders()\n",
    "len(borders[20])\n",
    "\n",
    "for i in range(len(borders)):\n",
    "    print i \n",
    "    print borders[i][:,0]\n",
    "    gbx.append(borders[i][:,0])\n",
    "    gbx.append( np.nan)\n",
    "    gby.append( borders[i][:,1])\n",
    "    gby.append( np.nan)\n",
    "\n",
    "GX = np.array(gbx)\n",
    "GY = np.array(gby)\n",
    "GX = np.concatenate((GX),axis=None)\n",
    "GY = np.concatenate((GY),axis=None)\n",
    "\n",
    "Gxy = np.vstack((GX.ravel().T, GY.ravel().T))\n",
    "\n",
    "print(Gxy)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "np.save('/automount/ftp/velibor/data/german_boarders.npy',Gxy)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "Gxy.shape\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "cmap = plt.cm.jet\n",
    "# extract all colors from the .jet map\n",
    "cmaplist = [cmap(i) for i in range(cmap.N)]\n",
    "# force the first color entry to be grey\n",
    "cmaplist[0] = (.5,.5,.5,1.0)\n",
    "# create the new map\n",
    "cmap = cmap.from_list('Custom cmap', cmaplist, cmap.N)\n",
    "\n",
    "# define the bins and normalize\n",
    "#bounds = np.linspace(0,5000,10)\n",
    "bounds = np.arange(0,6000,400)\n",
    "norm = mpl.colors.BoundaryNorm(bounds, cmap.N)\n",
    "\n",
    "\n",
    "\n",
    "cc=0.6\n",
    "######################################################################################################################\n",
    "\n",
    "\n",
    "fig = plt.figure(figsize=(16,16))\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(223, aspect='equal')\n",
    "colmap = 'jet'\n",
    "plt.pcolormesh(x,y,np.ma.masked_invalid(h_on_radolan_max)-np.ma.masked_invalid(val_dem[::-1, :]),\n",
    "               vmin=0, vmax=6000, cmap=cmap, norm=norm)\n",
    "#plt.title('c) Maximum beam height', fontsize=20)\n",
    "plot_borders(ax1)\n",
    "plot_all_cband2(ax1)\n",
    "cb = plt.colorbar(shrink=cc)\n",
    "cb.set_label(\"Height [m]\", fontsize=20)\n",
    "cb.ax.tick_params(labelsize=20)\n",
    "plt.grid()\n",
    "#plt.axis('off')\n",
    "\n",
    "plt.tick_params(axis='both', which='both', bottom='off', top='off',\n",
    "                labelbottom='off',right='off', left='off', labelleft='off')\n",
    "plt.xlim(-450,350)\n",
    "plt.ylim(-4650, -3770)\n",
    "\n",
    "######################################################################################################################\n",
    "ax2 = fig.add_subplot(222, aspect='equal')\n",
    "plt.pcolormesh(x,y,np.ma.masked_invalid(h_on_radolan_min-np.ma.masked_invalid(val_dem[::-1, :])), vmin=0,vmax=6000\n",
    "               , cmap=cmap, norm=norm )\n",
    "#plt.title('b) Minimum beam height', fontsize=20)\n",
    "plot_borders(ax2)\n",
    "plot_all_cband2(ax2)\n",
    "\n",
    "plt.xlim(-450,350)\n",
    "plt.ylim(-4650, -3770)\n",
    "cb = plt.colorbar(shrink=cc)\n",
    "cb.set_label(\"Height [m]\", fontsize=20)\n",
    "cb.ax.tick_params(labelsize=20)\n",
    "\n",
    "plt.grid()\n",
    "#plt.axis('off')\n",
    "plt.tick_params(axis='both', which='both', bottom='off', top='off',\n",
    "                labelbottom='off',right='off', left='off', labelleft='off')\n",
    "\n",
    "######################################################################################################################\n",
    "\n",
    "cmap2 = plt.cm.rainbow\n",
    "# extract all colors from the .jet map\n",
    "cmaplist2 = [cmap2(i) for i in range(cmap2.N)]\n",
    "# force the first color entry to be grey\n",
    "cmaplist2[0] = (.5,.5,.5,1.0)\n",
    "# create the new map\n",
    "cmap2 = cmap2.from_list('Custom cmap', cmaplist2, cmap2.N)\n",
    "\n",
    "# define the bins and normalize\n",
    "#bounds = np.linspace(0,5000,10)\n",
    "bounds2 = np.arange(0,4000,300)\n",
    "norm2 = mpl.colors.BoundaryNorm(bounds2, cmap2.N)\n",
    "\n",
    "ax3 = fig.add_subplot(224, aspect='equal')\n",
    "plt.pcolormesh(x,y,np.ma.masked_invalid(h_on_radolan_max-h_on_radolan_min), vmin=0,vmax=4000,\n",
    "               cmap=cmap2, norm=norm2 )\n",
    "#plt.title('d) Maximum beam height range', fontsize=20)\n",
    "\n",
    "plot_borders(ax3)\n",
    "plot_all_cband2(ax3)\n",
    "cb = plt.colorbar(shrink=cc)\n",
    "cb.set_label(\"Height [m]\", fontsize=20)\n",
    "cb.ax.tick_params(labelsize=20)\n",
    "\n",
    "plt.grid()\n",
    "plt.tick_params(axis='both', which='both', bottom='off', top='off',\n",
    "                labelbottom='off',right='off', left='off', labelleft='off')\n",
    "#plt.axis('off')\n",
    "\n",
    "plt.xlim(-450,350)\n",
    "plt.ylim(-4650, -3770)\n",
    "\n",
    "\n",
    "######################################################################################################################\n",
    "ax4 = fig.add_subplot(221, aspect='equal')\n",
    "\n",
    "cmap3 = get_miub_cmap()\n",
    "# extract all colors from the .jet map\n",
    "cmaplist3 = [cmap3(i) for i in range(cmap3.N)]\n",
    "# force the first color entry to be grey\n",
    "cmaplist3[0] = (.5,.5,.5,1.0)\n",
    "# create the new map\n",
    "cmap3 = cmap3.from_list('Custom cmap', cmaplist3, cmap3.N)\n",
    "\n",
    "# define the bins and normalize\n",
    "#bounds = np.linspace(0,5000,10)\n",
    "bounds3 = np.array([0,50,500,1200,1600,2400,2800,3400])\n",
    "norm3 = mpl.colors.BoundaryNorm(bounds3, cmap3.N)\n",
    "\n",
    "\n",
    "filename = wrl.util.get_wradlib_data_file('geo/radolan_900x900_cr_500.tif')\n",
    "# pixel_spacing is in output units (lonlat)\n",
    "ds = wrl.io.open_raster(filename)\n",
    "# pixel_spacing is in output units (lonlat)\n",
    "from matplotlib.colors import LogNorm\n",
    "ds = wrl.georef.reproject_raster_dataset(ds, spacing=1)\n",
    "rastervalues, rastercoords, proj = wrl.georef.extract_raster_dataset(ds)\n",
    "# specify kwargs for plotting, using terrain colormap and LogNorm\n",
    "dem = ax4.pcolormesh(rastercoords[..., 0], rastercoords[..., 1],\n",
    "                    val_dem0+8,  cmap=cmap3 ,\n",
    "                    vmin=0, vmax=1500)\n",
    "\n",
    "plot_all_cband2(ax4)\n",
    "plot_borders(ax4)\n",
    "\n",
    "cb = plt.colorbar(dem, shrink=cc)\n",
    "cb.set_label(\"Height [m]\", fontsize=20)\n",
    "cb.ax.tick_params(labelsize=20)\n",
    "\n",
    "plt.grid()\n",
    "plt.tick_params(axis='both', which='both', bottom='off', top='off',\n",
    "                labelbottom='off',right='off', left='off', labelleft='off')\n",
    "#plt.axis('off')\n",
    "plt.xlim(-450,350)\n",
    "plt.ylim(-4650, -3770)\n",
    "\n",
    "plt.tight_layout()\n",
    "# pad\n",
    "# W-pad\n",
    "# H_pad\n",
    "#plt.subplots_adjust(hspace=0.1,wspace = 0.1 )\n",
    "#plt.savefig(\"/automount/ftp/velibor/RY_H2_maxmindiff_RADOLANNETWORK.png\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "cmap = plt.cm.jet\n",
    "# extract all colors from the .jet map\n",
    "cmaplist = [cmap(i) for i in range(cmap.N)]\n",
    "# force the first color entry to be grey\n",
    "cmaplist[0] = (.5,.5,.5,1.0)\n",
    "# create the new map\n",
    "cmap = cmap.from_list('Custom cmap', cmaplist, cmap.N)\n",
    "\n",
    "# define the bins and normalize\n",
    "#bounds = np.linspace(0,5000,10)\n",
    "bounds = np.arange(0,6000,400)\n",
    "norm = mpl.colors.BoundaryNorm(bounds, cmap.N)\n",
    "\n",
    "\n",
    "\n",
    "cc=0.6\n",
    "######################################################################################################################\n",
    "\n",
    "\n",
    "fig = plt.figure(figsize=(16,16))\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(223, aspect='equal')\n",
    "colmap = 'jet'\n",
    "plt.pcolormesh(rhmax[0],rhmax[1],np.ma.masked_invalid(rhmax[2])-np.ma.masked_invalid(val_dem[::-1, :]),\n",
    "               vmin=0, vmax=6000, cmap='jet', norm=norm)\n",
    "#plt.title('c) Maximum beam height', fontsize=20)\n",
    "plot_borders(ax1)\n",
    "plot_all_cband2(ax1)\n",
    "cb = plt.colorbar(shrink=cc)\n",
    "cb.set_label(\"Height [m]\", fontsize=20)\n",
    "cb.ax.tick_params(labelsize=20)\n",
    "plt.grid()\n",
    "#plt.axis('off')\n",
    "\n",
    "plt.tick_params(axis='both', which='both', bottom='off', top='off',\n",
    "                labelbottom='off',right='off', left='off', labelleft='off')\n",
    "plt.xlim(-450,350)\n",
    "plt.ylim(-4650, -3770)\n",
    "\n",
    "######################################################################################################################\n",
    "ax2 = fig.add_subplot(222, aspect='equal')\n",
    "plt.pcolormesh(rhmin[0],rhmin[1],np.ma.masked_invalid(rhmin[2]-np.ma.masked_invalid(val_dem[::-1, :])), vmin=-100,vmax=6000\n",
    "               , cmap='jet', norm=norm )\n",
    "#plt.title('b) Minimum beam height', fontsize=20)\n",
    "plot_borders(ax2)\n",
    "plot_all_cband2(ax2)\n",
    "\n",
    "plt.xlim(-450,350)\n",
    "plt.ylim(-4650, -3770)\n",
    "cb = plt.colorbar(shrink=cc)\n",
    "cb.set_label(\"Height [m]\", fontsize=20)\n",
    "cb.ax.tick_params(labelsize=20)\n",
    "\n",
    "plt.grid()\n",
    "#plt.axis('off')\n",
    "plt.tick_params(axis='both', which='both', bottom='off', top='off',\n",
    "                labelbottom='off',right='off', left='off', labelleft='off')\n",
    "\n",
    "######################################################################################################################\n",
    "\n",
    "cmap2 = plt.cm.rainbow\n",
    "# extract all colors from the .jet map\n",
    "cmaplist2 = [cmap2(i) for i in range(cmap2.N)]\n",
    "# force the first color entry to be grey\n",
    "cmaplist2[0] = (.5,.5,.5,1.0)\n",
    "# create the new map\n",
    "cmap2 = cmap2.from_list('Custom cmap', cmaplist2, cmap2.N)\n",
    "\n",
    "# define the bins and normalize\n",
    "#bounds = np.linspace(0,5000,10)\n",
    "bounds2 = np.arange(0,4000,300)\n",
    "norm2 = mpl.colors.BoundaryNorm(bounds2, cmap2.N)\n",
    "\n",
    "ax3 = fig.add_subplot(224, aspect='equal')\n",
    "plt.pcolormesh(rhmin[0],rhmin[1],np.ma.masked_invalid(rhmax[2]-rhmin[2]), vmin=0,vmax=4000,\n",
    "               cmap='jet', norm=norm2 )\n",
    "#plt.title('d) Maximum beam height range', fontsize=20)\n",
    "\n",
    "plot_borders(ax3)\n",
    "plot_all_cband2(ax3)\n",
    "cb = plt.colorbar(shrink=cc)\n",
    "cb.set_label(\"Height [m]\", fontsize=20)\n",
    "cb.ax.tick_params(labelsize=20)\n",
    "\n",
    "plt.grid()\n",
    "plt.tick_params(axis='both', which='both', bottom='off', top='off',\n",
    "                labelbottom='off',right='off', left='off', labelleft='off')\n",
    "#plt.axis('off')\n",
    "\n",
    "plt.xlim(-450,350)\n",
    "plt.ylim(-4650, -3770)\n",
    "\n",
    "\n",
    "######################################################################################################################\n",
    "ax4 = fig.add_subplot(221, aspect='equal')\n",
    "\n",
    "cmap3 = get_miub_cmap()\n",
    "# extract all colors from the .jet map\n",
    "cmaplist3 = [cmap3(i) for i in range(cmap3.N)]\n",
    "# force the first color entry to be grey\n",
    "cmaplist3[0] = (.5,.5,.5,1.0)\n",
    "# create the new map\n",
    "cmap3 = cmap3.from_list('Custom cmap', cmaplist3, cmap3.N)\n",
    "\n",
    "# define the bins and normalize\n",
    "#bounds = np.linspace(0,5000,10)\n",
    "bounds3 = np.array([0,50,500,1200,1600,2400,2800,3400])\n",
    "norm3 = mpl.colors.BoundaryNorm(bounds3, cmap3.N)\n",
    "\n",
    "\n",
    "filename = wrl.util.get_wradlib_data_file('geo/radolan_900x900_cr_500.tif')\n",
    "# pixel_spacing is in output units (lonlat)\n",
    "ds = wrl.io.open_raster(filename)\n",
    "# pixel_spacing is in output units (lonlat)\n",
    "from matplotlib.colors import LogNorm\n",
    "ds = wrl.georef.reproject_raster_dataset(ds, spacing=1)\n",
    "rastervalues, rastercoords, proj = wrl.georef.extract_raster_dataset(ds)\n",
    "# specify kwargs for plotting, using terrain colormap and LogNorm\n",
    "dem = ax4.pcolormesh(rhmin[0],rhmin[1],\n",
    "                    np.ma.masked_invalid(val_dem[::-1, :])+8,  cmap='jet' ,\n",
    "                    vmin=0, vmax=1500)\n",
    "\n",
    "plot_all_cband2(ax4)\n",
    "plot_borders(ax4)\n",
    "\n",
    "cb = plt.colorbar(dem, shrink=cc)\n",
    "cb.set_label(\"Height [m]\", fontsize=20)\n",
    "cb.ax.tick_params(labelsize=20)\n",
    "\n",
    "plt.grid()\n",
    "plt.tick_params(axis='both', which='both', bottom='off', top='off',\n",
    "                labelbottom='off',right='off', left='off', labelleft='off')\n",
    "#plt.axis('off')\n",
    "plt.xlim(-450,350)\n",
    "plt.ylim(-4650, -3770)\n",
    "\n",
    "plt.tight_layout()\n",
    "# pad\n",
    "# W-pad\n",
    "# H_pad\n",
    "#plt.subplots_adjust(hspace=0.1,wspace = 0.1 )\n",
    "#plt.savefig(\"/automount/ftp/velibor/RY_H2_maxmindiff_RADOLANNETWORK.png\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "cmap = plt.cm.jet\n",
    "# extract all colors from the .jet map\n",
    "cmaplist = [cmap(i) for i in range(cmap.N)]\n",
    "# force the first color entry to be grey\n",
    "cmaplist[0] = (.5,.5,.5,1.0)\n",
    "# create the new map\n",
    "cmap = cmap.from_list('Custom cmap', cmaplist, cmap.N)\n",
    "\n",
    "# define the bins and normalize\n",
    "#bounds = np.linspace(0,5000,10)\n",
    "bounds = np.arange(0,6000,400)\n",
    "norm = mpl.colors.BoundaryNorm(bounds, cmap.N)\n",
    "\n",
    "\n",
    "\n",
    "cc=0.6\n",
    "######################################################################################################################\n",
    "\n",
    "\n",
    "fig = plt.figure(figsize=(16,16))\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(223, aspect='equal')\n",
    "colmap = 'jet'\n",
    "plt.pcolormesh(rhmax[0],rhmax[1],np.ma.masked_invalid(h_on_radolan_max)-np.ma.masked_invalid(val_dem[::-1, :]),\n",
    "               vmin=0, vmax=6000, cmap=cmap, norm=norm)\n",
    "#plt.title('c) Maximum beam height', fontsize=20)\n",
    "plot_borders(ax1)\n",
    "plot_all_cband2(ax1)\n",
    "cb = plt.colorbar(shrink=cc)\n",
    "cb.set_label(\"Height [m]\", fontsize=20)\n",
    "cb.ax.tick_params(labelsize=20)\n",
    "plt.grid()\n",
    "#plt.axis('off')\n",
    "\n",
    "plt.tick_params(axis='both', which='both', bottom='off', top='off',\n",
    "                labelbottom='off',right='off', left='off', labelleft='off')\n",
    "plt.xlim(-450,350)\n",
    "plt.ylim(-4650, -3770)\n",
    "\n",
    "######################################################################################################################\n",
    "ax2 = fig.add_subplot(222, aspect='equal')\n",
    "plt.pcolormesh(rhmin[0],rhmin[1],np.ma.masked_invalid(h_on_radolan_min-np.ma.masked_invalid(val_dem[::-1, :])), vmin=0,vmax=6000\n",
    "               , cmap=cmap, norm=norm )\n",
    "#plt.title('b) Minimum beam height', fontsize=20)\n",
    "plot_borders(ax2)\n",
    "plot_all_cband2(ax2)\n",
    "\n",
    "plt.xlim(-450,350)\n",
    "plt.ylim(-4650, -3770)\n",
    "cb = plt.colorbar(shrink=cc)\n",
    "cb.set_label(\"Height [m]\", fontsize=20)\n",
    "cb.ax.tick_params(labelsize=20)\n",
    "\n",
    "plt.grid()\n",
    "#plt.axis('off')\n",
    "plt.tick_params(axis='both', which='both', bottom='off', top='off',\n",
    "                labelbottom='off',right='off', left='off', labelleft='off')\n",
    "\n",
    "######################################################################################################################\n",
    "\n",
    "cmap2 = plt.cm.rainbow\n",
    "# extract all colors from the .jet map\n",
    "cmaplist2 = [cmap2(i) for i in range(cmap2.N)]\n",
    "# force the first color entry to be grey\n",
    "cmaplist2[0] = (.5,.5,.5,1.0)\n",
    "# create the new map\n",
    "cmap2 = cmap2.from_list('Custom cmap', cmaplist2, cmap2.N)\n",
    "\n",
    "# define the bins and normalize\n",
    "#bounds = np.linspace(0,5000,10)\n",
    "bounds2 = np.arange(0,4000,300)\n",
    "norm2 = mpl.colors.BoundaryNorm(bounds2, cmap2.N)\n",
    "\n",
    "ax3 = fig.add_subplot(224, aspect='equal')\n",
    "plt.pcolormesh(rhmin[0],rhmin[1],np.ma.masked_invalid(rhmax[2]-rhmin[2]), vmin=0,vmax=4000,\n",
    "               cmap=cmap2, norm=norm2 )\n",
    "#plt.title('d) Maximum beam height range', fontsize=20)\n",
    "\n",
    "plot_borders(ax3)\n",
    "plot_all_cband2(ax3)\n",
    "cb = plt.colorbar(shrink=cc)\n",
    "cb.set_label(\"Height [m]\", fontsize=20)\n",
    "cb.ax.tick_params(labelsize=20)\n",
    "\n",
    "plt.grid()\n",
    "plt.tick_params(axis='both', which='both', bottom='off', top='off',\n",
    "                labelbottom='off',right='off', left='off', labelleft='off')\n",
    "#plt.axis('off')\n",
    "\n",
    "plt.xlim(-450,350)\n",
    "plt.ylim(-4650, -3770)\n",
    "\n",
    "\n",
    "######################################################################################################################\n",
    "ax4 = fig.add_subplot(221, aspect='equal')\n",
    "\n",
    "cmap3 = get_miub_cmap()\n",
    "# extract all colors from the .jet map\n",
    "cmaplist3 = [cmap3(i) for i in range(cmap3.N)]\n",
    "# force the first color entry to be grey\n",
    "cmaplist3[0] = (.5,.5,.5,1.0)\n",
    "# create the new map\n",
    "cmap3 = cmap3.from_list('Custom cmap', cmaplist3, cmap3.N)\n",
    "\n",
    "# define the bins and normalize\n",
    "#bounds = np.linspace(0,5000,10)\n",
    "bounds3 = np.array([0,50,500,1200,1600,2400,2800,3400])\n",
    "norm3 = mpl.colors.BoundaryNorm(bounds3, cmap3.N)\n",
    "\n",
    "\n",
    "filename = wrl.util.get_wradlib_data_file('geo/radolan_900x900_cr_500.tif')\n",
    "# pixel_spacing is in output units (lonlat)\n",
    "ds = wrl.io.open_raster(filename)\n",
    "# pixel_spacing is in output units (lonlat)\n",
    "from matplotlib.colors import LogNorm\n",
    "ds = wrl.georef.reproject_raster_dataset(ds, spacing=1)\n",
    "rastervalues, rastercoords, proj = wrl.georef.extract_raster_dataset(ds)\n",
    "# specify kwargs for plotting, using terrain colormap and LogNorm\n",
    "dem = ax4.pcolormesh(rastercoords[..., 0], rastercoords[..., 1],\n",
    "                    val_dem0+8,  cmap=cmap3 ,\n",
    "                    vmin=0, vmax=1500)\n",
    "\n",
    "plot_all_cband2(ax4)\n",
    "plot_borders(ax4)\n",
    "\n",
    "cb = plt.colorbar(dem, shrink=cc)\n",
    "cb.set_label(\"Height [m]\", fontsize=20)\n",
    "cb.ax.tick_params(labelsize=20)\n",
    "\n",
    "plt.grid()\n",
    "plt.tick_params(axis='both', which='both', bottom='off', top='off',\n",
    "                labelbottom='off',right='off', left='off', labelleft='off')\n",
    "#plt.axis('off')\n",
    "plt.xlim(-450,350)\n",
    "plt.ylim(-4650, -3770)\n",
    "\n",
    "plt.tight_layout()\n",
    "# pad\n",
    "# W-pad\n",
    "# H_pad\n",
    "#plt.subplots_adjust(hspace=0.1,wspace = 0.1 )\n",
    "#plt.savefig(\"/automount/ftp/velibor/RY_H2_maxmindiff_RADOLANNETWORK.png\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "0.05/0.02"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from functools import reduce\n",
    "import re\n",
    "import scipy\n",
    "from scipy.spatial import cKDTree\n",
    "from scipy.interpolate import LinearNDInterpolator\n",
    "from scipy.ndimage.interpolation import map_coordinates\n",
    "from scipy.interpolate import griddata\n",
    "import numpy as np\n",
    "import warnings\n",
    "\n",
    "#from . import util as util\n",
    "\n",
    "class IpolBase:\n",
    "    \"\"\"\n",
    "    IpolBase(src, trg)\n",
    "    The base class for interpolation in N dimensions.\n",
    "    Provides the basic interface for all other classes.\n",
    "    Parameters\n",
    "    ----------\n",
    "    src : ndarray of floats, shape (npoints, ndims)\n",
    "        Data point coordinates of the source points.\n",
    "    trg : ndarray of floats, shape (npoints, ndims)\n",
    "        Data point coordinates of the target points.\n",
    "    \"\"\"\n",
    "\n",
    "    def __init__(self, src, trg, **kwargs):\n",
    "        src = self._make_coord_arrays(src)\n",
    "        trg = self._make_coord_arrays(trg)\n",
    "        self.numsources = len(src)\n",
    "        self.numtargets = len(trg)\n",
    "\n",
    "    def __call__(self, vals):\n",
    "        \"\"\"\n",
    "        Evaluate interpolator for values given at the source points.\n",
    "        Parameters\n",
    "        ----------\n",
    "        vals : ndarray of float, shape (numsources, ...)\n",
    "            Values at the source points which to interpolate\n",
    "        Returns\n",
    "        -------\n",
    "        output : None\n",
    "        \"\"\"\n",
    "        self._check_shape(vals)\n",
    "        return None\n",
    "\n",
    "    def _check_shape(self, vals):\n",
    "        \"\"\"\n",
    "        Checks whether the values correspond to the source points\n",
    "        Parameters\n",
    "        ----------\n",
    "        vals : ndarray of float\n",
    "        \"\"\"\n",
    "        assert len(vals) == self.numsources, \\\n",
    "            ('Length of value array %d does not correspond to number '\n",
    "             'of source points %d' % (len(vals), self.numsources))\n",
    "        self.valsshape = vals.shape\n",
    "        self.valsndim = vals.ndim\n",
    "\n",
    "    def _make_coord_arrays(self, x):\n",
    "        \"\"\"\n",
    "        Make sure that the coordinates are provided as ndarray\n",
    "        of shape (numpoints, ndim)\n",
    "        Parameters\n",
    "        ----------\n",
    "        x : ndarray of float with shape (numpoints, ndim)\n",
    "            OR a sequence of ndarrays of float with len(sequence)==ndim and\n",
    "            the length of the ndarray corresponding to the number of points\n",
    "        \"\"\"\n",
    "        if type(x) in [list, tuple]:\n",
    "            x = [item.ravel() for item in x]\n",
    "            x = np.array(x).transpose()\n",
    "        elif type(x) == np.ndarray:\n",
    "            if x.ndim == 1:\n",
    "                x = x.reshape(-1, 1)\n",
    "            elif x.ndim == 2:\n",
    "                pass\n",
    "            else:\n",
    "                raise Exception('Cannot deal wih 3-d arrays, yet.')\n",
    "        return x\n",
    "\n",
    "    def _make_2d(self, vals):\n",
    "        \"\"\"Reshape increase number of dimensions of vals if smaller than 2,\n",
    "        appending additional dimensions (as opposed to the atleast_nd methods\n",
    "        of numpy).\n",
    "        Parameters\n",
    "        ----------\n",
    "        vals : ndarray\n",
    "               values who are to be reshaped to the right shape\n",
    "        Returns\n",
    "        -------\n",
    "        output : ndarray\n",
    "                 if vals.shape==() [a scalar] output.shape will be (1,1)\n",
    "                 if vals.shape==(npt,) output.shape will be (npt,1)\n",
    "                 if vals.ndim > 1 vals will be returned as is\n",
    "        \"\"\"\n",
    "        if vals.ndim < 2:\n",
    "            # ndmin might be 0 so we get it to 1-d first\n",
    "            # then we add an axis as we assume that\n",
    "            return np.atleast_1d(vals)[:, np.newaxis]\n",
    "        else:\n",
    "            return vals\n",
    "\n",
    "        \n",
    "        \n",
    "class Nearest(IpolBase):\n",
    "    \"\"\"\n",
    "    Nearest(src, trg)\n",
    "    Nearest-neighbour interpolation in N dimensions.\n",
    "    Parameters\n",
    "    ----------\n",
    "    src : ndarray of floats, shape (npoints, ndims) or cKDTree object\n",
    "        Data point coordinates of the source points.\n",
    "    trg : ndarray of floats, shape (npoints, ndims)\n",
    "        Data point coordinates of the target points.\n",
    "    remove_missing : int\n",
    "        Number of neighbours to consider in the presence of NaN, defaults to 0.\n",
    "    Keyword Arguments\n",
    "    -----------------\n",
    "    **kwargs : keyword arguments of ipclass (see class documentation)\n",
    "    Examples\n",
    "    --------\n",
    "    See :ref:`/notebooks/interpolation/wradlib_ipol_example.ipynb`.\n",
    "    Note\n",
    "    ----\n",
    "    Uses :class:`scipy:scipy.spatial.cKDTree`\n",
    "    \"\"\"\n",
    "\n",
    "    def __init__(self, src, trg, remove_missing=0, **kwargs):\n",
    "        if isinstance(src, cKDTree):\n",
    "            self.tree = src\n",
    "        else:\n",
    "            src = self._make_coord_arrays(src)\n",
    "            if len(src) == 0:\n",
    "                raise MissingSourcesError\n",
    "            # plant a tree, use unbalanced tree as default\n",
    "            kwargs.update(balanced_tree=kwargs.pop('balanced_tree', False))\n",
    "            self.tree = cKDTree(src, **kwargs)\n",
    "\n",
    "        self.numsources = self.tree.n\n",
    "\n",
    "        trg = self._make_coord_arrays(trg)\n",
    "        self.numtargets = len(trg)\n",
    "        if self.numtargets == 0:\n",
    "            raise MissingTargetsError\n",
    "\n",
    "        self.nnearest = remove_missing + 1\n",
    "\n",
    "        # query tree\n",
    "        self.dists, self.ix = self.tree.query(trg, k=self.nnearest)\n",
    "        # avoid bug, if there is only one neighbor at all\n",
    "        if self.dists.ndim == 1:\n",
    "            self.dists = self.dists[:, np.newaxis]\n",
    "            self.ix = self.ix[:, np.newaxis]\n",
    "\n",
    "    def __call__(self, vals, maxdist=None):\n",
    "        \"\"\"\n",
    "        Evaluate interpolator for values given at the source points.\n",
    "        You can interpolate multiple datasets of source values (``vals``) at\n",
    "        once: the ``vals`` array should have the shape (number of source\n",
    "        points, number of source datasets). If you want to interpolate only one\n",
    "        set of source values, ``vals`` can have the shape (number of source\n",
    "        points, 1) or just (number of source points,) - which is a flat/1-D\n",
    "        array. The output will have the same number of dimensions as ``vals``,\n",
    "        i.e. it will be a flat 1-D array in case ``vals`` is a 1-D array.\n",
    "        Parameters\n",
    "        ----------\n",
    "        vals : ndarray of float, shape (numsourcepoints, ...)\n",
    "            Values at the source points which to interpolate\n",
    "        maxdist : the maximum distance up to which an interpolated values is\n",
    "            assigned - if maxdist is exceeded, np.nan will be assigned\n",
    "            If maxdist==None, values will be assigned everywhere\n",
    "        Returns\n",
    "        -------\n",
    "        output : ndarray of float with shape (numtargetpoints,...)\n",
    "        \"\"\"\n",
    "        self._check_shape(vals)\n",
    "\n",
    "        # get first neighbour\n",
    "        trgvals = vals[self.ix[:, 0]]\n",
    "        dists = self.dists[..., 0].copy()\n",
    "\n",
    "        # iteratively fill NaN with next neighbours\n",
    "        isnan = np.isnan(trgvals)\n",
    "        nanidx = np.argwhere(isnan)[..., 0]\n",
    "        if self.nnearest > 1 & np.count_nonzero(isnan):\n",
    "            for i in range(self.nnearest - 1):\n",
    "                trgvals[isnan] = vals[self.ix[:, i + 1]][isnan]\n",
    "                dists[nanidx] = self.dists[..., i + 1][nanidx]\n",
    "                isnan = np.isnan(trgvals)\n",
    "                nanidx = np.argwhere(isnan)[..., 0]\n",
    "                if not np.count_nonzero(isnan):\n",
    "                    break\n",
    "\n",
    "        if maxdist is None:\n",
    "            return trgvals\n",
    "        else:\n",
    "            return np.where(dists > maxdist, np.nan, trgvals)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "nn(vals)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Synthetic observations and source coordinates\n",
    "src = np.vstack( (np.array([4, 7, 3, 15, -4,4]), np.array([8, 18, 17, 3, -8,8]))).transpose()\n",
    "vals = np.arange(len(src))\n",
    "vals = vals/vals\n",
    "# Target coordinates\n",
    "xtrg = np.linspace(-20, 20, 40)\n",
    "ytrg = np.linspace(-20, 20, 40)\n",
    "trg = np.meshgrid(xtrg, ytrg)\n",
    "trg = np.vstack( (trg[0].ravel(), trg[1].ravel()) ).T\n",
    "\n",
    "# Interpolation objects\n",
    "nn = Nearest(src, trg)\n",
    "\n",
    "# Subplot layout\n",
    "def gridplot(interpolated, title=\"\"):\n",
    "    pm = ax.pcolormesh(xtrg, ytrg, interpolated.reshape( (len(xtrg), len(ytrg)) ) )\n",
    "    plt.axis(\"tight\")\n",
    "    ax.scatter(src[:, 0], src[:, 1], facecolor=\"None\", s=50, marker='s')\n",
    "    plt.title(title)\n",
    "    plt.xlabel(\"x coordinate\")\n",
    "    plt.ylabel(\"y coordinate\")\n",
    "\n",
    "# Plot results\n",
    "fig = plt.figure(figsize=(8,8))\n",
    "\n",
    "ax = fig.add_subplot(111, aspect=\"equal\")\n",
    "gridplot(nn(vals), \"Nearest Neighbours\")\n",
    "\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "data = np.array([1,1,1,1,1,1,1,1,1,1,np.nan, np.nan])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "np.count_nonzero(~np.isnan(data))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "print('Distance: ')\n",
    "print(150* np.cos(np.deg2rad(0.5)), 150* np.cos(np.deg2rad(1.8)))\n",
    "print('Height: ')\n",
    "print(150* np.sin(np.deg2rad(0.5)), 150* np.sin(np.deg2rad(1.8)))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "H = np.array([np.nan, 1000, 2000, 1000])\n",
    "D = np.array([np.nan, 1, 2, 1])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plt.hist(all_dis.ravel(), bins=100)\n",
    "print('')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "all_dis.ravel()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "np.nanargmin(D)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "wrl_test [conda env:wrl_test]",
   "language": "python",
   "name": "conda-env-wrl_test-wrl_test"
  },
  "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.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}