vanice-histogram.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Mapping of registered lightning arroud wind-turbines"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "from datetime import datetime, timedelta\n",
    "from datetime import timezone\n",
    "#import datetime as dt\n",
    "import json\n",
    "from IPython.display import Image as ImageDisp\n",
    "from pandas import DataFrame\n",
    "import string\n",
    "import os\n",
    "import glob\n",
    "import matplotlib.pyplot as plt\n",
    "#%pylab inline --no-import-all\n",
    "from pandas.plotting import register_matplotlib_converters\n",
    "register_matplotlib_converters()\n",
    "import time\n",
    "os.environ['TZ'] = 'UTC'\n",
    "time.tzset()\n",
    "import matplotlib.colors as colors\n",
    "import json\n",
    "import matplotlib\n",
    "import statistics\n",
    "\n",
    "hdr_t = np.dtype([\n",
    "    ('magic', 'u4'),\n",
    "    ('id', 'u4'),\n",
    "    ('timestamp', 'u4'),\n",
    "    ('gps_flags', 'u1'),\n",
    "    ('gps_sat', 'u1'),\n",
    "    ('gps_lat', 'f4'),\n",
    "    ('gps_long', 'f4'),\n",
    "    ('gps_alt', 'u2'),\n",
    "    ('gps_hour', 'u1'),\n",
    "    ('gps_min', 'u1'),\n",
    "    ('gps_sec', 'u1'),\n",
    "    ('num_energy_data', 'u2'),\n",
    "], align=True)\n",
    "\n",
    "energy_data_t = np.dtype([\n",
    "    ('info', 'u1'),\n",
    "    ('timestamp_us', 'u4'),\n",
    "    ('energy', 'f4'),\n",
    "], align=True)\n",
    "\n",
    "def read_rt_dat_file(filename):\n",
    "    ret_dtype = np.dtype([\n",
    "        ('time', 'u4'),\n",
    "        ('time_us', 'u4'),\n",
    "        ('energy', 'f4'),\n",
    "        ('info', 'u1'),\n",
    "        ('gps_lat', 'f4'),\n",
    "        ('gps_long', 'f4'),\n",
    "        ('gps_alt', 'u2')\n",
    "    ])\n",
    "    \n",
    "    ret = np.zeros(0, dtype=ret_dtype)\n",
    "    pieces = []\n",
    "\n",
    "    with open(filename, 'r') as f:\n",
    "        while True:\n",
    "            try:\n",
    "                hdr = np.fromfile(f, dtype=hdr_t, count=1)[0]\n",
    "            except IndexError:\n",
    "                break\n",
    "            if hdr['magic'] != 0xaaaabbbb:\n",
    "                raise Exception('corrupted data, magic mismatch')\n",
    "                \n",
    "            energy_data = np.fromfile(f, dtype=energy_data_t, count=hdr['num_energy_data'])\n",
    "            piece = np.zeros(len(energy_data), dtype=ret_dtype)\n",
    "            piece['time'] = hdr['timestamp']\n",
    "            piece['time_us'] = energy_data['timestamp_us']\n",
    "            piece['energy'] = energy_data['energy']\n",
    "            piece['info'] = energy_data['info']\n",
    "            piece['gps_lat'] = hdr['gps_lat']\n",
    "            piece['gps_long'] = hdr['gps_long']\n",
    "            piece['gps_alt'] = hdr['gps_alt']\n",
    "            pieces.append(piece)\n",
    "    return np.concatenate(pieces)\n",
    "\n",
    "# Haversine formula for calculating the distance between two points on a sphere\n",
    "def haversine(lat1, lon1, lat2, lon2):\n",
    "    R = 6371.0  # Earth's radius in km\n",
    "    dlat = lat2 - lat1\n",
    "    dlon = lon2 - lon1\n",
    "    a = np.sin(dlat / 2)**2 + np.cos(lat1) * np.cos(lat2) * np.sin(dlon / 2)**2\n",
    "    c = 2 * np.arctan2(np.sqrt(a), np.sqrt(1 - a))\n",
    "    return R * c\n",
    "\n",
    "# Vectorizing the Haversine formula\n",
    "vec_haversine = np.vectorize(haversine)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 90,
   "metadata": {},
   "outputs": [],
   "source": [
    "#lat=51.7656581\n",
    "#lon=10.6106400\n",
    "#start=pd.Timestamp('2023-02-01 12:00') # Winter s dirou\n",
    "#end=pd.Timestamp('2023-02-02 2:00')\n",
    "\n",
    "lat=50.1006197\n",
    "lon=10.7797819\n",
    "\n",
    "\n",
    "start=pd.Timestamp('2023-11-8 2:00') # Only one stroke\n",
    "end=pd.Timestamp('2023-11-8 3:00')\n",
    "\n",
    "#start=pd.Timestamp('2023-6-5 13:00') \n",
    "#end=pd.Timestamp('2023-6-5 18:00')\n",
    "\n",
    "# Please comment not relevant time intervals\n",
    "start=pd.Timestamp('2023-01-13 12:00') # Winter\n",
    "end=pd.Timestamp('2023-01-13 14:00')\n",
    "\n",
    "#start=pd.Timestamp('2023-02-01 20:00') # Winter\n",
    "#end=pd.Timestamp('2023-02-02 2:00')\n",
    "\n",
    "#start=pd.Timestamp('2023-02-25 08:00') # Winter\n",
    "#end=pd.Timestamp('2023-02-25 15:00')\n",
    "\n",
    "#start=pd.Timestamp('2023-03-10 10:00') # Winter\n",
    "#end=pd.Timestamp('2023-03-10 22:00')\n",
    "\n",
    "#start=pd.Timestamp('2023-03-13 15:00') # Winter\n",
    "#end=pd.Timestamp('2023-03-13 22:00')\n",
    "\n",
    "#start=pd.Timestamp('2023-03-30 11:00') # Winter\n",
    "#end=pd.Timestamp('2023-03-30 19:00')\n",
    "\n",
    "start=pd.Timestamp('2023-03-30 11:00') \n",
    "end=pd.Timestamp('2023-03-30 19:00')\n",
    "\n",
    "#start=pd.Timestamp('2023-6-8 2:00') \n",
    "#end=pd.Timestamp('2023-6-8 20:00')\n",
    "\n",
    "#start=pd.Timestamp('2023-11-24 07:00') # Winter\n",
    "#end=pd.Timestamp('2023-11-24 18:00')\n",
    "\n",
    "#start=pd.Timestamp('2023-12-21 17:00') # Winter\n",
    "#end=pd.Timestamp('2023-12-21 23:00')\n",
    "\n",
    "#start=pd.Timestamp('2023-12-22 19:00') # Winter\n",
    "#end=pd.Timestamp('2023-12-22 23:50')\n",
    "\n",
    "\n",
    "sitename = str(start)\n",
    "radius_error = 1.4\n",
    "#radius_error = 5.6\n",
    "\n",
    "BLITZ_DIST = 150\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 91,
   "metadata": {},
   "outputs": [],
   "source": [
    "# set random azimuth angles (seed for reproducablity)\n",
    "npoints = 1000\n",
    "np.random.seed(1235)\n",
    "rand_azimuths_deg = np.random.random(npoints)*360\n",
    "\n",
    "import cartopy.geodesic as cgeo\n",
    "\n",
    "# set random distances (seed for reproducablity)\n",
    "np.random.seed(6341)\n",
    "rand_distances = BLITZ_DIST*1000*np.sqrt((np.random.random(npoints))) # np.random.uniform(low=0, high=radius, size=npoints)\n",
    "\n",
    "blitz = pd.DataFrame(columns=['lat', 'lon'])\n",
    "blitz.lon = cgeo.Geodesic().direct((lon,lat),rand_azimuths_deg,rand_distances)[:,0]\n",
    "blitz.lat = cgeo.Geodesic().direct((lon,lat),rand_azimuths_deg,rand_distances)[:,1]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 92,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load of lightning data\n",
    "\n",
    "blitz = pd.read_csv('./data/blitz_'+sitename+'bourka.csv', sep=',', low_memory=False, on_bad_lines='skip')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 93,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Reduction of distant lightnings\n",
    "\n",
    "blitz['lonRad']=np.radians(blitz.lon)\n",
    "blitz['latRad']=np.radians(blitz.lat)\n",
    "lonRad=np.radians(lon)\n",
    "latRad=np.radians(lat)\n",
    "blitz['DlonRad']=lonRad-blitz.lonRad\n",
    "blitz['DlatRad']=latRad-blitz.latRad\n",
    "blitz['a']=np.sin(blitz.DlatRad/2.)**2. + np.cos(blitz.latRad) * np.cos(latRad) * np.sin(blitz.DlonRad/2.)**2.\n",
    "blitz['km'] = 2. * np.arcsin(np.sqrt(blitz.a))*6378.\n",
    "blitz = blitz.loc[blitz['km'] <= BLITZ_DIST]\n",
    "blitz.reset_index(inplace=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "https://overpass-turbo.eu/#\n",
    "\n",
    "// @name tower\n",
    "\n",
    "/*\n",
    "This has been generated by the overpass-turbo wizard.\n",
    "The original search was:\n",
    "“ man_made = tower”\n",
    "*/\n",
    "[out:json][timeout:25];\n",
    "// gather results\n",
    "nwr[\"man_made\"=\"tower\"](48.7022575, 7.0337444,51.5979428, 13.5376508);\n",
    "// print results\n",
    "out center;\n",
    "\n",
    "\n",
    "// @name wind_turbine\n",
    "\n",
    "/*\n",
    "This has been generated by the overpass-turbo wizard.\n",
    "The original search was:\n",
    "“generator:method=wind_turbine”\n",
    "*/\n",
    "[out:json][timeout:25];\n",
    "// gather results\n",
    "nwr[\"generator:source\"=\"wind\"](48.7022575, 7.0337444,51.5979428, 13.5376508);\n",
    "// print results\n",
    "out center;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 94,
   "metadata": {},
   "outputs": [],
   "source": [
    "json_file = '/home/kakl/vanice/data/windturbine.geojson'\n",
    "\n",
    "with open(json_file) as f:\n",
    "    data = json.load(f)\n",
    "windturbines = pd.json_normalize(data,'features')\n",
    "\n",
    "windturbines['lon'] = windturbines['geometry.coordinates'].apply(lambda x: x[0][0] if isinstance(x[0], (list, np.ndarray)) else x).apply(lambda x: x[0])\n",
    "windturbines['lat'] = windturbines['geometry.coordinates'].apply(lambda x: x[0][0] if isinstance(x[0], (list, np.ndarray)) else x).apply(lambda x: x[1])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 95,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Reduction of wind-turbines points\n",
    "\n",
    "windturbines['lonRad']=np.radians(windturbines.lon)\n",
    "windturbines['latRad']=np.radians(windturbines.lat)\n",
    "lonRad=np.radians(lon)\n",
    "latRad=np.radians(lat)\n",
    "windturbines['DlonRad']=lonRad-windturbines.lonRad\n",
    "windturbines['DlatRad']=latRad-windturbines.latRad\n",
    "windturbines['a']=np.sin(windturbines.DlatRad/2.)**2. + np.cos(windturbines.latRad) * np.cos(latRad) * np.sin(windturbines.DlonRad/2.)**2.\n",
    "windturbines['km'] = 2. * np.arcsin(np.sqrt(windturbines.a))*6378.\n",
    "windturbines = windturbines.loc[windturbines['km'] <= BLITZ_DIST]\n",
    "windturbines.reset_index(inplace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 96,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Identify lightning with nearby structures\n",
    "\n",
    "# Converting coordinates to radians\n",
    "blitz['latRad'] = np.radians(blitz['lat'])\n",
    "blitz['lonRad'] = np.radians(blitz['lon'])\n",
    "windturbines['latRad'] = np.radians(windturbines['lat'])\n",
    "windturbines['lonRad'] = np.radians(windturbines['lon'])\n",
    "\n",
    "# Converting pandas.Series to numpy.array\n",
    "latRad_blitz = blitz['latRad'].to_numpy()\n",
    "lonRad_blitz = blitz['lonRad'].to_numpy()\n",
    "latRad_windturbines = windturbines['latRad'].to_numpy()\n",
    "lonRad_windturbines = windturbines['lonRad'].to_numpy()\n",
    "\n",
    "# Calculating distance between each point in 'blitz' and all points in 'structures'\n",
    "distances_wt = vec_haversine(latRad_blitz[:, np.newaxis], lonRad_blitz[:, np.newaxis],\n",
    "                          latRad_windturbines, lonRad_windturbines)\n",
    "distances_wt2 = vec_haversine(latRad_windturbines[:, np.newaxis], lonRad_windturbines[:, np.newaxis],\n",
    "                          latRad_blitz, lonRad_blitz)\n",
    "\n",
    "# Marking points in 'blitz' that are within 1.4 km of any point in 'structures'\n",
    "blitz['is_close'] = False\n",
    "blitz['is_close'] = np.any(distances_wt <= (radius_error), axis=1)\n",
    "windturbines['is_close'] = np.any(distances_wt2 <= (radius_error), axis=1)\n",
    "vetrniky = len(blitz.loc[blitz['is_close']==True])\n",
    "blitz['vetrnik'] = blitz['is_close']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 97,
   "metadata": {},
   "outputs": [],
   "source": [
    "json_file = '/home/kakl/vanice/data/tower.geojson'\n",
    "\n",
    "with open(json_file) as f:\n",
    "    data = json.load(f)\n",
    "towers = pd.json_normalize(data,'features')\n",
    "\n",
    "towers['lon'] = towers['geometry.coordinates'].apply(lambda x: x[0][0] if isinstance(x[0], (list, np.ndarray)) else x).apply(lambda x: x[0])\n",
    "towers['lat'] = towers['geometry.coordinates'].apply(lambda x: x[0][0] if isinstance(x[0], (list, np.ndarray)) else x).apply(lambda x: x[1])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 98,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Reduction of distant tall structures\n",
    "\n",
    "towers['lonRad']=np.radians(towers.lon)\n",
    "towers['latRad']=np.radians(towers.lat)\n",
    "lonRad=np.radians(lon)\n",
    "latRad=np.radians(lat)\n",
    "towers['DlonRad']=lonRad-towers.lonRad\n",
    "towers['DlatRad']=latRad-towers.latRad\n",
    "towers['a']=np.sin(towers.DlatRad/2.)**2. + np.cos(towers.latRad) * np.cos(latRad) * np.sin(towers.DlonRad/2.)**2.\n",
    "towers['km'] = 2. * np.arcsin(np.sqrt(towers.a))*6378.\n",
    "towers = towers.loc[towers['km'] <= BLITZ_DIST]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 99,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Identify lightning with nearby structures\n",
    "\n",
    "# Converting coordinates to radians\n",
    "blitz['latRad'] = np.radians(blitz['lat'])\n",
    "blitz['lonRad'] = np.radians(blitz['lon'])\n",
    "towers['latRad'] = np.radians(towers['lat'])\n",
    "towers['lonRad'] = np.radians(towers['lon'])\n",
    "\n",
    "# Converting pandas.Series to numpy.array\n",
    "latRad_blitz = blitz['latRad'].to_numpy()\n",
    "lonRad_blitz = blitz['lonRad'].to_numpy()\n",
    "latRad_towers = towers['latRad'].to_numpy()\n",
    "lonRad_towers = towers['lonRad'].to_numpy()\n",
    "\n",
    "# Calculating distance between each point in 'blitz' and all points in 'structures'\n",
    "distances_tw = vec_haversine(latRad_blitz[:, np.newaxis], lonRad_blitz[:, np.newaxis],\n",
    "                          latRad_towers, lonRad_towers)\n",
    "\n",
    "# Marking points in 'blitz' that are within 1.4 km of any point in 'structures'\n",
    "blitz['is_close'] = False\n",
    "blitz['is_close'] = np.any(distances_tw <= (radius_error), axis=1)\n",
    "veze = len(blitz.loc[blitz['is_close']==True])\n",
    "blitz['both'] = blitz['vetrnik'] | blitz['is_close']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 100,
   "metadata": {},
   "outputs": [],
   "source": [
    "json_file = '/home/kakl/vanice/data/mast.geojson'\n",
    "\n",
    "with open(json_file) as f:\n",
    "    data = json.load(f)\n",
    "masts = pd.json_normalize(data,'features')\n",
    "\n",
    "masts['lon'] = masts['geometry.coordinates'].apply(lambda x: x[0][0] if isinstance(x[0], (list, np.ndarray)) else x).apply(lambda x: x[0])\n",
    "masts['lat'] = masts['geometry.coordinates'].apply(lambda x: x[0][0] if isinstance(x[0], (list, np.ndarray)) else x).apply(lambda x: x[1])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 101,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Reduction of distant tall structures\n",
    "\n",
    "masts['lonRad']=np.radians(masts.lon)\n",
    "masts['latRad']=np.radians(masts.lat)\n",
    "lonRad=np.radians(lon)\n",
    "latRad=np.radians(lat)\n",
    "masts['DlonRad']=lonRad-masts.lonRad\n",
    "masts['DlatRad']=latRad-masts.latRad\n",
    "masts['a']=np.sin(masts.DlatRad/2.)**2. + np.cos(masts.latRad) * np.cos(latRad) * np.sin(masts.DlonRad/2.)**2.\n",
    "masts['km'] = 2. * np.arcsin(np.sqrt(masts.a))*6378.\n",
    "masts = masts.loc[masts['km'] <= BLITZ_DIST]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 102,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Identify lightning with nearby structures\n",
    "\n",
    "# Haversine formula for calculating the distance between two points on a sphere\n",
    "def haversine(lat1, lon1, lat2, lon2):\n",
    "    R = 6371.0  # Earth's radius in km\n",
    "    dlat = lat2 - lat1\n",
    "    dlon = lon2 - lon1\n",
    "    a = np.sin(dlat / 2)**2 + np.cos(lat1) * np.cos(lat2) * np.sin(dlon / 2)**2\n",
    "    c = 2 * np.arctan2(np.sqrt(a), np.sqrt(1 - a))\n",
    "    return R * c\n",
    "\n",
    "# Vectorizing the Haversine formula\n",
    "vec_haversine = np.vectorize(haversine)\n",
    "\n",
    "# Converting coordinates to radians\n",
    "blitz['latRad'] = np.radians(blitz['lat'])\n",
    "blitz['lonRad'] = np.radians(blitz['lon'])\n",
    "masts['latRad'] = np.radians(masts['lat'])\n",
    "masts['lonRad'] = np.radians(masts['lon'])\n",
    "\n",
    "# Converting pandas.Series to numpy.array\n",
    "latRad_blitz = blitz['latRad'].to_numpy()\n",
    "lonRad_blitz = blitz['lonRad'].to_numpy()\n",
    "latRad_masts = masts['latRad'].to_numpy()\n",
    "lonRad_masts = masts['lonRad'].to_numpy()\n",
    "\n",
    "# Calculating distance between each point in 'blitz' and all points in 'structures'\n",
    "distances_mt = vec_haversine(latRad_blitz[:, np.newaxis], lonRad_blitz[:, np.newaxis],\n",
    "                          latRad_masts, lonRad_masts)\n",
    "\n",
    "# Marking points in 'blitz' that are within 1.4 km of any point in 'structures'\n",
    "blitz['is_close'] = False\n",
    "blitz['is_close'] = np.any(distances_mt <= (radius_error), axis=1)\n",
    "stozary = len(blitz.loc[blitz['is_close']==True])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 103,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "json_file = '/home/kakl/vanice/data/chimney.geojson'\n",
    "\n",
    "with open(json_file) as f:\n",
    "    data = json.load(f)\n",
    "chimneys = pd.json_normalize(data,'features')\n",
    "\n",
    "chimneys['lon'] = chimneys['geometry.coordinates'].apply(lambda x: x[0][0] if isinstance(x[0], (list, np.ndarray)) else x).apply(lambda x: x[0])\n",
    "chimneys['lat'] = chimneys['geometry.coordinates'].apply(lambda x: x[0][0] if isinstance(x[0], (list, np.ndarray)) else x).apply(lambda x: x[1])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 104,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Reduction of distant tall structures\n",
    "\n",
    "chimneys['lonRad']=np.radians(chimneys.lon)\n",
    "chimneys['latRad']=np.radians(chimneys.lat)\n",
    "lonRad=np.radians(lon)\n",
    "latRad=np.radians(lat)\n",
    "chimneys['DlonRad']=lonRad-chimneys.lonRad\n",
    "chimneys['DlatRad']=latRad-chimneys.latRad\n",
    "chimneys['a']=np.sin(chimneys.DlatRad/2.)**2. + np.cos(chimneys.latRad) * np.cos(latRad) * np.sin(chimneys.DlonRad/2.)**2.\n",
    "chimneys['km'] = 2. * np.arcsin(np.sqrt(chimneys.a))*6378.\n",
    "chimneys = chimneys.loc[chimneys['km'] <= BLITZ_DIST]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 105,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Identify lightning with nearby structures\n",
    "\n",
    "# Haversine formula for calculating the distance between two points on a sphere\n",
    "def haversine(lat1, lon1, lat2, lon2):\n",
    "    R = 6371.0  # Earth's radius in km\n",
    "    dlat = lat2 - lat1\n",
    "    dlon = lon2 - lon1\n",
    "    a = np.sin(dlat / 2)**2 + np.cos(lat1) * np.cos(lat2) * np.sin(dlon / 2)**2\n",
    "    c = 2 * np.arctan2(np.sqrt(a), np.sqrt(1 - a))\n",
    "    return R * c\n",
    "\n",
    "# Vectorizing the Haversine formula\n",
    "vec_haversine = np.vectorize(haversine)\n",
    "\n",
    "# Converting coordinates to radians\n",
    "blitz['latRad'] = np.radians(blitz['lat'])\n",
    "blitz['lonRad'] = np.radians(blitz['lon'])\n",
    "chimneys['latRad'] = np.radians(chimneys['lat'])\n",
    "chimneys['lonRad'] = np.radians(chimneys['lon'])\n",
    "\n",
    "# Converting pandas.Series to numpy.array\n",
    "latRad_blitz = blitz['latRad'].to_numpy()\n",
    "lonRad_blitz = blitz['lonRad'].to_numpy()\n",
    "latRad_chimneys = chimneys['latRad'].to_numpy()\n",
    "lonRad_chimneys = chimneys['lonRad'].to_numpy()\n",
    "\n",
    "# Calculating distance between each point in 'blitz' and all points in 'structures'\n",
    "distances_ch = vec_haversine(latRad_blitz[:, np.newaxis], lonRad_blitz[:, np.newaxis],\n",
    "                          latRad_chimneys, lonRad_chimneys)\n",
    "\n",
    "# Marking points in 'blitz' that are within 1.4 km of any point in 'structures'\n",
    "blitz['is_close'] = False\n",
    "blitz['is_close'] = np.any(distances_ch <= (radius_error), axis=1)\n",
    "kominy = len(blitz.loc[blitz['is_close']==True])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 106,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2023-01-13 12:00:00 - 2023-01-13 14:00:00\n",
      "Nothing found\t 6 \t 16.0 %\n",
      "\n",
      "Storm date\tStrokes\t\tMedian/Mean of distance to neares tall structure\n",
      "\t\t\t\tWind-turbines\tTowers\t\tMasts\t\tChimneys\n",
      "2023-01-13 \t 38 \t\t 0.5 / 1.3 \t 1.7 / 1.8 \t 4.5 / 4.0 \t 8.0 / 7.8\n",
      "\n",
      "Storm date\tStrokes\t\tStrokes/percent of tall structures up to 1.4 km\n",
      "\t\t\t\tWind-turbines\tTowers\t\tMasts\t\tChimneys\n",
      "2023-01-13 \t 38 \t\t 31 / 82.0 %\t 7 / 18.0 %\t 3 / 8.0 %\t 0 / 0.0 %\n"
     ]
    }
   ],
   "source": [
    "celkem = len(blitz)\n",
    "bez = len(blitz.loc[blitz['both']==False])\n",
    "print(start, '-', end)\n",
    "median_wt = round(statistics.median(np.amin(distances_wt, axis=1)),1)\n",
    "median_tw = round(statistics.median(np.amin(distances_tw, axis=1)),1)\n",
    "median_mt = round(statistics.median(np.amin(distances_mt, axis=1)),1)\n",
    "median_ch = round(statistics.median(np.amin(distances_ch, axis=1)),1)\n",
    "mean_wt = round(statistics.mean(np.amin(distances_wt, axis=1)),1)\n",
    "mean_tw = round(statistics.mean(np.amin(distances_tw, axis=1)),1)\n",
    "mean_mt = round(statistics.mean(np.amin(distances_mt, axis=1)),1)\n",
    "mean_ch = round(statistics.mean(np.amin(distances_ch, axis=1)),1)\n",
    "\n",
    "print('Nothing found\\t', bez, '\\t', round(bez/(celkem/100),0),\"%\")\n",
    "\n",
    "print()\n",
    "print('Storm date\\tStrokes\\t\\tMedian/Mean of distance to neares tall structure')\n",
    "print('\\t\\t\\t\\tWind-turbines\\tTowers\\t\\tMasts\\t\\tChimneys')\n",
    "print(sitename.split()[0],'\\t',celkem,'\\t\\t', median_wt,'/', mean_wt,'\\t', median_tw,'/', mean_tw,'\\t', median_mt,'/', mean_mt,'\\t', median_ch,'/', mean_ch)\n",
    "\n",
    "print()\n",
    "print('Storm date\\tStrokes\\t\\tStrokes/percent of tall structures up to 1.4 km')\n",
    "print('\\t\\t\\t\\tWind-turbines\\tTowers\\t\\tMasts\\t\\tChimneys')\n",
    "print(sitename.split()[0],'\\t',celkem,'\\t\\t', vetrniky,'/', round(vetrniky/(celkem/100),0),'%\\t', veze,'/', round(veze/(celkem/100),0),'%\\t', stozary,'/',round(stozary/(celkem/100),0),'%\\t', kominy,'/', round(kominy/(celkem/100),0),'%')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 107,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAncAAAFbCAYAAAC3TKiyAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABYyklEQVR4nO3deXxU1fnH8c/Dvm8ioCCiguK+EBUREBA3cCktAmqr/FCpddcCVnEBlxa3qq1LxapolSJSrSBaERRUwAXcN1xRUQSRNRAQwvP749wJk8kkmSyThMn3/XrNK5lzz733uXNnJk/Oufccc3dEREREJDPUqOwARERERKT8KLkTERERySBK7kREREQyiJI7ERERkQyi5E5EREQkgyi5ExEREckgSu5EREREMoiSOxGpcsxsTzO73sxeN7OfzGydmb1rZqPNrGGS+nuZ2X/NbJWZrTezV82sTzls949mNtvMlprZpujny2Y2oBTHlFKMUd0rzexJM/vKzNzMFpd0fyXdTnkeq4hULtMgxiJS1ZjZOOACYCrwOrAZ6A0MAt4Hurp7TlR3D+BNYAtwJ7AGOBfYDzjB3WeWZrtR/SeAHOBjYAXQAjgVOAy41t1vSPF4Uo4xqu/ASuBtoAuw1t07pLKv0m6nvI5VRCqfkjsRqXLMLAv43N3XJJTfCIwGLnL3u6OyycBvgC7u/m5U1gj4CNgIdPboi64k2y0itlrAQmB3oJm756ZwPCnHGC3b3d2/in7/EGhUyuSuTNspzbGKSOVTt6yIVDnuviAxAYs8Ef3cDyDqSj0ZmB1LmqL1s4F/AnsCh5Z0u8XEtgX4HmgI1C6ufkljjJZ9Vdx2U1HW7ZT0WEWkaqhV2QGIiJRAu+jnsujnAUBdYH6Suq9HPw8ldImWZLv5mFkLoCbQktBVeTzwsrtvTCHm8oqxQpTxWEWkClByJyLbBTOrCVxLuG5tYlS8c/Tz+ySrxMralmK7iT4Ddoh+3wL8Bzg/pcDLIcYKVpZjFZEqQMmdiGwv7gS6Ale5+6KorEH0c1OS+hsT6pRku4l+DdQjJGGnAvWBJsBPKcRdHjFWpLIcq4hUAUruRKTKM7MbgAuB8e7+l7hFG6KfdZOsVi+hTkm2m4+7vxL39GEz+zfwmpnt4+6rota/HRNWy4mu7ytTjEUxs6aE5CveT2W58aG4Yy3tdkWk4uiGChGp0sxsDHA18DBwXsLiH6Kfybo1Y2XJukOL225xHgHaEFq5AHYBliY87iprjCm4K8l+dynltgqTeKwiUsWp5U5Eqiwzuw64DngUOCd+uJDIB4TuziOSrN41+rmgFNstTqy1rEX080fgmIQ6saSuVDGm6BbgsYSyH0u5rcIkHquIVHFK7kSkSjKza4ExwL+A/3P3rYl13D3bzKYBvzazA939vWjdRsA5wOck3IWaynajeg0JY4FmJ5TXJAyEDNHdrtGdpDNJojQxpsrdPyYMOlwmJTlWEan6NIixiFQ5ZnYBcDfwLXANkJiALXP3F6O6HQnJ0WbgDmAtYfaH/YH+7v5CKbd7EDAHmAIsIsz00BY4DdgLeMTdh6Z4PCnHGNX/HbBr9PQioA5we/T8G3f/V4r7TWk75XmsIlL5lNyJSJVjZhOAs4qoMsfde8XV3xsYBxxFSGDeBsYkmdYr5e2aWUtgLNCdcB1bY8K0Ye8AE4CJJenOTTXGqO7sqF6RMaawz5S2U97HKiKVq8KTu2g6mxHA2UB7wu31T7r7ZXF1DLgS+ANhIM23gIvjR3c3s86EC307A88TrpvJjlveE/g3sFdiV4OIiIhIpqqMu2UfBi4GbgOOBf5EmKw63p8IXSY3AycB2cBMM2sTV2cC8AVhwu99gKtiC8ysBmHsqiuV2ImIiEh1UqEtd2Z2PDANODC6EDhZnXqEKYBud/fro7KGwGLgfne/OroQeR3Qyt1/MrPBwAh3PzSqfy7hQuWu6koQERGR6qSiW+6GAS8VlthFuhFGQ58cK3D39YSk8ISoqE70M9bityFWZmZNgBuAS5TYiYiISHVT0cnd4cBnZna3ma01sw1m9pSZ7RxXpzOQSxgeIN4n0TLcfSXwNXBRNMn1cLaNE3UNMNPdddu+iIiIVDsVndy1AYYCBwFDgP8DugBPRzdRADQHspNMn7MKaGBmsVa7C4DRwM+EW/XHRsMNnE24Zk9ERESk2qnoQYwtepzi7j8DmNlSwvhKfYBZUb1k3akWv8zdnzezVkA74Et3zzWzqcAd7r4kGs/qimidce5+b9KAzIYTWv6oX79+l112Ke+Ze/LbunUrNWpo1reqRuelatJ5qZp0XqoenZOqKd3n5bPPPlvh7onzWld4crcK+CqW2EVeA34h3PE6K6rT2MxqJrTeNQM2uPvmWIG7bwA+AzCzvsCBwGAzO5Bw3V23qOp8M3vN3d9PDMjdxwPjAbKysnzBgtLOApSa2bNn06tXr7TuQ0pO56Vq0nmpmnReqh6dk6op3efFzL5JVl7Raf4nhZQb20aK/xSoCXRMqNM5WlZw5TB23p3AKHfPAXoRbtz41N0/JSSNhQ3kKSIiIpIxKjq5exY4IBoNPaYnUBt4L3o+jzA1z6mxCmbWgDDe3fOFbPc8YJW7PxFX1iDu94Zs69YVERERyVgV3S07njCA8TQz+zNhipubCXe3vgZhAm4zGwdcY2arCK11lxMS0b8nbtDMmgPXAcfFFb8C3GJmw6LnfdBNFiIiIlINVGhy5+5rzawP8DdgEuFau2eAyxKqjiMkc1cCOxCGOTnG3Zcl2exYYKq7vx23n3fMbBRwU1Q0wt3fS7KuiIiISEap6JY73P0LoF8xdZyQmN1UVL2o7sWFlN8F3FWaGEVERES2V7pvWkRERCSDVHjLnYiIyPZk7dq1LF++nM2bNxdfuZI0bdqUTz4pbEAKqSylPS+1a9emVatWNGnSpFT7VXInIiJSiLVr17Js2TLatm1L/fr12TaZUtWybt06GjduXNlhSILSnBd3Jycnh++//x6gVAmeumVFREQKsXz5ctq2bUuDBg2qbGInmcXMaNCgAW3btmX58uWl2oaSOxERkUJs3ryZ+vXrV3YYUg3Vr1+/1JcCKLkTEREpglrspDKU5X2n5E5EREQkgyi5ExERyWAPPvggZsaSJUvylV9xxRWYGY899li+8hdffBEzY968eZgZzz77bLnEkZWVxdChQ4uss3z5csaMGcPixYvLZZ8xY8aMoWXLlsXW69ChAyNGjCjXfVcGJXciIiIlZGaV8iiNbt26ATBv3rx85fPmzaNBgwZJy+vWrUuXLl2YP38+3bt3L92LVArLly9n7Nix5Z7cperpp5/m4ouTzo2wXVFyJyIiksE6d+5MixYt8iVxmzdvZuHChZx55plJk7usrCzq1q1L165dadasWQVHXH42b95Mbm5uyvUPPvhg2rdvn8aIKoaSOxERkQxmZhxxxBH5krh33nkHd+f888/ngw8+YN26dQBs3bqVN954g27durF48eIC3bKxbss77riDdu3a0bx5c4YMGcLq1avz7fPDDz/kyCOPpF69euy9995MnTq12DgXL17M/vvvD0Dv3r3ztVZOmDABMyM7OzvfOondqL169WLgwIGMHz+ePfbYg3r16vHDDz/kLZ87dy6HHHII9erV46CDDuK1114rcntDhw4lKyuLF198kQMOOICGDRvSvXt3Pvroo3zrbd26lXHjxtGxY0fq1q3LnnvuySOPPJKvzmuvvUaPHj1o0qQJTZo04aCDDuLJJ58s9nUpDSV3IiIiGa5bt268++675OTkADB//ny6dOnCfvvtR7NmzXjjjTcA+Oijj1izZg1HHnlkoduaPHkys2bNYvz48dx88808++yzXHXVVXnLc3JyOO6448jOzmbixIlcffXVXHrppXz77bdFxrjTTjvx+OOPA3DPPfcwf/585s+fX+JjnTt3Lvfddx8333wz06ZNo2nTpgBs2LCB3/72t5x33nk8+eSTNGvWjBNOOIEff/yxyO19++23jBw5ktGjR/Pvf/+b5cuXM2jQINw9r85FF13EjTfeyPDhw5k+fToDBgxg2LBhPP/880AYDPvEE09k99135z//+Q9Tpkzhd7/7XYGkuLxohgoREZEM161bNzZv3sxbb71Fz549mTdvHkcccQRmRteuXZk3bx59+/bNa93r1q0b69evT7qt2rVr89///pdatUIK8fHHHzNp0iTuvfdeAB5++GGWL1/OG2+8Qbt27YDQIlbctXt169blgAMOAGCfffaha9eupTrW1atX884779CmTZt85Tk5Odx0002cfvrpQGgdbN++PXfeeSfjxo0rdHsrV65k7ty5dOrUCQitdAMGDGDRokV07tyZL774gvvuu4+HH36Ys846C4C+ffuydOlSxo0bx6BBg/jss89Ys2YNd999d96MFccee2ypji8VarkTERHJcIcddhi1atXKS95iyR2Ql9zFyjt16sSOO+5Y6LZ69+6dl9hBSMSWL1/OL7/8AsCbb75Jly5d8hI7gCOPPJJWrVrlPXd3tmzZkvcoyXVxxenSpUuBxC5mwIABeb83atSIY445hjfffLPI7XXo0CEvsYNwvEDe3cezZs2iRo0aDBgwIN8xHX300XzwwQfk5uayxx570KhRI04//XSeeeaZtLXYxSi5ExERyXANGjTgoIMOYt68eSxZsoQlS5bkJXdHHHEEr7/+Olu3bmXevHlFdskCBW6wqFOnDu6el9z9+OOP+RK5mPiyOXPmULt27bzH0UcfXcYj3KZ169ZJyxs1alRgtpFWrVqxdOnSIreX7HgBNm7cCMCKFSvIzc2ladOm+Y5p6NChbNmyhaVLl9K8eXNmzJjB5s2bGTRoEDvuuCP9+/fnq6++KuVRFk3dsiIiItVAt27dmDhxIvPmzaNDhw7stNNOABx++OGsW7eOOXPm8MUXX3DFFVeUaT9t2rTh008/LVAeP09qly5deOutt/Kex7oqC1OvXj2AvAQyZtWqVQXqFjZkTHZ2Njk5OfkSvOXLl+e9DqXVokULatWqxdy5c6lRI3+b2fr16/OS2iOOOIL//e9/5OTkMHPmTC6//HJOP/10Xn/99TLtPxm13ImIiFQDRx55JCtWrOCRRx7Ja7WDkFjtu+++3HbbbcC2cfFK69BDD2XhwoX5Bk2eO3duvuSucePGZGVl5T322msvoGCrWEysi/eTTz7JK3vjjTdYu3ZtiWJ7+umn837Pzs7mxRdf5LDDDivRNhL16dOH3Nxc1qxZk++YsrKyOOSQQ/KOKaZ+/fqcdNJJDBs2jI8//rhM+y6MWu5ERESqgVh36/PPP89dd92Vb9kRRxzBAw88QPPmzdl7773LtJ//+7//48Ybb6R///6MGTOGnJwcrrnmmpRmiGjfvj3169fnkUceyevmzMrK4rDDDqNt27ZcfPHF3HDDDaxcuZJbbrmFJk2apBxX/fr1GT16NNnZ2ey8887cdttt/PLLL1xyySVlOVz22msvzjvvPIYMGcKoUaPIyspi48aNfPTRR3z00Uc88sgjTJ8+nYceeohf/epXtG/fnu+//57777+fPn36lGnfhVHLnYiISDXQtm1b2rdvj7vna7mDkNzFyssyYT2E6/teeOEFGjZsyJAhQxg7diy33347u+66a7Hr1qtXjwceeICFCxdy1FFHceihhwKhRe/pp5+mRo0aDBw4kNtvv5377ruP5s2blyiuRx99lHvvvZff/OY3rFq1iueee67M3bIQhm655pprePTRR+nXrx9Dhw5l+vTpea2gHTt2xMy46qqrOPbYYxk1ahTHH388Dz30UJn3nYzFj9NS3WVlZfmCBQvSuo/Zs2fTq1evtO5DSk7npWrSeamaqtN5+eSTT8rcklUR1q1bV+x1a1Lxynpeinv/mdlCd89KLFfLnYiIiEgGUXInIiIikkGU3ImIiIhkECV3IiIiIhlEyZ2IiIhIBlFyJyIiIpJBlNyJiIiIZBAldyIiIiIZRMmdiIiISAZRciciIpLhxowZg5nRqVOnpMtj02ONGTOm3PY5Y8YM7rzzznLbnqROyZ2IiEg1UK9ePb7++msSp9l86623+Oabb6hXr1657k/JXeWpVdkBiIiIbHcmWuXs9/TSzwffsGFDDjnkECZNmkRW1rbpSCdNmkSfPn1YuHBheUQoVUCFttyZ2VAz8ySP8+LqmJldZWbfmVmOmb1iZgclbKezmb1hZmvMbJKZNUpY3tPMvk8sFxERqc6GDBnC5MmTcQ9JorszefJkhgwZkq/e/PnzOfnkk9l5551p2LAhBx10EI8//ni+OqtXr+acc85h5513pl69erRv355zzz0XCN3At99+O9988w1mhpkxdOhQAD766COOP/54WrRoQcOGDdl7772555570n/w1Uhltdz1AXLinn8V9/ufgGuAkcCnwOXATDPbz91/jOpMAL4ArgVuBa6KHphZDeBO4Ep3z07fIYiIiGxffv3rX/OHP/yB1157jR49evDqq6/y008/MWDAAEaOHJlX75tvvuHII4/kvPPOo169esydO5f/+7//o0aNGpx22mkAXH755cybN4877riDNm3a8N133/HKK68AcM455/D555/z0ksv8fTTTwOw4447AnDyySfTuXNnHnvsMerWrcuiRYtYu3ZtBb8Sma2ykru3kiVeZlaPkNz9xd3vjsrmA4uBC4Gro9a4w4GT3P0nM2sGjCBK7oCzgc3Av9J9ECIiItuTZs2acfzxxzNp0iR69OjBpEmTOP7442nWrFm+evEtee5Oz549WbJkCQ888EBecvfmm29ywQUXMHjw4Ly6v/3tbwFo164dO+20E3Xr1qVr1655y1esWMFXX33Ff//7X/bff38Ajj766HQdbrVV1W6o6AY0ASbHCtx9PTANOCEqqhP9jLX8bYiVmVkT4AbgEo+1OYuIiEieIUOGMGXKFDZt2sSUKVMKdMkCrFq1iosvvphdd92V2rVrU7t2bcaPH89nn32WV+eggw7i1ltv5d57781XXpQWLVqwyy67cN555/HEE0+wfPnycjsu2aaykrsvzWyLmS0ys9/HlXcGcoHPE+p/Ei3D3VcCXwMXmVkLYDgQu/XnGmCmu7+e1uhFRES2UyeffDLZ2dmMHj2a9evXc9JJJxWoM3ToUJ544glGjhzJjBkzeOuttxg2bBgbN27Mq3P33Xfzq1/9iuuvv5699tqLTp06MWnSpCL3XaNGDWbMmEGbNm0YNmwYbdq0oUePHrzzzjvlfpzVWUUnd0sJCdjvgJOAN4B/mNll0fLmQLa75yastwpoYGaxVrsLgNHAz8BewFgz60jokv1Teg9BRERk+9WwYUNOPPFE7rjjDk466SQaNmyYb/nGjRuZPn06Y8eO5cILL6RPnz5kZWWxdevWfPWaNWvG3/72N3788Ufee+89Dj/8cM444ww+/vjjIvffuXNn/vOf/7B69WpmzpzJxo0b6d+/f4HtS+lV6DV37v4C8EJc0fNmVpdwLd1dsWpJVrX4Ze7+vJm1AtoBX7p7rplNBe5w9yVmdgFwRbTOOHe/t7CYzGw4ofWP1q1bM3v27FIeXWqys7PTvo9qL9nt/F26FLmKzkvVVKnnZWWS91GLot9H1UV1+rw0bdqUdevWFShvXAmxAEljAcjNzS10GcCmTZtw97w6Z511FuvXr+fMM8/MK3N3Nm3axIoVK8jNzc1Xf926dTzzzDOYWdL97Lbbblx33XU8/vjjvPPOO+yyyy4A5OTkFBnXoYceyh/+8AfOPvtsvvvuO1q0aJHaC7GdKO68FGfjxo2l+qxVhXHupgCDgA6EFrrGZlYzofWuGbDB3TfHCtx9A/AZgJn1BQ4EBpvZgYTr7rpFVeeb2Wvu/n6ynbv7eGA8QFZWlvfq1av8jiyJ2bNnk+59VHu9excsK+YSTJ2XqqlSz8vEJO+jXrqUF6rX5+WTTz6hcePKSuUKKiyWdevWFRln3bp1MbO8Ov369aNfv3756pgZdevWpV27dhx66KHccssttGrViho1ajBu3DiaNWvG2rVr87bRvXt3BgwYwH777YeZ8cADD9CwYUOOOuooGjduzAEHHMDy5cv5z3/+w3777UfLli1Zu3YtI0aMYPDgwey+++6sWrWKv/3tbxx44IHsuuuu5fQqVR3FnZfi1KtXj4MPPrjE61WF5C7GCUOf1AQ6AovilnWOlhVgZrUIQ5+McvccM+sFvOTun0bLZwFHAUmTOxEREclv4sSJDB8+nDPPPJMddtiBCy+8kA0bNnD33Xfn1TniiCOYMGECixcvpmbNmhx88ME8//zztGvXDoBBgwbx8ssvM2rUKH766SfOOussbrnlFlq3bs1NN93EDz/8QLNmzejduzc333xzZR1qRqoKyd1vgBXAN4Rr8tYCpwI3AphZA8L1eeMLWf88YJW7PxFX1iDu94Zs69YVEREpuzLMFFEZxowZU+y8sStWrMj7vWPHjrz00ktJtxNz6623cuuttxa6vXr16vHwww8XKP/XvzRSWbpVaHJnZv8B3iS0otUEBkePi919K7DRzMYB15jZKrYNYlwD+HuS7TUHrgOOiyt+BbjFzIZFz/ugmyxERESkmqjolrtFwDBgF0Jr2sfAme4en8aPIyRzVwI7EIY5OcbdlyXZ3lhgqru/HStw93fMbBRwU1Q0wt3fK/cjEREREamCKvpu2bxpwoqo44TE7Kai6kV1Ly6k/C7grmTLRERERDJZVZuhQkRERETKQMmdiIiISAZRciciIiKSQZTciYiIiGQQJXciIiIiGUTJXQVbuHAhZpbvISIiIlJelNyJiIiIZBAldyIiIhkqsaco2WP27NmVHaaUs6owt6yIiMj2pbIuqfGSzWk7f/78vN9zcnLo06cPV199Nf37988r32effcotPKkalNyJiIhkqK5du+b9np2dDcAee+yRr7yqyMnJoX79+pUdRkZQt6yIiEg1lZuby5gxY2jfvj1169Zl3333ZeLEiXnLX3rpJcyMH374Ia/siCOOoGbNmqxevTqvbP/992f06NF5z7/99luGDBlCixYtaNCgAccddxyLFi3KW7548WLMjMcff5wzzzyTZs2acdJJJwEwdepUunTpQsOGDWnevDmHH344c+bMSeOrkHlSTu7MrK6Z/cHMHjSzGWbWKSofbGZ7py9EERERSYdrr72Wm266ieHDhzN16lSOPPJIzjjjDP79738DIZGrXbs2r776KgAbNmxg4cKF1KlTh7lz5wKwcuVKPvroI3r06JH3vHv37ixatIh//OMfTJ48mfXr19O3b19ycnLy7X/EiBE0btyYJ598kquuuoovv/ySgQMH0qdPH6ZNm8bjjz/OiSeeyMqVKyvwVdn+pdQta2Z7Ai8CTYGFQC+gcbS4B9AfODMN8YmIiEgarFy5kjvvvJOrr76aq6++GoDjjjuOJUuWMGbMGE477TTq169Ply5dePXVVxk8eDCvv/46TZo0oW/fvrz66qv079+f1157DTOjW7duANxxxx2sX7+ed999lxYtWgBw5JFH0qFDBx566CEuuOCCvBi6du3KPffck/d8ypQpNG7cmFtvvTWvrF+/fhXxcmSUVFvu/gZ8C3QAjgPirySdA3Qv37BEREQknT788EM2bNjAqaeemq988ODBfPbZZyxfvhyAHj165LXcvfLKK/To0YOjjjoqX9mBBx5IkyZNAJg5cybHHHMMTZo0YcuWLWzZsoXGjRvTpUsXFixYkG9f8Td2QOjeXbNmDWeddRYzZsxg/fr1aTn2TJdqctcD+Iu7rwYSb9VZBuxUnkGJiIhIei1duhSA1q1b5yuPPV+1ahUAPXv25MMPP2T16tW8+uqr9OjRgx49erBgwQI2btyYVxazYsUKnnjiCWrXrp3v8fLLL/Pdd98l3VfMXnvtxTPPPMNXX31Fv379aNmyJaeffjo//fRTuR9/Jks1udsIFHYLS1tgdblEIyIiIhVip51Cu0yshS5m2bJlAPm6VN2d2bNn8/rrr9OzZ0/23XdfGjVqxKxZs3j77bfzJXctWrTg5JNP5q233irwiO+CBZLO0tS/f39effVVfv75Zx588EFmzpzJRRddVK7HnulSHQrlReAqM5sJZEdlbmZ1gYuA59IRnIiIiKTHfvvtR4MGDXjyySe59tpr88onT57MnnvuyY477ghA8+bN2W+//bjjjjuoWbMmBx98MGZG9+7dueWWW9iyZQvdu2+7Ouvoo49m8uTJ7LvvvmUa2qRp06acfvrpzJkzJ994fVK8VJO7kcBc4AtCoufAtcC+QB3g12mJTkRERNKiRYsWXHrppdx4443UqlWLrKwsnnrqKZ577rm8u2VjevbsyT333MNxxx1HzZo1gXAt3siRI+nUqRNt2rTJq3v55Zfz2GOP0adPHy666CLatm3LsmXLmDNnDt27d+e0004rNKb777+f+fPnc/zxx7Pzzjvz+eef8+STT3LmmbpnsyRSSu7c/TszOxC4HDga+JJwnd2TwF/d/ef0hSgiIiLpcP3111OrVi3uu+8+li1bRseOHXnssccYMmRIvno9evTgnnvuoWfPnvnKgHytdgAtW7bk9ddfZ/To0Vx22WWsXr2anXbaie7du3PAAQcUGc8BBxzA1KlTufzyy1m5ciU77bQT5557Ltdff305HXH1kPIMFe6+CrgmeoiIiFRfJZwGrCpo1KgRnhB3zZo1GTt2LGPHji1y3cGDBzN48OB8ZYcffniB7cXsvPPOPPzww4Vur0OHDknXPeKII5g+fXqRsUjxUrqhwswONLOkA82YWT8zKzoVFxEREZEKkerdsncAhxey7NBouYiIiIhUslSTu0MIN1QkMx84uHzCEREREZGySDW5qwk0LGRZQ8IdsyIiIiJSyVJN7t4ChheybDiwoJBlIiIiIlKBUr1bdgww08zeAB4BfiQMhXImcCBwTFqiExEREZESSXWcu1fM7FjgL8DfAQO2Am8Ax7j7q+kLUURERERSVZJx7mYDR5hZA6A5sMrdN6QrMBEREREpuZSTu5gooVNSJyIiIlIFpXpDBWaWZWZ/NrNHzWxywuOJdAYpIiIiZffUU0/Rp08fmjVrRt26ddlzzz25+uqrWbFiBYsXL8bMePbZZ4vcxpgxY2jZsmUFRSylkeoMFX8gXF93DrAHsGPCo1W6AhQREZGy++Mf/8ipp57K7rvvzr/+9S9mzJjBZZddxrRp0zj33HNT3s4555zDCy+8kMZIpaxS7ZYdATwMnOfuW9IYj4iISJX3wfm7Vsp+97/3m1KtN23aNP7617/y4IMPMmzYsLzyo446iuHDhzNjxoyUt9WuXTvatWtXqjikYqTaLdsK+LcSOxERke3PHXfcwSGHHJIvsYupWbMmJ5xwQt7zDRs28Pvf/56mTZvSrl07rrvuOrZu3Zq3PLFbdvbs2ZgZs2bN4pRTTqFhw4Z06tSJGTNmkJuby8iRI2nZsiVt27blr3/9a759Dx06lKysLF588UUOOOAAGjZsSPfu3fnoo4/y1du6dSvjxo2jY8eOed3JjzzySN7ye+65h8aNG5OdnZ1vvZdffhkz4/333wdg6tSpdOnShYYNG9K8eXMOP/xw5syZU4pXtGpLNbl7nsLnli01M2trZtlm5mbWKK7czOwqM/vOzHLM7BUzOyhh3c5m9oaZrTGzSfHrR8t7mtn3ieUiIiLVyebNm5k3bx7HH398SvVHjRpFo0aNmDJlCr/97W+5/vrrmTJlSrHr/f73v6d79+48/fTT7LrrrgwcOJALL7yQdevWMXHiRAYOHMgf//hHXn/99Xzrffvtt4wcOZLRo0fz73//m+XLlzNo0CDcPa/ORRddxI033sjw4cOZPn06AwYMYNiwYXnXB55xxhls2bKlQJwTJkzgkEMO4YADDuDLL79k4MCB9OnTh2nTpvH4449z4oknsnLlypRel+1Jqt2y9wDjzaw28CKwOrGCu39civ3fCmRTcGqzPwHXACOBT4HLCYMo7+fuP0Z1JgBfANdG27kqemBmNYA7gSvdPX8aLyIiUo38/PPPbNq0ifbt26dUv2fPntx+++0AHHPMMfzvf//jqaeeYtCgQUWu97vf/Y6RI0cCoet23333ZdGiRbz00ksA9O3blyeeeIKnn36arl275q23cuVK5s6dS6dOnYDQSjdgwAAWLVpE586d+eKLL7jvvvt4+OGHOeuss/K2tXTpUsaOHcuJJ55Is2bN+M1vfsPDDz/M0KFDAcjOzuY///kP48aNA+Cdd96hcePG3HrrrXn77tevX0qvyfYm1Za7l4FOwHXAq8AHcY8Po58lYmY9gOOB2xLK6xGSu7+4+93uPhM4FXDgwqhOI0JL4qXu/gJwE/lnyTgb2Az8q6RxiYiIZCIzS6nesccem+/5Pvvsw5IlS4pd7+ijj877vWPHjgD06dMnr6xGjRrsvvvufP/99/nW69ChQ15iF9sfkLfPWbNmUaNGDQYMGMCWLVvyHkcffTTvvvsuubm5AJx99tm8+uqrfPXVVwBMnjyZLVu2cPrppwOw//77s2bNGs466yxmzJjB+vXrU3o9tkepJne94x59Eh6xspSZWU3CTBfXAysSFncDmgCTYwXuvh6YBsQuCqgT/cyJfm6IlZlZE+AG4BKPb9MVERGphnbYYQfq1q3Lt99+m1L9Zs2a5Xtep04dNm7cWKL16tSpk/K2ktUB8uqtWLGC3NxcmjZtSu3atfMeQ4cOZcuWLSxduhSAXr16sfvuuzNhwgQAHn74YU455RRatGgBwF577cUzzzzDV199Rb9+/WjZsiWnn346P/30U7HHtr1Jdfqx8r7a8DygHqG794yEZZ2BXODzhPJPgMFRPCvN7GvgIjO7HxgOLIjqXQPMdPfXERERqeZq167NkUceyQsvvMCNN95Y2eGUWIsWLahVqxZz586lRo2CbVKtWoXR2MyMYcOGMX78eH73u9/x2muv8fzzz+er279/f/r378+aNWuYPn06l156KRdddBGTJk2qkGOpKCkPYhxjZjXMrEHiowTr70BoWbvc3TcnqdIcyHb33ITyVUADM4u12l0AjAZ+BvYCxppZR0KX7J9KeFgiIiIZ69JLL2XBggX57jCN2bp1K//73/8qIarU9OnTh9zcXNasWUNWVlaBR6ylD8Ldt0uWLGHYsGG0bduWY445Juk2mzZtyumnn86AAQP4+OPS3DJQtaXUcmeho34UcC6wWyHVaqa4z5uAN9z9uSLqJOtOtfhl7v68mbUC2gFfunuumU0F7nD3JWZ2AXBFtM44d7832Y7MbDih5Y/WrVsze/bsFA+jdNq1a8dtt+W7zDDt+6x2El5fAIp5jbOzs3UeqqBKPS/1Sv4+qi6q0+eladOmrFu3rrLDyFNYLLm5uUXG2atXLy688ELOPvtsZs+eTb9+/WjUqBGfffYZDz30EO3bt+cvf/kLEIZCid/W5s2b2bp1a17Zpk2bcPe85xs2hBlJ169fXyCGjRs35ivLzc1ly5YteWWJ2wbyhjOJxbHzzjszbNgwBg8ezKWXXsrBBx/Mxo0b+fTTT/niiy+4++6789Zt3Lgxffv25YUXXuCPf/xjXmwADz30EG+++SZ9+/Zlp5124ssvv2Ty5MmcdtppaTvHxZ2X4mzcuLFUn7VU75a9mNAadgshObuR0HU6hHCt259T2YiZ7QsMA3qaWbOoONbq19TMcgktdI3NrGZC610zYEN8a180z+1n0bb7AgcCg83sQELrYLeo6nwze83d30+Myd3HA+MBsrKyvFevXqkcSqndfvvtjBgxIjGGtO6z2undu2BZMa/x7NmzSfe5l5Kr1PMyMcn7qJc+q1C9Pi+ffPIJjRs3ruww8hQWy7p164qN8+9//zu9evXi7rvv5pxzziEnJ4cOHTpw8sknM2LEiLxr3Bo0aJBvW7Vr16ZGjRp5ZXXr1sXM8p43aBD+jDds2LBADPXq1ctXVrNmTWrVqpVXlrhtgEaNGhWI44EHHmC//fbjgQce4KabbqJJkybss88+nH322QX2OXDgQF544QV+//vf51t22GGHMWPGDEaPHs3KlSvZaaedGD58ONdffz1169Yt8rUrrVTOS1Hq1avHwQcfXOL1LJXEwsw+JCRA9xDuQs1y97ejIUemAR+4e7FdoWb2K+DpIqo8CEwEZgGd3X1R3LoPAge5e5ck260FvAvc4O5PmNklQA93HxgtnwLMcfe/FxVfVlaWL1iwoKgqZabkrgIkuyNMyd12qXKTuyTvo9P1WYXq9Xn55JNP2HvvvSs7jGKVNYnIJIMGDWLp0qW8+uqrlR1Kmc9Lce8/M1vo7lmJ5am23O0GvBt1fW4mtKLh7lvN7F7gn6R2ndtrhLtr4x1P6D7tB3wFfAOsJQx/cmMUfAPgJKIWtiTOA1a5+xNxZfHXATZkW7euiIiIZJgPPviABQsW8NRTT2XcDRIllWpy9zMQm+nhW+Bg4KXoeXOgfiobcfcVwOz4MjPrEP36amzAYTMbB1xjZqvYNohxDcLwKSSs35ww/t5xccWvALeYWWyelT7oJgsREZGMddJJJ7FixQrOP/98Bg4cWNnhVKpUk7u5wKHAc4Ru0zFm1gL4hXDX6qxyjmscIZm7EtiBMMzJMe6+LEndscBUd387VuDu75jZKML1gQAj3P29co5RREREqojFixdXdghVRqrJ3RigbfT7nwndskMJLXYvAheVNgB3n0CYSiy+zAmJ2U1JVklc/+JCyu8C7iptXCIiIiLbo1QHMV4ELIp+3wRcEj1EREREpAop8SDGIiIi1YlGNJDKUJb3XaEtd2b2UmHLCgmiRPPLioiIVHW1a9cmJycnbyw3kYqSk5ND7dq1S7VuUS13Pyc89gR6EIYYyY5+dgc6AStKtXcREZEqrFWrVnz//fds2LBBLXhSIdydDRs28P333+fNm1tShbbcufupsd/N7GzC/K3d3P3buPL2wLOEmypEREQySpMmTQD44Ycf2Lw52XToVcPGjRupV69eZYchCUp7XmrXrk3r1q3z3n8llerdsqOBy+MTOwB3/9bMrgP+CjxQqghERESqsCZNmpT6j2xFmT17dqmmqZL0qqzzkuoNFW2AwiZeqwuUrt1QRERERMpVqsndbOBmM8s3f5mZHQrcDMwp57hEREREpBRSTe6GAyuBN8zsBzN718x+AF6PyoenK0ARERERSV2qgxgvAQ4xs36EacjaAD8Cb7n7c2mMT0RERERKINUbKgCIEjklcyIiIiJVVImSOzOrS5hjtsB9ve7+cXkFJSIiIiKlk1JyZ2Y7A+OBE5ItBhyoWY5xiYiIiEgppNpy90/gEOBy4GPgl7RFJCIiIiKllmpydyRwrrtPTmcwIiIiIlI2qQ6FshzISWcgIiIiIlJ2qSZ31wJXmFnVnn9FREREpJpLtVv210B74BszewtYnbDc3X1weQYmIiIiIiWXanLXEvgy+r02sGN6whERERGRskh1hore6Q5EREREKt8H5+9aoGz/e7+phEiktFK95k5EREREtgOpDmJ8S3F13H1U2cMRERERkbJI9Zq7U5OUNQeaAGuAVYCSOxEREZFKluo1d7slKzezwwnTkp1XnkGJiIiISOmU6Zo7d38DuBW4u3zCEREREZGyKI8bKn4G9iqH7YiIiIhIGaV6Q0WDJMV1gL2B64GPyjMoERERESmdVG+oyAY8SbkB3wO/Kq+ARERERKT0Uk3uhlEwudsILAHedPfN5RqViIiIiJRKqnfLTkhzHCIiIiJSDlJtuQPAzPYCDgV2ApYCC9z903QEJiIiIiIll+oNFU2AB4DfEO6wzQYaAVvN7CngHHdfm7YopUTMrECZe7JLJkVERCTTpDoUyr3AscCZQAN3bwI0AM4CjomWF8vMBprZPDP72cw2mtkiM7vazOrE1TEzu8rMvjOzHDN7xcwOSthOZzN7w8zWmNkkM2uUsLynmX2fWC4iIiKS6VJN7k4BRrr7RHffCODuG939ccK0Y6ekuJ0dgJeBc4ATgIeA0cBf4+r8CbgGuBk4idBKONPM2sTVmQB8AQwC9gGuii0wsxrAncCV7p6dYlwiIiIiGaEkQ6EsLWTZD8D6VDbi7vcnFL0cdfleYGYXAXUJyd1f3P1uADObDywGLgSujlrjDgdOcvefzKwZMIJtCd7ZwGbgX6kdmoiIiEjmSLXl7h5ghJnVjy+MBjceQYrdsoX4mTAgMkA3oAkwObbQ3dcD0wgtfcTVzYl+boiVRYniDcAlrovMREREpBoqtOXOzG5JKOoEfGdmLwLLgVaE6+1ygAUl2amZ1SS00h0CXAzc5+5uZp2BXODzhFU+AQYDuPtKM/sauMjM7geGx+3/GmCmu79eknhEREREMkVR3bKnJjzfHD26xpWti37+BhhZgv2uJyR3AI/GrdscyHb33IT6q4AGZlbH3X8BLgCeBP5MSAQvMLOOhC7ZA0oQh4iIiEhGscrovTSzQwh32x4GXAtMdPfzzWw0MMLdmyfUPxcYD9SJzYYRdQm3A75091wzmwq85e43mNkFwBXR6uPcvdBuYzMbTmj9o3Xr1l0mTZpUrseaaNmyZSxZsiRfWZcuXcp1HwsXLixQVt77qNKSHD/FHH92djaNGunm6qqmUs/LyiTvoxbV6HNUBH1eqp7yPCc5335QoKx++/3LZdvVTbo/K717917o7lmJ5ZWS3OULwOxM4BGgI3Ac8DegbnzrnZmNBMa4e8NCttEXeBDoDOxJuCO3W7R4PnCUu79fXCxZWVm+YEGJephL7Pbbb2fEiBH5ysr7HFT7ce6SHD/FHP/s2bPp1atXeuKRUqvU8zIxyfvo9Gr0OSqCPi9VT3mekw/O37VA2f73flMu265u0v1ZMbOkyV3KM1RYyBiOJCRP9RKXF9U6Voy3o5+7AZ8CNQmJ3qK4Op2jZcniqkUY+mSUu+eYWS/gpdjMGWY2CzgKKDa5ExEREdnepTpDRWtgFmFMOQdi/9LG/xtb2uTuyOjn18D3wFrC9X43RvtuQBjvbnwh658HrHL3J+LKGsT93jAuXhEREZGMlmrL3e3AGmAX4DvCOHPLgN8SZq3on8pGzOx/wEzgI8JdsUcCfwSecPcvozrjgGvMbBWhte5ywpAtf0+yvebAdYTu3JhXgFvMbFj0vA9h7DwRERGRjJdqcncUcAnbBjI2d/8W+HM0I8S95E+wCvMWMBToAGwBvgKuBP4RV2ccIZm7kjCjxQLgGHdflmR7Y4Gp7h7r2sXd3zGzUcBNUdEId38vhdhEREREtnupJnfNgJ/cfauZrSWMcRczj213phbJ3a8hjEVXVB0nJGY3FVUvqntxIeV3AXelEpOIiIhIJkl1hoqvgZ2i3z8CzohbdhKwsjyDEhEREZHSSbXlbjpwLGFasBuBZ8xsCWFQ4/ak2HInIiIiIumVUnLn7lfG/f68mXUDBgD1gRfd/fk0xSciIiIiJZDyOHfx3H0BJZxPVkRERETSL9Vr7kRERERkO6DkTkRERCSDKLkTERERySBK7kREREQyiJI7ERERkQyScnJnZgeY2RNm9qWZbTKzQ6Lym8zshPSFKCIiIiKpSmkolCh5m0qYauxR4Lq4xZuAiwCNdSciFWOiFSw73Ss+DpHIB+fvWqBs/3u/qYRIRFJvufsLMMHdj6LgnK/vAgeVY0wiIiIiUkqpJnedgSei3xP/PV4LtCi3iERERESk1FJN7pYDuxeybF/g2/IJR0RERETKItXkbhJwvZl1jytzM9sTuAJ4vNwjExEREZESS3Vu2WuAfYA5wI9R2TNAG2AG8OfyD01ERERESiql5M7dNwEnmtnRwNFAS2AlMMvdX0xjfCIiIiJSAqm23AHg7rOAWWmKRURERETKKKVr7sxssZndbGYHpzsgERERESm9VG+omAIMBhaY2Wdmdr2Z7ZfGuERERESkFFJK7tx9hLt3ALoDzwHDgPfM7EMzu8bMOqUxRhERERFJUcpzywK4+3x3vxTYBegFzCZMPfZJeQcmIiIiIiVXouQuTkOgPbAr0JQwv6yIiIiIVLKUkzszq29mg8zsP4QZKx4kTEU2DGidpvhEREREpARSGgrFzJ4A+gN1gZeAC4Cn3X11+kITERERkZJKdZy71sAIYIq7r0hjPCIiIiJSBqnOUNErzXGIiIiISDkoNLkzs32AL919U/R7kdz943KNTERERERKrKiWuw+BrsCb0e9eSD2LltUs39BEREREpKSKSu56A7HWuD4UntyJiIiISBVRaHLn7nPifp9dIdGIiIiISJmkNM6dmeWa2WGFLOtiZrnlG5aIiIiIlEaqgxhbEctqA1tS2ojZqWY21cy+N7NsM1toZqcl1DEzu8rMvjOzHDN7xcwOSqjT2czeMLM1ZjbJzBolLO8Z7SNfuYiIiEimK+pu2fZAh7iig82sXkK1esBZwNcp7u/yqO5lwAqgHzDRzFq6+9+jOn8CrgFGAp9G68w0s/3c/ceozgTgC+Ba4FbgquiBmdUA7gSudPfsFOMSERERyQhF3VDxf8B1hBspHLivkHo5wDkp7u+khEGQXzKznQkJ3N+j5PFPwF/c/W4AM5sPLAYuBK6OWuMOj7b1k5k1IwywfFW0zbOBzcC/UoxJREREJGMUldzdC0whdMm+D5wR/Yz3C/Ctu29KZWeFzG7xDnBK9Hs3oAkwOW6d9WY2DTgBuBqoEy3KiX5uiJWZWRPgBuBX7q67e0VERKTaKepu2Z+AnwDMbDdgqbv/koYYurFtyJXOQC7weUKdT4DBUVwrzexr4CIzux8YDiyI6l0DzHT319MQp4iIiEiVl+r0Y98AmFktoD3hWrvEOiWeocLMjia02g2LipoD2e6eePftKqCBmdWJEswLgCeBPxMSwQvMrCOhS/aAksYhIiIikiksld5LM6sN/I1w80TdZHXcvUQzVJhZB+ANYJ67D4jKRgMj3L15Qt1zgfFAHXffHJU1ANoRpkjLNbOpwFvufoOZXQBcEa0+zt3vLSKO4YTWP1q3bt1l0qRJJTmMElu2bBlLlizJV9alS5dy3cfChQsLlJX3Pqq0JMdPMcefnZ1No0a6ubqqKfS8rExyjluU83u8IvaxndLnpaCcbz8oUFa//f4Vtv/yPCeVfSyZJN2fld69ey9096zE8lSTuxuAocAo4HFCy9l64LfAHsBF7v5cqsGYWQtgLpAN9HL39VH5+YQksm58652ZjQTGuHvDQrbXF3iQ0K27J/AyobsXYD5wlLsnXi9YQFZWli9YsKC4amVy++23M2LEiHxl5X15oFnBkWuq1SWISY6fYo5/9uzZ9OrVKz3xSKkVel4mJjnHp5fze7wi9rGd0ueloA/O37VA2f73flNh+y/Pc1LZx5JJ0v1ZMbOkyV2q49wNAsaw7UaHN939UXc/FniNbTdEpBJIA+BZwk0Q/WOJXeRTwhy1HRNW6xwtS7a9WoShT0a5ew7QC3jJ3T9190+BWcBRqcYnIiIisj1LNbnbBfgsak3bSLg2LuZx4DepbCRKxJ4EOgEnuPvyhCrzgLXAqXHrNABOAp4vZLPnAavc/Ym4sgZxvzek6EGYRURERDJGSjdUAEuBZtHvXwM9gZnR8z1KsL97CQMXXwK0MLOuccvecfeNZjYOuMbMVrFtEOMawN8TN2ZmzQlj8R0XV/wKcIuZxW7S6EMYO09EREQk46Wa3M0GegDTgAeA26K7UzcRhij5d4rbOTb6eVeSZbsRBiseR0jmrgR2IAxzcoy7L0uyzlhgqru/HStw93fMbBRwU1Q0wt3fSzE+ERERke1aqsndaKAlgLvfaeGK/YFAfUKL2vWpbMTdO6RQxwmJ2U0p1L24kPK7SJ5AioiIiGS0VMe5+xH4Me75HcAd6QpKREREREon1RsqRERERGQ7UGjLnZm9BaQ8qJO7H1YuEYmIiIhIqRXVLfsRJUjuRERERKTyFZrcufvQCoxDRERERMpBia65s2AXM+tmZkmnAhMRERGRypNychfN+/o98A3wKrBXVP6UmV2aluhEREREpERSSu7MbCTwV8IAxn3IP53XbMJAxiIiIiJSyVIdxPgC4Fp3v8XMaiYsWwTsWb5hiYiIiEhppNot2wZYWMiyrUC98glHRERERMoi1eTuC+CoQpb1BD4un3BEREREpCxS7Za9E7jXzH4BpkRlrczsbOBy4Nw0xCYiIiIiJZTq3LL/NLPmwLXA2Kj4OWADMMbdJ6YpPpGqyaxgmWvMbxERqXypttzh7rea2T+AI4CWwEpgvruvSVdwIiIiIlIyKSd3AO6+DpgRX2ZmvYFR7n5CeQYmIiIiIiVXZHJnZs2A44FdgK+BZ9x9c7TsVOAK4BDgs/SGKSIiIiKpKDS5M7P9Ca10reOK3zaz3wATCd2zHwFnAE+kM0gRERERSU1RQ6H8GVhLSOIaAHsTrrN7C9gPONPd93f3f7v71rRHKiIiIiLFKqpbNgu4xN3fiJ4vMrM/AJ8Dw939sbRHJyIiIiIlUlTLXWtgcUJZ7Pl76QhGRERERMqmuBkqChu4a0t5ByIiIiIiZVfcUCgvmFmyRG5WYrm7tyq/sERERESkNIpK7sYWsUxEREREqqBCkzt3V3InIiIisp0p0QwVImWi+VhFRErtg/N3LVC2/73fVEIkUtUVd0OFiIiIiGxHlNyJiIiIZBAldyIiIiIZRMmdiIiISAZRciciIiKSQZTciYiIiGQQJXciIiIiGUTJnYiIiEgGqfDkzsw6mtn9ZvaemeWa2ewkdczMrjKz78wsx8xeMbODEup0NrM3zGyNmU0ys0YJy3ua2feJ5SIiIiKZrDJa7vYF+gGfRY9k/gRcA9wMnARkAzPNrE1cnQnAF8AgYB/gqtgCM6sB3Alc6e7Z5Ru+iIiISNVVGcndNHffxd1PBT5KXGhm9QjJ3V/c/W53nwmcCjhwYVSnEXA4cKm7vwDcBBwTt5mzgc3Av9J6JCIiIiJVTIUnd+6+tZgq3YAmwOS4ddYD04AToqI60c+c6OeGWJmZNQFuAC5x18SlKTMr+BAREZHtTlW8oaIzkAt8nlD+SbQMd18JfA1cZGYtgOHAgqjeNcBMd3+9YsIVERERqTpqVXYASTQHst09N6F8FdDAzOq4+y/ABcCTwJ8JieAFZtaR0CV7QEUGLCIiIlJVWGX2XJrZFKClu/eKKxsNjHD35gl1zwXGA3XcfXNU1gBoB3zp7rlmNhV4y91vMLMLgCui1ce5+72FxDCc0PJH69atu0yaNKlcjzHRsmXLWLJkSb6yLl26lOs+Fi5cWKCs2H0kWYdyjqtC9lHK/WRnZ9OoUQlurK6oY6nmCj0vK5O8/i3K+fWviH1sp0r8eakGcr79oEBZ/fb7V9g+yvOcVMSxVBfp/qz07t17obtnJZZXxeTufOBvQN341jszGwmMcfeGhWyrL/Agoet2T+BlwvV7APOBo9z9/aLiycrK8gULFhRVpcxuv/12RowYka+svM+BJblerth9JLvGrrzfGxWxj1LuZ/bs2fTq1Sut+5CSK/S8TEzy+p9ezq9/RexjO1Xiz0s18MH5uxYo2//ebypsH+V5TiriWKqLdH9WzCxpclcVr7n7FKgJdEwo7xwtK8DMahGGPhnl7jlAL+Ald//U3T8FZgFHpStgERERkaqiKiZ384C1hOFPgLzu15OA5wtZ5zxglbs/EVfWIO73hoBu/xQREZGMV+E3VESJWr/oaVugiZkNjJ4/5+4bzGwccI2ZrSK01l1OSET/nmR7zYHrgOPiil8BbjGzYdHzPoSx80REREQyWmXcLduKcJdrvNjz3YDFwDhCMnclsANhmJNj3H1Zku2NBaa6+9uxAnd/x8xGEQY3hnCDxnvldgQiIiIiVVSFJ3fuvphiukijwYdvYltyVlTdiwspvwu4qxQhioiIiGy3quI1dyIiIiJSSkruRERERDJIVZyhQqTCLVy4kN69e+cry7ipiQsbt01j9omIlFpVHBdQLXciIiIiGUTJnYiIiEgGUXInIiIikkGU3ImIiIhkECV3IiIiIhlEyZ2IiIhIBlFyJyIiIpJBlNyJiIiIZBAldyIiIiIZRMmdiIiISAZRciciIiKSQTS3rJS/wuYwFRERkbRTy52IiIhIBlFyJyIiIpJBlNyJiIiIZBAldyIiIiIZRMmdiIiISAZRciciIiKSQZTciYiIiGQQJXciIiIiGUTJnYiIiEgGUXInIiIikkGU3ImIiIhkEM0tW92U47yvZgW35V6+c8gWuo9kxwFFHssH5+9aoGz/e78pOoCKmCe3sH0kOXYA3Et3LCIpyqT3V1U4lor4rly4cCG9e/dO6z7KU1U4L5lMLXciIiIiGUTJnYiIiEgGUXInIiIikkGU3ImIiIhkECV3IiIiIhmkyiZ3ZraPmc0ysw1m9oOZXW9mNRPqXG9mP5nZV2Z2UpJtzDKzyyouahEREZHKVSWHQjGz5sBM4GPgFGAP4HZCMnp1VOd44GJgeLT8cTPbzd1/jpYPANoBd1f4AYiIiIhUkiqZ3AHnAfWBX7v7WuBFM2sCjDGzW6KyvsDj7j4ZwMzOBLoC082sLnAbcLG7b66cQxARERGpeFW1W/YE4IUoiYuZREj4joqe1wFy4pZviMoALgM+d/fp6Q5UREREpCqpqsldZ+DT+AJ3/5aQwHWOihYCvzaz3czsaGA/4F0zawOMIiR4IiIiItVKVe2WbQ6sTlK+KloGMBE4DfgKcOAad//azB4G/uXun1REoCIiIiJViVXFuefMbDMwwt3vSij/Hpjg7qPjyjoA2e6+wsy6AM8DewENgQeBQ4EFwFB3/yHJvoYTbsogWm9R+R9RPi2BFWneh5SczkvVpPNSNem8VD06J1VTus/Lru6+Y2JhVW25WwU0S1LelIQWPXdfHPf0LmCMu68yswcJidoA4Bbgb8DAxA26+3hgfHkEnQozW+DuWRW1P0mNzkvVpPNSNem8VD06J1VTZZ2Xqprcfcq2a+sAMLNdCK1xnyZbwcxOIyR/90dFvYEe7r7BzP4BzElfuCIiIiJVQ1W9oeJ54DgzaxxXNphwd2yBJM3M6gM3A5e5e27cogbRz4aApSlWERERkSqjqiZ3/wA2AU+ZWd/ourgxwF8ThkeJGQW87e4z48rmADeY2THAjcDs9IacsgrrApYS0XmpmnReqiadl6pH56RqqpTzUiVvqIAw/RhhdokjCNfZ/ZNwPV1uQr12wAfAoe7+RVx5e+ARIAt4CzjT3ZdUTPQiIiIilaPKJnciIiIiUnJVtVs2o5jZPmY2y8w2mNkPZna9mdWs7LiqMzM71cymmtn3ZpZtZgujm3KkCjGzttH5cTNrVNnxVGdmVsvM/mRmn5vZJjNbYmZ3VHZc1Z2ZDTGzt6PPyfdm9qiZ7VzZcVUXZtbRzO43s/fMLNfMZiepY2Z2lZl9Z2Y5ZvaKmR2UzriU3KWZmTUHZhIGWj4FuB74IzC2MuMSLgeyCTOZnAy8DEw0s4sqNSpJdCvhPEnlexi4mDBv97HAn8g/BaRUMDM7Gfg3MI/w9+UKoCfwrJnp73vF2BfoB3wWPZL5E3AN4cbPkwjfaTOjGbXSQt2yaWZmVxJu+Ng1djOImY0i3CDSppAbRCTNzKylu69IKJsIHOHuu1VSWBLHzHoAzwB/JiR5jd1diV4lMLPjgWnAge7+cWXHI4GZTQI6uXuXuLKTCZ+bfTRTU/qZWQ133xr9PgVo6e694pbXA5YBt7v79VFZQ2AxcL+7X52OuJTZp98JwAsJSdwkoD5wVOWEJImJXeQdoFVFxyIFRZct/J3Q0q1R9yvfMOAlJXZVTm1gTULZ6uinhv+qALHErgjdgCbA5Lh11hP+WTohXXEpuUu/ziQMvOzu3wIbSBioWSpdN0B/vKqG84B6wD2VHYgAcDjwmZndbWZro+uHn9K1XZXuIaCHmZ1pZk3MbE/C0F8vKxGvMjoDucDnCeWfkMYcQMld+jUnYcq0yKpomVQBZnY04ZoVJROVzMx2AG4ALnf3zZUdjwDQBhgKHAQMAf4P6AI8bWZqIaok7j6dcF7GE1rwFgE1gV9XYliSX3MgO3EYN0IO0MDM6qRjp1V1+rFMk+zCRiukXCqYmXUAJgLPuPuEyo1GgJuAN9z9ucoORPJY9DjF3X8GMLOlhMHi+wCzKjG2asvMehMG/b+LMLNTa8L13E+bWd8kCYVUjsJygMKWlZmSu/RbBTRLUt6U5C16UoHMrAXhS/Fb4LeVHE61Z2b7Eq7v6mlmzaLi2DSCTc0s1911h2bFWwV8FUvsIq8BvwD7oOSustwOTHX3K2IFZvYu4VKgU4CnKiku2WYV0NjMaiYk282ADenqnVC3bPp9SkK/upntQpjv9tOka0iFMLMGwLNAHaB/dJGrVK5OhIvE5xO+FFexrat8CeEmC6l4hd11aUBxF5RL+nQG3o0vcPdFhCFq9qiMgKSATwld5R0Tygtcj1+elNyl3/PAcWbWOK5sMOHDN6dyQhIzqwU8SUgmTnD35ZUckgSvAb0THjdHy/oRhkSRivcscICZtYwr60lIxN+rnJAE+AY4JL7AzPYmjMawuDICkgLmAWuBU2MFUcPCSYT8IC3ULZt+/yAM/PmUmd0M7E64JuKvGuOuUt1LSBYuAVqYWde4Ze+4+6bKCat6i4aomR1fFl0TCfCqxrmrNOMJ32PTzOzPQGNC0j3T3V+r1Miqt38Ad5jZD2y75u5aQmKna1YrQJSo9YuetgWamNnA6Plz7r7BzMYB15jZKkJr3eWExrW09URoEOMKYGb7AHcDRxCus/snMEYXu1YeM1sM7FrI4t3cfXHFRSNFMbOhhNkRNIhxJTKzjsDfCONz/kIYKPcyd19VqYFVY9GdyucBfyB0w64mtH5f6e5fVWJo1Ub0z+fXhSzezd0XR+fpKsJ52gFYAFzs7u+kLS4ldyIiIiKZQ9fciYiIiGQQJXciIiIiGUTJnYiIiEgGUXInIiIikkGU3ImIiIhkECV3IiIiIhlEyZ3kY2ZjzMyjx1YzW2Vmb5nZTWbWJqFuh6jeiSluu060/YPSEnwalXfsZtYq2l6HhPJe0Wu6X3nsJ5OZ2bFmdmmKdUeZWa9S7me2mU2Jez7GzFaUZltx2xgUjd9XocxsuJn9qqL3WxgzO9fMvjazLWY2u7LjiWdmE8xsQdzzodFns1EJt5PS+6WM79FSxVZacd9Tbmark5Sn5fsr+izG9nthOvaRKZTcSTJrCAMudwOGECaf/h3wgZl1iau3NKqX6gj1dYDrgIPKLdKKU96xt4q216GctlcdHQtcmmLdUUCvtEVScoOAoZWw3+HAryphvwVE/yzeRxgM+Sjg/MqNqFjTCd93G9K0/bK8R9MdW2HOAPpW4P7OJxynFEPTj0kyW9z99bjnL5jZfcArwBNmtpe750ZTdL2efBMiBUUjtdd1942VHcv2wMxqA1ur6mw2ZlavDOeyI2FC9Yfc/f0i9lHf3XNKuY9y4+4/AT9VdhzJVGJs77v7hxW1M3f/GCB8jUhR1HInKXH31YT/LPcAjoHk3bJmdrKZLTSz9VGX7htmdlS0eF308+G4pvUO0XrjzOwDM8s2syVm9niSbuDFZnabmV0W1VllZpPMrFlCvR3M7H4zW2pmG81sUXz3nZnVMLM/mdkXZrbJzD4zs7OKeQmKir2lmT1iZj+b2Yao6yCrsA1F630QPX05tr2Eai3N7Mno9fjKzAq0aphZdzObE+3zZzN7wMwaF3UQsa4mMzvGzN6PztNrZrZvQr1iXyMz629mL5rZcjNba2avm9mxCXXGmNmKKNa3gI1EE2gXF7+ZNTOzf5rZD9F5/NbMHohtF/gjsGvc+ZhQyDEvJkz5c11c3V7Rsj9auOxgjZktM7NpFqbZKhMza2dmk6PXJsfMvjSzG6JlE4DfAEfFxTMmWjbbzKZY6D79Mnq9draELsKobrLPX00zuzI6X5uiz8mE2LaBLsBZcfsdGi0r0M1lCd2Jtq3r77AozhxgZLRsPzObbmbroseTlvD5Tdw28Gr09L1YLHHHdIaZPWqhy29atM5uZvbf6L22Ltm5ita9zMxuj95TK8xsRLTsrOiztNrMHjKzesWcxsSYC3R9mll7M3s+OsdfR3WmWJIuZjM7OPqMbDCzd8ysR9yyxRT+HnUzu8TM/mxmP0XvqXvMrG5hscW9joMsfBeuid4LY82sRkJcp5rZ59ExvBzFmffeKCszG2Jmv5jZedHz2HfC4Ra+i3IsfAftZuFylf9a+N77xMz6lEcM1ZFa7qQkXga2AF2B/yUuNLM9gCnAXYQv/XqEPyYtoip9gJeAGwndCBC6diF0U/4Z+AHYkfCH+yUz2z+h1WIQ8D6he6kd8NdovfOjGOoTJp5vBYwlTNLcMXrE/B04C7geeJuQrD5kZj+7+7OFHHtRsf832v4IYEV07C+b2cHu/kWSbS0ldGc8DlwQxZDoAeARwoTtpwH3mNkCd38zOs4jgVnRvgcS/jCMA5pHz4vSHrgVuAnIAW4DJpvZfr5tPsJUXqPdCH94bwO2AicAz5tZT3efG7e/BtGx3AJ8BvyQYvx/JVwacBnwI7AL0DNa9k+gE+G8DIjKCmu5GEB4706J1gP4OPrZjjDv8zdAE8I8nXPNbE93X1PI9lLxKFCf8D5dDewOdI6W3UA4B83Y1hW5JG7dIwn/RF1B6GYrSRz3A2cSXus5hM9e7PU8H/gP8FUUA8CXJdh2zL8J3aljgdVRgjWXMF/m7witcTcA08zsME8+x+U/geXAPYTPwldRLA2j5bcRLgc5FciNEplZwGbgXML30FhgTvQdsTJu238kfEZPA04EbjWzVsChwMWE1/4OwntxXCmOH8hrhZ5KOI/DCIn4NYTvr8TXNfYZuIPwXr4OeNrM2rv7Bop+j8aO6SXgt8ABwF8I79lbignzFsI5HwgcDVwLfARMjo4hC5gU7fciYG/gidRegeJFCeJ4YLi7T4hb1CAqvwVYT5iz+F/AJuB54F5CY8KTZrZL9BpJSbi7HnrkPYAxwIoili8F7ot+7wA4cGL0fCDwcxHrNorqDy0mhppA26huz7jyxYQvzVpxZXcCP8Y9/z0h0TiokG13jJaflVD+KPBWSWMHjo/Kj4ora0hINO4vYnv7Rev1SijvFZVfH1dWO9reuLiyV4GXE9btE627XxH7nUD4w9gpruxX0XqdS/saEXoBagEvELrZ4t9PDpySUL/Y+IEPgYuKOJbbgMUpvq9XAGNSeN/VJ7TSnhlXPhuYkupnJKqTDZxUxPIpwOwk5bMJCXebJOdtQUJZB/J//jpHzy8uYr8LgAlJyh24MKEs33ESrhF04JKEev8CFgF14so6AblA/yJiib3X90tyTE8n1D0vet/uHlfWDvgFuDLhOF6Oe16D8J21CmgSVz4ZeKOYc5jvNY87/kbR8/7R88Pi6rQlJKCz48pin4E+cWUHRWXHF/cejeq9klD2X+D1ImKLvY6PJqz3LjAp7vmThM+ZxZWNopjv6WTnLrE8OmebgCFJ3leJ35nnR2XXxpXtE5WdkMr7VY/8D3XLSkkVdbHDB0BTC12Ux5pZwyLq5t+o2QlmNs/M1hC+xGMtGXsmVH3Z3bfEPf8YaGVmdaLnfYB33P3dQnZ1NCFxedrMasUehFaBg8ysZqoxRw4DfnL3ObECd18PPAt0L+G24s2I295m4HPCHzPMrAHhouLJCcfwGuEPS5ck24u32N0/j3se34oFKb5GFroeHzGz7wnnbDPhJofEc+aE/8YpYfzvAiPN7HwzS9xmuTCzrha6ln+OjmEDIZEv6/7eBf4SdZe1L+G6C939x1Lss3f0c0Ip1i2J6QnP+wJPA1vjzuXXhH/GCr08oYT7OAx4292/ihW4+xJCi2Hi52xWXJ2tUSwL3X1tXJ0vCIlYWRxK+Mfyzbj9fQ8sTFJ3MyFxj0n8zBVnRsLzj1Nct7j1DgWmeZQxRaamGFNRLib84z3E3SclWf4L27rlIZwPCK2TiWVlPU/VkpI7SVl0jcoOwLJky919EXAKoQvqOWCFmU00sx2L2e6hhC+UJYRunSMIXb8QunbjrU54/gsh4Ywldzuwrbs0mZaEFpo1hC/c2GMCoeVpp6JiTWInkr8ey9jWHV0aqxOe/8K216I54RjuJf8xbCK08u1Sim0Tt/1iX6Poup2phG7TawmJxaGEJC7xnK1y91/inqca/4WEFoprgUXRdUFDijm2lEVJ1wzC++f3hO7QQwndhSW6HiuJwYRWsjuAb8zsXTM7OsV1k36+UrADsD4hiUmHxPhaErqQNyc8dqf492Kq+yjJ52x1wvNfCikr6zluQ/JLAZKVrY0STQDiPg+pxrA64Xmq8Re3XrJjKI8bM35DSM5mFrJ8XfzrwbbvoNWxglK8RhJH19xJSfQmvGfmF1bB3acD082sKaHb4k7C9VtF/VEeQPhCGRz7D9LMdi1ljD+T//q6RCsJLTRHElqnEi0v4f6WEq7vS9Q62lc6rCa0ho0hJNGJfijj9lN5jToCBxO6TPKuv4yueUyUeM3ValKI38NNPBcDF5vZAYTuosfN7H2P7poro+MJ1/6cErW2ErU6lSUpB/JacIZGSfBhhGOdGl1j9XNxqycp28i2f2BiEuP8GWhoZk1KkeBtSmH7hcW3ktBy988kdUs7HmDiPpYC+yapl87PWXF+JFxfl2hHwvnaHiQ7hiL/GU/RGYTrMqeZ2QleBe52rm7UcicpsXBH6s0U/d9YHndf4+4TCV/6+0TFhf0nVh/YnNA1cEYpQ50FHBwlA8m8RGg1auruC5I8filkvcJif4PQLRy70D/W7difosf/K/V/pVEi8jqwVyHHUNbkLpXXKJbEbYqtFCXkR6Yjfg9DZYwkfGfFbkwoSetLsrr1CclrfDf/IMrxn1533+phWKGxhEQy9k9LSVuOlgAdLP8dnsck1Il1aZ1ZxHYK2+8SwsX0QLhbmnCJQypmEa6xWpjkXC5OcRvFeQPoYma7xcXYltBynOo4m+XtLaCNmR2WEFNxl0UUpjxaE0vqLeCk6OaQmJPLYbtLCJd3dAKmWBjSRyqQWu4kmVpmFusWbUz4svoD4Y/T8V7ImFtm9ntCl+r/CK0vnQh3uz0KoZndzL4GBpnZh4T/bt8HXgQuNbM7CXdfdiPcFVYajxLuQJ1hYbiFRYS7Ovd09z+5+yIz+wcwycxuIXSd1SO0Cuzp7uck22hhsbv7C2Y2lzD+358IrScjCInDrUXE+S3hwvmzousMN7v7giLqJxoFzDKzrYSL89cR7gLsD4x2989KsK18UnyNPiV8gd9uZtcQ3idjge/LK34ze43wz8GHhJaccwl31sWucfoUaB3dkfch4eL/xYXs71Ogv5n9j3CzwyK2JbEPm9mD0fGNoGBXVolErdYvEN6LnwF1CXc7/gh8EhfPKRZmi1gC/FBMUv5fwp3L/7QwtMnBwP/FV4jO23jCOWlFGJeyGTDQ3WMt558Cx5nZcYT36tdRS+LTwAVm9g7hztVzCHcPp2IM4ZxMN7OHCK11bQnJ5wR3n53idooygdD1+7yZXUu4WWNMtK/7y2H7pfEc8B7h2tErCZ/n6whdxclavItT4D3q7uuKWaesbiYkzpPM7GFCgn9utKw0x5DH3b8ys76E9+FjZnZaQlespJFa7iSZpoSu13mEu6kGAo8B+7t7souFY94nNOn/lXAt09WEIT2uiKtzHuEanZmE/xp3dvfnojq/IVzHdRRhCIMS8zCgah9Ckng94RqwUeTvqryAMFTDmYQv6AmEpOKVYjZfIPaofAAhQb2T8HoZ4c64ZMOgxMd5LiFxnhNtL2Xu/hphWJAdCXcrTiMc53eU/pqteEW+Rh4GsP41odVrSlT3L4RjKa/45xPuApxCuLuxJaEbOHazzeQorlsIr9+YInY5kpAYTo/qdnH3DwgJ0uGEG2BOJ/wzUpYhUCAk/h8AlxDez48QbtQ4Nq576l7CZ+ShKJ7hRW3Qw0Cxwwj/PMU+I8OSVD2fkGT/lnDe7iQkHTE3EhLMydF+T4rKxxLeuzcSXtN3o9iKFf0j0TU6xvGEz9xYQqtuoZ+Bkojeb30JCdCDhNf0G8Ld5pXSLRv1NJwSxfQwYQio+wg3LZTmuscC79HyibRw0T+Up0X7+i/hO/gP0eIyX7vp7p8QbrI6DnggoYVQ0sjy94SJiIhIaUSttl8Bd7v7dZUdT2mY2W8J/3Dt7u5fF1KnF2FcvoOADwvrzUlDbDUJ/zxvJgyTdHdF7Hd7pG5ZERGRUrAw68JWwlBFOwKXE7rhU2r1rAosTC35ImEswEMIPS7TC0vsErxLaOlulq74EswitFpLMZTciYiIlM4mwiUl7QnXhb4J9HX3byo1qpLZgXCZwA6E6zCfIFwiUZSFhGGDIP8NSen2e8L1vRC65aUQ6pYVERERySC6oUJEREQkgyi5ExEREckgSu5EREREMoiSOxEREZEMouROREREJIMouRMRERHJIP8PZkO69hf8v8wAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 720x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "#plt.figure(figsize=(10,5), facecolor='lightyellow')\n",
    "plt.figure(figsize=(10,5))\n",
    "matplotlib.rcParams.update({'font.size': 15})\n",
    "plt.gca().yaxis.set_major_formatter(matplotlib.ticker.PercentFormatter(1))\n",
    "\n",
    "BINS = 20\n",
    "\n",
    "#plt.hist(np.amin(distances_wt,axis=1), bins=50)\n",
    "#plt.hist((np.amin(distances_wt, axis=1),np.amin(distances_mt, axis=1),np.amin(distances_tw, axis=1)),color=('orange','purple','red'),edgecolor='white',label=('wind-turbines','masts','towers'),bins=np.linspace(0, 10, 40))\n",
    "hh1 = np.histogram(np.amin(distances_wt, axis=1), bins= np.linspace(0, 10, BINS+1))\n",
    "hh2 = np.histogram(np.amin(distances_mt, axis=1), bins= np.linspace(0, 10, BINS+1))\n",
    "hh3 = np.histogram(np.amin(distances_tw, axis=1), bins= np.linspace(0, 10, BINS+1))\n",
    "hh4 = np.histogram(np.amin(distances_ch, axis=1), bins= np.linspace(0, 10, BINS+1))\n",
    "#plt.hist((np.amin(distances_wt, axis=1),np.amin(distances_mt, axis=1),np.amin(distances_tw, axis=1)),color=('black','orange','red'),edgecolor='white',label=('wind-turbines','masts','towers'),bins=np.linspace(0, 10, 40),density=True)\n",
    "plt.ylim(0,0.6)\n",
    "plt.bar(hh1[1][:-1]+0.05,hh1[0]/celkem, width=0.1, label='Wind-turbines', color='black')\n",
    "plt.bar(hh2[1][:-1]+0.05+0.12,hh2[0]/celkem, width=0.1, label='Masts', color='orange')\n",
    "plt.bar(hh3[1][:-1]+0.05+2*0.12,hh3[0]/celkem, width=0.1, label='Towers', color='red')\n",
    "plt.bar(hh4[1][:-1]+0.05+3*0.12,hh4[0]/celkem, width=0.1, label='Chimneys', color='chocolate')\n",
    "plt.legend()\n",
    "plt.xlabel('Distance to the nearest tall structure from lightning [km]')\n",
    "plt.ylabel('Relative abundance')\n",
    "plt.title(sitename.split()[0])\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "ipython_py3"
  },
  "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.8.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}