{ "cells": [ { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", "## You are using the Python ARM Radar Toolkit (Py-ART), an open source\n", "## library for working with weather radar data. Py-ART is partly\n", "## supported by the U.S. Department of Energy as part of the Atmospheric\n", "## Radiation Measurement (ARM) Climate Research Facility, an Office of\n", "## Science user facility.\n", "##\n", "## If you use this software to prepare a publication, please cite:\n", "##\n", "## JJ Helmus and SM Collis, JORS 2016, doi: 10.5334/jors.119\n", "\n" ] } ], "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 ccrs\n", "import cartopy.feature as cfeature\n", "from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter\n", "import matplotlib.ticker as mticker\n", "import matplotlib\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 matplotlib.lines as mlines\n", "import matplotlib.patheffects as pe\n", "\n", "\n", "import wrf\n", "from wrf import (getvar, interplevel, to_np, latlon_coords, get_cartopy,\n", " cartopy_xlim, cartopy_ylim)\n", "\n", "import pyart" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "month='06'" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "############## Read in NEXRAD ################\n", "\n", "directory= f'/pscratch/sd/d/dbrooks/acc2017_analysis/model_evaluation/nexrad_data/{month}'\n", "nexrad_files = sorted(glob.glob(os.path.join(directory, f\"nexrad_3d_data_2017*.nc\")))\n", "\n", "datasets=[]\n", "for file in nexrad_files:\n", " ds = xr.open_dataset(file)\n", " ds = ds['Reflectivity'].max(dim='Altitude')\n", " datasets.append(ds)\n", "\n", "nexrad_ds = xr.concat(datasets, dim='time')\n" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "########## Read in WRF Reflectivity ############\n", "\n", "def decode_times(times_char_array):\n", " \"\"\"\n", " Convert the Times character array into datetime64 objects, handling the\n", " WRF datetime format with underscores.\n", " \"\"\"\n", " # Decode the character array into strings\n", " times_str = [''.join(t.astype(str)) for t in times_char_array]\n", " \n", " # Replace underscore with space to make it compatible with datetime format\n", " times_str = [t.replace('_', ' ') for t in times_str]\n", " \n", " # Convert to numpy datetime64 array\n", " return np.array(times_str, dtype=\"datetime64[ns]\")\n", "\n", "# Use to read in reflectivity data and grid it to the proper lat lon grid\n", "def read_in_monthly_data(month, hour_interval, climate_state, var_name):\n", " \"\"\"\n", " Reads in all NetCDF files for a given month and combines them into one dataset.\n", "\n", " Parameters:\n", " month (str): month you want data for, e.g., '04'\n", " hour_interval (str): '1hr' or '3hr'\n", " climate_state (str): 'current', 'future', or 'future_urban'\n", " var_name (str): Variable name you want to load in, e.g., \"RAINNC\" for accumulated precipitation\n", "\n", " Returns:\n", " xarray.Dataset of full month of data\n", " \"\"\"\n", "\n", " # Determine the file path based on the input parameters\n", " if hour_interval == '3hr':\n", " file_path = f'/pscratch/sd/y/yuwei/Climate_Impact/long-term/data/{climate_state}/{hour_interval}/wrfout_d01_2017-{month}*'\n", " elif hour_interval == '1hr':\n", " file_path = f'/pscratch/sd/y/yuwei/Climate_Impact/long-term/data/{climate_state}/{hour_interval}/wrfout_hourly_d01_2017-{month}*'\n", "\n", " # also need to grab xlat and xlon from the 3hr files\n", " file_path2 = f'/pscratch/sd/y/yuwei/Climate_Impact/long-term/data/{climate_state}/3hr/wrfout_d01_2017-04-01_00:00:00'\n", " ncfile2 = netCDF4.Dataset(file_path2, 'r') \n", " # Extract latitude and longitude\n", " lats = ncfile2.variables['XLAT'][:] # Latitude\n", " lons = ncfile2.variables['XLONG'][:] # Longitude\n", "\n", " # Use glob to find all matching files for the month\n", " file_list = sorted(glob.glob(file_path))\n", "\n", " array_list = []\n", " for file in file_list:\n", " print(file)\n", " # Read file\n", " ncfile = netCDF4.Dataset(file, 'r') \n", " # Check if REFL_10CM is the variable we want\n", " if var_name == 'REFL_10CM':\n", " # Extract REFL_10CM variable directly\n", " ref_data = ncfile.variables['REFL_10CM'][:]\n", " times_char_array = ncfile.variables['Times'][:]\n", "\n", " ref_data = ref_data.max(axis=1) # Get max dbz for vertical column\n", "\n", " # Decode times (convert to string or datetime64)\n", " times = decode_times(times_char_array)\n", " \n", " # Convert to xarray DataArray with time, lat, lon as dimensions\n", " ref_da = xr.DataArray(\n", " ref_data, \n", " dims=[\"Time\", \"south_north\", \"west_east\"], \n", " coords={\"Time\": times, \"XLAT\": ([\"south_north\", \"west_east\"], lats[0]), \"XLONG\": ([\"south_north\", \"west_east\"], lons[0])},\n", " name=\"REFL_10CM\",\n", " attrs={\"Description\": \"Composite ref for vertical columns (dbz)\"}\n", " )\n", " \n", " \n", " # Convert to xarray Dataset\n", " data = ref_da.to_dataset(name=\"REFL_10CM\")\n", " #print(data)\n", " else:\n", " # For other variables, use wrf-python getvar\n", " data = getvar(ncfile, var_name).to_dataset(name=var_name)\n", " \n", " array_list.append(data)\n", " ncfile.close()\n", "\n", " print('done')\n", " \n", " # Combine all datasets along the Time dimension\n", " combined_ds = xr.concat(array_list, dim='Time')\n", " \n", " return combined_ds\n", "\n", "#wrf_dbz = read_in_monthly_data(month, hour_interval='1hr', climate_state='current', var_name='REFL_10CM')\n", "\n", "filename = f'/pscratch/sd/d/dbrooks/acc2017_analysis/radar_data/current/composite_ref_month{month}.nc'\n", "wrf_dbz = xr.open_dataset(filename)" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "#wrf_dbz.to_netcdf(f'/pscratch/sd/d/dbrooks/acc2017_analysis/radar_data/current/composite_ref_month{month}.nc')" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " Size: 4GB\n", "Dimensions: (Time: 720, south_north: 1080, west_east: 1210)\n", "Coordinates:\n", " * Time (Time) datetime64[ns] 6kB 2017-06-01 ... 2017-06-30T23:00:00\n", " XLAT (south_north, west_east) float32 5MB ...\n", " XLONG (south_north, west_east) float32 5MB ...\n", "Dimensions without coordinates: south_north, west_east\n", "Data variables:\n", " REFL_10CM (Time, south_north, west_east) float32 4GB ...\n" ] } ], "source": [ "print(wrf_dbz)" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " Size: 4GB\n", "array([[[ nan, nan, nan, ..., nan, nan,\n", " nan],\n", " [ nan, nan, nan, ..., nan, nan,\n", " nan],\n", " [ nan, nan, nan, ..., nan, nan,\n", " nan],\n", " ...,\n", " [ nan, nan, 12.892142, ..., nan, nan,\n", " nan],\n", " [18.415363, 14.983327, nan, ..., nan, nan,\n", " nan],\n", " [21.09337 , 16.295795, nan, ..., nan, nan,\n", " nan]],\n", "\n", " [[ nan, nan, nan, ..., nan, nan,\n", " nan],\n", " [ nan, nan, nan, ..., nan, nan,\n", " nan],\n", " [ nan, nan, nan, ..., nan, nan,\n", " nan],\n", "...\n", " [ nan, nan, nan, ..., nan, nan,\n", " nan],\n", " [ nan, nan, nan, ..., nan, nan,\n", " nan],\n", " [ nan, nan, nan, ..., nan, nan,\n", " nan]],\n", "\n", " [[ nan, nan, nan, ..., nan, nan,\n", " nan],\n", " [ nan, nan, nan, ..., nan, nan,\n", " nan],\n", " [ nan, nan, nan, ..., nan, nan,\n", " nan],\n", " ...,\n", " [ nan, nan, nan, ..., nan, nan,\n", " nan],\n", " [ nan, nan, nan, ..., nan, nan,\n", " nan],\n", " [ nan, nan, nan, ..., nan, nan,\n", " nan]]], dtype=float32)\n", "Coordinates:\n", " * Latitude (Latitude) float32 4kB 26.34 26.36 26.38 ... 46.36 46.38 46.4\n", " * Longitude (Longitude) float32 6kB -113.0 -113.0 -112.9 ... -82.04 -82.02\n", "Dimensions without coordinates: time\n" ] } ], "source": [ "print(nexrad_ds)" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import xarray as xr\n", "from shapely.geometry import Point, MultiPoint\n", "from shapely.prepared import prep\n", "\n", "#----------- Mask NEXRAD data to keep data only within WRF domain ----------------\n", "# Step 1: Boolean mask for ds1 footprint (first time slice)\n", "mask1_bool = ~np.isnan(wrf_dbz['REFL_10CM'].isel(Time=0))\n", "\n", "# Step 2: Extract lat/lon arrays from ds1\n", "lat1 = wrf_dbz.XLAT.data[mask1_bool.data]\n", "lon1 = wrf_dbz.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: Build mask for ds2\n", "lat2, lon2 = np.meshgrid(nexrad_ds.Latitude.data, nexrad_ds.Longitude.data, indexing='ij')\n", "mask2 = np.vectorize(lambda lo, la: prep_poly.contains(Point(lo, la)))(lon2, lat2)\n", "\n", "# Step 5: Apply mask to ds2\n", "# mask2 has shape (len(latitude), len(longitude)), broadcast to ds2\n", "nexrad_ds = nexrad_ds.where(mask2)" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "'\\n# Step 1: Boolean mask for ds1 footprint (first time slice)\\nmask1_bool = ~np.isnan(nexrad_ds.isel(time=0))\\n\\n# Step 2: Extract lat/lon arrays from ds1\\nlat1 = nexrad_ds.Latitude.data\\nlon1 = nexrad_ds.Longitude.data\\n\\n# Step 3: Build polygon from ds1 valid points\\npoints = MultiPoint(list(zip(lon1, lat1)))\\npolygon = points.convex_hull # or cascaded_union for exact pixel shapes\\nprep_poly = prep(polygon) # speeds up \"contains\" checks\\n\\n# Step 4: Build mask for ds2\\nlat2, lon2 = np.meshgrid(wrf_dbz.XLAT.data, wrf_dbz.XLONG.data, indexing=\\'ij\\')\\nmask2 = np.vectorize(lambda lo, la: prep_poly.contains(Point(lo, la)))(lon2, lat2)\\n\\n# Step 5: Apply mask to ds2\\n# mask2 has shape (len(latitude), len(longitude)), broadcast to ds2\\nwrf_dbz = wrf_dbz.where(mask2)\\n'" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "'''\n", "# Step 1: Boolean mask for ds1 footprint (first time slice)\n", "mask1_bool = ~np.isnan(nexrad_ds.isel(time=0))\n", "\n", "# Step 2: Extract lat/lon arrays from ds1\n", "lat1 = nexrad_ds.Latitude.data\n", "lon1 = nexrad_ds.Longitude.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: Build mask for ds2\n", "lat2, lon2 = np.meshgrid(wrf_dbz.XLAT.data, wrf_dbz.XLONG.data, indexing='ij')\n", "mask2 = np.vectorize(lambda lo, la: prep_poly.contains(Point(lo, la)))(lon2, lat2)\n", "\n", "# Step 5: Apply mask to ds2\n", "# mask2 has shape (len(latitude), len(longitude)), broadcast to ds2\n", "wrf_dbz = wrf_dbz.where(mask2)\n", "'''" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [], "source": [ "import cmweather\n", "wind_cmap = cmweather.cm_colorblind.ChaseSpectral\n", "\n", "p_clevs=np.arange(10,61,5)\n", "p_cmap = mcolors.ListedColormap(wind_cmap(np.linspace(0.15,0.9,len(p_clevs))))\n", "p_norm = mcolors.BoundaryNorm(p_clevs, len(p_clevs))\n", "#norm = mcolors.BoundaryNorm(clevs, cmap.N)\n", "p_cmap.set_over('#fdd8fb')\n", "p_cmap.set_under('white')\n", "\n", "def plot_nexrad(ds):\n", " # Create the plot\n", " fig, ax_wind = plt.subplots(\n", " figsize=(8, 8),\n", " subplot_kw={\"projection\": ccrs.PlateCarree()}\n", " )\n", " #wind_speed_current = ds['wspd_wdir10'].sel(wspd_wdir='wspd') \n", " lat=ds['Latitude']\n", " lon=ds['Longitude']\n", " reflectivity=ds.max(dim='time')\n", " #print(max_wind_current)\n", " #time = datetime(time)\n", "\n", " #ax_wind.set_extent([-97.953608,-88.879741,28,32.545423])\n", "\n", " #max_precip_threshold = max_wind_current >= 25\n", " \n", " #ax_wind = axs[0]\n", " # Plot the reflectivity data\n", " pb = ax_wind.pcolormesh(\n", " lon,\n", " lat,\n", " reflectivity,\n", " cmap=p_cmap,\n", " norm=p_norm,\n", " transform=ccrs.PlateCarree()\n", " )\n", " \n", " # Add features and title\n", " ax_wind.add_feature(cfeature.STATES, edgecolor=\"black\", linewidths=0.85, alpha=1, zorder=1)\n", " #ax_wind.add_feature(cfeature.STATES, edgecolor=\"white\", linewidths=0.85, alpha=1, zorder=1)\n", " ax_wind.add_feature(cfeature.COASTLINE, edgecolor=\"black\", linewidths=0.85, alpha=1)\n", " #ax_wind.add_feature(USCOUNTIES, linewidth=0.35, alpha=0.3)\n", " gl = ax_wind.gridlines(draw_labels=True, linewidth=1, 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 = False\n", " ax_wind.set_title(f'Hourly Composite Reflectivity', loc='left', fontsize=12)\n", " ax_wind.set_title(f'(NEXRAD)', 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.5, extend='both', drawedges=True)\n", " cbar.set_label('dbz')\n", "\n", " plt.show()\n", "\n", "\n", "def plot_wrf(ds):\n", " # Create the plot\n", " fig, ax_wind = plt.subplots(\n", " figsize=(8, 8),\n", " subplot_kw={\"projection\": ccrs.PlateCarree()}\n", " )\n", " #wind_speed_current = ds['wspd_wdir10'].sel(wspd_wdir='wspd') \n", " lat=ds['XLAT']\n", " lon=ds['XLONG']\n", " reflectivity=ds['REFL_10CM'].max(dim='Time')\n", " #print(max_wind_current)\n", " #time = datetime(time)\n", "\n", " #ax_wind.set_extent([-97.953608,-88.879741,28,32.545423])\n", "\n", " #max_precip_threshold = max_wind_current >= 25\n", " \n", " #ax_wind = axs[0]\n", " # Plot the reflectivity data\n", " pb = ax_wind.pcolormesh(\n", " lon,\n", " lat,\n", " reflectivity,\n", " cmap=p_cmap,\n", " norm=p_norm,\n", " transform=ccrs.PlateCarree()\n", " )\n", " \n", " # Add features and title\n", " ax_wind.add_feature(cfeature.STATES, edgecolor=\"black\", linewidths=0.85, alpha=1, zorder=1)\n", " #ax_wind.add_feature(cfeature.STATES, edgecolor=\"white\", linewidths=0.85, alpha=1, zorder=1)\n", " ax_wind.add_feature(cfeature.COASTLINE, edgecolor=\"black\", linewidths=0.85, alpha=1)\n", " #ax_wind.add_feature(USCOUNTIES, linewidth=0.35, alpha=0.3)\n", " gl = ax_wind.gridlines(draw_labels=True, linewidth=1, 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 = False\n", " ax_wind.set_title(f'Hourly Composite Reflectivity (Month:{month})', loc='left', fontsize=12)\n", " ax_wind.set_title(f'(Current)', 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.5, extend='both', drawedges=True)\n", " cbar.set_label('dbz')\n", "\n", " plt.show()\n", "\n", "#plot_nexrad(nexrad_ds)\n", "#plot_wrf(wrf_dbz)" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYYAAAGGCAYAAAB/gCblAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAABIl0lEQVR4nO3deVxO6f8/8NfdvkdoI4USsmcLg0ZkGfs+ljBj37OMPoOGsTOGmTGMscQMM9YxjH0ixi60UPbsJaJSllTX7w+/ztdR0cm9pF7Px+N+PNzXOfd1vy+ZXnPOdc65VEIIASIiov9PT9cFEBFRwcJgICIiGQYDERHJMBiIiEiGwUBERDIMBiIikmEwEBGRDIOBiIhkGAxERCTDYCAiIhkGA5FCQUFBUKlUCA0N1XUpRBrBYCAiIhkGAxERyTAYiD5Qs2bN0KxZs2zt/fv3h4uLi/T+5s2bUKlUWLhwIVasWIEKFSrA2NgYdevWxZkzZ7J9/tKlS+jatStsbGxgYmKCOnXqYMeOHRocCdFrBrougKio2bBhA54+fYohQ4ZApVJh/vz56Ny5M27cuAFDQ0MAwMWLF9GoUSOULl0akydPhrm5OTZt2oSOHTti69at6NSpk45HQYUZg4FIy27fvo2rV6+iePHiAAB3d3d06NAB+/btw2effQYAGDNmDMqWLYszZ87A2NgYADB8+HA0btwYX331FYOBNIqnkoi0rEePHlIoAMAnn3wCALhx4wYA4PHjxzh48CC6d++Op0+f4tGjR3j06BESEhLg6+uLq1ev4t69ezqpnYoGHjEQaVnZsmVl77NC4smTJwCAa9euQQiBqVOnYurUqTn2ER8fj9KlS2u2UCqyGAxEH0ilUiGnFXIzMjJy3F9fXz/H9qw+MjMzAQATJkyAr69vjvu6urrmp1SiPGEwEH2g4sWLS6eB3nTr1q189Ve+fHkAgKGhIXx8fD6oNqL84BwD0QeqUKECLl26hIcPH0pt4eHhOHbsWL76s7W1RbNmzfDLL78gNjY22/Y3v4dIE3jEQPSBBg4ciEWLFsHX1xdffPEF4uPjsXz5cnh4eCA5OTlffS5duhSNGzdGtWrVMGjQIJQvXx4PHjzAiRMncPfuXYSHh6t5FET/h0cMRAplzQVkzRVUrlwZ69atQ1JSEvz9/bFjxw789ttvqF27dr6/o0qVKggNDUXbtm0RFBSEESNGYPny5dDT08O0adPUMg6i3KhETrNmRJSrH374AWPGjMG1a9dQoUIFXZdDpHY8YiBS6MyZMzA3N4ezs7OuSyHSCM4xEOXR1q1bERISgvXr1+PLL7+EgQH/86HCiaeSiPKoXLlyePr0KTp16oTFixfD3Nxc1yURaQSDgYiIZDjHQEREMgwGIiKSKfKzZ5mZmbh//z4sLS2hUql0XQ4RkUYIIfD06VM4OjpCT+/dxwRFPhju378PJycnXZdBRKQVd+7cQZkyZd65T5EPBktLSwCv/7KsrKx0XA0RkWYkJyfDyclJ+p33LkU+GLJOH1lZWTEYiKjQy8spc04+ExGRDIOBiIhkGAxERCRT5OcYiChvMjIy8OrVK12XQbkwNDTMddlYpRgMRPROQgjExcUhMTFR16XQexQrVgz29vYffE8Wg4GI3ikrFGxtbWFmZsYbQQsgIQSePXuG+Ph4AICDg8MH9cdgIKJcZWRkSKFQokQJXZdD72BqagoAiI+Ph62t7QedVuLkMxHlKmtOwczMTMeVUF5k/Zw+dC6IwUBE78XTRx8Hdf2cGAxERCTDYCAiIhlOPhNRvkyfPl1r3xUYGKj4M/3798fatWsxZ84cTJ48WWrfvn07OnXqBCEEQkJC4O3tnePnY2NjYW9vjx49eiAmJgYnTpyQJnRfvXqFBg0aoFKlSli/fj0A+WkcS0tLuLu7Y8qUKejQoUO2vufMmYMpU6Zg7ty5mDhxomxbUFAQBgwYAADQ09ODlZUVKlasiLZt22LMmDGwtrZW/HehFI8YiKjQMjExwbx58/DkyZN37nf58mXExsbKXra2tgCAn3/+Gbdv38bcuXOl/b/99lvExsbip59+kvWzZs0axMbGIjQ0FI0aNULXrl0RGRmZ7ftWr16NSZMmYfXq1TnWY2VlhdjYWNy9exfHjx/H4MGDsW7dOtSsWRP3799X+tegGIOBiAotHx8f2NvbY86cOe/cz9bWFvb29rJX1mI2JUqUwIoVKzBjxgxEREQgNDQUc+bMwcqVK1G8eHFZP1k3mFWsWBHffvst0tPTcejQIdk+hw8fxvPnzzFjxgwkJyfj+PHj2epRqVSwt7eHg4MDKleujC+++ALHjx9HSkoKJk2a9IF/K+/HYCCiQktfXx+zZ8/Gjz/+iLt37+a7n/bt26Nnz57o168f/Pz84OfnhzZt2uS6f3p6OlatWgUAMDIykm1btWoVevXqBUNDQ/Tq1Uva731sbW3Ru3dv7NixAxkZGfkeS14wGIioUOvUqRNq1qz5znmKMmXKwMLCQnp5eHhk22fx4sW4cuUKEhISsGjRohz76dWrFywsLGBsbIxx48bBxcUF3bt3l7YnJydjy5Yt6NOnDwCgT58+2LRpE1JSUvI0lkqVKuHp06dISEjI0/75xWAgokJv3rx5WLt2LaKjo3Pc/t9//yEsLEx67d69O9s+f/zxB1QqFR49eoRLly7l2M/333+PsLAw7NmzB1WqVMHKlSthY2Mj66NChQqoUaMGAKBmzZpwdnbGxo0b8zQOIQQAzd9XwmAgokKvSZMm8PX1RUBAQI7by5UrB1dXV+nl7Ows237jxg1MmjQJy5YtQ9++fdG/f3+8fPkyWz/29vZwdXVFy5YtsWbNGvTo0UN6fhHw+jTSxYsXYWBgIL2ioqJynYR+W3R0NKysrDT+eBIGAxEVCXPnzsXOnTtx4sQJRZ/LzMxE//790bx5c/Tr1w+LFy/G06dPMW3atHd+rl69evD09MSsWbMAAJGRkQgNDUVISIjs6CQkJAQnTpzI9SgkS3x8PDZs2ICOHTtKE+OawvsYiKhIqFatGnr37o0ffvgh27b4+Hi8ePFC1laiRAkYGhpiyZIluHjxIi5evAgAsLa2xsqVK/HZZ5+hS5cuqFevXq7fOXbsWHTq1AmTJk3CqlWrUK9ePTRp0iTbfnXr1sWqVauwYMECAP/3qHMhBBITE3HixAnMnj0b1tbWsstmNYVHDERUZMyYMQOZmZnZ2t3d3eHg4CB7nT17FleuXMHXX3+NH3/8Efb29tL+vr6+GDBgQK6nlLK0atUK5cqVw6xZs/D777+jS5cuOe7XpUsXrFu3Tnr4XXJyMhwcHFC6dGl4eXnhl19+gZ+fH86fP//Bj9TOC5XIms0oopKTk2FtbY2kpCRYWVnpuhyiAuXFixeIiYlBuXLlYGJiouty6D3e9fNS8ruORwxERCTDYCAiIhkGAxERyXxQMLxr0oWIiD5OioJhz5498PPzQ/ny5WFoaAgzMzNYWVmhadOmmDVrllae+kdERJqVp2D466+/ULFiRQwcOBAGBgb46quvsG3bNuzbtw8rV65E06ZN8e+//6J8+fIYOnQoHj58qOm6iYhIQ/J0g9v8+fPx/fffo3Xr1jnecZf1kKh79+7hxx9/xO+//45x48apt1IiItKKPAVDXm8hL126tFbuyiMiIs354KuSUlNTkZycrI5aiIioAMh3MERFRaFOnTqwtLRE8eLFUa1aNYSGhqqzNiIi0oF8B8OQIUMwcuRIpKSkICEhAZ07d4afn586ayOiAkylUmntlV9xcXEYNWoUypcvD2NjYzg5OaFdu3YIDg5W49+EeqlUKmzfvl2nNeQ5GDp06IB79+5J7x8+fIj27dvDzMwMxYoVQ5s2bfDgwQONFElEpNTNmzfh6emJgwcPYsGCBYiMjMTevXvh7e2NESNG5KtPIQTS09OztaelpX1ouQVKnoOhT58++PTTT/HDDz9ACIGRI0fCw8MDPXv2RJcuXdCqVSuMHTtWg6USEeXd8OHDoVKpcPr0aXTp0gUVK1aEh4cH/P39cfLkSdy8eRMqlQphYWHSZxITE6FSqRASEgIACAkJgUqlwp49e+Dp6QljY2McPXoUzZo1w8iRIzF27FiULFkSvr6+AIALFy6gdevWsLCwgJ2dHfr27YtHjx5J/Tdr1gyjR4/GpEmTYGNjA3t7e3zzzTfSdhcXFwCvlyNVqVTSe23LczB069YNp0+fRlRUFBo0aIBGjRph//79aNSoET755BPs378fU6ZM0WStRER58vjxY+zduxcjRoyAubl5tu3FihVT1N/kyZMxd+5cREdHo3r16gCAtWvXwsjICMeOHcPy5cuRmJiITz/9FLVq1UJoaCj27t2LBw8eyNZ8zvqcubk5Tp06hfnz52PGjBk4cOAAAODMmTMAgDVr1iA2NlZ6r22KFuqxtrbG8uXLcfToUfj5+aFFixb49ttvYWZmpqn6iIgUu3btGoQQqFSpklr6mzFjBlq0aCFrc3Nzw/z586X3M2fORK1atTB79mypbfXq1XBycsKVK1dQsWJFAED16tURGBgo9fHTTz8hODgYLVq0QKlSpQC8Dq4313/QNkWTz48fP8bZs2dRrVo1nD17FlZWVqhVq1aOC2cTEemKupeZqVOnTrY2T09P2fvw8HAcOnQIFhYW0isrmK5fvy7tl3XEkcXBwUG2LnRBkOcjhg0bNuDLL7+ElZUVXrx4gXXr1iEwMBA9evTA0KFDERQUhB9//BF2dnaarJeI6L3c3NygUqneuY5y1lMc3gyRrBXU3pbT6ai321JSUtCuXTvMmzcv275vrrpmaGgo26ZSqXJcVU6X8nzEEBAQgNWrVyMuLg7BwcGYOnUqAKBSpUoICQlBixYt4OXlpbFCc3Pnzh00a9YMVapUQfXq1bF582at10BEBYuNjQ18fX2xdOlSpKamZtuemJgonbaJjY2V2t+ciFaqdu3auHjxIlxcXODq6ip75RQsuTE0NERGRka+61CHPAdDSkoK3N3dAQAVKlTAs2fPZNsHDRqEkydPqre6PDAwMMDixYsRFRWF/fv3Y+zYsTn+QyCiomXp0qXIyMhAvXr1sHXrVly9ehXR0dH44Ycf4OXlBVNTUzRo0ECaVD58+PAHXUAzYsQIPH78GL169cKZM2dw/fp17Nu3DwMGDFD0i97FxQXBwcGIi4vDkydP8l3Ph8hzMPj5+aFt27b4/PPPUa9ePfTt2zfbPra2tmotLi8cHBxQs2ZNAIC9vT1KliyJx48fa70OIipYypcvj3PnzsHb2xvjx49H1apV0aJFCwQHB2PZsmUAXk8Op6enw9PTE2PHjsXMmTPz/X2Ojo44duwYMjIy0LJlS1SrVg1jx45FsWLFcnz4aG6+++47HDhwAE5OTqhVq1a+6/kgQoEdO3aI+fPni3379in52DsdPnxYfPbZZ8LBwUEAEH/99Ve2fX766Sfh7OwsjI2NRb169cSpU6dy7Cs0NFR4eHgo+v6kpCQBQCQlJeWnfKJC7fnz5yIqKko8f/5c16VQHrzr56Xkd52iq5LatWuHiRMnomXLlmoLptTUVNSoUQNLly7NcfvGjRvh7++PwMBAnDt3DjVq1ICvr2+2WfzHjx+jX79+WLFihdpqIyIqivIUDH/++WeeO7xz5w6OHTuW5/1bt26NmTNnolOnTjluX7RoEQYNGoQBAwagSpUqWL58OczMzLB69Wppn5cvX6Jjx46YPHkyGjZs+M7ve/nyJZKTk2UvIiL6P3kKhmXLlqFy5cqYP38+oqOjs21PSkrC7t278fnnn6N27dpISEhQS3FpaWk4e/YsfHx8/q9gPT34+PhIa0QIIdC/f398+umnOc57vG3OnDmwtraWXk5OTmqplYiosMhTMBw+fBjz5s3DgQMHULVqVVhZWcHNzQ3VqlVDmTJlUKJECQwcOBBly5bFhQsX0L59e7UU9+jRI2RkZGS7N8LOzg5xcXEAgGPHjmHjxo3Yvn07atasiZo1ayIyMjLXPgMCApCUlCS97ty5o5ZaiYgKizzf4Na+fXu0b98ejx49wtGjR3Hr1i08f/4cJUuWRK1atVCrVi1FM+/q0rhxY0U3hxgbG8PY2FiDFRERfdwUPSsJAEqWLImOHTtqoJScv0tfXz/b47wfPHig0+eIEBU1Be3OXMqZun5OioNBm4yMjODp6Yng4GApjDIzMxEcHIyRI0fqtjiiIsDIyAh6enq4f/8+SpUqBSMjow9aOIc0QwiBtLQ0PHz4EHp6ejAyMvqg/nQeDCkpKbh27Zr0PiYmBmFhYbCxsUHZsmXh7+8PPz8/1KlTB/Xq1cPixYuRmpqKAQMG6LBqoqJBT08P5cqVQ2xsLO7fv6/rcug9zMzMULZs2Q8+ra/zYAgNDYW3t7f03t/fH8DrO62DgoLQo0cPPHz4ENOmTUNcXBxq1qyJvXv38mF9RFpiZGSEsmXLIj09XefP8KHc6evrw8DAQC1HdCoh1Px82o9McnIyrK2tkZSUBCsrK12XQ0SkEUp+133wZUQZGRkICwvT2cOeiIhIvRQHw9ixY7Fq1SoAr0OhadOmqF27NpycnKR1UomI6OOlOBi2bNmCGjVqAAB27tyJmJgYXLp0CePGjcPXX3+t9gKJiEi7FAfDo0ePpHsIdu/ejW7duqFixYoYOHDgO+84JiKij4PiYLCzs0NUVBQyMjKwd+9eaYHsZ8+eQV9fX+0FEhGRdim+XHXAgAHo3r07HBwcoFKppAfcnTp1Slr4moiIPl6Kg+Gbb75B1apVcefOHXTr1k167pC+vj4mT56s9gKJiEi7Pug+hhcvXsDExESd9Wgd72MgoqJAo/cxZGRk4Ntvv0Xp0qVhYWGBGzduAACmTp0qXcZKREQfL8XBMGvWLAQFBWH+/PmyBzVVrVoVK1euVGtxRESkfYqDYd26dVixYgV69+4tuwqpRo0auHTpklqLIyIi7VMcDPfu3YOrq2u29szMTLx69UotRRERke4oDoYqVargv//+y9a+ZcsW1KpVSy1FERGR7ii+XHXatGnw8/PDvXv3kJmZiW3btuHy5ctYt24d/vnnH03USEREWqT4iKFDhw7YuXMn/v33X5ibm2PatGmIjo7Gzp07pbugiYjo48X1GHgfAxEVARq9j+HMmTM4depUtvZTp04hNDRUaXdERFTAKA6GESNG4M6dO9na7927hxEjRqilKCIi0h3FwRAVFYXatWtna69VqxaioqLUUhQREemO4mAwNjbGgwcPsrXHxsbCwEDxRU5ERFTAKA6Gli1bIiAgAElJSVJbYmIi/ve///GqJCKiQkDx/+IvXLgQTZo0gbOzs3RDW1hYGOzs7PDbb7+pvUAiItIuxcFQunRpREREYP369QgPD4epqSkGDBiAXr16wdDQUBM1EhGRFuVrUsDc3ByDBw9Wdy1ERFQA5CsYrl69ikOHDiE+Ph6ZmZmybdOmTVNLYUREpBuKg+HXX3/FsGHDULJkSdjb20OlUknbVCoVg4GI6COnOBhmzpyJWbNm4auvvtJEPUREpGOKL1d98uQJunXrpolaiIioAFB8xNCtWzfs378fQ4cO1UQ9Rdb06dPV0k9gYKBa+iGioktxMLi6umLq1Kk4efIkqlWrlu0S1dGjR6utOCIi0j7FwbBixQpYWFjg8OHDOHz4sGybSqViMBARfeQUB0NMTIwm6iAiogJC8eRzlrS0NFy+fBnp6enqrOejolKp1PYiIiooFAfDs2fP8MUXX8DMzAweHh64ffs2AGDUqFGYO3eu2gskIiLtUhwMAQEBCA8PR0hICExMTKR2Hx8fbNy4Ua3FERGR9imeY9i+fTs2btyIBg0ayE6BeHh44Pr162otjoiItE/xEcPDhw9ha2ubrT01NZXnyomICgHFRwx16tTBrl27MGrUKACQwmDlypXw8vJSb3VU6Knrxj6AN/cRqYviYJg9ezZat26NqKgopKenY8mSJYiKisLx48ez3ddAREQfH8Wnkho3bozw8HCkp6ejWrVq2L9/P2xtbXHixAl4enpqokYiItIiRUcMr169wpAhQzB16lT8+uuvmqqJiIh0SNERg6GhIbZu3aqpWoiIqABQfCqpY8eO2L59uwZKISKigkDx5LObmxtmzJiBY8eOwdPTE+bm5rLtfIgeEdHHTXEwrFq1CsWKFcPZs2dx9uxZ2TY+XZWI6OOnKBiEEAgJCYGtrS1MTU01VRMREemQojkGIQTc3Nxw9+5dTdVDREQ6pigY9PT04ObmhoSEBE3VQ0REOqb4qqS5c+di4sSJuHDhgibqISIiHVMcDP369cPp06dRo0YNmJqawsbGRvaiwo8LFBEVboqvSlq8eLEGyiAiooJCcTD4+flpog4iIiogFAdD1lKeuSlbtmy+iyEiIt1THAwuLi7vPDeckZHxQQUREZFuKQ6G8+fPy96/evUK58+fx6JFizBr1iy1FUZERLqhOBhq1KiRra1OnTpwdHTEggUL0LlzZ7UURkREuqH4ctXcuLu748yZM+rqjoiIdETxEUNycrLsvRACsbGx+Oabb+Dm5qa2woiISDcUB0OxYsWyTT4LIeDk5IQ///xTbYUREZFuKA6GgwcPyoJBT08PpUqVgqurKwwMFHdHREQFjOLf5M2aNdNAGUREVFAonnyeM2cOVq9ena199erVmDdvnlqKIiIi3VEcDL/88gsqVaqUrd3DwwPLly9XS1FERKQ7ioMhLi4ODg4O2dpLlSqF2NhYtRRFRES6ozgYnJyccOzYsWztx44dg6Ojo1qKIiIi3VE8+Txo0CCMHTsWr169wqeffgoACA4OxqRJkzB+/Hi1F0hERNqlOBgmTpyIhIQEDB8+HGlpaQAAExMTfPXVV5g8ebLaCyQiIu1SHAwqlQrz5s3D1KlTER0dDVNTU7i5ucHY2FgT9REVKOpadU4IoZZ+iDRBcTAkJSUhIyMDNjY2qFu3rtT++PFjGBgYwMrKSq0FEhGRdimefO7Zs2eOj77YtGkTevbsqZaiiIhIdxQHw6lTp+Dt7Z2tvVmzZjh16pRaiiIiIt1RHAwvX75Eenp6tvZXr17h+fPnaimKiIh0R3Ew1KtXDytWrMjWvnz5cnh6eqqlKCIi0h3Fk88zZ86Ej48PwsPD0bx5cwCv72M4c+YM9u/fr/YCiQqj6dOnq62vwMBAtfVFBOTjiKFRo0Y4ceIEypQpg02bNmHnzp1wdXVFREQEPvnkE03USEREWpSvBRRq1qyJDRs2qLsWIiIqABQHw71797B161ZcuXIFwOu1nrt06cLnJBERFRKKguHnn3+Gv78/0tLSpBvZkpOTMXHiRCxatAjDhw/XSJFERKQ9eZ5j2LVrF0aPHo2RI0fi3r17SExMRGJiIu7du4fhw4djzJgx2L17tyZrJSIiLcjzEcOCBQswefJkzJw5U9bu4OCARYsWwczMDPPnz0ebNm3UXiQREWlPno8Yzp07h759++a6vW/fvjh37pxaiiIiIt3JczBkZGTA0NAw1+2GhobIyMhQS1FERKQ7eQ4GDw8P/P3337lu3759Ozw8PNRSFBER6U6e5xhGjBiBYcOGwdjYGIMHD4aBweuPpqen45dffsGUKVPw888/a6xQIiLSjjwHg5+fHyIjIzFy5EgEBASgQoUKEELgxo0bSElJwejRo9G/f38NlkpERNqg6D6GhQsXomvXrvjjjz9w9epVAEDTpk3Rs2dPNGjQQCMFEhGRdim+87lBgwYMASKiQkzxQ/QKok6dOqF48eLo2rWrrkshIvroFYpgGDNmDNatW6frMoiICoVCEQzNmjWDpaWlrssgIioUdB4MR44cQbt27eDo6AiVSoXt27dn22fp0qVwcXGBiYkJ6tevj9OnT2u/UCKiIkLnwZCamooaNWpg6dKlOW7fuHEj/P39ERgYiHPnzqFGjRrw9fVFfHy8lislKrpUKpXaXlTw5emqpFq1auX5B6r0eUmtW7dG69atc92+aNEiDBo0CAMGDADwem3pXbt2YfXq1Zg8ebKi7wKAly9f4uXLl9L75ORkxX0QERVmeQqGjh07ariMnKWlpeHs2bMICAiQ2vT09ODj44MTJ07kq885c+aodb1dIqLCJk/BoKvFxh89eoSMjAzY2dnJ2u3s7HDp0iXpvY+PD8LDw5GamooyZcpg8+bN8PLyyrHPgIAA+Pv7S++Tk5Ph5OSkmQEQEX2E8rXmc0Hz77//5nlfY2NjGBsba7AaIqKPm+JgyMjIwPfff49Nmzbh9u3bSEtLk21//Pix2oorWbIk9PX18eDBA1n7gwcPYG9vr7bvISKi/6P4qqTp06dj0aJF6NGjB5KSkuDv74/OnTtDT08P33zzjVqLMzIygqenJ4KDg6W2zMxMBAcH53qqiIiIPoziI4b169fj119/Rdu2bfHNN9+gV69eqFChAqpXr46TJ09i9OjRivpLSUnBtWvXpPcxMTEICwuDjY0NypYtC39/f/j5+aFOnTqoV68eFi9ejNTUVOkqJSIiUi/FwRAXF4dq1aoBACwsLJCUlAQA+OyzzzB16lTFBYSGhsLb21t6nzUx7Ofnh6CgIPTo0QMPHz7EtGnTEBcXh5o1a2Lv3r3ZJqSJiEg9FAdDmTJlEBsbi7Jly6JChQrYv38/ateujTNnzuRrUrdZs2YQQrxzn5EjR2LkyJGK+yYiIuUUzzF06tRJOuc/atQoTJ06FW5ubujXrx8GDhyo9gKJiEi7FB8xzJ07V/pzjx494OzsjOPHj8PNzQ3t2rVTa3FERKR9ioPhyJEjaNiwobTmc9bCPenp6Thy5AiaNGmi9iKJiEh7FJ9K8vb2zvFehaSkJNkkMhERfZwUB4MQIscH6iUkJMDc3FwtRRERke7k+VRS586dAbx+/G7//v1lVyBlZGQgIiICDRs2VH+FRESkVXkOBmtrawCvjxgsLS1hamoqbTMyMkKDBg0waNAg9VdIRERaledgWLNmDQDAxcUFEyZM4GkjIqJCSvFVSVmP4H748CEuX74MAHB3d0epUqXUWxkREemE4snnZ8+eYeDAgXBwcECTJk3QpEkTODo64osvvsCzZ880USMREWmR4mAYN24cDh8+jJ07dyIxMRGJiYn4+++/cfjwYYwfP14TNRIRkRYpPpW0detWbNmyBc2aNZPa2rRpA1NTU3Tv3h3Lli1TZ31ERKRl+TqVlNOTTW1tbXkqiYioEFAcDF5eXggMDMSLFy+ktufPn2P69OlcPIeIqBDI86kkfX19xMbGYvHixWjVqhXKlCmDGjVqAADCw8NhYmKCffv2aaxQIiLSjjwHQ9aaCdWqVcPVq1exfv16XLp0CQDQq1cv9O7dW3bTGxERfZwUTz4DgJmZGe9yJiIqpBQFw8qVK2FhYfHOfZSu+UxERAWLomBYvnw59PX1c92uUqkYDEREHzlFwRAaGgpbW1tN1UJERAVAni9XzWkNBiIiKnzyHAxZVyUREVHhludgCAwMfO/EMxERffzyPMeQ9bhtIiIq3BQ/EoOIiAo3BgMREcnkKRh27NiBV69eaboWIiIqAPIUDJ06dUJiYiKA1w/Ti4+P12RNRESkQ3kKhlKlSuHkyZMAXl+2ynsaiIgKrzxdlTR06FB06NABKpUKKpUK9vb2ue6bkZGhtuKIiEj78hQM33zzDXr27Ilr166hffv2WLNmDYoVK6bh0oiISBfyfB9DpUqVUKlSJQQGBqJbt24wMzPTZF1ERKQjitdjyLrR7eHDh7h8+TIAwN3dHaVKlVJvZUREpBOK72N49uwZBg4cCEdHRzRp0gRNmjSBo6MjvvjiCzx79kwTNRIRkRYpDoZx48bh8OHD2LFjBxITE5GYmIi///4bhw8fxvjx4zVRIxGRTmRdcKOO18dE8amkrVu3YsuWLWjWrJnU1qZNG5iamqJ79+5YtmyZOusjIiIty9epJDs7u2zttra2PJVERFQIKA4GLy8vBAYG4sWLF1Lb8+fPMX36dHh5eam1OCIi0j7Fp5KWLFkCX19flClTBjVq1AAAhIeHw8TEBPv27VN7gUREpF2Kg6Fq1aq4evUq1q9fj0uXLgEAevXqhd69e8PU1FTtBRIRkXYpDgYAMDMzw6BBg9RdCxERFQD5CgYiovyaPn26WvrhqpKaw4V6iIhIhsFAREQyDAYiIpLJVzAkJiZi5cqVCAgIwOPHjwEA586dw71799RaHBERaZ/iyeeIiAj4+PjA2toaN2/exKBBg2BjY4Nt27bh9u3bWLdunSbqJCIiLVF8xODv74/+/fvj6tWrMDExkdrbtGmDI0eOqLU4IiLSPsXBcObMGQwZMiRbe+nSpREXF6eWooiISHcUB4OxsTGSk5OztV+5coWL9RARFQKKg6F9+/aYMWMGXr16BeD188pv376Nr776Cl26dFF7gUREpF2Kg+G7775DSkoKbG1t8fz5czRt2hSurq6wtLTErFmzNFEjERFpkeKrkqytrXHgwAEcPXoUERERSElJQe3ateHj46OJ+oiISMvy/aykxo0bo3HjxuqshYiICgDFwfDDDz/k2K5SqWBiYgJXV1c0adIE+vr6H1wcERFpn+Jg+P777/Hw4UM8e/YMxYsXBwA8efIEZmZmsLCwQHx8PMqXL49Dhw7ByclJ7QUTEZFmKZ58nj17NurWrYurV68iISEBCQkJuHLlCurXr48lS5bg9u3bsLe3x7hx4zRRLxERaZjiI4YpU6Zg69atqFChgtTm6uqKhQsXokuXLrhx4wbmz5/PS1eJiD5Sio8YYmNjkZ6enq09PT1duvPZ0dERT58+/fDqiIhI6xQHg7e3N4YMGYLz589LbefPn8ewYcPw6aefAgAiIyNRrlw59VVJRERaozgYVq1aBRsbG3h6esLY2BjGxsaoU6cObGxssGrVKgCAhYUFvvvuO7UXS0REmqd4jsHe3h4HDhzApUuXcOXKFQCAu7s73N3dpX28vb3VVyEREWlVvm9wq1SpEipVqqTOWoiIqADIVzDcvXsXO3bswO3bt5GWlibbtmjRIrUURkREuqE4GIKDg9G+fXuUL18ely5dQtWqVXHz5k0IIVC7dm1N1EhERFqkePI5ICAAEyZMQGRkJExMTLB161bcuXMHTZs2Rbdu3TRRIxERaZHiYIiOjka/fv0AAAYGBnj+/DksLCwwY8YMzJs3T+0FEhGRdikOBnNzc2lewcHBAdevX5e2PXr0SH2VERGRTiieY2jQoAGOHj2KypUro02bNhg/fjwiIyOxbds2NGjQQBM1EhGRFikOhkWLFiElJQUAMH36dKSkpGDjxo1wc3PjFUlERIWA4mAoX7689Gdzc3MsX75crQUREZFuKZ5jKF++PBISErK1JyYmykKDiIg+ToqD4ebNm8jIyMjW/vLlS9y7d08tRRERke7k+VTSjh07pD/v27cP1tbW0vuMjAwEBwfDxcVFrcUREZH25TkYOnbsCOD12s5+fn6ybYaGhnBxceETVYmICoE8B0NmZiYAoFy5cjhz5gxKliypsaKIiEh3FF+VFBMTo4k6iIiogMjX01WDg4MRHByM+Ph46Ugiy+rVq9VSGBER6YbiYJg+fTpmzJiBOnXqwMHBASqVShN1ERGRjigOhuXLlyMoKAh9+/bVRD1ERKRjiu9jSEtLQ8OGDTVRCxERFQCKg+HLL7/Ehg0bNFELEREVAIpPJb148QIrVqzAv//+i+rVq8PQ0FC2nQ/SIyL6uCkOhoiICNSsWRMAcOHCBdk2TkQTEX38FAfDoUOHNFEHEREVEIrnGLJcu3YN+/btw/PnzwEAQgi1FUVERLqjOBgSEhLQvHlzVKxYEW3atEFsbCwA4IsvvsD48ePVXiAREWmX4mAYN24cDA0Ncfv2bZiZmUntPXr0wN69e9VaXF79888/cHd3h5ubG1auXKmTGoiICgvFcwz79+/Hvn37UKZMGVm7m5sbbt26pbbC8io9PR3+/v44dOgQrK2t4enpiU6dOqFEiRJar4WIqDBQfMSQmpoqO1LI8vjxYxgbG6ulKCVOnz4NDw8PlC5dGhYWFmjdujX279+v9TqIiAoLxcHwySefYN26ddJ7lUqFzMxMzJ8/H97e3ooLOHLkCNq1awdHR0eoVCps37492z5Lly6Fi4sLTExMUL9+fZw+fVradv/+fZQuXVp6X7p0aa4kR0T0ARSfSpo/fz6aN2+O0NBQpKWlYdKkSbh48SIeP36MY8eOKS4gNTUVNWrUwMCBA9G5c+ds2zdu3Ah/f38sX74c9evXx+LFi+Hr64vLly/D1tZW8fcREdG7KT5iqFq1Kq5cuYLGjRujQ4cOSE1NRefOnXH+/HlUqFBBcQGtW7fGzJkz0alTpxy3L1q0CIMGDcKAAQNQpUoVLF++HGZmZtLjvR0dHWVHCPfu3YOjo2Ou3/fy5UskJyfLXkRE9H/ytR6DtbU1vv76a3XXkk1aWhrOnj2LgIAAqU1PTw8+Pj44ceIEAKBevXq4cOEC7t27B2tra+zZswdTp07Ntc85c+Zg+vTpGq+diOhjpfiIYc2aNdi8eXO29s2bN2Pt2rVqKSrLo0ePkJGRATs7O1m7nZ0d4uLiAAAGBgb47rvv4O3tjZo1a2L8+PHvvCIpICAASUlJ0uvOnTtqrZmI6GOnOBjmzJmT43rPtra2mD17tlqKUqp9+/a4cuUKrl27hsGDB79zX2NjY1hZWcleRET0fxQHw+3bt1GuXLls7c7Ozrh9+7ZaispSsmRJ6Ovr48GDB7L2Bw8ewN7eXq3fRURErykOBltbW0RERGRrDw8PV/tNZUZGRvD09ERwcLDUlpmZieDgYHh5ean1u4iI6DXFk8+9evXC6NGjYWlpiSZNmgAADh8+jDFjxqBnz56KC0hJScG1a9ek9zExMQgLC4ONjQ3Kli0Lf39/+Pn5oU6dOqhXrx4WL16M1NRUDBgwQPF3ERHR+ykOhm+//RY3b95E8+bNYWDw+uOZmZno169fvuYYQkNDZTfG+fv7AwD8/PwQFBSEHj164OHDh5g2bRri4uJQs2ZN7N27N9uENBERqYeiYBBCIC4uDkFBQZg5cybCwsJgamqKatWqwdnZOV8FNGvW7L2P7B45ciRGjhyZr/6JiEgZxcHg6uqKixcvws3NDW5ubpqqi4iIdETR5LOenh7c3NyQkJCgqXqIiEjHFF+VNHfuXEycODHbes9ERFQ4KJ587tevH549e4YaNWrAyMgIpqamsu2PHz9WW3FERKR9ioNh8eLFGiiDiIgKCsXB4Ofnp4k6iIiogFA8xwAA169fx5QpU9CrVy/Ex8cDAPbs2YOLFy+qtTgiItI+xcFw+PBhVKtWDadOncK2bduQkpIC4PUjMQIDA9VeIBERaZfiYJg8eTJmzpyJAwcOwMjISGr/9NNPcfLkSbUWR0RE2qc4GCIjI3Ncbc3W1haPHj1SS1FERKQ7ioOhWLFiiI2NzdZ+/vx5lC5dWi1FERGR7igOhp49e+Krr75CXFwcVCoVMjMzcezYMUyYMAH9+vXTRI1ERKRFioNh9uzZqFSpEpycnJCSkoIqVaqgSZMmaNiwIaZMmaKJGomISIsU38dgZGSEX3/9FdOmTUNkZCRSUlJQq1YtPlCPiKiQyHMwZGZmYsGCBdixYwfS0tLQvHlzBAYGZnskBhERfdzyfCpp1qxZ+N///gcLCwuULl0aS5YswYgRIzRZGxER6UCeg2HdunX4+eefsW/fPmzfvh07d+7E+vXrkZmZqcn6iIhIy/IcDLdv30abNm2k9z4+PlCpVLh//75GCiMiIt3IczCkp6fDxMRE1mZoaIhXr16pvSgiItKdPE8+CyHQv39/GBsbS20vXrzA0KFDYW5uLrVt27ZNvRUSEZFW5TkYcnrcdp8+fdRaDBER6V6eg2HNmjWarIOIiAqIfK3HQEREhReDgYiIZBgMREQkw2AgIiIZBgMREckwGIiISIbBQEREMgwGIiKSYTAQEZEMg4GIiGQYDEREJMNgICIiGQYDERHJMBiIiEgmz4/dJiKi/Js+fbra+goMDFRbXznhEQMREckwGIiISIbBQEREMgwGIiKSYTAQEZEMg4GIiGQYDEREJMNgICIiGQYDERHJFPk7n4UQAIDk5GSd1vHixQu19KPrcSilrnEDH9fYi+q4Af5bV4f8jD3rM1m/895FJfKyVyF29+5dODk56boMIiKtuHPnDsqUKfPOfYp8MGRmZuL+/fuwtLSESqXSSQ3JyclwcnLCnTt3YGVlpZMadIHjLlrjBoru2AvCuIUQePr0KRwdHaGn9+5ZhCJ/KklPT++96aktVlZWReo/liwcd9FTVMeu63FbW1vnaT9OPhMRkQyDgYiIZBgMBYCxsTECAwNhbGys61K0iuMuWuMGiu7YP7ZxF/nJZyIikuMRAxERyTAYiIhIhsFAREQyDIYPcOTIEbRr1w6Ojo5QqVTYvn27bLsQAtOmTYODgwNMTU3h4+ODq1evvrPP8PBw9OrVC05OTjA1NUXlypWxZMmSbPuFhISgdu3aMDY2hqurK4KCgtQ4snebM2cO6tatC0tLS9ja2qJjx464fPmybJ8XL15gxIgRKFGiBCwsLNClSxc8ePDgnf1evnwZ3t7esLOzg4mJCcqXL48pU6bg1atXsv02b96MSpUqwcTEBNWqVcPu3bvVPsacLFu2DNWrV5euRffy8sKePXuk7fkZ85uuXbsGS0tLFCtWLNs2XY05J3PnzoVKpcLYsWOltvyM/ebNm1CpVNleJ0+elO1XkMaeHy4uLtnGOHfuXNk+ERER+OSTT2BiYgInJyfMnz9fR9X+f4Lybffu3eLrr78W27ZtEwDEX3/9Jds+d+5cYW1tLbZv3y7Cw8NF+/btRbly5cTz589z7XPVqlVi9OjRIiQkRFy/fl389ttvwtTUVPz444/SPjdu3BBmZmbC399fREVFiR9//FHo6+uLvXv3amqoMr6+vmLNmjXiwoULIiwsTLRp00aULVtWpKSkSPsMHTpUODk5ieDgYBEaGioaNGggGjZs+M5+r1+/LlavXi3CwsLEzZs3xd9//y1sbW1FQECAtM+xY8eEvr6+mD9/voiKihJTpkwRhoaGIjIyUmPjzbJjxw6xa9cuceXKFXH58mXxv//9TxgaGooLFy4IIfI35ixpaWmiTp06onXr1sLa2lq2TZdjftvp06eFi4uLqF69uhgzZozUnp+xx8TECADi33//FbGxsdIrLS1N2qcgjf1Nqamped7X2dlZzJgxQzbGN/9bSUpKEnZ2dqJ3797iwoUL4o8//hCmpqbil19+0UTpecJgUJO3gyEzM1PY29uLBQsWSG2JiYnC2NhY/PHHH4r6Hj58uPD29pbeT5o0SXh4eMj26dGjh/D19c1f8R8oPj5eABCHDx8WQrwep6Ghodi8ebO0T3R0tAAgTpw4oajvcePGicaNG0vvu3fvLtq2bSvbp379+mLIkCEfMIL8K168uFi5cuUHj3nSpEmiT58+Ys2aNdmCoaCM+enTp8LNzU0cOHBANG3aVAqG/I49KxjOnz+f6z4FZexv++qrr4Srq6sYPXq02Lt3r3jx4kWu+zo7O4vvv/8+1+0///yzKF68uHj58qWsf3d3d3WWrAhPJWlITEwM4uLi4OPjI7VZW1ujfv36OHHihKK+kpKSYGNjI70/ceKErF8A8PX1VdyvuiQlJQGAVOPZs2fx6tUrWY2VKlVC2bJlFdV47do17N27F02bNpXaCsrYMzIy8OeffyI1NRVeXl4fNOaDBw9i8+bNWLp0aY7bC8qYR4wYgbZt22ar5UN/3u3bt4etrS0aN26MHTt2yLYVlLG/7auvvsKMGTOQkJCA3r17w8bGBu3bt8fy5ctx+/btbPvPnTsXJUqUQK1atbBgwQKkp6dL206cOIEmTZrAyMhIavP19cXly5fx5MkTrYznbUX+WUmaEhcXBwCws7OTtdvZ2Unb8uL48ePYuHEjdu3aJes7p36Tk5Px/PlzmJqafkDlymRmZmLs2LFo1KgRqlatKtVnZGSU7Vx5XsfesGFDnDt3Di9fvsTgwYMxY8YMaVtuY1fyd/ohIiMj4eXlhRcvXsDCwgJ//fUXqlSpgrCwsHyNOSEhAf3798fvv/+e6zN0dD1mAPjzzz9x7tw5nDlzJtu2/P68LSws8N1336FRo0bQ09PD1q1b0bFjR2zfvh3t27eX+tb12HNSvHhx9OrVC7169UJmZiZOnjyJXbt2YdmyZRg2bBg8PDywf/9+ODo6YvTo0ahduzZsbGxw/PhxBAQEIDY2FosWLQLweozlypWT9Z815ri4OBQvXlzr42Mw6FDr1q3x33//AQCcnZ1x8eJF2fYLFy6gQ4cOCAwMRMuWLXVR4nuNGDECFy5cwNGjRxV9zsPDA7du3QIAfPLJJ7JJ3I0bN+Lp06cIDw/HxIkTsXDhQkyaNEmtdeeXu7s7wsLCkJSUhC1btsDPzw+HDx/O02dzGvOgQYPw+eefo0mTJpos+4PcuXMHY8aMwYEDB2BiYpKvPnIae8mSJeHv7y/tU7duXdy/fx8LFiyQguFj8PTpU9y/fx+xsbF4+PAhTE1N4ezsDENDQwCQjbF69eowMjLCkCFDMGfOnAJ7JzSDQUPs7e0BAA8ePICDg4PU/uDBA9SsWRMAsHLlSjx//hwApH9EWaKiotC8eXMMHjwYU6ZMydb321d8PHjwAFZWVlo9Whg5ciT++ecfHDlyRPaEWnt7e6SlpSExMVH2f5EPHjyQ/l52794tXW30ds1Z62NUqVIFGRkZGDx4MMaPHw99ff1cx57Vr6YZGRnB1dUVAODp6YkzZ85gyZIl6NGjR77GfPDgQezYsQMLFy4E8PpKtszMTBgYGGDFihUYOHCgzsd89uxZxMfHo3bt2lJbRkYGjhw5gp9++gn79u37oJ/3m+rXr48DBw5I73U99tzExMRg06ZN2L17N44fPw5nZ2e0bt0aq1atgre39zsDtH79+khPT8fNmzfh7u6e6xgB6G6cOpvdKGSQy+TzwoULpbakpKQ8TT5fuHBB2NraiokTJ+a4fdKkSaJq1aqytl69emlt8jkzM1OMGDFCODo6iitXrmTbnjUZuWXLFqnt0qVL+Zp8Xrt2rTAwMJCuVOnevbv47LPPZPt4eXnpbDLS29tb+Pn55XvMUVFRIjIyUnrNnDlTWFpaisjISPH48WMhhO7HnJycLKsxMjJS1KlTR/Tp00dERkaq9ef95Zdfilq1aknvdT323EybNk00b95cfPfddyI6OlrRZ3///Xehp6cn/XyzJp/fvBorICBAp5PPDIYP8PTpU3H+/Hlx/vx5AUAsWrRInD9/Xty6dUsI8fpy1WLFiom///5bREREiA4dOrz3ctXIyEhRqlQp0adPH9nlbfHx8dI+WZerTpw4UURHR4ulS5dq9XLVYcOGCWtraxESEiKr8dmzZ9I+Q4cOFWXLlhUHDx4UoaGhwsvLS3h5eb2z399//11s3LhRREVFievXr4uNGzcKR0dH0bt3b2mfY8eOCQMDA7Fw4UIRHR0tAgMDtXb54uTJk8Xhw4dFTEyMiIiIEJMnTxYqlUrs379fCJG/Mb8tp6uSdDnm3Lx5VZIQ+Rt7UFCQ2LBhg4iOjhbR0dFi1qxZQk9PT6xevVrapyCOXYjXV+Jl1Z3b6+XLl+L48ePi+++/F2FhYeL69evi999/F6VKlRL9+vWT+kpMTBR2dnaib9++4sKFC+LPP/8UZmZmvFz1Y3Xo0CEBINvLz89PCPH6/6ynTp0q7OzshLGxsWjevLm4fPnyO/sMDAzMsU9nZ+ds312zZk1hZGQkypcvL9asWaOZQeYgp/oAyGp4/vy5GD58uChevLgwMzMTnTp1ErGxse/s988//xS1a9cWFhYWwtzcXFSpUkXMnj07W5Bu2rRJVKxYURgZGQkPDw+xa9cuTQwzm4EDBwpnZ2dhZGQkSpUqJZo3by6FghD5G/PbcgoGIXQ35ty8HQz5GXtQUJCoXLmyMDMzE1ZWVqJevXqyS16zFLSxC/H6ctLc/jvIekVHR4uzZ8+K+vXrC2tra2FiYiIqV64sZs+ene3y1vDwcNG4cWNhbGwsSpcuLebOnaujkb3Gp6sSEZEM72MgIiIZBgMREckwGIiISIbBQEREMgwGIiKSYTAQEZEMg4GIiGQYDEREJMNgoELlm2++kR5SWBAJITB48GDY2NhApVIhLCwMzZo1ky2TqQkhISFQqVRITEzM0/5Zy26GhYVppJ5Vq1a994nB/fv3R8eOHTXy/ZMnT8aoUaM00ndhwGAoguLi4jBq1CiUL18exsbGcHJyQrt27RAcHKzr0j7YhAkTZONQ1y+Xt9cntrGxQdOmTaXHpufV3r17ERQUhH/++QexsbHSGhbqlFPQNGzYELGxsbC2ts5TH05OTrL6lAbLu7x48QJTp05FYGCgos/1799f9jMoUaIEWrVqhYiIiFz3efu1du1aAK//naxduxY3btz44PEURgyGIubmzZvw9PTEwYMHsWDBAkRGRmLv3r3w9vbGiBEjdF3eB7OwsECJEiU01v+///6L2NhYHDlyBI6Ojvjss8/eu+j9m65fvw4HBwc0bNgQ9vb2MDDQzpPvjYyMYG9vD5VKlaf9sx5xron6tmzZAisrKzRq1EjxZ1u1aoXY2FjExsYiODgYBgYG+Oyzz6TtS5Yskba/+fLx8YGLiwvatm0LAChZsiR8fX2xbNkytY2rUNHpk5pI61q3bi1Kly4tW4w8y5MnT6Q/37p1S7Rv316Ym5sLS0tL0a1bNxEXFydtDwwMFDVq1BCrVq0STk5OwtzcXAwbNkykp6eLefPmCTs7O1GqVCkxc+ZM2XcAED///LNo1aqVMDExEeXKlcv24LSIiAjh7e0tTExMhI2NjRg0aJB4+vSptP3QoUOibt26wszMTFhbW4uGDRuKmzdvyurK+jPeerDZoUOHhBBC3L59W3Tr1k1YW1uL4sWLi/bt24uYmJhc/95yWp84IiJCABB///231BYZGSlatWolzM3Nha2trejTp494+PChEEIIPz+/HB+M+PYD6V68eCHGjx8vHB0dhZmZmahXr55Ud5ajR4+Kpk2bClNTU1GsWDHRsmVL8fjx42zfAUDExMRID3x88uSJSEpKEiYmJmL37t2yPrdt2yYsLCxEamqqbLxZf37z5efnJ9auXStsbGyyPRCuQ4cOok+fPrn+XbZt21ZMmDBB1paeni7GjRsnrK2thY2NjZg4caLo16+f6NChg7SPn5+f7L0QQvz3338CgOzpw2+bOXOmMDc3F2FhYbL2tWvXijJlyuT6uaKMwVCEJCQkCJVKJWbPnv3O/TIyMkTNmjVF48aNRWhoqDh58qTw9PQUTZs2lfYJDAwUFhYWomvXruLixYtix44dwsjISPj6+opRo0aJS5cuidWrVwsA4uTJk9LnAIgSJUqIX3/9VVy+fFlMmTJF6Ovri6ioKCGEECkpKcLBwUF07txZREZGiuDgYFGuXDnpibWvXr0S1tbWYsKECeLatWsiKipKBAUFSY86fzMYnj59Krp37y5atWolPRr85cuXIi0tTVSuXFkMHDhQREREiKioKPH5558Ld3d32YLsb3o7GJ49eyYmTJggAIg9e/YIIV4Ha6lSpURAQICIjo4W586dEy1atBDe3t5CiNePV54xY4YoU6aM7FHqbwfDl19+KRo2bCiOHDkirl27JhYsWCCMjY2ltS/Onz8vjI2NxbBhw0RYWJi4cOGC+PHHH8XDhw9FYmKi8PLyEoMGDZLGnJ6eLgsGIYTo2rVrtl/eXbp0kdreHG96errYunWrACAuX74sYmNjRWJionj27JmwtrYWmzZtkvp48OCBMDAwEAcPHsz135e1tbX4888/ZW3z5s0TxYsXF1u3bhVRUVHiiy++EJaWlu8MhqdPn4ohQ4YIV1dXkZGRkeN37dy5U+jp6clqzBIdHS0FJ8kxGIqQU6dOCQBi27Zt79xv//79Ql9fX9y+fVtqu3jxogAgTp8+LYR4/QvYzMxMJCcnS/v4+voKFxcX2X+k7u7uYs6cOdJ7AGLo0KGy76tfv74YNmyYEEKIFStWiOLFi8uOaHbt2iX09PREXFycSEhIEABESEhIjrW/GQxC5Px/mb/99ptwd3cXmZmZUtvLly+Fqamp2LdvX479Zv2iNDU1Febm5kKlUgkAwtPTU1pg5dtvvxUtW7aUfe7OnTvSL1QhhPj++++zPUL9zWC4deuW0NfXF/fu3ZPt07x5cxEQECCEeL0oU6NGjXKs8+3+srwdDH/99Zd0dCCEkI4iskLu7SB8+/NZhg0bJlq3bi29/+6770T58uVlf7dvevLkiQAgjhw5Imt3cHAQ8+fPl96/evVKlClTJlsw6OvrC3Nzc2Fubi4ACAcHB3H27Nkcvys6OlpYWVmJr7/+OsftSUlJ7/y3VJRxjqEIEXl8wnp0dDScnJykJTaB18tsFitWDNHR0VKbi4sLLC0tpfd2dnaoUqUK9PT0ZG3x8fGy/r28vLK9z+o3OjoaNWrUgLm5ubS9UaNGyMzMxOXLl2FjY4P+/fvD19cX7dq1k84pKxEeHo5r167B0tISFhYWsLCwgI2NDV68eIHr16+/87MbN27E+fPnsXXrVri6uiIoKEhaljU8PByHDh2S+rSwsEClSpUA4L39ZomMjERGRgYqVqwo6+fw4cNSH2FhYWjevLmiMb+tTZs2MDQ0xI4dOwAAW7duhZWVFXx8fBT1M2jQIOzfvx/37t0DAAQFBUkTwDnJWsr2zaUvk5KSEBsbi/r160ttBgYGqFOnTrbPe3t7IywsDGFhYTh9+jR8fX3RunVraT3pN/vs2LEjmjZtim+//TbHWrKWGH327JmCERcNXPO5CHFzc4NKpcKlS5fU0t/b61SrVKoc2zIzM9XyfVnWrFmD0aNHY+/evdi4cSOmTJmCAwcOoEGDBnn6fEpKCjw9PbF+/fps20qVKvXOzzo5OcHNzQ1ubm5IT09Hp06dcOHCBRgbGyMlJQXt2rXDvHnzsn3uzXW/31ebvr4+zp49C319fdk2CwsLAO9eMzmvjIyM0LVrV2zYsAE9e/bEhg0b0KNHD8WTzbVq1UKNGjWwbt06tGzZEhcvXsSuXbty3b9EiRJQqVR48uRJvuo2NzeX1twGXq+bbm1tjV9//RUzZ84EAGRmZuLzzz+Hnp4e1q9fn2tIPX78GMD7f+ZFEY8YihAbGxv4+vpi6dKlSE1NzbY961LEypUr486dO7hz5460LSoqComJiahSpcoH13Hy5Mls7ytXrix9d3h4uKy+Y8eOQU9PD+7u7lJbrVq1EBAQgOPHj6Nq1arYsGFDjt9lZGSEjIwMWVvt2rVx9epV2NrawtXVVfbK6+WcANC1a1cYGBjg559/lvq9ePEiXFxcsvX75hHQu9SqVQsZGRmIj4/P1kfWwvDVq1d/56XFOY05J71798bevXtx8eJFHDx4EL17935nnwBy7PfLL79EUFAQ1qxZAx8fH9mRZk79VKlSBVFRUVKbtbU1HBwccOrUKaktPT0dZ8+efe8YVCoV9PT0pCMRAJgyZQqOHz+Ov//+W3ZE+7YLFy7A0NAQHh4e7/2eoobBUMQsXboUGRkZqFevHrZu3YqrV68iOjoaP/zwg3SKx8fHB9WqVUPv3r1x7tw5nD59Gv369UPTpk1zPLxXavPmzVi9ejWuXLmCwMBAnD59GiNHjgTw+peViYkJ/Pz8cOHCBRw6dAijRo1C3759YWdnh5iYGAQEBODEiRO4desW9u/fj6tXr0rB8jYXFxdERETg8uXLePToEV69eoXevXujZMmS6NChA/777z/ExMQgJCQEo0ePxt27d/M8DpVKhdGjR2Pu3Ll49uwZRowYgcePH6NXr144c+YMrl+/jn379mHAgAF5+kUNABUrVkTv3r3Rr18/bNu2DTExMTh9+jTmzJkj/Z94QEAAzpw5g+HDhyMiIgKXLl3CsmXL8OjRI2nMp06dws2bN/Ho0aNcj9iaNGkCe3t79O7dG+XKlZOdynmbs7MzVCoV/vnnHzx8+BApKSnSts8//xx3797Fr7/+ioEDB753jL6+vjh69KisbcyYMZg7dy62b9+OS5cuYfjw4TneM/Hy5UvExcUhLi4O0dHRGDVqlHSkBgCbNm3C3LlzsXjxYlhaWkr7Zr3erPu///7DJ598opYjsEJH15McpH33798XI0aMkNYvLl26tGjfvr3sksi8Xq76ppwmet+eCAUgli5dKlq0aCGMjY2Fi4uL2Lhxo+wz77pcNS4uTnTs2FE4ODgIIyMj4ezsLKZNmyZNeL9dV3x8vGjRooWwsLCQXa4aGxsr+vXrJ0qWLCmMjY1F+fLlxaBBg0RSUlKOf2c5Xa4qhBCpqamiePHiYt68eUIIIa5cuSI6deokihUrJkxNTUWlSpXE2LFjpcnY900+CyFEWlqamDZtmnBxcRGGhobCwcFBdOrUSUREREj7hISEiIYNGwpjY2NRrFgx4evrK00MX758WTRo0ECYmprmeLnqmyZNmiQAiGnTpr13vDNmzBD29vZCpVJJV4ll6du3b46Xrubk4sWLwtTUVCQmJkptr169EmPGjBFWVlaiWLFiwt/fP8fLVfHGJbOWlpaibt26YsuWLdI+zZo1e+c6zIGBgdK+7u7u4o8//nhvvUUR13wmrVKpVPjrr7809qgD0o3mzZvDw8MDP/zwQ57279atG2rXro2AgAANV5azPXv2YPz48YiIiNDaTYYfE55KIqJ8e/LkCf766y+EhIQounN+wYIF0mS6LqSmpmLNmjUMhVzwb4WI8q1WrVp48uQJ5s2bJ7s44H1cXFx0+hC7rl276uy7PwY8lURERDI8lURERDIMBiIikmEwEBGRDIOBiIhkGAxERCTDYCAiIhkGAxERyTAYiIhIhsFAREQy/w//MJ7ut+h9bAAAAABJRU5ErkJggg==", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "\n", "# Constants\n", "PRECIP_BINS = [(10, 20), (20, 30), (30, 40), (40, 50), (50, np.inf)]\n", "BIN_LABELS = ['10-20', '20-30', '30-40', '40-50', '≥50']\n", "\n", "# Function to calculate **normalized** frequency per 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(nexrad, current):\n", " return ((current / nexrad) - 1) * 100\n", "\n", "def plot_precip_intensity_percentage_bar(month=month):\n", " # Step 1: Use all occurrences of precipitation values for each dataset\n", " all_radar_current = wrf_dbz['REFL_10CM'].where(wrf_dbz['REFL_10CM'] >= 10)\n", " all_radar_nexrad = nexrad_ds.where(nexrad_ds >= 10)\n", "\n", " # Step 2: Calculate normalized percentage for each bin\n", " percent_current = calculate_percentage_per_bin(all_radar_current, PRECIP_BINS)\n", " percent_mrms = calculate_percentage_per_bin(all_radar_nexrad, PRECIP_BINS)\n", "\n", " relative_change_future = calculate_relative_change(np.array(percent_mrms), np.array(percent_current))\n", "\n", " # Step 3: 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(PRECIP_BINS))\n", "\n", " # Plot bars for each dataset\n", " ax1.bar(x + bar_width / 2, percent_mrms, width=bar_width, color='gray', label='NEXRAD', zorder=2)\n", " ax1.bar(x - bar_width / 2, percent_current, width=bar_width, color='black', label='Current', zorder=2)\n", "\n", " # Labels and formatting\n", " ax1.set_xlabel('Composite Reflectivity (dBZ)')\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, 300) # Since we're dealing with percentages\n", "\n", " # Add legend\n", " ax1.legend(loc='upper right')\n", "\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(-40, 1500) # Adjust based on expected range of relative changes\n", " ax2.legend(loc='upper left')\n", " '''\n", " # Convert month code to full name\n", " month_dict = {'04': 'April', '05': 'May', '06': 'June'}\n", " month_name = month_dict.get(month, month)\n", "\n", " # Title and layout adjustments\n", " ax1.set_title(f'{month_name}')\n", " fig.tight_layout()\n", "\n", " plt.show()\n", "\n", "\n", "plot_precip_intensity_percentage_bar()\n", "\n", "######## Need to normalize the frequency because the nexrad is higher resolution" ] } ], "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.10" } }, "nbformat": 4, "nbformat_minor": 2 }