{ "cells": [ { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "from netCDF4 import Dataset\n", "import h5py\n", "import matplotlib.pyplot as plt\n", "import matplotlib.colors as mcolors\n", "import matplotlib.colors as Normalize\n", "import cartopy.crs as crs\n", "import cartopy.feature as cfeature\n", "from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter\n", "import matplotlib.ticker as mticker\n", "import matplotlib as mpl\n", "import xarray as xr\n", "import netCDF4\n", "import numpy as np\n", "import pandas as pd\n", "import glob\n", "import dask\n", "import os\n", "import re\n", "import warnings\n", "import gc\n", "from datetime import datetime\n", "import seaborn as sns\n", "from matplotlib.colors import LinearSegmentedColormap, TwoSlopeNorm\n", "import requests\n", "from io import StringIO\n", "from datetime import timedelta\n" ] }, { "cell_type": "code", "execution_count": 26, "metadata": {}, "outputs": [], "source": [ "####### Read in ASOS data ######\n", "month='06'\n", "directory= f'/pscratch/sd/d/dbrooks/acc2017_analysis/model_evaluation/asos_data/{month}'\n", "asos_files = sorted(glob.glob(os.path.join(directory, f\"*wind_mean_{month}.csv\")))\n", "\n", "datasets=[]\n", "for file in asos_files:\n", " ds = pd.read_csv(file, parse_dates=[\"anchor_time\"])\n", " datasets.append(ds)\n", "\n", "asos_df = pd.concat(datasets, ignore_index=True)\n", "asos_df = asos_df.drop([0])\n", "\n", "######## Read in WRF data ########\n", "filename = f'/pscratch/sd/d/dbrooks/acc2017_analysis/wind_data/Current/wind_speed_data_month{month}.nc'\n", "ds = xr.open_dataset(filename)" ] }, { "cell_type": "code", "execution_count": 27, "metadata": {}, "outputs": [], "source": [ "def find_max_asos_wind_gusts(df):\n", "\n", " # Convert necessary columns to appropriate data types\n", " df[\"sped\"] = pd.to_numeric(df[\"sped\"], errors=\"coerce\") # Convert gust to numeric, set invalid to NaN\n", " df[\"lon\"] = pd.to_numeric(df[\"lon\"], errors=\"coerce\")\n", " df[\"lat\"] = pd.to_numeric(df[\"lat\"], errors=\"coerce\")\n", " df[\"valid\"] = pd.to_datetime(df[\"valid\"], errors=\"coerce\") # Convert valid column to datetime\n", "\n", " # Drop rows with invalid or missing values in relevant columns\n", " df = df.dropna(subset=[\"sped\", \"lon\", \"lat\", \"valid\"])\n", "\n", " # Find the maximum wind gust for each station\n", " max_gusts = df.loc[df.groupby(\"station\")[\"sped\"].idxmax()]\n", "\n", " # Select and rename the columns for the final output\n", " result = max_gusts[[\"station\", \"sped\", \"lat\", \"lon\", \"valid\"]].rename(\n", " columns={\"sped\": \"max_speed\"}\n", " )\n", "\n", " # Reset the index for a clean output\n", " result = result.reset_index(drop=True)\n", "\n", " return result" ] }, { "cell_type": "code", "execution_count": 28, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
stationlonlat
005U-116.005139.6042
106D-99.620848.8844
208D-102.406448.3008
30A9-82.173436.3712
40CO-105.763939.7939
............
5068SCMK-73.7389-43.8950
5069SCON-73.6333-43.1167
5070SCRD-71.5830-33.0500
5071SCSN-71.6144-33.6547
5072TNCE-62.979417.4965
\n", "

5073 rows × 3 columns

\n", "
" ], "text/plain": [ " station lon lat\n", "0 05U -116.0051 39.6042\n", "1 06D -99.6208 48.8844\n", "2 08D -102.4064 48.3008\n", "3 0A9 -82.1734 36.3712\n", "4 0CO -105.7639 39.7939\n", "... ... ... ...\n", "5068 SCMK -73.7389 -43.8950\n", "5069 SCON -73.6333 -43.1167\n", "5070 SCRD -71.5830 -33.0500\n", "5071 SCSN -71.6144 -33.6547\n", "5072 TNCE -62.9794 17.4965\n", "\n", "[5073 rows x 3 columns]" ] }, "execution_count": 28, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# keep only unique station-coordinate pairs\n", "df_stations_unique = asos_df[[\"station\", \"lon\", \"lat\"]].drop_duplicates().reset_index(drop=True)\n", "df_stations_unique" ] }, { "cell_type": "code", "execution_count": 29, "metadata": {}, "outputs": [], "source": [ "#max_df = find_max_asos_wind_gusts(asos_df)\n", "\n", "from shapely.geometry import Point, MultiPoint\n", "from shapely.prepared import prep\n", "\n", "# Step 1: Boolean mask for ds1 footprint (first time slice)\n", "mask1_bool = ~np.isnan(ds['wspd_wdir10'].sel(wspd_wdir='wspd').isel(Time=0))\n", "\n", "# Step 2: Extract lat/lon arrays from ds1\n", "lat1 = ds.XLAT.data[mask1_bool.data]\n", "lon1 = ds.XLONG.data[mask1_bool.data]\n", "\n", "# Step 3: Build polygon from ds1 valid points\n", "points = MultiPoint(list(zip(lon1, lat1)))\n", "polygon = points.convex_hull # or cascaded_union for exact pixel shapes\n", "prep_poly = prep(polygon) # speeds up \"contains\" checks\n", "\n", "# Step 4: Filter DataFrame rows inside polygon\n", "df_stations_unique = df_stations_unique[df_stations_unique.apply(lambda row: prep_poly.contains(Point(row['lon'], row['lat'])), axis=1)]" ] }, { "cell_type": "code", "execution_count": 30, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " Unnamed: 0 station anchor_time lon lat sped\n", "0 23 0CO 2017-06-01 03:00:00 -105.7639 39.7939 3.084576\n", "1 24 0CO 2017-06-01 06:00:00 -105.7639 39.7939 7.197344\n", "2 25 0CO 2017-06-01 09:00:00 -105.7639 39.7939 3.598672\n", "3 26 0CO 2017-06-01 12:00:00 -105.7639 39.7939 5.655056\n", "4 27 0CO 2017-06-01 15:00:00 -105.7639 39.7939 2.056384\n", "... ... ... ... ... ... ...\n", "4885090 744752 YKN 2017-06-30 12:00:00 -97.3859 42.9167 2.570480\n", "4885091 744753 YKN 2017-06-30 15:00:00 -97.3859 42.9167 4.626864\n", "4885092 744754 YKN 2017-06-30 18:00:00 -97.3859 42.9167 5.655056\n", "4885093 744755 YKN 2017-06-30 21:00:00 -97.3859 42.9167 7.711440\n", "4885094 744756 YKN 2017-07-01 00:00:00 -97.3859 42.9167 6.683248\n", "\n", "[4885095 rows x 6 columns]\n" ] } ], "source": [ "# Filter original df to keep only stations within model domain\n", "# df_stations_unique contains the stations to keep\n", "df_filtered = asos_df[asos_df[\"station\"].isin(df_stations_unique['station'])].reset_index(drop=True)\n", "print(df_filtered)" ] }, { "cell_type": "code", "execution_count": 31, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " Size: 3GB\n", "Dimensions: (Time: 240, wspd_wdir: 2, south_north: 1080, west_east: 1210)\n", "Coordinates:\n", " * Time (Time) datetime64[ns] 2kB 2017-06-01 ... 2017-06-30T21:00:00\n", " * wspd_wdir (wspd_wdir) = 25\n", " df1 = df1.groupby(by='station').agg('mean') # average by station\n", " \n", " if station_data == True:\n", " wind_speed_current1 = df1.where(df1['sped'] < 17)\n", " wind_speed_current1 = wind_speed_current1.dropna()\n", " wind_speed_current2 = df1.where(df1['sped'] >=17)\n", " wind_speed_current2 = wind_speed_current2.dropna()\n", "\n", " lats1 = wind_speed_current1['lat']\n", " lons1 = wind_speed_current1['lon']\n", "\n", " lats2 = wind_speed_current2['lat']\n", " lons2 = wind_speed_current2['lon']\n", "\n", " mean = df1['sped'].mean()\n", " max = df1['sped'].max()\n", "\n", " pb = ax_wind.scatter(lons1, lats1, c=wind_speed_current1['sped'], cmap=p_cmap, norm=p_norm, linewidths=0.6,edgecolors='black',s=30,transform=crs.PlateCarree(), zorder=3, alpha=0.9)\n", " ax_wind.scatter(lons2, lats2, c=wind_speed_current2['sped'], cmap=p_cmap, norm=p_norm, linewidths=0.6,edgecolors='black',s=30,transform=crs.PlateCarree(), zorder=4, alpha=0.9)\n", " sim='ASOS Stations'\n", " elif station_data == False:\n", " wind_speed_current1 = df1['wrf_sped']\n", " wind_speed_current1 = wind_speed_current1.dropna()\n", " #wind_speed_current2 = df1.where(df1['wrf_sped'] >=17)\n", " #wind_speed_current2 = wind_speed_current2.dropna()\n", "\n", " lats1 = df1['lat']\n", " lons1 = df1['lon']\n", "\n", " #lats2 = wind_speed_current2['lat']\n", " #lons2 = wind_speed_current2['lon']\n", "\n", " mean = df1['wrf_sped'].mean()\n", " max = df1['wrf_sped'].max()\n", "\n", " pb = ax_wind.scatter(lons1, lats1, c=wind_speed_current1, cmap=p_cmap, norm=p_norm, linewidths=0.6,edgecolors='black',s=30,transform=crs.PlateCarree(), zorder=3, alpha=0.9)\n", " #ax_wind.scatter(lons2, lats2, c=wind_speed_current2['wrf_sped'], cmap=p_cmap, norm=p_norm, linewidths=0.6,edgecolors='black',s=30,transform=crs.PlateCarree(), zorder=4, alpha=0.9)\n", " sim='Current'\n", " # Add features and title\n", " ax_wind.add_feature(cfeature.STATES, edgecolor=\"black\", linewidths=0.35, alpha=0.8)\n", " ax_wind.add_feature(cfeature.COASTLINE, edgecolor=\"black\", linewidths=0.35, alpha=0.8)\n", " gl = ax_wind.gridlines(draw_labels=True, linewidth=0.8, color='gray', alpha=0.5, linestyle='--', zorder=2)\n", " gl.top_labels = False\n", " gl.right_labels = False\n", " gl.left_labels = True\n", " gl.bottom_labels = True\n", "\n", " # Set month name based on the list provided\n", " if month == '04':\n", " month = 'April'\n", " elif month == '05':\n", " month = 'May'\n", " elif month == '06':\n", " month = 'June'\n", "\n", " ax_wind.set_title(f'Max Wind Speed (m s$^{{-1}}$)', loc='left', fontsize=12)\n", " ax_wind.set_title(f'({sim}, {month})', loc='right', fontsize=12, weight='bold')\n", "\n", " text_str = (\n", " f\"Mean: {mean:.2f} m s$^{{-1}}$\\n\"\n", " f\"Max: {max:.2f} m s$^{{-1}}$\"\n", " )\n", " ax_wind.text(\n", " -113.5, 26.7, text_str, # Adjust coordinates as needed\n", " fontsize=8, color=\"black\", weight='bold',transform=crs.PlateCarree(),\n", " bbox=dict(facecolor=\"white\", alpha=0.9, boxstyle=\"round,pad=0.3\")\n", " )\n", " \n", " # Colorbar for the left plot\n", " cbar = plt.colorbar(pb, ax=ax_wind, orientation='vertical', fraction=0.05, pad=0.01, shrink=0.6, extend='both')\n", " #cbar.set_label('m $^{-1}$')\n", "\n", " plt.show()\n", "\n", "#plot_max_wind_asos(df_winds, station_data=False)\n", "#plot_max_wind_asos(df_winds, station_data=True)" ] }, { "cell_type": "code", "execution_count": 35, "metadata": {}, "outputs": [], "source": [ "\n", "def plot_max_wind(ds, month=month):\n", " # Create the plot\n", " fig, ax_wind = plt.subplots(\n", " figsize=(6, 6),\n", " subplot_kw={\"projection\": crs.PlateCarree()}\n", " )\n", " #wind_speed_current = ds['wspd_wdir10'].sel(wspd_wdir='wspd') \n", " wind_speed_current = ds['wspd_wdir10'].sel(wspd_wdir='wspd') \n", " #wind_speed_current = wind_speed_current.where(wind_speed_current >= 25)\n", " max_wind_current = wind_speed_current.mean(dim='Time') # Maximum wind speed over time for each grid point\n", " lats = ds['XLAT']\n", " lons = ds['XLONG']\n", "\n", " #ax_wind.set_extent([-97.953608,-88.879741,28,32.545423])\n", " #ax_wind.set_extent([-93,-90,29,30.5])\n", "\n", " #max_precip_threshold = max_wind_current >= 25\n", " \n", " #ax_wind = axs[0]\n", " #pb = ax_wind.pcolormesh(lons, lats, max_wind_current, cmap=p_cmap,norm=p_norm, transform=ccrs.PlateCarree(), zorder=1)\n", " pb = ax_wind.contourf(lons, lats, max_wind_current, cmap=p_cmap, norm=p_norm, levels=p_clevs, transform=crs.PlateCarree(), zorder=1, extend='both')\n", " ax_wind.contour(lons, lats, max_wind_current, cmap=p_cmap, norm=p_norm, levels=p_clevs, transform=crs.PlateCarree(), zorder=1, extend='both')\n", " \n", " # Add features and title\n", " ax_wind.add_feature(cfeature.STATES, edgecolor=\"black\", linewidths=0.65, alpha=0.7)\n", " ax_wind.add_feature(cfeature.COASTLINE, edgecolor=\"black\", linewidths=0.65, alpha=0.7)\n", " #ax_wind.add_feature(USCOUNTIES, linewidth=0.35, alpha=0.3)\n", " gl = ax_wind.gridlines(draw_labels=True, linewidth=0.8, color='gray', alpha=0.35, linestyle='--', zorder=2)\n", " gl.top_labels = False\n", " gl.right_labels = False\n", " gl.left_labels = True\n", " gl.bottom_labels = True\n", "\n", " if month == '04':\n", " month = 'April'\n", " elif month == '05':\n", " month = 'May'\n", " elif month == '06':\n", " month = 'June'\n", "\n", " ax_wind.set_title(f'Max 10m Wind Speed (m s$^{{-1}}$)', loc='left', fontsize=12)\n", " ax_wind.set_title(f'(Current, {month})', loc='right', fontsize=12, weight='bold')\n", " \n", " # Colorbar for the left plot\n", " cbar = plt.colorbar(pb, ax=ax_wind, orientation='vertical', fraction=0.05, pad=0.01, shrink=0.6, extend='both')\n", " #cbar.set_label('m $s^{-1}$')\n", "\n", " plt.show()\n", "\n", "#plot_max_wind(ds)" ] }, { "cell_type": "code", "execution_count": 36, "metadata": {}, "outputs": [], "source": [ "def plot_max_wind_diff(df1, month=month):\n", " fig, ax_wind = plt.subplots(\n", " figsize=(6, 6),\n", " subplot_kw={\"projection\": crs.PlateCarree()}\n", " )\n", " \n", "\n", " df1 = df1.groupby(by='station').agg('mean') # average over time by station\n", " \n", "\n", " wind_speed_asos = df1['sped'] \n", " wind_speed_wrf = df1['wrf_sped']\n", "\n", " lats1 = df1['lat']\n", " lons1 = df1['lon']\n", " \n", " wind_diff = wind_speed_wrf - wind_speed_asos\n", "\n", " mean = wind_diff.mean()\n", " max = wind_diff.max()\n", "\n", " pb = ax_wind.scatter(lons1, lats1, c=wind_diff, cmap='coolwarm', vmin=-3,vmax=3, linewidths=0.6,edgecolors='black',s=30,transform=crs.PlateCarree(), zorder=3, alpha=1)\n", " sim='Current - ASOS'\n", "\n", " # Add features and title\n", " ax_wind.add_feature(cfeature.STATES, edgecolor=\"black\", linewidths=0.35, alpha=0.8)\n", " ax_wind.add_feature(cfeature.COASTLINE, edgecolor=\"black\", linewidths=0.35, alpha=0.8)\n", " gl = ax_wind.gridlines(draw_labels=True, linewidth=0.8, color='gray', alpha=0.5, linestyle='--', zorder=2)\n", " gl.top_labels = False\n", " gl.right_labels = False\n", " gl.left_labels = True\n", " gl.bottom_labels = True\n", "\n", " # Set month name based on the list provided\n", " if month == '04':\n", " month = 'April'\n", " elif month == '05':\n", " month = 'May'\n", " elif month == '06':\n", " month = 'June'\n", "\n", " ax_wind.set_title(f'Mean Wind Speed Difference (m s$^{{-1}}$)\\n{sim}', loc='left', fontsize=12)\n", " ax_wind.set_title(f'({month})', loc='right', fontsize=12, weight='bold')\n", "\n", " text_str = (\n", " f\"Mean: {mean:.2f} m s$^{{-1}}$\"\n", " )\n", " ax_wind.text(\n", " -113.5, 26.7, text_str, # Adjust coordinates as needed\n", " fontsize=8, color=\"black\", weight='bold',transform=crs.PlateCarree(),\n", " bbox=dict(facecolor=\"white\", alpha=0.9, boxstyle=\"round,pad=0.3\")\n", " )\n", " \n", " # Colorbar for the left plot\n", " cbar = plt.colorbar(pb, ax=ax_wind, orientation='vertical', fraction=0.05, pad=0.01, shrink=0.6, extend='both')\n", " #cbar.set_label('m $^{-1}$')\n", "\n", " plt.show()\n", "\n", "#plot_max_wind_diff(df_winds)" ] }, { "cell_type": "code", "execution_count": 37, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/tmp/ipykernel_470328/2331693398.py:54: UserWarning: Attempt to set non-positive ylim on a log-scaled axis will be ignored.\n", " ax1.set_ylim(0, 500) # Since we're dealing with percentages\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYYAAAGGCAYAAAB/gCblAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAAA/NklEQVR4nO3dd1QU5+I+8GfpVSxIFbGABUXsXQFFCSZ2jXqjYonXGLChXjVRsUaD0WCMVxK7uZqrsRuNogiaWFBUsGHHYJAmCEhRYJnfH37ZX/aCuAOzuyDP5xzOyb4zzD4L0ceZd4pMEAQBRERE/0dH2wGIiKhyYTEQEZESFgMRESlhMRARkRIWAxERKWExEBGREhYDEREpYTEQEZESFgMRESlhMRARkRIWA1EFbd++HTKZDFFRUdqOQiQJFgMRESlhMRARkRIWA5HEPDw84OHhUWJ83LhxaNCggeL1kydPIJPJ8M033+DHH39E48aNYWhoiA4dOuDKlSslvv/u3bsYNmwYateuDSMjI7Rv3x5HjhxR4yeh6kpP2wGIqrvdu3fj5cuXmDx5MmQyGYKCgjBkyBA8fvwY+vr6AIDbt2+jW7dusLe3x7x582Bqaoq9e/di0KBB2L9/PwYPHqzlT0HvExYDkZbFx8fjwYMHqFWrFgCgadOmGDhwIE6ePImPPvoIADB9+nTUr18fV65cgaGhIQDg888/R/fu3TF37lwWA0mKh5KItGzEiBGKUgCAHj16AAAeP34MAEhPT8eZM2fw8ccf4+XLl3j+/DmeP3+OtLQ0eHt748GDB0hISNBKdno/cY+BSMvq16+v9Lq4JF68eAEAePjwIQRBwMKFC7Fw4cJSt5GSkgJ7e3v1BqVqg8VAJDGZTIbSnpgrl8tLXV9XV7fU8eJtFBUVAQBmz54Nb2/vUtd1cnIqT1SiUrEYiCRWq1YtxWGgv/vzzz/Ltb1GjRoBAPT19eHl5VWhbESq4BwDkcQaN26Mu3fvIjU1VTEWExOD8+fPl2t7VlZW8PDwwA8//IDExMQSy//+PkRS4B4DkcQmTJiAtWvXwtvbGxMnTkRKSgpCQkLQokULZGVllWubGzZsQPfu3eHq6opJkyahUaNGSE5OxsWLF/HXX38hJiZG4k9B1Rn3GIgqqHguoHiuoHnz5ti5cycyMzMREBCAI0eO4KeffkLbtm3L/R4uLi6IiorChx9+iO3bt8PPzw8hISHQ0dHBokWLJPkcRMVkQmmzZESksu+++w7Tp0/Hw4cP0bhxY23HIaow7jEQVdCVK1dgamoKR0dHbUchkgTnGIjKaf/+/YiIiMCuXbvw6aefQk+Pf5zo/cBDSUTl1LBhQ7x8+RKDBw9GcHAwTE1NtR2JSBIsBiIiUsI5BiIiUsJiICIiJZwtK0NRURGePXsGc3NzyGQybcchIio3QRDw8uVL2NnZQUen7H0CFkMZnj17BgcHB23HICKSzNOnT1GvXr0y12ExlMHc3BzAmx9kjRo1tJyGiKj8srKy4ODgoPh7rSwshjIUHz6qUaMGi4GI3guqHBbn5DMRESlhMRARkRIWAxERKeEcAxFViFwuR0FBgbZjVHv6+vpvfUysWCwGIioXQRCQlJSEjIwMbUeh/1OzZk3Y2NhU+LorFgMRlUtxKVhZWcHExIQXgWqRIAjIzc1FSkoKAMDW1rZC22MxEJFocrlcUQp16tTRdhwCYGxsDABISUmBlZVVhQ4rcfKZiEQrnlMwMTHRchL6u+LfR0XnfFgMRFRuPHxUuUj1+2AxEBGREhYDEREp4eQzEUlmyZIlGn2/wMDAcn3fxYsX0b17d3zwwQc4duyY0rKDBw/i66+/RmxsLIqKilC/fn306dMHwcHBinXy8vKwatUq/Pzzz/jzzz9hbm4OT09PLF68GC1atFCsl5ubi2XLlmHv3r1ISEiAubk5XFxcEBAQgIEDB5YruyZwj4GIqp0tW7Zg6tSpOHfuHJ49e6YYDwsLw4gRIzB06FBcvnwZV69exYoVK5Qmc1+/fg0vLy9s3boVy5cvx/3793H8+HEUFhaiU6dOuHTpkmLdzz77DAcOHMD69etx9+5dnDhxAsOGDUNaWppGP69Y3GMgomolOzsbe/bsQVRUFJKSkrB9+3Z88cUXAICjR4+iW7dumDNnjmL9Jk2aYNCgQYrXwcHBuHjxIq5fvw43NzcAgKOjI/bv349OnTph4sSJuHXrFmQyGY4cOYJ169ahX79+AIAGDRqgXbt2mvuw5cQ9BiKqVvbu3YtmzZqhadOmGD16NLZu3QpBEAAANjY2uH37Nm7duvXW79+9ezf69OmjKIViOjo6mDlzJu7cuYOYmBjF9o4fP46XL1+q7wOpAYuBiKqVLVu2YPTo0QCADz74AJmZmTh79iwAYOrUqejQoQNcXV3RoEEDjBw5Elu3bsXr168V33///n00b9681G0Xj9+/fx8A8OOPP+LChQuoU6cOOnTogJkzZ+L8+fPq/HiSYDEQUbVx7949XL58GaNGjQIA6OnpYcSIEdiyZQsAwNTUFMeOHcPDhw+xYMECmJmZYdasWejYsSNyc3MV2ynew3iXnj174vHjxwgLC8OwYcNw+/Zt9OjRA8uWLZP+w0mIxUBE1caWLVtQWFgIOzs76OnpQU9PDxs3bsT+/fuRmZmpWK9x48b49NNPsXnzZly7dg137tzBnj17ALyZc4iNjS11+8XjTZo0UYzp6+ujR48emDt3LkJDQ7F06VIsW7YM+fn5avykFcNiIKJqobCwEDt37sSaNWsQHR2t+IqJiYGdnR1+/vnnUr+vQYMGMDExQU5ODgBg5MiROH36tGIeoVhRURG+/fZbuLi4lJh/+DsXFxcUFhbi1atX0n04ifGsJCKqFn799Ve8ePECEydOhIWFhdKyoUOHYsuWLUhKSkJubi769esHR0dHZGRk4LvvvkNBQQH69OkDAJg5cyYOHz6M/v37Y82aNejUqROSk5Px1VdfITY2FqdPn1bcmsLDwwOjRo1C+/btUadOHdy5cwdffPEFPD09K/Vz5LnHQETVwpYtW+Dl5VWiFIA3xRAVFYVatWrh8ePHGDt2LJo1awYfHx8kJSUhNDQUTZs2BQAYGRnhzJkzGDt2LL744gs4OTnhgw8+gK6uLi5duoTOnTsrtuvt7Y0dO3agb9++aN68OaZOnQpvb2/s3btXY5+7PGSCqrMopXj9+jUMDQ2lzFOpZGVlwcLCApmZmZW63Yk07dWrV4iLi0PDhg1hZGSk7Tj0f8r6vYj5+0zUHsNvv/0GX19fNGrUCPr6+jAxMUGNGjXg7u6OFStWKF1BSEREVZNKxXDw4EE0adIEEyZMgJ6eHubOnYsDBw7g5MmT2Lx5M9zd3XH69Gk0atQIn332GVJTU9Wdm4iI1ESlyeegoCB8++238PHxgY5OyS75+OOPAQAJCQlYv349/vOf/2DmzJnSJiUiIo1QqRguXryo0sbs7e2xatWqCgUiIiLtqvBZSTk5OcjKypIiCxERVQLlLoY7d+6gffv2MDc3R61ateDq6oqoqCgpsxERkRaUuxgmT54Mf39/ZGdnIy0tDUOGDIGvr6+U2YiISAtULoaBAwciISFB8To1NRUDBgyAiYkJatasiX79+iE5OVktIYmISHNUviXG6NGj0atXL/j5+WHq1Knw9/dHixYt4O7ujoKCApw5cwazZs1SZ1YiItIAlfcYhg8fjsuXL+POnTvo3LkzunXrhtDQUHTr1g09evRAaGgoFixYoM6sRESkAaJuomdhYYGQkBD88ccf8PX1RZ8+fbBs2TKYmJioKx8RVSHFN4/TlPLe0ScpKQkrVqzAsWPHkJCQACsrK7Ru3RozZsxA7969JU4pDZlMhoMHDyo9ZlRdRE0+p6en4+rVq3B1dcXVq1dRo0YNtGnTBsePH1dXPiIiST158gTt2rXDmTNnsHr1aty8eRMnTpyAp6cn/Pz8yrVNQRBQWFhYYrwyP3OhTIKKdu3aJRgbGwvW1taChYWFcPjwYUEQBCE2NlZwd3cXhg8fLiQlJam6uSohMzNTACBkZmZqOwpRpZKXlyfcuXNHyMvLUxoHoNGv8vDx8RHs7e2F7OzsEstevHghxMXFCQCE69evK40DEMLDwwVBEITw8HABgHD8+HGhbdu2gr6+vhAeHi64u7sLfn5+wvTp04U6deoIHh4egiAIws2bN4UPPvhAMDU1FaysrITRo0cLqampiu27u7sLU6dOFebMmSPUqlVLsLa2FgIDAxXLHR0dlT63o6OjqN+LIIj7+0zlPYb58+dj69atSEpKQlhYGBYuXAgAaNasGSIiItCnTx906dKlXOVERKQJ6enpOHHiBPz8/GBqalpiec2aNUVtb968eVi1ahViY2PRqlUrAMCOHTtgYGCA8+fPIyQkBBkZGejVqxfatGmDqKgonDhxAsnJyYpbCRXbsWMHTE1NERkZiaCgICxduhSnTp0CAFy5cgUAsG3bNiQmJipeq4vKcwzZ2dmK+5E3btxY6fmnADBp0iQMHDhQ2nRERBJ6+PAhBEFAs2bNJNne0qVLFQ/wKebs7IygoCDF6+XLl6NNmzb46quvFGNbt26Fg4MD7t+/r3gMaKtWrRAYGKjYxvfff4+wsDD06dMHdevWBfCmuGxsbCTJXhaVi8HX1xcffvghPDw8EBUVhTFjxpRYx8rKStJwRERSEsr/+JlStW/fvsRYu3btlF7HxMQgPDwcZmZmJdZ99OiRUjH8na2tLVJSUiRMqzqVi2Ht2rXw9PTE3bt3MW7cOPTt21eduYiIJOfs7AyZTIa7d+++dZ3iO0j/vUQKCgpKXbe0w1H/O5adnY3+/fvj66+/LrGura2t4r/19fWVlslkMhQVFb01pzqJOiupf//+mDNnTpUqhadPn8LDwwMuLi5o1aoVfvnlF21HIiItqV27Nry9vbFhwwbk5OSUWJ6RkaE4bJOYmKgYj46OLvd7tm3bFrdv30aDBg3g5OSk9FVasbyNvr4+5HJ5uXOIoVIx/Pe//1V5g0+fPsX58+fLHUhqenp6CA4Oxp07dxAaGooZM2aU+j8EEVUPGzZsgFwuR8eOHbF//348ePAAsbGx+O6779ClSxcYGxujc+fOiknls2fPVujiXT8/P6Snp2PUqFG4cuUKHj16hJMnT2L8+PGi/qJv0KABwsLCkJSUhBcvXpQ7jypUKoaNGzeiefPmCAoKQmxsbInlmZmZOH78OP7xj3+gbdu2SEtLkzxoedna2qJ169YAABsbG1haWiI9PV27oYhIaxo1aoRr167B09MTs2bNQsuWLdGnTx+EhYVh48aNAN5MDhcWFqJdu3aYMWMGli9fXu73s7Ozw/nz5yGXy9G3b1+4urpixowZqFmzZqkPPnubNWvW4NSpU3BwcECbNm3KnUcl7zyh9f8cPnxY8PLyEnR0dARzc3PByclJaNmypWBvby/o6uoK1tbWwty5cyW/luHs2bPCRx99JNja2goAhIMHD5ZY5/vvvxccHR0FQ0NDoWPHjkJkZGSp24qKihJatGih8nvzOgai0pV1vjxpj1TXMag8+TxgwAAMGDAAz58/xx9//IE///wTeXl5sLS0RJs2bdCmTRtR7aeqnJwcuLm5YcKECRgyZEiJ5Xv27EFAQABCQkLQqVMnBAcHw9vbG/fu3VM6Syo9PR1jx47Fpk2bJM9IRPQ+kQmCxOdvqVFp9wrp1KkTOnTogO+//x4AUFRUBAcHB0ydOhXz5s0DALx+/Rp9+vTBpEmTSj3N9m2ysrJgYWGBzMxM1KhRQ9LPQlSVvXr1CnFxcWjYsCGMjIy0HYf+T1m/FzF/n0n/T3wNys/Px9WrV+Hl5aUY09HRgZeXl+I51YIgYNy4cejVq9c7S+H169fIyspS+iIiqm5E3V21snn+/Dnkcjmsra2Vxq2trRXnKZ8/fx579uxBq1atcOjQIQDATz/9BFdX1xLbW7lyJZYsWaL23Jqi7c9SfBUnEVUtVboYVNG9e3eVLxKZP38+AgICFK+zsrLg4OCgrmhERJVSlS4GS0tL6OrqlnikaHJycrnuJ2JoaAhDQ0Op4hG997R1ZS6VTqrfR4WLQS6X4+bNm3B0dEStWrWkyKQyAwMDtGvXDmFhYYoJ6aKiIoSFhcHf31+jWYiqEwMDA+jo6ODZs2eoW7cuDAwMNP6QHvr/BEFAfn4+UlNToaOjAwMDgwptT3QxzJgxA66urpg4cSLkcjnc3d1x4cIFmJiY4Ndff4WHh0eFAv2v7OxsPHz4UPE6Li4O0dHRqF27NurXr4+AgAD4+vqiffv26NixI4KDg5GTk4Px48dLmkOsyvCHZPHixdqOQO8pHR0dNGzYEImJiXj27Jm249D/MTExQf369St86YDoYti3bx9Gjx4NADh69Cji4uJw9+5d/PTTT/jyyy8lvx1GVFQUPD09Fa+L5wB8fX2xfft2jBgxAqmpqVi0aBGSkpLQunVrnDhxosSENBFJy8DAAPXr10dhYaHG7uFDb6erqws9PT1J/lEquhieP3+uOH5//PhxDB8+HE2aNMGECROwbt26Cgf6Xx4eHu+8Va6/vz8PHRFpgUwmg76+fok7g1LVJnp/w9raGnfu3IFcLseJEycUD6nIzc2Frq6u5AGJiEizRO8xjB8/Hh9//DFsbW0hk8kUF5dFRkZK9lQkIiLSHtHFsHjxYrRs2RJPnz7F8OHDFad36urqKm5BQUREVVe5TlcdNmwYgDf35Sjm6+srTSIiItIq0XMMcrkcy5Ytg729PczMzPD48WMAwMKFC7FlyxbJAxIRkWaJLoYVK1Zg+/btCAoKUrqIomXLlti8ebOk4YiISPNEF8POnTvx448/4pNPPlE6C8nNza3MB2wTEVHVILoYEhIS4OTkVGK8qKgIBQUFkoQiIiLtEV0MLi4u+P3330uM79u3T/3PISUiIrUTfVbSokWL4Ovri4SEBBQVFeHAgQO4d+8edu7ciV9//VUdGYmISINE7zEMHDgQR48exenTp2FqaopFixYhNjYWR48eVVwFTUREVVe5rmPo0aMHTp06JXUWIiKqBETvMVy5cgWRkZElxiMjIxEVFSVJKCIi0h7RxeDn54enT5+WGE9ISICfn58koYiISHtEF8OdO3fQtm3bEuNt2rTBnTt3JAlFRETaI7oYDA0NSzxjGQASExOhp1elHyFNREQoRzH07dsX8+fPR2ZmpmIsIyMDX3zxBc9KIiJ6D4j+J/4333yDnj17wtHRUXFBW3R0NKytrfHTTz9JHpCIiDRLdDHY29vjxo0b2LVrF2JiYmBsbIzx48dj1KhRfLwfEdF7oFyTAqampvjnP/8pdRYiIqoEylUMDx48QHh4OFJSUlBUVKS0bNGiRZIEIyIi7RBdDJs2bcKUKVNgaWkJGxsbyGQyxTKZTMZiICKq4kQXw/Lly7FixQrMnTtXHXmIiEjLRJ+u+uLFCwwfPlwdWYiIqBIQXQzDhw9HaGioOrIQEVElIPpQkpOTExYuXIhLly7B1dW1xCmq06ZNkywcERFpnuhi+PHHH2FmZoazZ8/i7NmzSstkMhmLgYioihNdDHFxcerIQURElYToOYZi+fn5uHfvHgoLC6XMQ0REWia6GHJzczFx4kSYmJigRYsWiI+PBwBMnToVq1atkjwgERFpluhimD9/PmJiYhAREQEjIyPFuJeXF/bs2SNpOCIi0jzRcwyHDh3Cnj170LlzZ6Wrnlu0aIFHjx5JGo6IiDRP9B5DamoqrKysSozn5OQoFQUREVVNoouhffv2OHbsmOJ1cRls3rwZXbp0kS4ZERFphehDSV999RV8fHxw584dFBYWYt26dbhz5w4uXLhQ4roGIiKqekTvMXTv3h0xMTEoLCyEq6srQkNDYWVlhYsXL6Jdu3bqyEhERBokao+hoKAAkydPxsKFC7Fp0yZ1ZSIiIi0Stcegr6+P/fv3qysLERFVAqIPJQ0aNAiHDh1SQxQiIqoMRE8+Ozs7Y+nSpTh//jzatWsHU1NTpeW8iR4RUdUmuhi2bNmCmjVr4urVq7h69arSMt5dlYio6hNVDIIgICIiAlZWVjA2NlZXJiIi0iJRcwyCIMDZ2Rl//fWXuvIQEZGWiSoGHR0dODs7Iy0tTV15iIhIy0SflbRq1SrMmTMHt27dUkceIiLSMtGTz2PHjkVubi7c3NxgYGBQYq4hPT1dsnBERKR5ooshODhYDTGIiKiyEF0Mvr6+6shBRESVhOhiKH6U59vUr1+/3GGIiEj7RBdDgwYNynwgj1wur1AgIiLSLtHFcP36daXXBQUFuH79OtauXYsVK1ZIFoyIiLRDdDG4ubmVGGvfvj3s7OywevVqDBkyRJJgRESkHaKvY3ibpk2b4sqVK1JtjoiItET0HkNWVpbSa0EQkJiYiMWLF8PZ2VmyYEREpB2ii6FmzZolJp8FQYCDgwP++9//ShaMiIi0Q3QxnDlzRqkYdHR0ULduXTg5OUFPT/TmNGLw4MGIiIhA7969sW/fPm3HISKq1ET/Te7h4aGGGOo1ffp0TJgwATt27NB2FCKiSk/05PPKlSuxdevWEuNbt27F119/LUkoqXl4eMDc3FzbMYiIqgTRxfDDDz+gWbNmJcZbtGiBkJAQSUL93blz59C/f3/Y2dlBJpOV+rzpDRs2oEGDBjAyMkKnTp1w+fJlyXMQEVUXooshKSkJtra2Jcbr1q2LxMRESUL9XU5ODtzc3LBhw4ZSl+/ZswcBAQEIDAzEtWvX4ObmBm9vb6SkpEiehYioOhBdDA4ODjh//nyJ8fPnz8POzk6SUH/n4+OD5cuXY/DgwaUuX7t2LSZNmoTx48fDxcUFISEhMDExKfVwFxERvZvoyedJkyZhxowZKCgoQK9evQAAYWFh+Ne//oVZs2ZJHrAs+fn5uHr1KubPn68Y09HRgZeXFy5evCh6e69fv8br168Vr//3mg0ioupAdDHMmTMHaWlp+Pzzz5Gfnw8AMDIywty5czFv3jzJA5bl+fPnkMvlsLa2Vhq3trbG3bt3Fa+9vLwQExODnJwc1KtXD7/88gu6dOlSYnsrV67EkiVL1J6biKgyE10MMpkMX3/9NRYuXIjY2FgYGxvD2dkZhoaG6sgnidOnT6u03vz58xEQEKB4nZWVBQcHB3XFIiKqlEQXQ2ZmJuRyOWrXro0OHTooxtPT06Gnp4caNWpIGrAslpaW0NXVRXJystJ4cnIybGxsRG/P0NCwUhccEZEmiJ58HjlyZKm3vti7dy9GjhwpSShVGRgYoF27dggLC1OMFRUVISwsrNRDRURE9G6iiyEyMhKenp4lxj08PBAZGSlJqL/Lzs5GdHQ0oqOjAQBxcXGIjo5WPEkuICAAmzZtwo4dOxAbG4spU6YgJycH48ePlzwLEVF1IPpQ0uvXr1FYWFhivKCgAHl5eZKE+ruoqCilIiqeA/D19cX27dsxYsQIpKamYtGiRUhKSkLr1q1x4sSJEhPSRESkGtHF0LFjR/z4449Yv3690nhISAjatWsnWbBiHh4eEAShzHX8/f3h7+8v+XsTEVVHooth+fLlitM/e/fuDeDNdQxXrlxBaGio5AGJiEizRM8xdOvWDRcvXkS9evWwd+9eHD16FE5OTrhx4wZ69OihjoxERKRB5XqAQuvWrbF7926psxARUSUguhgSEhKwf/9+3L9/H8CbZz0PHTpULfdJIiIizRNVDP/+978REBCA/Px8xYVsWVlZmDNnDtauXYvPP/9cLSGJiEhzVJ5jOHbsGKZNmwZ/f38kJCQgIyMDGRkZSEhIwOeff47p06fj+PHj6sxKREQaoPIew+rVqzFv3jwsX75cadzW1hZr166FiYkJgoKC0K9fP8lDEhGR5qi8x3Dt2jWMGTPmrcvHjBmDa9euSRKKiIi0R+VikMvl0NfXf+tyfX19yOVySUIREZH2qFwMLVq0wOHDh9+6/NChQ2jRooUkoYiISHtUnmPw8/PDlClTYGhoiH/+85/Q03vzrYWFhfjhhx+wYMEC/Pvf/1ZbUCIi0gyVi8HX1xc3b96Ev78/5s+fj8aNG0MQBDx+/BjZ2dmYNm0axo0bp8aoRESkCaKuY/jmm28wbNgw/Pzzz3jw4AEAwN3dHSNHjkTnzp3VEpCIiDRL9JXPnTt3ZgkQEb3HRN9Ej4iI3m8sBiIiUsJiICIiJSwGIiJSwmIgIiIlKp2V1KZNG8hkMpU2yPslERFVbSoVw6BBg9Qcg4iIKguViiEwMFDdOYiIqJLgHAMRESkRfeWzXC7Ht99+i7179yI+Ph75+flKy9PT0yULR0REmid6j2HJkiVYu3YtRowYgczMTAQEBGDIkCHQ0dHB4sWL1RCRiIg0SXQx7Nq1C5s2bcKsWbOgp6eHUaNGYfPmzVi0aBEuXbqkjoxERKRBooshKSkJrq6uAAAzMzNkZmYCAD766CMcO3ZM2nRERKRxoouhXr16SExMBAA0btwYoaGhAIArV67A0NBQ2nRERKRxooth8ODBCAsLAwBMnToVCxcuhLOzM8aOHYsJEyZIHpCIiDRL9FlJq1atUvz3iBEj4OjoiAsXLsDZ2Rn9+/eXNBwREWme6GI4d+4cunbtqnjmc/GDewoLC3Hu3Dn07NlT8pBERKQ5og8leXp6lnqtQmZmJjw9PSUJRURE2iO6GARBKPWGemlpaTA1NZUkFBERaY/Kh5KGDBkCAJDJZBg3bpzSGUhyuRw3btxA165dpU9IREQapXIxWFhYAHizx2Bubg5jY2PFMgMDA3Tu3BmTJk2SPiEREWmUysWwbds2AECDBg0we/ZsHjYiInpPiT4rqfgW3Kmpqbh37x4AoGnTpqhbt660yYiISCtETz7n5uZiwoQJsLW1Rc+ePdGzZ0/Y2dlh4sSJyM3NVUdGIiLSINHFMHPmTJw9exZHjx5FRkYGMjIycPjwYZw9exazZs1SR0YitZPJZFr9IqpMRB9K2r9/P/bt2wcPDw/FWL9+/WBsbIyPP/4YGzdulDIfUbWwZMkSrb4/n9JIf1euQ0nW1tYlxq2srHgoiYjoPSC6GLp06YLAwEC8evVKMZaXl4clS5agS5cukoYjIiLNU/lQkq6uLhITExEcHIwPPvgA9erVg5ubGwAgJiYGRkZGOHnypNqCEhGRZqhcDIIgAABcXV3x4MED7Nq1C3fv3gUAjBo1Cp988onSRW9ERFQ1iZ58BgATExNe5UxE9J4SVQybN2+GmZlZmetMmzatQoGIiEi7RBVDSEgIdHV137pcJpOxGIiIqjhRxRAVFQUrKyt1ZSEiokpA5dNVeXUmEVH1oHIxFJ+VRERE7zeViyEwMPCdE89ERFT1qTzHwHupEBFVD6JviUFERO83FgMRESlRqRiOHDmCgoICdWchIqJKQKViGDx4MDIyMgC8uZleSkqKOjNJ6tdff0XTpk3h7OyMzZs3azsOEVGlp1Ix1K1bF5cuXQLw5rTVqnJNQ2FhIQICAnDmzBlcv34dq1evRlpamrZjERFVaioVw2effYaBAwdCV1cXMpkMNjY20NXVLfWrMrl8+TJatGgBe3t7mJmZwcfHB6GhodqORURUqal0uurixYsxcuRIPHz4EAMGDMC2bdtQs2ZNNUcDzp07h9WrV+Pq1atITEzEwYMHMWjQIKV1NmzYgNWrVyMpKQlubm5Yv349OnbsCAB49uwZ7O3tFeva29sjISFB7bmJiKoyla9jaNasGZo1a4bAwEAMHz4cJiYm6swFAMjJyYGbmxsmTJiAIUOGlFi+Z88eBAQEICQkBJ06dUJwcDC8vb1x79493tOJiKicRD+PofhCt9TUVNy7dw8A0LRpU9StW1faZAB8fHzg4+Pz1uVr167FpEmTMH78eABv7v567NgxbN26FfPmzYOdnZ3SHkJCQoJib4KIiEon+jqG3NxcTJgwAXZ2dujZsyd69uwJOzs7TJw4Ebm5uerIWKr8/HxcvXoVXl5eijEdHR14eXnh4sWLAICOHTvi1q1bSEhIQHZ2Nn777Td4e3u/dZuvX79GVlaW0hcRUXUjuhhmzpyJs2fP4siRI8jIyEBGRgYOHz6Ms2fPYtasWerIWKrnz59DLpfD2tpaadza2hpJSUkAAD09PaxZswaenp5o3bo1Zs2ahTp16rx1mytXroSFhYXiy8HBQa2fgYioMhJ9KGn//v3Yt28fPDw8FGP9+vWDsbExPv74Y2zcuFHKfBU2YMAADBgwQKV158+fj4CAAMXrrKwslgMRVTuiiyE3N7fEv9IBwMrKSqOHkiwtLaGrq4vk5GSl8eTkZNjY2JRrm4aGhjA0NJQiHhFRlSX6UFKXLl0QGBiIV69eKcby8vKwZMkSdOnSRdJwZTEwMEC7du0QFhamGCsqKkJYWJhGcxARvW9E7zGsW7cO3t7eqFevHtzc3AAAMTExMDIywsmTJyUNl52djYcPHypex8XFITo6GrVr10b9+vUREBAAX19ftG/fHh07dkRwcDBycnIUZykREZF4oouhZcuWePDgAXbt2oW7d+8CAEaNGoVPPvkExsbGkoaLioqCp6en4nXx8X9fX19s374dI0aMQGpqKhYtWoSkpCS0bt0aJ06cKPVQFxERqUZ0MQCAiYkJJk2aJHWWEjw8PN75SFF/f3/4+/urPQsRUXXB5zEQEZESFgMRESlhMRARkRIWAxERKSlXMWRkZGDz5s2YP38+0tPTAQDXrl3jLa2JiN4Dos9KunHjBry8vGBhYYEnT55g0qRJqF27Ng4cOID4+Hjs3LlTHTmJiEhDRO8xBAQEYNy4cXjw4AGMjIwU4/369cO5c+ckDUdERJonuhiuXLmCyZMnlxi3t7dX3NWUiIiqLtHFYGhoWOpzCu7fv6+Wh/UQEZFmiS6GAQMGYOnSpSgoKAAAyGQyxMfHY+7cuRg6dKjkAYmISLNEF8OaNWuQnZ0NKysr5OXlwd3dHU5OTjA3N8eKFSvUkZGIiDRI9FlJFhYWOHXqFP744w/cuHED2dnZaNu2rdIjNomIqOoq1030AKB79+7o3r27lFmIiKgSEF0M3333XanjMpkMRkZGcHJyQs+ePaGrq1vhcEREpHmii+Hbb79FamoqcnNzUatWLQDAixcvYGJiAjMzM6SkpKBRo0YIDw/n85KJiKog0ZPPX331FTp06IAHDx4gLS0NaWlpuH//Pjp16oR169YhPj4eNjY2mDlzpjryEhGRmoneY1iwYAH279+Pxo0bK8acnJzwzTffYOjQoXj8+DGCgoJ46ioRURUleo8hMTERhYWFJcYLCwsVVz7b2dnh5cuXFU9HREQaJ7oYPD09MXnyZFy/fl0xdv36dUyZMgW9evUCANy8eRMNGzaULiUREWmM6GLYsmULateujXbt2sHQ0BCGhoZo3749ateujS1btgAAzMzMsGbNGsnDEhGR+omeY7CxscGpU6dw9+5d3L9/HwDQtGlTNG3aVLGOp6endAmJiEijyn2BW7NmzdCsWTMpsxARUSVQrmL466+/cOTIEcTHxyM/P19p2dq1ayUJRkRE2iG6GMLCwjBgwAA0atQId+/eRcuWLfHkyRMIgoC2bduqIyMREWmQ6Mnn+fPnY/bs2bh58yaMjIywf/9+PH36FO7u7hg+fLg6MhIRkQaJLobY2FiMHTsWAKCnp4e8vDyYmZlh6dKl+PrrryUPSEREmiW6GExNTRXzCra2tnj06JFi2fPnz6VLRkREWiF6jqFz5874448/0Lx5c/Tr1w+zZs3CzZs3ceDAAXTu3FkdGYmISINEF8PatWuRnZ0NAFiyZAmys7OxZ88eODs784wkIqL3gOhiaNSokeK/TU1NERISImkgIiLSLtFzDI0aNUJaWlqJ8YyMDKXSICKiqkl0MTx58gRyubzE+OvXr5GQkCBJKCIi0h6VDyUdOXJE8d8nT56EhYWF4rVcLkdYWBgaNGggaTgiItI8lYth0KBBAN4829nX11dpmb6+Pho0aMA7qhIRvQdULoaioiIAQMOGDXHlyhVYWlqqLRQREWmP6LOS4uLi1JGDiIgqiXLdXTUsLAxhYWFISUlR7EkU27p1qyTBiIhIO0QXw5IlS7B06VK0b98etra2kMlk6shFRERaIroYQkJCsH37dowZM0YdeYiISMtEX8eQn5+Prl27qiMLERFVAqKL4dNPP8Xu3bvVkYWIiCoB0YeSXr16hR9//BGnT59Gq1atoK+vr7ScN9IjIqraRBfDjRs30Lp1awDArVu3lJZxIpqIqOoTXQzh4eHqyEFERJWE6DmGYg8fPsTJkyeRl5cHABAEQbJQRESkPaKLIS0tDb1790aTJk3Qr18/JCYmAgAmTpyIWbNmSR6QiIg0S3QxzJw5E/r6+oiPj4eJiYlifMSIEThx4oSk4YiISPNEzzGEhobi5MmTqFevntK4s7Mz/vzzT8mCERGRdojeY8jJyVHaUyiWnp4OQ0NDSUIREZH2iC6GHj16YOfOnYrXMpkMRUVFCAoKgqenp6ThiIhI80QfSgoKCkLv3r0RFRWF/Px8/Otf/8Lt27eRnp6O8+fPqyMjERFpkOg9hpYtW+L+/fvo3r07Bg4ciJycHAwZMgTXr19H48aN1ZGRiIg0qFzPY7CwsMCXX34pdRYiIqoERO8xbNu2Db/88kuJ8V9++QU7duyQJBQREWmP6GJYuXJlqc97trKywldffSVJKCIi0h7RxRAfH4+GDRuWGHd0dER8fLwkoYiISHtEF4OVlRVu3LhRYjwmJgZ16tSRJJTUBg8ejFq1amHYsGHajkJEVOmJLoZRo0Zh2rRpCA8Ph1wuh1wux5kzZzB9+nSMHDlSHRkrbPr06UrXXhAR0duJPitp2bJlePLkCXr37g09vTffXlRUhLFjx1baOQYPDw9ERERoOwYRUZUgao9BEAQkJSVh+/btuHfvHnbt2oUDBw7g0aNH2Lp1KwwMDEQHOHfuHPr37w87OzvIZDIcOnSoxDobNmxAgwYNYGRkhE6dOuHy5cui34eIiFQjao9BEAQ4OTnh9u3bcHZ2hrOzc4UD5OTkwM3NDRMmTMCQIUNKLN+zZw8CAgIQEhKCTp06ITg4GN7e3rh37x6srKwAAK1bt0ZhYWGJ7w0NDYWdnV2FMxIRVSeiikFHRwfOzs5IS0uTpBQAwMfHBz4+Pm9dvnbtWkyaNAnjx48HAISEhODYsWPYunUr5s2bBwCIjo6WJMvr16/x+vVrxeusrCxJtktEVJWInnxetWoV5syZU+J5z+qQn5+Pq1evwsvLSzGmo6MDLy8vXLx4UfL3W7lyJSwsLBRfDg4Okr8HEVFlJ3ryeezYscjNzYWbmxsMDAxgbGystDw9PV2ycM+fP4dcLoe1tbXSuLW1Ne7evavydry8vBATE4OcnBzUq1cPv/zyC7p06VJivfnz5yMgIEDxOisri+VARNWO6GIIDg5WQwz1On36tErrGRoa8pkSRFTtiS4GX19fdeQolaWlJXR1dZGcnKw0npycDBsbG43lICKqTkTPMQDAo0ePsGDBAowaNQopKSkAgN9++w23b9+WNJyBgQHatWuHsLAwxVhRURHCwsJKPRREREQVJ7oYzp49C1dXV0RGRuLAgQPIzs4G8OaWGIGBgaIDZGdnIzo6WnFmUVxcHKKjoxX3XQoICMCmTZuwY8cOxMbGYsqUKcjJyVGcpURERNISfShp3rx5WL58OQICAmBubq4Y79WrF77//nvRAaKiopQeCVo8+evr64vt27djxIgRSE1NxaJFi5CUlITWrVvjxIkTJSakiYhIGqKL4ebNm9i9e3eJcSsrKzx//lx0AA8PDwiCUOY6/v7+8Pf3F71tIiIST/ShpJo1ayIxMbHE+PXr12Fvby9JKCIi0h7RxTBy5EjMnTsXSUlJkMlkKCoqwvnz5zF79myMHTtWHRmJiEiDRBfDV199hWbNmsHBwQHZ2dlwcXFBz5490bVrVyxYsEAdGYmISINEzzEYGBhg06ZNWLRoEW7evIns7Gy0adNGsnsnERGRdqlcDEVFRVi9ejWOHDmC/Px89O7dG4GBgSVuiUFERFWbyoeSVqxYgS+++AJmZmawt7fHunXr4Ofnp85sRESkBSoXw86dO/Hvf/8bJ0+exKFDh3D06FHs2rULRUVF6sxHREQapnIxxMfHo1+/forXXl5ekMlkePbsmVqCERGRdqhcDIWFhTAyMlIa09fXR0FBgeShiIhIe1SefBYEAePGjVO6LfWrV6/w2WefwdTUVDF24MABaRMSEZFGqVwMpd1ue/To0ZKGISIi7VO5GLZt26bOHEREVEmU63kMRET0/mIxEBGREhYDEREpYTEQEZES0TfRIyLSNplMptX3f9fDxao67jEQEZESFgMRESlhMRARkRIWAxERKeHkMxGRSEuWLNHq+wcGBqp1+9xjICIiJSwGIiJSwmIgIiIlLAYiIlLCYiAiIiUsBiIiUsJiICIiJSwGIiJSwmIgIiIlLAYiIlLCW2KUofie61lZWVpOUj6vXr3S6vtX1Z+bNvB3VbVUxd9X8feo8iwJmfC+P3GiAv766y84ODhoOwYRkWSePn2KevXqlbkOi6EMRUVFePbsGczNzbX+xCixsrKy4ODggKdPn6JGjRrajkNl4O+qaqmqvy9BEPDy5UvY2dlBR6fsWQQeSiqDjo7OO5u1sqtRo0aV+p+3OuPvqmqpir8vCwsLldbj5DMRESlhMRARkRIWw3vK0NAQgYGBMDQ01HYUegf+rqqW6vD74uQzEREp4R4DEREpYTEQEZESFgMRESlhMRARkRIWQyW3cuVKdOjQAebm5rCyssKgQYNw7969d37fihUr0LVrV5iYmKBmzZqlrhMfH48PP/wQJiYmsLKywpw5c1BYWCjxJ3i/nTt3Dv3794ednR1kMhkOHTqktFwQBCxatAi2trYwNjaGl5cXHjx48M7tTps2De3atYOhoSFat25dYvnixYshk8lKfJmamkr0yag07/pztX379lJ/LzKZDCkpKZoPXE4shkru7Nmz8PPzw6VLl3Dq1CkUFBSgb9++yMnJKfP78vPzMXz4cEyZMqXU5XK5HB9++CHy8/Nx4cIF7NixA9u3b8eiRYvU8THeWzk5OXBzc8OGDRtKXR4UFITvvvsOISEhiIyMhKmpKby9vVW6CduECRMwYsSIUpfNnj0biYmJSl8uLi4YPnx4hT5PdZSbm6vyuu/6czVixIgSvxdvb2+4u7vDyspKqsjqJ1CVkpKSIgAQzp49q9L627ZtEywsLEqMHz9+XNDR0RGSkpIUYxs3bhRq1KghvH79Wqq41QoA4eDBg4rXRUVFgo2NjbB69WrFWEZGhmBoaCj8/PPPKm0zMDBQcHNze+d60dHRAgDh3LlzYmNXe3PnzhWcnJyEadOmCSdOnBBevXr1zu9525+r/5WSkiLo6+sLO3fulCCp5nCPoYrJzMwEANSuXbtC27l48SJcXV1hbW2tGPP29kZWVhZu375doW3TG3FxcUhKSoKXl5dizMLCAp06dcLFixclfa/NmzejSZMm6NGjh6TbrQ7mzp2LpUuXIi0tDZ988glq166NAQMGICQkBPHx8RXa9s6dO2FiYoJhw4ZJlFYzWAxVSFFREWbMmIFu3bqhZcuWFdpWUlKSUikAULxOSkqq0LbpjeKfY2k/Zyl/xq9evcKuXbswceJEybZZndSqVQujRo3Cf/7zH6SkpODUqVNwdXXFxo0b4ejoiJYtW+LZs2fl2vaWLVvwj3/8A8bGxhKnVi8WQxXi5+eHW7du4b///a9i7LPPPoOZmZnii6oWHx8fxe+uRYsW5drGwYMH8fLlS/j6+kqcrvp5+fIlnj17hsTERKSmpsLY2BiOjo7Q19cXva2LFy8iNja2ShY2b7tdRfj7++PXX3/FuXPnlG4FvnTpUsyePVv09mxsbHD58mWlseTkZMUyqrjin2NycjJsbW0V48nJyYozjTZv3oy8vDwAKNdfPsXb+Oijj0rsmZBq4uLisHfvXhw/fhwXLlyAo6MjfHx8sGXLFnh6esLIyKhc2928eTNat26Ndu3aSZxY/VgMlZwgCJg6dSoOHjyIiIgINGzYUGm5lZVVuc526NKlC1asWIGUlBTF9586dQo1atSAi4uLJNmru4YNG8LGxgZhYWGKIsjKykJkZKTirBZ7e/sKvUdcXBzCw8Nx5MiRisattrZv347z589j4MCB+OGHH9CsWbMKbzM7Oxt79+7FypUrJUioeSyGSs7Pzw+7d+/G4cOHYW5urjg2bWFhUeZxy/j4eKSnpyM+Ph5yuRzR0dEAACcnJ5iZmaFv375wcXHBmDFjEBQUhKSkJCxYsAB+fn7v9V0jpZadnY2HDx8qXsfFxSE6Ohq1a9dG/fr1MWPGDCxfvhzOzs5o2LAhFi5cCDs7OwwaNKjM7T58+BDZ2dlISkpCXl6e4vfn4uICAwMDxXpbt26Fra0tfHx81PHxqgV/f3+MGjVK8fru3bsl1mnUqBEMDAze+eeq2J49e1BYWIjRo0erPb9aaPu0KCobgFK/tm3bVub3+fr6lvp94eHhinWePHki+Pj4CMbGxoKlpaUwa9YsoaCgQL0f6D0THh5e6s/Z19dXEIQ3p6wuXLhQsLa2FgwNDYXevXsL9+7de+d23d3dS91uXFycYh25XC7Uq1dP+OKLL9T06aqHuXPnvvXPWfFXbGysIAiq/bkSBEHo0qWL8I9//EMLn0YavO02EREp4VlJRESkhMVARERKWAxERKSExUBEREpYDEREpITFQERESlgMRESkhMVARERKWAxERKSExUBEGjV48GDUqlWryj28pjphMRCRRk2fPh07d+7UdgwqA4uBqhQPDw/MmDGj0mxHW1TJn5aWBisrKzx58kQjmVTl4eEBc3PzEuMjR47EmjVrtJCI/heLgbQiJCQE5ubmKCwsVIxlZ2dDX18fHh4eSutGRERAJpPh0aNHOHDgAJYtW6b2fKmpqZgyZQrq168PQ0ND2NjYwNvbG+fPn1f7e0tlxYoVGDhwIBo0aKDtKCpZsGABVqxYoXiuOWkPn8dAWuHp6Yns7GxERUWhc+fOAIDff/8dNjY2iIyMxKtXrxRPzgoPD0f9+vXRuHFjjeUbOnQo8vPzsWPHDjRq1AjJyckICwtDWlqaxjJURG5uLrZs2YKTJ09q/L1bt26tVPjFQkNDYWdn99bva9myJRo3boz//Oc/8PPzU2dEegfuMZBWNG3aFLa2toiIiFCMRUREYODAgWjYsCEuXbqkNO7p6Qmg5CEUDw8PTJs2Df/6179Qu3Zt2NjYYPHixUrvlZOTg7Fjx8LMzAy2trbvPFyRkZGB33//HV9//TU8PT3h6OiIjh07Yv78+RgwYIDSe/v7+8Pf3x8WFhawtLTEwoUL8fc72RcVFWHlypVo2LAhjI2N4ebmhn379qm8vDz5AeD48eMwNDRUlG5x3qlTp2LGjBmoVasWrK2tsWnTJuTk5GD8+PEwNzeHk5MTfvvttzK3vW/fPri6usLY2Bh16tSBl5cXcnJyFMujo6Nx69atEl9llUKx/v37Kz3TnLSDxUBa4+npifDwcMXr8PBweHh4wN3dXTGel5eHyMhIRTGUZseOHTA1NUVkZCSCgoKwdOlSnDp1SrF8zpw5OHv2LA4fPozQ0FBERETg2rVrb92emZkZzMzMcOjQIbx+/brMz7Bjxw7o6enh8uXLWLduHdauXYvNmzcrlq9cuRI7d+5ESEgIbt++jZkzZ2L06NE4e/asSsvLkx94s/dV2rOGd+zYAUtLS1y+fBlTp07FlClTMHz4cHTt2hXXrl1D3759MWbMGOTm5pa63cTERIwaNQoTJkxAbGwsIiIiMGTIEEj1WJeOHTvi8uXL7/y5k5pp9zlBVJ1t2rRJMDU1FQoKCoSsrCxBT09PSElJEXbv3i307NlTEARBCAsLEwAIf/75pyAIb55sNn36dMU23N3dhe7duyttt0OHDsLcuXMFQRCEly9fCgYGBsLevXsVy9PS0gRjY2Ol7fyvffv2CbVq1RKMjIyErl27CvPnzxdiYmKU1nF3dxeaN28uFBUVKcbmzp0rNG/eXBAEQXj16pVgYmIiXLhwQen7Jk6cKIwaNeqdyyuSf+DAgcKECRNK5P37z6qwsFAwNTUVxowZoxhLTEwUAAgXL14sdbtXr14VAAhPnjx563u/S+/evQVLS0vB2NhYsLe3V/r8MTExFd4+VRznGEhrPDw8kJOTgytXruDFixdo0qQJ6tatC3d3d4wfPx6vXr1CREQEGjVqhPr16791O61atVJ6bWtri5SUFADAo0ePkJ+fj06dOimW165dG02bNi0z29ChQ/Hhhx/i999/x6VLl/Dbb78hKCgImzdvxrhx4xTrde7cGTKZTPG6S5cuWLNmDeRyOR4+fIjc3Fz06dNHadv5+flo06bNO5dXJH9eXp5ijubv/v6z0tXVRZ06deDq6qoYs7a2BgDFz+9/ubm5oXfv3nB1dYW3tzf69u2LYcOGoVatWmXm+bvTp0+/dVnxc8zftsdCmsFiIK1xcnJCvXr1EB4ejhcvXsDd3R0AYGdnBwcHB1y4cAHh4eHo1atXmdvR19dXei2TyVBUVFThfEZGRujTpw/69OmDhQsX4tNPP0VgYKBSMZQlOzsbAHDs2DHY29srLTM0NMSzZ8/KXF4RlpaWePHiRYnx0n5Wfx8rLrm3/fx0dXVx6tQpXLhwAaGhoVi/fj2+/PJLREZGomHDhhXKDADp6ekAgLp161Z4W1R+nGMgrfL09ERERAQiIiKUTlPt2bMnfvvtN1y+fLnM+YV3ady4MfT19REZGakYe/HiBe7fvy96Wy4uLkqTrACUtgsAly5dgrOzM3R1deHi4gJDQ0PEx8fDyclJ6cvBweGdyyuSv02bNrhz547oz6gKmUyGbt26YcmSJbh+/ToMDAxw8OBBSbZ969Yt1KtXD5aWlpJsj8qHewykVZ6envDz80NBQYFijwEA3N3d4e/vj/z8/AoVg5mZGSZOnIg5c+agTp06sLKywpdffgkdnbf/mygtLQ3Dhw/HhAkT0KpVK5ibmyMqKgpBQUEYOHCg0rrx8fEICAjA5MmTce3aNaxfv15x1pC5uTlmz56NmTNnoqioCN27d0dmZibOnz+PGjVqwNfX953Ly5MfALy9vTF//ny8ePFC1GGed4mMjERYWBj69u0LKysrREZGIjU1Fc2bN5dk+7///jv69u0rybao/FgMpFWenp7Iy8tDs2bNFMe3gTfF8PLlS8VprRWxevVqZGdno3///jA3N8esWbPKvIjKzMwMnTp1wrfffotHjx6hoKAADg4OmDRpEr744guldceOHYu8vDx07NgRurq6mD59Ov75z38qli9btgx169bFypUr8fjxY9SsWRNt27ZVbOddy8uTHwBcXV3Rtm1b7N27F5MnTy7Pj61UNWrUwLlz5xAcHIysrCw4OjpizZo18PHxqfC2X716hUOHDuHEiRMSJKWKkAmCROeZEVUzHh4eaN26NYKDg7UdpVTHjh3DnDlzcOvWrXfuYVQGGzduxMGDBxEaGqrtKNUe9xiI3lMffvghHjx4gISEBMWcRWWmr6+P9evXazsGgcVA9F6rSjcK/PTTT7Udgf4PDyUREZGSyn/gkYiINIrFQERESlgMRESkhMVARERKWAxERKSExUBEREpYDEREpITFQERESlgMRESkhMVARERKWAxERKSExUBEREpYDEREpOT/AeV9TSp/5OSvAAAAAElFTkSuQmCC", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "############## Wind PDF ################\n", "# Bins\n", "WIND_BINS = [(2, 10), (10, 17), (17, np.inf)]\n", "BIN_LABELS = ['2-10', '10-17', '≥17']\n", "\n", "# Function to calculate frequency (number of occurrences) per intensity bin\n", "def calculate_percentage_per_bin(precip_data, bins):\n", " total_occurrences = precip_data.count().item() # Total number of grid points across all time steps\n", " bin_percentages = []\n", "\n", " for lower, upper in bins:\n", " bin_mask = (precip_data >= lower) & (precip_data < upper)\n", " bin_frequency = bin_mask.sum().item() # Count occurrences in this bin\n", " bin_percentage = (bin_frequency / total_occurrences) * 100 # Normalize to percentage\n", " bin_percentages.append(bin_percentage)\n", "\n", " return bin_percentages\n", "\n", "# Functions to calculate relative changes\n", "def calculate_relative_change(asos, current):\n", " return ((current / asos) - 1) * 100\n", "\n", "def plot_precip_intensity_frequency_bar(month=month):\n", " # Step 1: Use all occurrences of precipitation values for each simulation across all time steps and grid cells\n", " all_winds_current = df_winds['wrf_sped']\n", " asos_winds = df_winds['sped']\n", "\n", " # Step 2: Calculate frequency for each bin for each simulation\n", " freq_current = calculate_percentage_per_bin(all_winds_current, WIND_BINS)\n", " freq_mrms = calculate_percentage_per_bin(asos_winds, WIND_BINS)\n", "\n", " relative_change_future = calculate_relative_change(np.array(freq_mrms), np.array(freq_current))\n", "\n", " # Step 4: Plotting\n", " fig, ax1 = plt.subplots(figsize=(4, 4))\n", " #ax2 = ax1.twinx()\n", " #ax2.grid(True, axis='y', alpha=0.5, zorder=0)\n", "\n", " # Bar width and positions\n", " bar_width = 0.3\n", " x = np.arange(len(WIND_BINS))\n", "\n", " # Plot bars for each simulation\n", " ax1.bar(x+bar_width/2, freq_mrms, width=bar_width, color='gray', label='ASOS',zorder=2)\n", " ax1.bar(x - bar_width/2, freq_current, width=bar_width, color='black', label='Current',zorder=2)\n", "\n", " # Primary Y-axis (left)\n", " ax1.set_xlabel('Wind Speed (m s$^{-1}$)')\n", " ax1.set_ylabel('Percentage of Total Occurrences (%)')\n", " ax1.set_yscale('log')\n", " ax1.set_xticks(x)\n", " ax1.set_xticklabels(BIN_LABELS)\n", " ax1.minorticks_off()\n", " ax1.set_ylim(0, 500) # Since we're dealing with percentages\n", " #ax1.set_ylim(10e1, 10e8)\n", "\n", " # Add legend for secondary y-axis lines\n", " ax1.legend(loc='upper right')\n", " '''\n", " ax2.plot(x, relative_change_future, color='blue', marker='o', linestyle='--', label='% Change', linewidth=2)\n", " ax2.set_ylabel('Relative change (%)', color='black')\n", " ax2.set_ylim(-20, 40) # Adjust based on expected range of relative changes\n", " #ax2.set_yscale('symlog', linthresh=10)\n", " #custom_ticks = [-100,-10,0, 10,100,1000]\n", " #ax2.minorticks_off()\n", " #ax2.set_yticks(custom_ticks)\n", " ax2.legend(loc='upper left')\n", " '''\n", " if month == '04':\n", " month = 'April'\n", " elif month == '05':\n", " month = 'May'\n", " elif month == '06':\n", " month = 'June'\n", "\n", " # Title and layout adjustments\n", " ax1.set_title(f'{month}')\n", " fig.tight_layout()\n", "\n", " plt.show()\n", "\n", "plot_precip_intensity_frequency_bar()" ] } ], "metadata": { "kernelspec": { "display_name": "myenv", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.11.6" } }, "nbformat": 4, "nbformat_minor": 2 }