{ "cells": [ { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Running on local URL: http://127.0.0.1:7861\n", "\n", "To create a public link, set `share=True` in `launch()`.\n" ] }, { "data": { "text/html": [ "
" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/plain": [] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Gradio app that takes seismic waveform as input and marks 2 phases on the waveform as output.\n", "\n", "import gradio as gr\n", "import numpy as np\n", "import pandas as pd\n", "from phasehunter.data_preparation import prepare_waveform\n", "import torch\n", "import io\n", "\n", "from scipy.stats import gaussian_kde\n", "from scipy.signal import resample\n", "from scipy.interpolate import interp1d\n", "\n", "from bmi_topography import Topography\n", "import earthpy.spatial as es\n", "\n", "import obspy\n", "from obspy.clients.fdsn import Client\n", "from obspy.clients.fdsn.header import FDSNNoDataException, FDSNTimeoutException, FDSNInternalServerException\n", "from obspy.geodetics.base import locations2degrees\n", "from obspy.taup import TauPyModel\n", "from obspy.taup.helper_classes import SlownessModelError\n", "\n", "from obspy.clients.fdsn.header import URL_MAPPINGS\n", "\n", "import matplotlib.pyplot as plt\n", "import matplotlib.dates as mdates\n", "from mpl_toolkits.axes_grid1 import ImageGrid\n", "\n", "from glob import glob\n", "\n", "\n", "def resample_waveform(waveform, original_freq, target_freq):\n", " \"\"\"\n", " Resample a waveform from original frequency to target frequency using SciPy's resample function.\n", " \n", " Args:\n", " waveform (numpy.ndarray): The input waveform as a 1D array.\n", " original_freq (float): The original sampling frequency of the waveform.\n", " target_freq (float): The target sampling frequency of the waveform.\n", " \n", " Returns:\n", " resampled_waveform (numpy.ndarray): The resampled waveform as a 1D array.\n", " \"\"\"\n", " # Calculate the resampling ratio\n", " resampling_ratio = target_freq / original_freq\n", " # Calculate the new length of the resampled waveform\n", " resampled_length = int(waveform.shape[-1] * resampling_ratio)\n", " # Resample the waveform using SciPy's resample function\n", " resampled_waveform = resample(waveform, resampled_length, axis=-1)\n", " \n", " return resampled_waveform\n", "\n", "def sort_channels_to_ZNE(waveform, channels):\n", " # Input:\n", " # waveform: a 2D numpy array with shape (3, n), where n is the number of samples\n", " # channels: a list or tuple of 3 strings representing the channel order, e.g. ('N', 'Z', 'E')\n", " channels = list(channels)\n", "\n", " if len(channels) != 3 or set(channels) != {'Z', 'N', 'E'}:\n", " raise ValueError(\"Invalid channel input. It should be a permutation of 'Z', 'N', and 'E'.\")\n", "\n", " # Find the indices of the Z, N, and E channels\n", " z_index = channels.index('Z')\n", " n_index = channels.index('N')\n", " e_index = channels.index('E')\n", " \n", " print(z_index, n_index, e_index)\n", " # Sort the channels to ZNE\n", " sorted_waveform = waveform[[z_index, n_index, e_index], :]\n", " \n", " return sorted_waveform\n", "\n", "def make_prediction(waveform, sampling_rate, order):\n", " waveform = np.load(waveform)\n", " print('Loaded', waveform.shape)\n", "\n", " if len(waveform.shape) == 1:\n", " waveform = waveform.reshape(1, waveform.shape[0])\n", "\n", " elif waveform.shape[0] == 3:\n", " waveform = sort_channels_to_ZNE(waveform, order)\n", "\n", " if sampling_rate != 100:\n", " waveform = resample_waveform(waveform, sampling_rate, 100)\n", " print('Resampled', waveform.shape)\n", "\n", "\n", " orig_waveform = waveform[:, :6000].copy()\n", " processed_input = prepare_waveform(waveform)\n", "\n", " # Make prediction\n", " with torch.inference_mode():\n", " output = model(processed_input)\n", "\n", " p_phase = output[:, 0]\n", " s_phase = output[:, 1]\n", "\n", " return processed_input, p_phase, s_phase, orig_waveform\n", "\n", "\n", "def mark_phases(waveform, uploaded_file, p_thres, s_thres, sampling_rate, order):\n", "\n", " if uploaded_file is not None:\n", " waveform = uploaded_file.name\n", "\n", " processed_input, p_phase, s_phase, orig_waveform = make_prediction(waveform, sampling_rate, order)\n", "\n", " # Create a plot of the waveform with the phases marked\n", " if sum(processed_input[0][2] == 0): #if input is 1C\n", " fig, ax = plt.subplots(nrows=2, figsize=(10, 2), sharex=True)\n", "\n", " ax[0].plot(orig_waveform[0], color='black', lw=1)\n", " ax[0].set_ylabel('Norm. Ampl.')\n", "\n", " else: #if input is 3C\n", " fig, ax = plt.subplots(nrows=4, figsize=(10, 6), sharex=True)\n", " ax[0].plot(orig_waveform[0], color='black', lw=1)\n", " ax[1].plot(orig_waveform[1], color='black', lw=1)\n", " ax[2].plot(orig_waveform[2], color='black', lw=1)\n", "\n", " ax[0].set_ylabel('Z')\n", " ax[1].set_ylabel('N')\n", " ax[2].set_ylabel('E')\n", "\n", "\n", " do_we_have_p = (p_phase.std().item()*60 < p_thres)\n", " if do_we_have_p:\n", " p_phase_plot = p_phase*processed_input.shape[-1]\n", " p_kde = gaussian_kde(p_phase_plot)\n", " p_dist_space = np.linspace( min(p_phase_plot)-10, max(p_phase_plot)+10, 500 )\n", " ax[-1].plot( p_dist_space, p_kde(p_dist_space), color='r')\n", " else:\n", " ax[-1].text(0.5, 0.75, 'No P phase detected', horizontalalignment='center', verticalalignment='center', transform=ax[-1].transAxes)\n", "\n", " do_we_have_s = (s_phase.std().item()*60 < s_thres)\n", " if do_we_have_s:\n", " s_phase_plot = s_phase*processed_input.shape[-1]\n", " s_kde = gaussian_kde(s_phase_plot)\n", " s_dist_space = np.linspace( min(s_phase_plot)-10, max(s_phase_plot)+10, 500 )\n", " ax[-1].plot( s_dist_space, s_kde(s_dist_space), color='b')\n", "\n", " for a in ax:\n", " a.axvline(p_phase.mean()*processed_input.shape[-1], color='r', linestyle='--', label='P', alpha=do_we_have_p)\n", " a.axvline(s_phase.mean()*processed_input.shape[-1], color='b', linestyle='--', label='S', alpha=do_we_have_s)\n", " else:\n", " ax[-1].text(0.5, 0.25, 'No S phase detected', horizontalalignment='center', verticalalignment='center', transform=ax[-1].transAxes)\n", "\n", " ax[-1].set_xlabel('Time, samples')\n", " ax[-1].set_ylabel('Uncert., samples')\n", " ax[-1].legend()\n", "\n", " plt.subplots_adjust(hspace=0., wspace=0.)\n", "\n", " # Convert the plot to an image and return it\n", " fig.canvas.draw()\n", " image = np.array(fig.canvas.renderer.buffer_rgba())\n", " plt.close(fig)\n", " return image\n", "\n", "def bin_distances(distances, bin_size=10):\n", " # Bin the distances into groups of `bin_size` kilometers\n", " binned_distances = {}\n", " for i, distance in enumerate(distances):\n", " bin_index = distance // bin_size\n", " if bin_index not in binned_distances:\n", " binned_distances[bin_index] = (distance, i)\n", " elif i < binned_distances[bin_index][1]:\n", " binned_distances[bin_index] = (distance, i)\n", "\n", " # Select the first distance in each bin and its index\n", " first_distances = []\n", " for bin_index in binned_distances:\n", " first_distance, first_distance_index = binned_distances[bin_index]\n", " first_distances.append(first_distance_index)\n", " \n", " return first_distances\n", "\n", "def variance_coefficient(residuals):\n", " # calculate the variance of the residuals\n", " var = residuals.var()\n", " # scale the variance to a coefficient between 0 and 1\n", " coeff = 1 - (var / (residuals.max() - residuals.min()))\n", " return coeff\n", "\n", "def predict_on_section(client_name, timestamp, eq_lat, eq_lon, radius_km, source_depth_km, velocity_model, max_waveforms, conf_thres_P, conf_thres_S):\n", " distances, t0s, st_lats, st_lons, waveforms, names = [], [], [], [], [], []\n", " \n", " taup_model = TauPyModel(model=velocity_model)\n", " client = Client(client_name)\n", "\n", " window = radius_km / 111.2\n", " max_waveforms = int(max_waveforms)\n", "\n", " assert eq_lat - window > -90 and eq_lat + window < 90, \"Latitude out of bounds\"\n", " assert eq_lon - window > -180 and eq_lon + window < 180, \"Longitude out of bounds\"\n", "\n", " starttime = obspy.UTCDateTime(timestamp)\n", " endtime = starttime + 120\n", "\n", " try:\n", " print('Starting to download inventory')\n", " inv = client.get_stations(network=\"*\", station=\"*\", location=\"*\", channel=\"*H*\", \n", " starttime=starttime, endtime=endtime, \n", " minlatitude=(eq_lat-window), maxlatitude=(eq_lat+window),\n", " minlongitude=(eq_lon-window), maxlongitude=(eq_lon+window), \n", " level='station')\n", " print('Finished downloading inventory')\n", " \n", " except (IndexError, FDSNNoDataException, FDSNTimeoutException, FDSNInternalServerException):\n", " fig, ax = plt.subplots()\n", " ax.text(0.5,0.5,'Something is wrong with the data provider, try another')\n", " fig.canvas.draw();\n", " image = np.array(fig.canvas.renderer.buffer_rgba())\n", " plt.close(fig)\n", " return image\n", " \n", " waveforms = []\n", " cached_waveforms = glob(\"data/cached/*.mseed\")\n", "\n", " for network in inv:\n", " if network.code == 'SY':\n", " continue\n", " for station in network:\n", " print(f\"Processing {network.code}.{station.code}...\")\n", " distance = locations2degrees(eq_lat, eq_lon, station.latitude, station.longitude)\n", "\n", " arrivals = taup_model.get_travel_times(source_depth_in_km=source_depth_km, \n", " distance_in_degree=distance, \n", " phase_list=[\"P\", \"S\"])\n", "\n", " if len(arrivals) > 0:\n", "\n", " starttime = obspy.UTCDateTime(timestamp) + arrivals[0].time - 15\n", " endtime = starttime + 60\n", " try:\n", " filename=f'{network.code}_{station.code}_{starttime}'\n", " if f\"data/cached/{filename}.mseed\" not in cached_waveforms:\n", " print(f'Downloading waveform for {filename}')\n", " waveform = client.get_waveforms(network=network.code, station=station.code, location=\"*\", channel=\"*\", \n", " starttime=starttime, endtime=endtime)\n", " waveform.write(f\"data/cached/{network.code}_{station.code}_{starttime}.mseed\", format=\"MSEED\")\n", " print('Finished downloading and caching waveform')\n", " else:\n", " print('Reading cached waveform')\n", " waveform = obspy.read(f\"data/cached/{network.code}_{station.code}_{starttime}.mseed\")\n", " \n", "\n", " except (IndexError, FDSNNoDataException, FDSNTimeoutException, FDSNInternalServerException):\n", " print(f'Skipping {network.code}_{station.code}_{starttime}')\n", " continue\n", " \n", " waveform = waveform.select(channel=\"H[BH][ZNE]\")\n", " waveform = waveform.merge(fill_value=0)\n", " waveform = waveform[:3].sort(keys=['channel'], reverse=True)\n", "\n", " len_check = [len(x.data) for x in waveform]\n", " if len(set(len_check)) > 1:\n", " continue\n", "\n", " if len(waveform) == 3:\n", " try:\n", " waveform = prepare_waveform(np.stack([x.data for x in waveform]))\n", "\n", " distances.append(distance)\n", " t0s.append(starttime)\n", " st_lats.append(station.latitude)\n", " st_lons.append(station.longitude)\n", " waveforms.append(waveform)\n", " names.append(f\"{network.code}.{station.code}\")\n", "\n", " print(f\"Added {network.code}.{station.code} to the list of waveforms\")\n", "\n", " except:\n", " continue\n", " \n", " \n", " # If there are no waveforms, return an empty plot\n", " if len(waveforms) == 0:\n", " print('No waveforms found')\n", " fig, ax = plt.subplots()\n", " # prints \"No waveforms found\" on the plot aligned at center and vertically\n", " ax.text(0.5,0.5,'No waveforms found', horizontalalignment='center', verticalalignment='center', transform=ax.transAxes)\n", " fig.canvas.draw();\n", " image = np.array(fig.canvas.renderer.buffer_rgba())\n", " plt.close(fig)\n", "\n", " output_picks = pd.DataFrame()\n", " output_picks.to_csv('data/picks.csv', index=False)\n", " output_csv = 'data/picks.csv'\n", " return image, output_picks, output_csv\n", " \n", "\n", " first_distances = bin_distances(distances, bin_size=10/111.2)\n", "\n", " # Edge case when there are way too many waveforms to process\n", " selection_indexes = np.random.choice(first_distances, \n", " np.min([len(first_distances), max_waveforms]),\n", " replace=False)\n", "\n", " waveforms = np.array(waveforms)[selection_indexes]\n", " distances = np.array(distances)[selection_indexes]\n", " t0s = np.array(t0s)[selection_indexes]\n", " st_lats = np.array(st_lats)[selection_indexes]\n", " st_lons = np.array(st_lons)[selection_indexes]\n", " names = np.array(names)[selection_indexes]\n", "\n", " waveforms = [torch.tensor(waveform) for waveform in waveforms]\n", "\n", " print('Starting to run predictions')\n", " with torch.no_grad():\n", " waveforms_torch = torch.vstack(waveforms)\n", " output = model(waveforms_torch)\n", "\n", " p_phases = output[:, 0]\n", " s_phases = output[:, 1]\n", "\n", " p_phases = p_phases.reshape(len(waveforms),-1)\n", " s_phases = s_phases.reshape(len(waveforms),-1)\n", "\n", " # Max confidence - min variance \n", " p_max_confidence = p_phases.std(axis=-1).min()\n", " s_max_confidence = s_phases.std(axis=-1).min()\n", "\n", " print(f\"Starting plotting {len(waveforms)} waveforms\")\n", " fig, ax = plt.subplots(ncols=3, figsize=(10, 3))\n", " \n", " # Plot topography\n", " print('Fetching topography')\n", " params = Topography.DEFAULT.copy()\n", " extra_window = 0.5\n", " params[\"south\"] = np.min([st_lats.min(), eq_lat])-extra_window\n", " params[\"north\"] = np.max([st_lats.max(), eq_lat])+extra_window\n", " params[\"west\"] = np.min([st_lons.min(), eq_lon])-extra_window\n", " params[\"east\"] = np.max([st_lons.max(), eq_lon])+extra_window\n", "\n", " topo_map = Topography(**params)\n", " topo_map.fetch()\n", " topo_map.load()\n", "\n", " print('Plotting topo')\n", " hillshade = es.hillshade(topo_map.da[0], altitude=10)\n", " \n", " topo_map.da.plot(ax = ax[1], cmap='Greys', add_colorbar=False, add_labels=False)\n", " topo_map.da.plot(ax = ax[2], cmap='Greys', add_colorbar=False, add_labels=False)\n", " ax[1].imshow(hillshade, cmap=\"Greys\", alpha=0.5)\n", "\n", " output_picks = pd.DataFrame({'station_name' : [], \n", " 'st_lat' : [], 'st_lon' : [],\n", " 'starttime' : [], \n", " 'p_phase, s' : [], 'p_uncertainty, s' : [], \n", " 's_phase, s' : [], 's_uncertainty, s' : [],\n", " 'velocity_p, km/s' : [], 'velocity_s, km/s' : []})\n", " \n", " for i in range(len(waveforms)):\n", " print(f\"Plotting waveform {i+1}/{len(waveforms)}\")\n", " current_P = p_phases[i]\n", " current_S = s_phases[i]\n", " \n", " x = [t0s[i] + pd.Timedelta(seconds=k/100) for k in np.linspace(0,6000,6000)]\n", " x = mdates.date2num(x)\n", "\n", " # Normalize confidence for the plot\n", " p_conf = 1/(current_P.std()/p_max_confidence).item()\n", " s_conf = 1/(current_S.std()/s_max_confidence).item()\n", "\n", " delta_t = t0s[i].timestamp - obspy.UTCDateTime(timestamp).timestamp\n", "\n", " ax[0].plot(x, waveforms[i][0, 0]*10+distances[i]*111.2, color='black', alpha=0.5, lw=1)\n", "\n", " if (current_P.std().item()*60 < conf_thres_P) or (current_S.std().item()*60 < conf_thres_S):\n", " ax[0].scatter(x[int(current_P.mean()*waveforms[i][0].shape[-1])], waveforms[i][0, 0].mean()+distances[i]*111.2, color='r', alpha=p_conf, marker='|')\n", " ax[0].scatter(x[int(current_S.mean()*waveforms[i][0].shape[-1])], waveforms[i][0, 0].mean()+distances[i]*111.2, color='b', alpha=s_conf, marker='|')\n", " \n", " velocity_p = (distances[i]*111.2)/(delta_t+current_P.mean()*60).item()\n", " velocity_s = (distances[i]*111.2)/(delta_t+current_S.mean()*60).item()\n", "\n", " # Generate an array from st_lat to eq_lat and from st_lon to eq_lon\n", " x = np.linspace(st_lons[i], eq_lon, 50)\n", " y = np.linspace(st_lats[i], eq_lat, 50)\n", " \n", " # Plot the array\n", " ax[1].scatter(x, y, c=np.zeros_like(x)+velocity_p, alpha=0.1, vmin=0, vmax=8)\n", " ax[2].scatter(x, y, c=np.zeros_like(x)+velocity_s, alpha=0.1, vmin=0, vmax=8)\n", "\n", " else:\n", " velocity_p = np.nan\n", " velocity_s = np.nan\n", " \n", " ax[0].set_ylabel('Z')\n", " print(f\"Station {st_lats[i]}, {st_lons[i]} has P velocity {velocity_p} and S velocity {velocity_s}\")\n", "\n", " output_picks = output_picks.append(pd.DataFrame({'station_name': [names[i]], \n", " 'st_lat' : [st_lats[i]], 'st_lon' : [st_lons[i]],\n", " 'starttime' : [str(t0s[i])], \n", " 'p_phase, s' : [(delta_t+current_P.mean()*60).item()], 'p_uncertainty, s' : [current_P.std().item()*60], \n", " 's_phase, s' : [(delta_t+current_S.mean()*60).item()], 's_uncertainty, s' : [current_S.std().item()*60],\n", " 'velocity_p, km/s' : [velocity_p], 'velocity_s, km/s' : [velocity_s]}))\n", " \n", " \n", " # Add legend\n", " ax[0].scatter(None, None, color='r', marker='|', label='P')\n", " ax[0].scatter(None, None, color='b', marker='|', label='S')\n", " ax[0].xaxis.set_major_formatter(mdates.DateFormatter('%H:%M:%S'))\n", " ax[0].xaxis.set_major_locator(mdates.SecondLocator(interval=20))\n", " ax[0].legend()\n", "\n", " print('Plotting stations')\n", " for i in range(1,3):\n", " ax[i].scatter(st_lons, st_lats, color='b', label='Stations')\n", " ax[i].scatter(eq_lon, eq_lat, color='r', marker='*', label='Earthquake')\n", " ax[i].set_aspect('equal')\n", " ax[i].set_xticklabels(ax[i].get_xticks(), rotation = 50)\n", "\n", " fig.subplots_adjust(bottom=0.1, top=0.9, left=0.1, right=0.8,\n", " wspace=0.02, hspace=0.02)\n", " \n", " cb_ax = fig.add_axes([0.83, 0.1, 0.02, 0.8])\n", " cbar = fig.colorbar(ax[2].scatter(None, None, c=velocity_p, alpha=0.5, vmin=0, vmax=8), cax=cb_ax)\n", "\n", " cbar.set_label('Velocity (km/s)')\n", " ax[1].set_title('P Velocity')\n", " ax[2].set_title('S Velocity')\n", "\n", " for a in ax:\n", " a.tick_params(axis='both', which='major', labelsize=8)\n", " \n", " plt.subplots_adjust(hspace=0., wspace=0.5)\n", " fig.canvas.draw();\n", " image = np.array(fig.canvas.renderer.buffer_rgba())\n", " plt.close(fig)\n", "\n", " output_csv = f'data/velocity/{eq_lat}_{eq_lon}_{source_depth_km}_{timestamp}_{len(waveforms)}.csv'\n", " output_picks.to_csv(output_csv, index=False)\n", " \n", " return image, output_picks, output_csv\n", "\n", "import numpy as np\n", "from matplotlib import colors, cm\n", "from scipy.interpolate import griddata\n", "\n", "def interpolate_vel_model(velocity_model, initial_velocity, lat_values, lon_values, depth_values, n_lat, n_lon, n_depth):\n", " # Create a mask for points with the initial velocity\n", " initial_velocity_mask = (velocity_model == initial_velocity)\n", "\n", " # Find the indices of points with non-initial velocities\n", " non_initial_velocity_indices = np.argwhere(~initial_velocity_mask)\n", "\n", " # Extract the coordinates and corresponding velocities of the known points\n", " known_points = np.column_stack([lat_values[non_initial_velocity_indices[:, 0]],\n", " lon_values[non_initial_velocity_indices[:, 1]],\n", " depth_values[non_initial_velocity_indices[:, 2]]])\n", " \n", " # Find the maximum depth in the known_points\n", " max_known_depth = np.max(known_points[:, 2])\n", "\n", " known_velocities = velocity_model[~initial_velocity_mask]\n", "\n", " # Create a grid of points for the entire volume\n", " grid_points = np.array(np.meshgrid(lat_values, lon_values, depth_values, indexing='ij')).reshape(3, -1).T\n", "\n", " # Create a mask for grid points that are deeper than the maximum known depth\n", " depth_mask = grid_points[:, 2] <= max_known_depth\n", "\n", " # Interpolate the velocities at the grid points\n", " interpolated_velocities = griddata(known_points, known_velocities, grid_points[depth_mask], method='linear')\n", "\n", " # Fill nan values with the nearest known velocities\n", " interpolated_velocities_filled = griddata(known_points, known_velocities, grid_points[depth_mask], method='nearest')\n", " interpolated_velocities[np.isnan(interpolated_velocities)] = interpolated_velocities_filled[np.isnan(interpolated_velocities)]\n", "\n", " # Initialize an array with the same length as grid_points and fill it with nan values\n", " interpolated_velocities_with_depth_limit = np.full(grid_points.shape[0], np.nan)\n", "\n", " # Update the array with the interpolated velocities for the masked grid points\n", " interpolated_velocities_with_depth_limit[depth_mask] = interpolated_velocities\n", "\n", " # Reshape the interpolated velocities to match the shape of the velocity_model\n", " interpolated_velocity_model = interpolated_velocities_with_depth_limit.reshape(n_lat, n_lon, n_depth)\n", "\n", " return interpolated_velocity_model\n", "\n", "\n", "# Function to find the closest index for a given value in an array\n", "def find_closest_index(array, value):\n", " return np.argmin(np.abs(array - value))\n", "\n", "def compute_velocity_model(azimuth, elevation, interpolate, n_lat, n_lon, n_depth):\n", " filename = list(output_csv.temp_files)[0]\n", " \n", " df = pd.read_csv(filename)\n", " filename = filename.split('/')[-1]\n", " \n", " # Current EQ location\n", " eq_lat = float(filename.split(\"_\")[0])\n", " eq_lon = float(filename.split(\"_\")[1])\n", " eq_depth = float(filename.split(\"_\")[2])\n", "\n", " # Define the region of interest (latitude, longitude, and depth ranges)\n", " lat_range = (np.min([df.st_lat.min(), eq_lat]), np.max([df.st_lat.max(), eq_lat]))\n", " lon_range = (np.min([df.st_lon.min(), eq_lon]), np.max([df.st_lon.max(), eq_lon]))\n", " depth_range = (0, 50)\n", "\n", " # Define the number of nodes in each dimension\n", " num_points = 100\n", "\n", " taup_model = TauPyModel(model='1066a')\n", "\n", " # Create the grid\n", " lat_values = np.linspace(lat_range[0], lat_range[1], n_lat)\n", " lon_values = np.linspace(lon_range[0], lon_range[1], n_lon)\n", " depth_values = np.linspace(depth_range[0], depth_range[1], n_depth)\n", "\n", " # Initialize the velocity model with constant values\n", " initial_velocity = 0 # km/s, this can be P-wave or S-wave velocity\n", " velocity_model = np.full((n_lat, n_lon, n_depth), initial_velocity, dtype=float)\n", "\n", " # Loop through the stations and update the velocity model\n", " for i in range(len(df)):\n", " if ~np.isnan(df['velocity_p, km/s'].iloc[i]):\n", "\n", " ray_path = taup_model.get_ray_paths_geo(source_depth_in_km=eq_depth,\n", " source_latitude_in_deg=eq_lat,\n", " source_longitude_in_deg=eq_lon,\n", " receiver_latitude_in_deg=df.st_lat.iloc[i],\n", " receiver_longitude_in_deg=df.st_lon.iloc[i],\n", " phase_list=['P', 'S'])\n", "\n", " # Create the interpolator objects for latitude, longitude, and depth\n", " interp_latitude = interp1d(np.linspace(0, ray_path[0].path['lat'].max(), len(ray_path[0].path['lat'])), ray_path[0].path['lat'])\n", " interp_longitude = interp1d(np.linspace(0, ray_path[0].path['lon'].max(), len(ray_path[0].path['lon'])), ray_path[0].path['lon'])\n", " interp_depth = interp1d(np.linspace(0, ray_path[0].path['depth'].max(), len(ray_path[0].path['depth'])), ray_path[0].path['depth'])\n", "\n", " # Resample the ray path to N points\n", " lat_values_interp = interp_latitude(np.linspace(0, ray_path[0].path['lat'].max(), num_points))\n", " lon_values_interp = interp_longitude(np.linspace(0, ray_path[0].path['lon'].max(), num_points))\n", " depth_values_interp = interp_depth(np.linspace(0, ray_path[0].path['depth'].max(), num_points))\n", "\n", " # Loop through the interpolated coordinates and update the grid cells with the average P-wave velocity\n", " for lat, lon, depth in zip(lat_values_interp, lon_values_interp, depth_values_interp):\n", " lat_index = find_closest_index(lat_values, lat)\n", " lon_index = find_closest_index(lon_values, lon)\n", " depth_index = find_closest_index(depth_values, depth)\n", " \n", " if velocity_model[lat_index, lon_index, depth_index] == initial_velocity:\n", " velocity_model[lat_index, lon_index, depth_index] = df['velocity_p, km/s'].iloc[i]\n", " else:\n", " velocity_model[lat_index, lon_index, depth_index] = (velocity_model[lat_index, lon_index, depth_index] +\n", " df['velocity_p, km/s'].iloc[i]) / 2\n", "\n", " # Create the figure and axis\n", " fig = plt.figure(figsize=(8, 8))\n", " ax = fig.add_subplot(111, projection='3d')\n", "\n", " # Set the plot limits\n", " ax.set_xlim3d(lat_range[0], lat_range[1])\n", " ax.set_ylim3d(lon_range[0], lon_range[1])\n", " ax.set_zlim3d(depth_range[1], depth_range[0])\n", "\n", " ax.set_xlabel('Latitude')\n", " ax.set_ylabel('Longitude')\n", " ax.set_zlabel('Depth (km)')\n", " ax.set_title('Velocity Model')\n", " \n", " # Create the meshgrid\n", " x, y, z = np.meshgrid(\n", " np.linspace(lat_range[0], lat_range[1], velocity_model.shape[0]+1),\n", " np.linspace(lon_range[0], lon_range[1], velocity_model.shape[1]+1),\n", " np.linspace(depth_range[0], depth_range[1], velocity_model.shape[2]+1),\n", " indexing='ij'\n", " )\n", "\n", " # Create the color array\n", " norm = plt.Normalize(vmin=2, vmax=8)\n", " colors_vel = plt.cm.plasma(norm(velocity_model)) \n", " \n", " # Plot the voxels\n", " if interpolate:\n", " interpolated_velocity_model = interpolate_vel_model(velocity_model, initial_velocity, lat_values, lon_values, depth_values, n_lat, n_lon, n_depth)\n", " colors_interp = plt.cm.plasma(norm(interpolated_velocity_model))\n", " ax.voxels(x, y, z, interpolated_velocity_model > 0, facecolors=colors_interp, alpha=0.5, edgecolor='k')\n", " \n", " ax.voxels(x, y, z, velocity_model > 0, facecolors=colors_vel, alpha=1, edgecolor='black')\n", "\n", " # Set the view angle\n", " ax.view_init(elev=elevation, azim=azimuth)\n", "\n", " m = cm.ScalarMappable(cmap=plt.cm.plasma, norm=norm)\n", " m.set_array([])\n", " plt.colorbar(m)\n", "\n", " # Show the plot\n", " fig.canvas.draw();\n", " image = np.array(fig.canvas.renderer.buffer_rgba())\n", " plt.close(fig)\n", "\n", " return image\n", "\n", "# model = torch.jit.load(\"model.pt\")\n", "model = torch.jit.load(\"model.pt\")\n", "\n", "model.eval()\n", "\n", "with gr.Blocks() as demo:\n", " gr.HTML(\"\"\"\n", "
\n", "

PhaseHunter 🏹\n", "\n", "

\n", " \n", "

Detect P and S seismic phases with uncertainty

\n", " \n", "
\n", "\"\"\")\n", " with gr.Tab(\"Try on a single station\"):\n", " with gr.Row(): \n", " # Define the input and output types for Gradio\n", " inputs = gr.Dropdown(\n", " [\"data/sample/sample_0.npy\", \n", " \"data/sample/sample_1.npy\", \n", " \"data/sample/sample_2.npy\"], \n", " label=\"Sample waveform\", \n", " info=\"Select one of the samples\",\n", " value = \"data/sample/sample_0.npy\"\n", " )\n", " with gr.Column(scale=1):\n", " P_thres_inputs = gr.Slider(minimum=0.01,\n", " maximum=1,\n", " value=0.1,\n", " label=\"P uncertainty threshold (s)\",\n", " step=0.01,\n", " info=\"Acceptable uncertainty for P picks expressed in std() seconds\",\n", " interactive=True,\n", " )\n", " \n", " S_thres_inputs = gr.Slider(minimum=0.01,\n", " maximum=1,\n", " value=0.2,\n", " label=\"S uncertainty threshold (s)\",\n", " step=0.01,\n", " info=\"Acceptable uncertainty for S picks expressed in std() seconds\",\n", " interactive=True,\n", " )\n", " with gr.Column(scale=1):\n", " upload = gr.File(label=\"Upload your waveform\")\n", " with gr.Row():\n", " sampling_rate_inputs = gr.Slider(minimum=10,\n", " maximum=1000,\n", " value=100,\n", " label=\"Samlping rate, Hz\",\n", " step=10,\n", " info=\"Sampling rate of the waveform\",\n", " interactive=True,\n", " )\n", " order_input = gr.Text(value='ZNE', \n", " label='Channel order', \n", " info='Order of the channels in the waveform file (e.g. ZNE)')\n", "\n", " button = gr.Button(\"Predict phases\")\n", " outputs = gr.Image(label='Waveform with Phases Marked', type='numpy', interactive=False)\n", " \n", " button.click(mark_phases, inputs=[inputs, upload, \n", " P_thres_inputs, S_thres_inputs,\n", " sampling_rate_inputs, order_input], \n", " outputs=outputs) \n", " with gr.Tab(\"Select earthquake from catalogue\"):\n", "\n", " gr.HTML(\"\"\"\n", "
\n", "

Using PhaseHunter to Analyze Seismic Waveforms

\n", "

Select an earthquake from the global earthquake catalogue (e.g. USGS) and the app will download the waveform from the FDSN client of your choice. The app will use a velocity model of your choice to select appropriate time windows for each station within a specified radius of the earthquake.

\n", "

The app will then analyze the waveforms and mark the detected phases on the waveform. Pick data for each waveform is reported in seconds from the start of the waveform.

\n", "

Velocities are derived from distance and travel time determined by PhaseHunter picks (v = distance/predicted_pick_time). The background of the velocity plot is colored by DEM.

\n", "
\n", " \"\"\")\n", " with gr.Row(): \n", " with gr.Column(scale=2):\n", " client_inputs = gr.Dropdown(\n", " choices = list(URL_MAPPINGS.keys()), \n", " label=\"FDSN Client\", \n", " info=\"Select one of the available FDSN clients\",\n", " value = \"IRIS\",\n", " interactive=True\n", " )\n", "\n", " velocity_inputs = gr.Dropdown(\n", " choices = ['1066a', '1066b', 'ak135', \n", " 'ak135f', 'herrin', 'iasp91', \n", " 'jb', 'prem', 'pwdk'], \n", " label=\"1D velocity model\", \n", " info=\"Velocity model for station selection\",\n", " value = \"1066a\",\n", " interactive=True\n", " )\n", "\n", " with gr.Column(scale=2):\n", " timestamp_inputs = gr.Textbox(value='2019-07-04T17:33:49-00',\n", " placeholder='YYYY-MM-DDTHH:MM:SS-TZ',\n", " label=\"Timestamp\",\n", " info=\"Timestamp of the earthquake\",\n", " max_lines=1,\n", " interactive=True)\n", " \n", " source_depth_inputs = gr.Number(value=10,\n", " label=\"Source depth (km)\",\n", " info=\"Depth of the earthquake\",\n", " interactive=True)\n", " \n", " with gr.Column(scale=2):\n", " eq_lat_inputs = gr.Number(value=35.766, \n", " label=\"Latitude\", \n", " info=\"Latitude of the earthquake\",\n", " interactive=True)\n", " \n", " eq_lon_inputs = gr.Number(value=-117.605,\n", " label=\"Longitude\",\n", " info=\"Longitude of the earthquake\",\n", " interactive=True)\n", " \n", " with gr.Column(scale=2):\n", " radius_inputs = gr.Slider(minimum=1, \n", " maximum=200, \n", " value=50, \n", " label=\"Radius (km)\", \n", " step=10,\n", " info=\"\"\"Select the radius around the earthquake to download data from.\\n \n", " Note that the larger the radius, the longer the app will take to run.\"\"\",\n", " interactive=True)\n", " \n", " max_waveforms_inputs = gr.Slider(minimum=1,\n", " maximum=100,\n", " value=10,\n", " label=\"Max waveforms per section\",\n", " step=1,\n", " info=\"Maximum number of waveforms to show per section\\n (to avoid long prediction times)\",\n", " interactive=True,\n", " )\n", " with gr.Column(scale=2):\n", " P_thres_inputs = gr.Slider(minimum=0.01,\n", " maximum=1,\n", " value=0.1,\n", " label=\"P uncertainty threshold, s\",\n", " step=0.01,\n", " info=\"Acceptable uncertainty for P picks expressed in std() seconds\",\n", " interactive=True,\n", " )\n", " S_thres_inputs = gr.Slider(minimum=0.01,\n", " maximum=1,\n", " value=0.2,\n", " label=\"S uncertainty threshold, s\",\n", " step=0.01,\n", " info=\"Acceptable uncertainty for S picks expressed in std() seconds\",\n", " interactive=True,\n", " )\n", " \n", " button = gr.Button(\"Predict phases\")\n", " output_image = gr.Image(label='Waveforms with Phases Marked', type='numpy', interactive=False)\n", " \n", " with gr.Row():\n", " with gr.Column(scale=2):\n", " azimuth_input = gr.Slider(minimum=-180, maximum=180, value=0, step=5, label=\"Azimuth\", interactive=True)\n", " elevation_input = gr.Slider(minimum=-90, maximum=90, value=30, step=5, label=\"Elevation\", interactive=True)\n", "\n", " with gr.Row():\n", " interpolate_input = gr.Checkbox(label=\"Interpolate\", info=\"Interpolate velocity model\")\n", " n_lat_input = gr.Slider(minimum=5, maximum=100, value=50, step=5, label=\"N lat\", info='Number of Lat grid points', interactive=True)\n", " n_lon_input = gr.Slider(minimum=5, maximum=100, value=50, step=5, label=\"N lon\", info='Number of Lon grid points', interactive=True)\n", " n_depth_input = gr.Slider(minimum=5, maximum=100, value=50, step=5, label=\"N depth\", info='Number of Depth grid points', interactive=True)\n", " \n", " button = gr.Button(\"Look at 3D Velocities\")\n", " outputs_vel_model = gr.Image(label=\"3D Velocity Model\")\n", "\n", " button.click(compute_velocity_model, \n", " inputs=[azimuth_input, elevation_input, \n", " interpolate_input, n_lat_input, \n", " n_lon_input, n_depth_input], \n", " outputs=[outputs_vel_model])\n", " \n", " with gr.Row():\n", " output_picks = gr.Dataframe(label='Pick data', \n", " type='pandas', \n", " interactive=False)\n", " output_csv = gr.File(label=\"Output File\", file_types=[\".csv\"])\n", "\n", " button.click(predict_on_section, \n", " inputs=[client_inputs, timestamp_inputs, \n", " eq_lat_inputs, eq_lon_inputs, \n", " radius_inputs, source_depth_inputs, \n", " velocity_inputs, max_waveforms_inputs,\n", " P_thres_inputs, S_thres_inputs],\n", " outputs=[output_image, output_picks, output_csv])\n", "\n", "demo.launch()" ] }, { "cell_type": "code", "execution_count": 89, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "['__class__',\n", " '__delattr__',\n", " '__dict__',\n", " '__dir__',\n", " '__doc__',\n", " '__eq__',\n", " '__format__',\n", " '__ge__',\n", " '__getattribute__',\n", " '__getstate__',\n", " '__gt__',\n", " '__hash__',\n", " '__init__',\n", " '__init_subclass__',\n", " '__le__',\n", " '__lt__',\n", " '__module__',\n", " '__ne__',\n", " '__new__',\n", " '__reduce__',\n", " '__reduce_ex__',\n", " '__repr__',\n", " '__setattr__',\n", " '__sizeof__',\n", " '__str__',\n", " '__subclasshook__',\n", " '__weakref__',\n", " 'azimuth',\n", " 'distance',\n", " 'incident_angle',\n", " 'name',\n", " 'path',\n", " 'phase',\n", " 'pierce',\n", " 'purist_dist',\n", " 'purist_distance',\n", " 'purist_name',\n", " 'ray_param',\n", " 'ray_param_index',\n", " 'ray_param_sec_degree',\n", " 'receiver_depth',\n", " 'source_depth',\n", " 'takeoff_angle',\n", " 'time']" ] }, "execution_count": 89, "metadata": {}, "output_type": "execute_result" } ], "source": [ "dir(ray_path[0])" ] }, { "cell_type": "code", "execution_count": 107, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 107, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlIAAAGdCAYAAADZiZ2PAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABeE0lEQVR4nO3deZhV1Z3v//fa+8w1MhRQyqDEIGJCyqYVISHodUAS02oSpLWNEDC5pk2rrd5u0PwM5mlDbtvaxr7dMWlAzA1JxCQYHGJrTKk3kZh2KIwDBLRMoVLIVOfUqTrzXr8/TnGkrIGqglPj5/U8+0nO3muvvXZtq86Xtdb+LmOttYiIiIhIrzkD3QARERGRoUqBlIiIiEgfKZASERER6SMFUiIiIiJ9pEBKREREpI8USImIiIj0kQIpERERkT5SICUiIiLSR76BbsBw53ke7733HmVlZRhjBro5IiIi0gPWWpqbmznuuONwnK77nRRIFdl7773HpEmTBroZIiIi0ge7du1i4sSJXR5XIFVkZWVlQP5BlJeXD3BrREREpCdisRiTJk0qfI93RYFUkR0azisvL1cgJSIiMsQcaVqOJpuLiIiI9FHRAqnbb7+duXPnEolEqKys7LTMtddey6xZswgGg9TU1HQ4vmrVKowxHbaSkpIjXn/9+vXMnDmTUCjEuHHjuOaaa9odf+WVV5g3bx6hUIhJkybxz//8zx3qePDBB5k+fTqhUIiPf/zjPPbYYz26dxERERkZihZIpdNpFi1axNe+9rVuyy1btozFixd3euymm25i9+7d7bYZM2awaNGibuu86667uOWWW1ixYgWvvfYav/71r1mwYEHheCwW4/zzz2fKlCm8+OKL3HHHHaxatYof/OAHhTLPPfccl112GcuXL+fll1/m4osv5uKLL+bVV1/txU9BREREhjNjrbXFvMD69eu5/vrraWpq6rLMqlWreOihh6irq+u2rq1bt1JTU8Ozzz7LvHnzOi1z8OBBjj/+eB5++GHOOeecTst873vf45ZbbqGxsZFAIADAihUreOihh9i2bRsAixcvpqWlhUceeaRw3plnnklNTQ333ntvt+08XCwWo6Kigmg0qjlSIiIiQ0RPv7+H1BypNWvWMG3atC6DKIAnn3wSz/N49913OeWUU5g4cSKXXnopu3btKpTZsmULn/70pwtBFMCCBQvYvn07Bw8eLJQ599xz29W9YMECtmzZ0m0bU6kUsVis3SYiIiLD05AJpJLJJBs2bGD58uXdlnvrrbfwPI9vf/vb3H333fzsZz/jwIEDnHfeeaTTaQAaGxsZP358u/MOfW5sbOy2zKHjXVm9ejUVFRWFTTmkREREhq9eBVIrVqzodPL34duhobFjbdOmTTQ3N7NkyZJuy3meRyaT4Z577mHBggWceeaZ/OQnP2HHjh3U1tYWpW2HW7lyJdFotLAd3hMmIiIiw0uv8kjdeOONLF26tNsyU6dOPZr2dGnNmjVceOGFHXqJPqy6uhqAGTNmFPZVVVUxduxYGhoaAJgwYQJ79uxpd96hzxMmTOi2zKHjXQkGgwSDwR7ckYiIiAx1vQqkqqqqqKqqKlZbulRfX09tbS2bN28+YtlPfvKTAGzfvr2Q0v3AgQPs27ePKVOmADBnzhxuueUWMpkMfr8fyM+tOvnkkxk1alShzFNPPcX1119fqPvJJ59kzpw5x/LW+iSd3cfBxP/DI0N5YCYlwekD3SQREZERqWhzpBoaGqirq6OhoYFcLkddXR11dXXE4/FCmZ07d1JXV0djYyOJRKJQ5tBcpkPWrVtHdXU1Cxcu7HCdTZs2MX36B4HEtGnTuOiii7juuut47rnnePXVV1myZAnTp0/n7LPPBuDyyy8nEAiwfPlyXnvtNR544AG++93vcsMNNxTque6663j88ce588472bZtG6tWreKFF17g61//+rH+UfVYzkvyp70r+MM7n+RP+/8XO/ffzEu7L+Tl975Aa/rtAWuXiIjIiGWLZMmSJRbosNXW1hbKzJ8/v9My9fX1hTK5XM5OnDjR3nzzzZ1e57777rMfvo1oNGqXLVtmKysr7ejRo+0ll1xiGxoa2pXZunWr/dSnPmWDwaA9/vjj7Xe+850OdW/cuNFOmzbNBgIBe+qpp9pHH3201z+HaDRqARuNRnt97uE8z7Ov7F5qn60/yT5b/5EPbR+1z/35L2w89eZRXUNERETyevr9XfQ8UiPdscojdTDxHK/uubKbEobRoU8xdfQ3CfmnHHFtIBEREenasMwjNZLtif8ccLspYWlK/oFY62O0JF4il0t3U1ZERESOhV5NNpeBk8ntA3LdlvFIkc28Q6uNYe1+wqG/xOeO7p8GioiIjEAKpIaIgDuOfAei12UZ1wQxBmymiWTu/+Hl3iccmovff6KG+kRERIpAQ3tDxPjSL9JdEAWGct8JuGRxjcEYi5d5jWT8EVLJF8jlkv3VVBERkRFDgdQQURE6g9Hhc4HOepYMPhOi0n8ijvUwJodjcxjrYL1GcsknySRryeX293ezRUREhjUFUkOEMYbpVfcwofQyDP52x8LOeI4PnYXfhDFYjLUY8sGU3+4mkHocX/M3sPu/iNd8N15WCymLiIgcC0p/UGTHKv3B4RKZ93g/upFk9g2CJkTALcVHDmMtDlmMAcdaIrkX8bEbi8FgOfSgjTMeRv8ExzfxmLRHRERkuFH6g2Es7D+OSaOv4fiKrxP2n4yDwVgHx7gY48fBJWT/hMtuAExbCGXaNuvthYN/i2JoERGRo6NAaohyHD+l4Y8zqvyvCAVmYhw/WHBwcSz4vfpOZ1MBGDzIbcOmX+jXNouIiAw3CqSGuIBvImUlCwkHP43PLccYi2taMGSOcKYD6d/3SxtFRESGKwVSw4DrlhGOzCdUciE+3+Que6I60tCeiIjI0VBCzmHCcXwEgx/D544m2zoGm336CL1SHgTO6Lf2iYiIDEfqkRpmXN9x+EoWYoMLsF32TTngfhQTmN2vbRMRERluFEgNQ65bilP+Txj/XwJ8KKAy4IyBUf+hZWNERESOkob2hinHjWBH/xCbeBiTeAByu8CUQ/hzEL4Mx60c6CaKiIgMeQqkhjFjXEzkYohc3KfzrfWALPnsUz71YImIiHyIAinpwFoLNo71YmDT+Z0mBE45xikZ2MaJiIgMIgqkpB1rLdbbD95BwAcmAFiwrdhcK1CFcY7NUjciIiJDnQIpac+2ghcFE8GYw/7zMH6sTWJz+8GEMCYwcG0UEREZJPTWnrRjbQtg2wdRbYwJARmwyX5vl4iIyGCkQErasxkwbjcFHKzN9ltzREREBjMFUvIhLlivm+MepttAS0REZORQICXtGKcU8NpSH7RnbYb8BPRQv7dLRERkMFIgJe2ZMJhSsM3YttQH1lqsTYFNgFOBMcEBbqSIiMjgoLf2pB1jXHCrsJ4fvGasTQAmnwbBGYtxKga6iSIiIoOGAinpwBgfxh2Ldcrzk8/bAqnO3uQTEREZyfTNKF0yJtCWkLN38sOAqbZPvra8UxpFFhGR4UeBlBwz1ubasqLH+WCNPsBEwB2ruVUiIjLsqJtAjokPlpZpAhPEOBX5pWRMCdgENre37a0/ERGR4UOBlBwjKfCawZS0m0tljNv2FmBrfhMRERlGihZI3X777cydO5dIJEJlZWWnZa699lpmzZpFMBikpqamw/FVq1ZhjOmwlZSUHPH669evZ+bMmYRCIcaNG8c111xTOPb0009z0UUXUV1dTUlJCTU1NWzYsKHD+R++biik/Eldsikg18XSMgbwY714vzdLRESkmIo2RyqdTrNo0SLmzJnD2rVruyy3bNkynn/+eV555ZUOx2666SauvvrqdvvOOeccTj/99G6vfdddd3HnnXdyxx13MHv2bFpaWnj77bcLx5977jlmzpzJP/7jPzJ+/HgeeeQRrrzySioqKrjwwgsL5crLy9m+fXvhcz4gkM5ZCnOiOmNMWxkREZHho2iB1G233Qbke3a6cs899wCwd+/eTgOp0tJSSktLC5+3bt3K66+/zr333ttlnQcPHuQb3/gGDz/8MOecc05h/8yZMwv//+abb253znXXXccTTzzBL37xi3aBlDGGCRMmdHktOVz+PyVrbecBp82CU9pxv4iIyBA2pOZIrVmzhmnTpjFv3rwuyzz55JN4nse7777LKaecwsSJE7n00kvZtWtXt3VHo1FGjx7dbl88HmfKlClMmjSJiy66iNdee+2IbUylUsRisXbbiGDC+aVjbEuHQ9amAAfjHHlIVkREZCgZMoFUMplkw4YNLF++vNtyb731Fp7n8e1vf5u7776bn/3sZxw4cIDzzjuPdDrd6TkbN27kv//7v/nyl79c2HfyySezbt06fvnLX/KjH/0Iz/OYO3cu77zzTrfXX716NRUVFYVt0qRJvb/ZIcgYF+OMAVysF8PaBNamsF5zPqmnMwpjwgPdTBERkWOqV4HUihUrOp38ffi2bdu2ojR006ZNNDc3s2TJkm7LeZ5HJpPhnnvuYcGCBZx55pn85Cc/YceOHdTW1nYoX1tby5e//GX+8z//k1NPPbWwf86cOVx55ZXU1NQwf/58fvGLX1BVVcX3v//9bq+/cuVKotFoYTtST9hwYpwIxncctAVUADjlGLcaxx01oG0TEREphl7NkbrxxhtZunRpt2WmTp16NO3p0po1a7jwwgsZP358t+Wqq6sBmDFjRmFfVVUVY8eOpaGhoV3ZZ555hs997nP867/+K1deeWW39fr9fk477TR27tzZbblgMEgwOHITTxoTwLijsXZU2+feT9Dvcp6ViIjIINOrQKqqqoqqqqpitaVL9fX11NbWsnnz5iOW/eQnPwnA9u3bmThxIgAHDhxg3759TJkypVDu6aef5sILL+R//+//zVe/+tUj1pvL5fjjH//IZz7zmT7excjS20DI2izY1vxQIDksvraEnhEtLyMiIoNW0b6hGhoaqKuro6GhgVwuR11dHXV1dcTjH+QS2rlzJ3V1dTQ2NpJIJAplPjyXad26dVRXV7Nw4cIO19m0aRPTp08vfJ42bRoXXXQR1113Hc899xyvvvoqS5YsYfr06Zx99tlAfjjvs5/9LNdeey1f+MIXaGxspLGxkQMHDhTq+da3vsUTTzzBW2+9xUsvvcQVV1zBn//8Z6666qpj/aMa8azNYHPvY3N7wLY9e5vA5t7Devuw1hvYBoqIiHTFFsmSJUss+cRB7bba2tpCmfnz53dapr6+vlAml8vZiRMn2ptvvrnT69x33332w7cRjUbtsmXLbGVlpR09erS95JJLbENDwxHbNn/+/EKZ66+/3k6ePNkGAgE7fvx4+5nPfMa+9NJLvf45RKNRC9hoNNrrc0eKXHaPzaX/ZHOZ3dbLNha2XOYdm0tvt16uaaCbKCIiI0xPv7+NtVZZEosoFotRUVFBNBqlvLx8oJsz6FibxmbfBePHGH8nx1sBP8Y9TkN8IiLSb3r6/a1vJhlYNgNkOg2i8gJtZbL92CgREZGeUSAlA8wAhq47Ru1h5URERAYXBVIysEwwv5Hs/LhNggl302MlIiIycBRIyYAyxgVTDjbTtpRMnrW2bX6Um0+DICIiMggVbdFikZ4yTgUWC14T1h62NqEJYNwxGCcycI0TERHphgIpGXDGGIw7CuuUgE0BHuCCCWGM/hMVEZHBS99SMmgYEwAT6PV5+YnqabBZMA4QVKoEERHpFwqkZEizNoXNHQTbCuQAAyYMTgXGKR3o5omIyDCnQEqGLGvTHywrY8IY48svJ2OT+f2gYEpERIpK4x8yZFmvGWwS45QV5lIZ47RNTnewXpPW6RMRkaJSICVDkrU5sHEwoc4LmFA+B5XtIj+ViIjIMaBASoYoD2zb232dyE82t/lyIiIiRaJASoYoF4xLV2vwWXto4nnngZaIiMixoEBKhiRjHDBlYFOdz4Oyrfm39+hi6E9EROQY0Ft7MmQZpwxrE2CbsQTJ/+fstc2L8mPc0RijxY5FRKR4FEjJkGWMD9xxWC8EtjmfBgHackiVY7qaiC4iInKMKJCSIc0YH8YdjbUVHErIaYx/oJslIiIjhAIpGRaMcenqDb6uWJsBm8DaJPkALNyW2FMT1EVEpGcUSMmIZL0WrLevbZFkF7BYmsBEwK3Kr/snIiJyBAqkZMTJLy2zF/AwTsVh+y3YeP6YW62Fj0VE5Ij0TSEjjvVagDTGKWm33xgDpiSfOsG2DkzjRERkSFEgJSOPTQCdT0g/1AtlD70BKCIi0g0FUiIiIiJ9pEBKRh4TBjKdHjqUJV2TzUVEpCcUSMmIk58bFWibK/WB/GTzlvybeyYyMI0TEZEhRW/tyYhjTADcKqy3D+tFOZT+ACyYCMat0ht7IiLSIwqkZEQyTgmYQFtCzgTgKCGniIj0mgIpGbGM8YPxYyjv9bnWWj6YZ+VTD5aIyAilQEqkl/JZ0WNgk/kdxg9OOZhSBVQiIiOMAimRXrBec1tWdAsmCBggjc3tAScFzth8Yk8RERkRivbP59tvv525c+cSiUSorKzstMy1117LrFmzCAaD1NTUdDi+atUqjDEdtpKSko6Vfcj69euZOXMmoVCIcePGcc011xSOvf32253W+/vf/75dHQ8++CDTp08nFArx8Y9/nMcee6xXPwMZXqzNYHP7wTgYpxRj/Bjjwxx6y89ryr/1JyIiI0bRAql0Os2iRYv42te+1m25ZcuWsXjx4k6P3XTTTezevbvdNmPGDBYtWtRtnXfddRe33HILK1as4LXXXuPXv/41CxYs6FDu17/+dbu6Z82aVTj23HPPcdlll7F8+XJefvllLr74Yi6++GJeffXVHty9DEs2CWTyk9I/xBgfYLBevN+bJSIiA8fY/KzZolm/fj3XX389TU1NXZZZtWoVDz30EHV1dd3WtXXrVmpqanj22WeZN29ep2UOHjzI8ccfz8MPP8w555zTaZm3336bE088kZdffrnTnjCAxYsX09LSwiOPPFLYd+aZZ1JTU8O9997bbTsPF4vFqKioIBqNUl7e+0nNMnhYrwmb24txOn+O1iYBB8c3qX8bJiIix1xPv7+H1MzYNWvWMG3atC6DKIAnn3wSz/N49913OeWUU5g4cSKXXnopu3bt6lD2r/7qrxg3bhyf+tSn2Lx5c7tjW7Zs4dxzz223b8GCBWzZsuXY3IwMQQ75fFNdsDnyOalERGSkGDKBVDKZZMOGDSxfvrzbcm+99Rae5/Htb3+bu+++m5/97GccOHCA8847j3Q6vxBtaWkpd955Jw8++CCPPvoon/rUp7j44ovbBVONjY2MHz++Xd3jx4+nsbGx2+unUilisVi7TYYJEwL8nS5onF9aJodxyvq9WSIiMnB6FUitWLGi00nah2/btm0rSkM3bdpEc3MzS5Ys6bac53lkMhnuueceFixYwJlnnslPfvITduzYQW1tLQBjx47lhhtuYPbs2Zx++ul85zvf4YorruCOO+446nauXr2aioqKwjZpkoZ5hgtjAuBUtCXxTHJoVNzaNNhmMCVaWkZEZITpVfqDG2+8kaVLl3ZbZurUqUfTni6tWbOGCy+8sEMv0YdVV1cDMGPGjMK+qqoqxo4dS0NDQ5fnzZ49myeffLLwecKECezZs6ddmT179jBhwoRur79y5UpuuOGGwudYLKZgahgxzigsDtgo2GbysZQPnFEYp1JZ0UVERpheBVJVVVVUVVUVqy1dqq+vp7a2tsM8ps588pOfBGD79u1MnDgRgAMHDrBv3z6mTJnS5Xl1dXWFIAxgzpw5PPXUU1x//fWFfU8++SRz5szp9vrBYJBgMHjEdsrQZIzBuJVYWwY21bbTl++tEhGREadoCTkbGho4cOAADQ0N5HK5wht5J510EqWlpQDs3LmTeDxOY2MjiUSiUGbGjBkEAh98Ma1bt47q6moWLlzY4TqbNm1i5cqVhSHFadOmcdFFF3Hdddfxgx/8gPLyclauXMn06dM5++yzAbj//vsJBAKcdtppAPziF79g3bp1rFmzplDvddddx/z587nzzjv57Gc/y09/+lNeeOEFfvCDHxzzn5UMPca4fRrGyw8DpgAPcMGE2lIniIjIkGSLZMmSJZb8K07tttra2kKZ+fPnd1qmvr6+UCaXy9mJEyfam2++udPr3HffffbDtxGNRu2yZctsZWWlHT16tL3kkktsQ0ND4fj69evtKaecYiORiC0vL7dnnHGGffDBBzvUvXHjRjtt2jQbCATsqaeeah999NFe/xyi0agFbDQa7fW5Mnx4nmdz2QM2l37L5tJ/+mDLvG29XPNAN09ERD6kp9/fRc8jNdIpj5QAeLkoeO+DCWJMfug3/6uXyK82447HOJqoLiIyWAzLPFIiQ5G1ufzkdOMvBFHQNt/KRIBcfhFkEREZchRIiRSbTbXNiwp1ftyE2lIqdMxPJSIig5sCKZGiy4+eG2O6ON7VfhERGewUSIkUm/EDPqzNdFEgUygjIiJDiwIpkSLLZ0QvaRu+a/9uR37+VBrjlGGMfh1FRIYa/RNYpB8YZxTWZtuyofvAuGAzQA6cSjBao09EZChSICXSD4zxgzsObCvWawZyYMIYpxxMRL1RIiJDlAIpkX5ijA9MOcYpx1rbzeTzrvX1PBERKQ4FUiIDoDfBkLVeW09WDMhicTFOaVtPlr94jRQRkSNSICUyiFnrYb394DUBTtvbfWlsbg+YEnCrtGCyiMgAUiAlMpjZeD6IMpF2ixtbbH7ieu4AuOM13CciMkA0w1VkkLLWtk1M97ULoqBtaNBEwLYCyoguIjJQFEiJDFpZsOm24byO8sFVri2NgoiIDAQFUiKDlmnbbKdHP0juqWE9EZGBokBKZJAyxgdOBGyyixIpMAEwwX5tl4iIfECBlMggZpwywIf1WtotL2NtCmwqn5fK6J0REZGBor/AIoOYMWFwx2G9A2Bj5GMpA/jAGYtxKge2gSIiI5wCKZFBzjglYEJgE0COfD6poPJHiYgMAgqkRIYAY1wwpX0619p025t9BkxAQ4EiIseQ/qKKDFPWZrFeE3jNQLZtrx/rVGCcCi2ULCJyDCiQEhmGrPWwuX1gY2DCGBNpm6yeBm8fFg+c0cqILiJylPRPUpHhyLa2BVFlhblUxhiMCYIJgxdFGdFFRI6eAimRYch6ccDtdPjOGD/5jOiJfm+XiMhwo0BKZFjKQbdzoAzW5vqtNSIiw5UCKZHhyASg20DJa+uZEhGRo6FASmQYMqYEyL+592HWJgF/PjeViIgcFb21JzIcmTA4o8Dbj7X+fA8VNr+sDGDcKiX0FBE5BhRIiQxDxhhwRoMJYL0Y2LY39EwY45RjnL4l9xQRkfYUSIkMU8YYMGVgSjg8IadyR4mIHDsKpESGuXwKhN4N41mbT49gbStgMSbUlthTw4EiIodTICUi7Vibxub25pN6YgAHSxQIgDtWw4IiIocp2lt7t99+O3PnziUSiVBZWdlpmWuvvZZZs2YRDAapqanpcHzVqlVt2ZjbbyUlJUe8/vr165k5cyahUIhx48ZxzTXX9Kre9evXdzgeCuktJxnePlhapjWfFd0pwzglGKcCjMXm9mLbJqyLiEgRe6TS6TSLFi1izpw5rF27tstyy5Yt4/nnn+eVV17pcOymm27i6quvbrfvnHPO4fTTT+/22nfddRd33nknd9xxB7Nnz6alpYW333671/WWl5ezffv2wmfNLZFhzybBtoAp6fDfe369vijWa8G4wQFqoIjI4FK0QOq2224D8j07XbnnnnsA2Lt3b6eBVGlpKaWlHwwjbN26lddff5177723yzoPHjzIN77xDR5++GHOOeecwv6ZM2f2ul5jDBMmTOjyWiLDT/7tPmPczg+bQNuQ3+j+a5KIyCA2pBJyrlmzhmnTpjFv3rwuyzz55JN4nse7777LKaecwsSJE7n00kvZtWtXr+uNx+NMmTKFSZMmcdFFF/Haa68dsY2pVIpYLNZuExERkeFpyARSyWSSDRs2sHz58m7LvfXWW3iex7e//W3uvvtufvazn3HgwAHOO+880umOq913Ve/JJ5/MunXr+OUvf8mPfvQjPM9j7ty5vPPOO91ef/Xq1VRUVBS2SZMm9f5mRQZM/q28LtfhsxkwkX5sj4jI4NarQGrFihWdTtI+fNu2bVtRGrpp0yaam5tZsmRJt+U8zyOTyXDPPfewYMECzjzzTH7yk5+wY8cOamtre1zvnDlzuPLKK6mpqWH+/Pn84he/oKqqiu9///vdXn/lypVEo9HC1l1PmMigY0L5vFO2BWu9does1wr4MM6RX/YQERkpejVH6sYbb2Tp0qXdlpk6derRtKdLa9as4cILL2T8+PHdlquurgZgxowZhX1VVVWMHTuWhoaGPtfr9/s57bTT2LlzZ7flgsEgwaAm4srQZIwD7lhszoJtxlqHfAqEHBDAuGMxRv99i4gc0qtAqqqqiqqqqmK1pUv19fXU1tayefPmI5b95Cc/CcD27duZOHEiAAcOHGDfvn1MmTKlz/Xmcjn++Mc/8pnPfKYPdyAydBgTAHdCW0LOJODlgycl5BQR6aBoc6QaGhqoq6ujoaGBXC5HXV0ddXV1xOPxQpmdO3dSV1dHY2MjiUSiUObDc5nWrVtHdXU1Cxcu7HCdTZs2MX369MLnadOmcdFFF3Hdddfx3HPP8eqrr7JkyRKmT5/O2Wef3eN6v/Wtb/HEE0/w1ltv8dJLL3HFFVfw5z//mauuuupofzQig54xLsYpxXHH4rjjME5Fj4Moay1e+o94yWfxMm8WuaUiIgOraOkPbr31Vu6///7C59NOOw2A2tpazjrrLACuuuoqnnnmmQ5l6uvrOeGEE4D8nKf169ezdOlSXLfjK9nRaLRdrieAH/7wh/z93/89n/3sZ3Ech/nz5/P444/j9/sLZY5U78GDB/nKV75CY2Mjo0aNYtasWTz33HPthgxFpD0v8TjE74DcB3MDPd9MKP8GTqBm4BomIlIkxlprB7oRw1ksFqOiooJoNEp5eflAN0ekaLzEQxD9B/Jzqg7/s+IAPhj9IwVTIjJk9PT7e8ikPxCRwcvzUhD7p7ZPH/63mQdkofnb/dwqEZHiUyAlIkcv9Wuw3SWf9SBTh5d5q9+aJCLSHxRIicjRy71Dj/6ceO8WvSkiIv1JgZSIHD1nNPkhvCMwWqNPRIYXBVIicvSCF3BoeZnOGXBPwPj11quIDC8KpETkqDluGZR8tZsSFspuwhjTb20SEekPCqRE5JgwpX8HJX8HHFpCpi1oMhVQ8S84ofMHqmkiIkVTtIScIjKyGGMwZX+HF1kKqSfAOwju8RD8HziO1ucTkeFJgZSIHFOOWwaRL/T6PGtzYBPkF0h2wAS1tp+IDHoKpERkwFkvjvUOgE2RT+hpAB/WqcA4ozS3SkQGLQVSIjKgrE1gc3vzH0xZIWiyNgXePiwG444awBaKiHRNk81FZEBZrxnIYpxIu54nY4JggmBjWJsduAaKiHRDgZSIDBhrs+C1gAl1USIINt025CciMvgokBKRAXRogePO/xR90EP14YWQRUQGBwVSIjKAfGAC+V6nTuSH9Fww/v5tlohIDymQEpEBY4zBOGVAtsM8KGst2FYwkfx8KRGRQUhv7YnIwDKl4KTAa8Jap633KQc2A6YE42qhYxEZvBRIiciAMsYBZwyYcNsbfBkgkA+gTASjYT0RGcQUSInIgDPGAVOKcUqx1vYpAWdfzxMRORoKpERkUOltMGRtAuu1gE3k3+0zJRinVMvLiEi/UCAlIkOW9WLY3D7Aa5tbZfPZ0L0YuOMwTmSgmygiw5ze2hORIcnaJDa3H4wP45RhTAhjwhinAvCw3j5lRBeRolMgJSJDkvVagGznqRFMJJ8N3Sb7vV0iMrIokBKRockmwXQ+OyE/z8pgu0j0KSJyrCiQEpEhyhxh5RibfxtQRKSI9FdGRIYk45SSz4jeMZqyNgcYQG/uiUhxKZASkaHJhPObjWOtV9htbRZsPJ8x3YQHsIEiMhIo/YGIDEnG+MGtyqc/sC1twZQBHDAVGHeMEnSKSNEpkBKRIcuYILjVYBNAW6oDEwSCCqJEpF8okBKRIS2/vExJn87Nz6VKg7Vg/FrXT0R6rahzpG6//Xbmzp1LJBKhsrKy0zLXXnsts2bNIhgMUlNT0+H4qlWrMMZ02EpKuv7DuX79+k7PMcbw/vvvF8o9/fTT/MVf/AXBYJCTTjqJ9evXd6jr3//93znhhBMIhULMnj2bP/zhD739MYjIIGOtbcuK/i42+27b/76Dl9uvJJ4i0itFDaTS6TSLFi3ia1/7Wrflli1bxuLFizs9dtNNN7F79+5224wZM1i0aFGX9S1evLjDOQsWLGD+/PmMGzcOgPr6ej772c9y9tlnU1dXx/XXX89VV13Ff/3XfxXqeeCBB7jhhhv45je/yUsvvcQnPvEJFixY0C4YE5Ghx3oHsbk95JeWKcU45WBc8PZjc3vbeqpERI7M2M7eHT7G1q9fz/XXX09TU1OXZVatWsVDDz1EXV1dt3Vt3bqVmpoann32WebNm9ej6+/du5fjjz+etWvX8qUvfQmAf/zHf+TRRx/l1VdfLZT767/+a5qamnj88ccBmD17Nqeffjr/5//8HwA8z2PSpEn83d/9HStWrOjRtWOxGBUVFUSjUcrLy3t0jogUj7VpbPZdMG6HrOjWemDjGLe6Lb2CiIxUPf3+HnLpD9asWcO0adN6HEQB/PCHPyQSifDFL36xsG/Lli2ce+657cotWLCALVu2APnetBdffLFdGcdxOPfccwtlOpNKpYjFYu02ERlEbBLIdLq0TD6Bp5Nf9FhEpAeGVCCVTCbZsGEDy5cv79V5a9eu5fLLLycc/iCnTGNjI+PHj29Xbvz48cRiMRKJBPv27SOXy3VaprGxsctrrV69moqKisI2adKkXrVVRIrtUJqELhi3rYyIyJH1OpBasWJFlxO5D23btm0rRlvZtGkTzc3NLFmypMfnbNmyhTfeeKPXwVdfrVy5kmg0Wth27drVL9cVkZ5yuz9ss+iFZhHpqV7/tbjxxhtZunRpt2WmTp3a1/Z0a82aNVx44YUdeomOdE5NTQ2zZs1qt3/ChAns2bOn3b49e/ZQXl5OOBzGdV1c1+20zIQJE7q8XjAYJBjsZDV6ERkcTBDwYW0SY0LtDuXf2LOaHyUiPdbrQKqqqoqqqqpitKVb9fX11NbWsnnz5h6fE4/H2bhxI6tXr+5wbM6cOTz22GPt9j355JPMmTMHgEAgwKxZs3jqqae4+OKLgfxk86eeeoqvf/3rfb8RERlQxgTAHYPNvZ9/O88EyQ/1pcGmwakEExngVorIUFHUOVINDQ3U1dXR0NBALpejrq6Ouro64vF4oczOnTupq6ujsbGRRCJRKJNOp9vVtW7dOqqrq1m4cGGH62zatInp06d32P/AAw+QzWa54oorOhy7+uqreeutt/iHf/gHtm3bxn/8x3+wceNG/v7v/75Q5oYbbuA///M/uf/++3njjTf42te+RktLC1/+8peP5sciIgPMOOUYdwKYUD54sknAxbjjMc7YtknnIiI9YItoyZIlFuiw1dbWFsrMnz+/0zL19fWFMrlczk6cONHefPPNnV7nvvvus53dypw5c+zll1/eZftqa2ttTU2NDQQCdurUqfa+++7rUObf/u3f7OTJk20gELBnnHGG/f3vf9/j+7fW2mg0agEbjUZ7dZ6IFJ/nedbzUm1bbqCbIyKDSE+/v/slj9RIpjxSIsNL/k9mEuslAC+/rIwJ54cMRWTY6On3t15NERHpIWs9rHcAvCiH0ihYPMAP7miMUzHALRSR/qZASkSkh6zXBN4BMCUY88GfT2tT2NxewNUbfyIjjGZUioj0gLUZ8GJgQu2CKKAtS7rJL4Ss2RIiI4oCKRGRnrBpIAN0MRfKBMGmgGw/NkpEBpoCKRGRHrMY09XyMqZQRkRGDgVSIiI9YfyAH2vTnR+3KTABwN+frRKRAaZASkSkB4wJgFMGNpHPiH4YazNALp/os8seKxEZjvTWnohIDxmnMr8en23GWkN+AeQs4IAzBoze2BMZaRRIiYj0kDE+cMeBLcPaFrBZMOUYJwKE1BslMgIpkBIR6QVjnHweKUr6dL61HpADHIxxj2nbRKT/KZASEekH1ubAxrFeM9gMGAdrSjFOmZaXERnCNNlcRKTIrM1hc+9jc3uATNsbgAa8/dhcI9amBrqJItJHCqRERIrMejGwzWDKMCaCMX6MCebX5rNpbO6AMqKLDFEKpEREishary2ICuTnV32YCYNNAOqVEhmKFEiJiBRVFmyOrhJ15tfty7WVEZGhRoGUiEhROeSXj+k8UMq/xWf4YIkZERlKFEiJiBSRMT5wStoWNO5MMr/gsQn1a7tE5NhQICUiUmTGKQf8WK+5sLyMtRZrW8HmME5l5/OnRGTQUx4pEZEiMyYIvvHY3AGwrVgsWMAEMe5YjFM20E0UkT5SICUi0g+MCYN7HJBsm1huwISU3VxkiFMgJSLST/Jr8YV7Pa88n2MqlV/bD9pSKSgbushgoEBKRGQQszaNze0H2wp4+SDM+rBOedvcKvVoiQwkBVIiIoOUtVlsbm8+iDKRtpxTYEmDdyA/18oZ09bTJSIDQa+JiIgMVjbRFkSVFYIoID+sZyLgxVBGdJGBpUBKRGSQsl4z4Hba45QPrLLd5KcSkf6gQEpEZNCy0G1+KQN4/dUYEemEAikRkcHKBMBmOj2Uf5MPNNVVZGApkBIRGaSMKQEcbGfBlE3kAy0T7vd2icgH9E8ZEZHByoTBGZV/Q8+m8oETFmwacDHOmHaT0EWk/+k3UERkkDLG5AMpE8B6sXwAZQw4lRinFKOFjkUGXFGH9m6//Xbmzp1LJBKhsrKy0zLXXnsts2bNIhgMUlNT0+H4qlWrMMZ02EpKSrq87vr16zs9xxjD+++/D8AvfvELzjvvPKqqqigvL2fOnDn813/91xGvPX369D7/PEREessYkw+a3GqMbxLGnYjjjlUQJTJIFDWQSqfTLFq0iK997Wvdllu2bBmLFy/u9NhNN93E7t27220zZsxg0aJFXda3ePHiDucsWLCA+fPnM27cOACeffZZzjvvPB577DFefPFFzj77bD73uc/x8ssvt6vr1FNPbVfPb3/7217+FEREjl7+H3NurzKZW5vBelG87Hv5Ldc2RCgix0xRh/Zuu+02IN9D1JV77rkHgL179/LKK690OF5aWkppaWnh89atW3n99de59957u6wzHA4TDn8wAXPv3r385je/Ye3atYV9d999d7tzvv3tb/PLX/6Shx9+mNNOO62w3+fzMWHChC6vJSIyGFmbbMuKngD8bZkSWrBeFNwqjFN6pCpEpAeG3Ft7a9asYdq0acybN6/H5/zwhz8kEonwxS9+scsynufR3NzM6NGj2+3fsWMHxx13HFOnTuVv/uZvaGho6PZaqVSKWCzWbhMR6U/W5trW50thnAqME8GYCMYpB2Owub3qmRI5RoZUIJVMJtmwYQPLly/v1Xlr167l8ssvb9dL9WH/8i//Qjwe59JLLy3smz17NuvXr+fxxx/ne9/7HvX19cybN4/m5uYu61m9ejUVFRWFbdKkSb1qq4jIUbMJsC1gOvY6GRMGslivpf/bJTIM9TqQWrFiRZcTuQ9t27ZtK0Zb2bRpE83NzSxZsqTH52zZsoU33nij2+Drxz/+MbfddhsbN24szKECWLhwIYsWLWLmzJksWLCAxx57jKamJjZu3NhlXStXriQajRa2Xbt29bitIiLHgrVpwHS9mLHxtQ35icjR6vUcqRtvvJGlS5d2W2bq1Kl9bU+31qxZw4UXXsj48eN7dU5NTQ2zZs3q9PhPf/pTrrrqKh588EHOPffcbuuqrKxk2rRp7Ny5s8sywWCQYDDY4/aJiBxrxjhYbNcFLPk0CiJy1HodSFVVVVFVVVWMtnSrvr6e2tpaNm/e3ONz4vE4GzduZPXq1Z0e/8lPfsKyZcv46U9/ymc/+9ke1ffmm2/ypS99qcdtEBHpfwHAYG2ui7f8MhhndCf7RaS3ijpHqqGhgbq6OhoaGsjlctTV1VFXV0c8Hi+U2blzJ3V1dTQ2NpJIJApl0ul0u7rWrVtHdXU1Cxcu7HCdTZs2dZrf6YEHHiCbzXLFFVd0OPbjH/+YK6+8kjvvvJPZs2fT2NhIY2Mj0Wi0UOamm27imWee4e233+a5557jkksuwXVdLrvssqP5sYiIFJcJgSkBG8faDxY1ttZivea241paRuRYKGr6g1tvvZX777+/8PlQWoHa2lrOOussAK666iqeeeaZDmXq6+s54YQTgPwbdevXr2fp0qW4bsd/XUWjUbZv395h/9q1a/n85z/faTLQH/zgB2SzWa655hquueaawv4lS5YU0jW88847XHbZZezfv5+qqio+9alP8fvf/35AeuRERHrKGAfcsdgcYFvaFjg2gAUTwrhjMSYwwK0UGR6M/WAJcSmCWCxGRUUF0WiU8vLygW6OiIwg1npgE4VUB6ZtkePeJPUUGal6+v2ttfZERIYpYxwwJRi6XlKrK/l/Y2fIz0x3MMZ/rJsnMiwokBIRkXasTWBz0bYUCR7gYJ0yjFOuIUGRD1EgJSIiBdZrxeb2ALm2CekOkAFvf36I0B2n3imRwwypzOYiIlI81npY7yDgYZwyjPFhjIMxQTAV+YnrXtcrO4iMRAqkRESkTSo/nNdJagRjDJgg2GaszQ1A20QGJwVSIiKSl8+X0M1bfT6wHvl5UyICCqREROQQ45LPiN5VoJQD46CvDpEP6LdBRETaBPPDep0saGytBZsEU6Y8VCKHUSAlIiJAPu+UcSoBg/XiWJvNLytjM2Bj+ZxUTulAN1NkUFH6AxERKTBOCTAO60XzPVDkAB84lRinUnmkRD5EgZSIiLRjnBIwESBNfmK5qwBKpAsKpEREpANjDBDs9Xn5YcBDPVmmbW0/BWEyfCmQEhGRY8J6zdjcfvJr9EF+nT4/1qnAOKPagjOR4UWBlIiIHDXrtWBze8E4GFP+wX6bBm8fFgfjVg5cA0WKRG/tiYjIUbHWYr0Y+WSe7bOiGxNoy4gexdrswDRQpIgUSImIyFFqmxdluppTFQSbAZvu11aJ9AcFUiIicox0PgcqPzfKtm0iw4sCKREROUo+MH7y6RI6sjZ9WBmR4UWBlIiIHJV8RvRysNkO86Cs9fJLzjhlSoMgw5Le2hMRkaNnSsFJg9eEtYDxgfXIL3Rc2rb0jMjwo0BKRESOmjEOOGPAhLFeHEiDCWCcMjARLXQsw5YCKREROSaMMW0LG5f0uQ5rvba6NPNEhgYFUiIiMqCstWBbsV4zkAIM1pRgnFJMlykVRAYHBVIiIjJg8sk8D4J3ELBgAvn/9fbnAyt3HMaJDHQzRbqkQEpERAaOTeSDKBNo/1afCeWXnfH2gwlqjpUMWhqEFhGRAWNtnPzSMp2kRjCRfMZ0m+j3don0lAIpEREZODbVZaLOfEZ0AK3RJ4OXAikRERlATlu+qa5Yulp6RmQwUCAlIiIDxjilQCb/5t6H5LOk+7pZDFlk4CmQEhGRgWMibXOh4libK+y2NgO2BZwyQIGUDF5FDaRuv/125s6dSyQSobKystMy1157LbNmzSIYDFJTU9Ph+KpVqzDGdNhKSrpO+LZ+/fpOzzHG8P777wPw9NNPd3q8sbGxXV3//u//zgknnEAoFGL27Nn84Q9/6PPPQ0RE2jPGj3Gr2oKpVqwXw3oxsBlwRmOc0YfNlRIZfIqa/iCdTrNo0SLmzJnD2rVruyy3bNkynn/+eV555ZUOx2666SauvvrqdvvOOeccTj/99C7rW7x4MRdccEG7fUuXLiWZTDJu3Lh2+7dv3055eXnh8+HHH3jgAW644QbuvfdeZs+ezd13382CBQvYvn17h3pERKRvjAmCW932dl4WMG0pD7TIsQx+RQ2kbrvtNiDfQ9SVe+65B4C9e/d2GkiVlpZSWlpa+Lx161Zef/117r333i7rDIfDhMPhwue9e/fym9/8ptNgbty4cV32lt1111185Stf4ctf/jIA9957L48++ijr1q1jxYoVXV5fRER6J7+8TN8Sb3qZNyD1DNgcBE7DBOaoF0v6zZCbI7VmzRqmTZvGvHnzenzOD3/4QyKRCF/84hc7HKupqaG6uprzzjuP3/3ud4X96XSaF198kXPPPbewz3Eczj33XLZs2dLltVKpFLFYrN0mIiLHnpfbj7f/b2D/Rdj4v2Jb/g0OLsXbdz5e+o2Bbp6MEEMqkEomk2zYsIHly5f36ry1a9dy+eWXt+ulqq6u5t577+XnP/85P//5z5k0aRJnnXUWL730EgD79u0jl8sxfvz4dnWNHz++wzyqw61evZqKiorCNmnSpF61VUREjszzUnDgSmzmRQAMFkNbGoVcA/bA35BNbx/AFspI0etAasWKFV1O5D60bdu2rRhtZdOmTTQ3N7NkyZIen7NlyxbeeOONDsHXySefzP/8n/+TWbNmMXfuXNatW8fcuXP513/916Nq48qVK4lGo4Vt165dR1WfiIh0IvUryO34IHg6jMECLdjm/002/VanqRVEjpVez5G68cYbWbp0abdlpk6d2tf2dGvNmjVceOGFHXqJjnROTU0Ns2bNOmLZM844g9/+9rcAjB07Ftd12bNnT7sye/bsYcKECV3WEQwGCQb1qq6ISFElHsbidBpIQVsPVeZF0olHcLwz8ftn4rqhfm6kjAS9DqSqqqqoqqoqRlu6VV9fT21tLZs3b+7xOfF4nI0bN7J69eoela+rq6O6uhqAQCDArFmzeOqpp7j44osB8DyPp556iq9//eu9br+IiBxDXlOXQdQhhgypTCPknsIL7icYPAPXN6afGigjRVHf2mtoaODAgQM0NDSQy+Woq6sD4KSTTiq8ibdz507i8TiNjY0kEolCmRkzZhAIfPDq67p166iurmbhwoUdrrNp0yZWrlzZYUjxgQceIJvNcsUVV3Q45+677+bEE0/k1FNPJZlMsmbNGn7zm9/wxBNPFMrccMMNLFmyhL/8y7/kjDPO4O6776alpaXwFp+IiAwQdzI2+1qXwZQFPCLgQNaLkm39Lens+4TDcwn4P4IxQ2qKsAxiRQ2kbr31Vu6///7C59NOOw2A2tpazjrrLACuuuoqnnnmmQ5l6uvrOeGEE4B8T9D69etZunQprut2uE40GmX79o6TCteuXcvnP//5TtMbpNNpbrzxRt59910ikQgzZ87k17/+NWeffXahzOLFi9m7dy+33norjY2N1NTU8Pjjj/dqaFFERIogcikm9Wi3RVLOFHJeFmscMJDMbCftNREOnkkk+AkN9ckxYaxm4RVVLBajoqKCaDTaLvGniIj0nbUW23QDNvVohyWNLYYcFTS7c/Hw4TkO1jp41pefeO6WEAl8gpLwbAK+sQPSfhn8evr9XdQeKRERkWIwxkDlv2CbJ2Jbf4ShBQCLS8ZMocX5GJ4xeNZibY4chqyXoiW3i1jyTXLNG3CdMsaVXMTEiqsJ+EYN8B3JUKUeqSJTj5SISHHlsgfJtj5MLvNHMp7BMwE8cuTw8GyOnLFkvAx7kltIewc/dLYh4FbxiQkbCfknDkj7ZXDq6fe3ZtuJiMiQ5vpG4S+9HF/JMlz/NHDBMwZrHazxYa3DgdRrpL2mTs62pHP7eGPvtf3dbBkmFEiJiMiQ5zg+AqFTCJd+Dr97Gg75t76t9ZHxcrTkdpF/l68zHvH0K+xv3ULWSyiBp/SKAikRERk2XN94wiUXEAyejXErsFjSthmOkHMKoCnxO1ozu0lm9+LZbPEbK8OCAikRERlWXLeESORTlEX+Cr97AnR4r6+L85wIPidM2ouRzn14LpVI5/TWnoiIDDuO4xIOzcDxjYaWcvYmn8XSfS9TeegvyS8u4xDPvI9nDX63FJ8J5t8SFOmEeqRERGTYCvomMLr0QsaEz6XrnimHsuAZ+N3jiGf3k8g20ZLZT3OmkebMHlqy+8lpqE+6oEBKRESGNZ9bwkfH3kFF8Iy2PYe++vKBVch3AieMupXWbBNZL4XfCRFwIwTdUvxOiLTXTCLbpEno0ikN7YmIyLDnc8N8fML/ZX/rr9nd/ADJbAOuKWNsyYWMjiwga9PkbAsBJ4LF4uDgGD+OcTH4iWf2AoagW4rPBDTUJwUKpEREZEQwxmFsyfmMLTmfrNdKa2YPjuPHNX6SmSiOcQBL1mvF75QAPloyUZJeK8lcnLTnEfaVE3QjlPgqcI2+QkVDeyIiMgL5nAgh3xiszZLOxch4CbJekoyXwO9ECLiVtOZiJHLN+IyPkBsh7EbwOwGS2TjxzEE8mxvo25BBQOG0iIiMSAG3HNcE2oKoHCkvTtg/Fp8JkbVZ0l6CgBPCGAc8i2N8uMZHwAkRz0YxuIR9Zfgd/0DfigwgBVIiIjJiuU4I1wnhcyI0Z/ZgcDHGJZNrzc+VMi5pm8B1ArgmSCLXSmu2hdZcnNZsgrLAGEJOmFJfGT5HX6kjkYb2RERkxPM5QSK+0Xg2RzIXJ+slyXppUl4LBoewW0HKSxJLR7FAyAkTckP48dOaixPNHFSKhBFK4bOIiAgQdEtxjI9MrpWclwNrCDsV+J0wGIfWXAyf48vPk8q14jN+fI4Px4vQnG3GNT71TI1AetoiIiJt/E4IvxMi6JbhOmEM4Dp+krkEWZsl4paQ9TIYDH43RNpL05JtJZ6JE8vEGR0YS9gXotRXorlTI4SG9kRERD7EdfyU+irxrEcyFyfjZch5OVK5BDmbJeKW4XmGg+koSS9JwA0QdP04xiGebeFgOkrW01DfSKBASkREpBNBN0JFoIqQW4oBLB5+J0SZfzQht5SWbAseHmE3DFhc4yfg+ClxIyRyCaLpGDmlSBj2NLQnIiLSBb8TxO8ECbvluE4ILAScEGkvQ8pLE3QDeNaSsx5lbpic9WjNJollWtlvm2n10pS6EUr8EQKaOzUsqUdKRETkCHyOn3J/JdZAItdKJpciZ7PkbJZkLkHYF8HBR1M6RjQTx3FcfMYFC7FMnP3Jg6S9zEDfhhSBAikREZEeCLlhRgVGE3IjeMaSJYvFUuGvpMxXTjKXIpFLU+IL4TMOPscl6AYo9ZeQ9jIcTDfjWW+gb0OOMfUzioiI9FDACRIIBCn1leI3+SG+sC9CzuZoyaUIuPmv1bSXptRXisGhNZuiOZNmT6qZRDZDhb+EEn9IQ33DhHqkREREesnn+KkIVOBzfLRmW8l4WTwvh7XQmksQcIKEnBDRdCv7UjFy1gMsOZvjYCbO+8kmUjkN9Q0HCqRERET6IOQGGR2oJOwLk8llSdkUGS9LiS9Chb+MrPWIZVsJuQFCro+A4yPkBijzRch4nob6hgn1K4qIiPRR0A0ScAKU+Urxu0Fas0kq/GVYa2nKxnBxcI1DIpci6OTTIyRzGZLZLI2tTbRkMowOllLiCxaGBWVo0VMTERE5CsYY/MbP6EA51nq0ZBP4HR8ZL4drIJFLAYZSf4R4JsWBdAuetXjWI+1l2Z+KE8skqQqVEfEFBvp2pJc0tCciInIMBN0Ao4OVhN0gKS9DIpekNZfG7/gYFSjFWIcD6RZcHEp9IYI+P2FfgPJAPv/U/lScrKcEnkONeqRERESOkZAbIOj4KfcyhJwgTekWRgfLMDg0pVrJWY+IL0jGZnGNQ9Dxk/VyeNbyTvwgqVyWqlAZYV+AgOMO9O1IDyiQEhEROYaMMW29U+XkrKUlkyLsC5LysviNj7SXIWNzVPrCZK3HvmQLqVyGlJclnknhWQg4LlXhMko01DfoaWhPRESkCAKOj7HBcsI+P4lcipZsgni2FQ+PUf4IEX+IfckWMl6W8kCIsBsg4gtQEQhjgb2JZtIa6hv0ihpI3X777cydO5dIJEJlZWWnZa699lpmzZpFMBikpqamw/FVq1ZhjOmwlZSUdHnd9evXd3qOMYb3338fgKVLl3Z6/NRTT+322tOnTz+qn4mIiIwcQdfP+NAoJoRHMblkLBWBEsYHKykLREjlsqRyGUp8AXLWwyEffHnWYjDsTcRpaD5ALJ0knVNANVgVdWgvnU6zaNEi5syZw9q1a7sst2zZMp5//nleeeWVDsduuukmrr766nb7zjnnHE4//fQu61u8eDEXXHBBu31Lly4lmUwybtw4AL773e/yne98p3A8m83yiU98gkWLFrU779RTT+XXv/514bPPp9FQERHpOWMMYTeAL+TiWZtfyNgJkfJyOMbk92XTlPqD+IzLvkQL8WyKZDZDxrPkLARdlzHBEsoCwYG+HfmQokYFt912G5DvIerKPffcA8DevXs7DaRKS0spLS0tfN66dSuvv/469957b5d1hsNhwuFw4fPevXv5zW9+0y6Yq6iooKKiovD5oYce4uDBg3z5y19uV5fP52PChAldXktERKQn/I7L2FAZ+1NxWrIp4ukU8XQSE4BSf5AKf4SmdIJ4JkXE78c1hpDrpzIQojWbZm8yjt9xCPn8A30rcpgh172yZs0apk2bxrx583p8zg9/+EMikQhf/OIXuyyzdu1azj33XKZMmdJu/44dOzjuuOMIhULMmTOH1atXM3ny5C7rSaVSpFKpwudYLNbjdoqIyPAWdP1MCFeQzGUIuwFc41AZCBP2BUh5OVqyacI+P65x2t7wC2Ct5ZX397D5zTdozWY5sWIUfzP9E5w6dvxA344wxAKpZDLJhg0bWLFiRa/OW7t2LZdffnm7XqrDvffee/zqV7/ixz/+cbv9s2fPZv369Zx88sns3r2b2267jXnz5vHqq69SVlbWaV2rV68u9MSJiIh8mGMcIr4gQdePtYbWbAoLZL0cOc/iug7xbIqg4yOX87j66c28vHc3DgYPy+9372LDtq1cPPUU7pq/ENdVmoSB1OvJ5itWrOhyIvehbdu2bcVoK5s2baK5uZklS5b0+JwtW7bwxhtvsHz58i7L3H///VRWVnLxxRe3279w4UIWLVrEzJkzWbBgAY899hhNTU1s3Lixy7pWrlxJNBotbLt27epxW0VEZORwjUNVqISIL0g8k6Q5kySRTdGSTRI0PsaESvmn55+mbu9uADxsu/996K03uLb2UVrS6QG7B+lDj9SNN97I0qVLuy0zderUvranW2vWrOHCCy9k/Pied2euWbOGmpoaZs2a1elxay3r1q3jS1/6EoFA9/k6KisrmTZtGjt37uyyTDAYJBjUZEARETmygOtjQqSMRDZENJUgm7NUBEJEfAF2x2M8++6fuz3/8T/vYNJ/l3H5jE8wedTofmq1HK7XgVRVVRVVVVXFaEu36uvrqa2tZfPmzT0+Jx6Ps3HjRlavXt1lmWeeeYadO3d222N1eH1vvvkmX/rSl3rcBhERke64xqHUH6TEF8AxDs3pFK4x/L/3ug+iALKex8/feJ0X3nmPvz51Jp+fcSqOoxSR/amoP+2Ghgbq6upoaGggl8tRV1dHXV0d8Xi8UGbnzp3U1dXR2NhIIpEolEl/qKty3bp1VFdXs3Dhwg7X2bRpU6f5nR544AGy2SxXXHFFl21cu3Yts2fP5mMf+1iHYzfddBPPPPMMb7/9Ns899xyXXHIJruty2WWX9ebHICIickTGGMaESoj4A0QzKeKZNKarwhZMyuC0OuxvSfFS4x7+4aknOf//3s+f9u3rz2aPeEWdbH7rrbdy//33Fz6fdtppANTW1nLWWWcBcNVVV/HMM890KFNfX88JJ5wAgOd5rF+/nqVLl3Y6qS4ajbJ9+/YO+9euXcvnP//5LpOBRqNRfv7zn/Pd73630+PvvPMOl112Gfv376eqqopPfepT/P73vx+QHjkRERn+Aq5LdaSM1myGmWMntM2G+pC2IMrkDOZDoVZ900EuffABHrn8CiYeluJHisdYazt9TnJsxGIxKioqiEajlJeXD3RzRERkiLDW8qkHfsC78Vj7gCoHbrLrN/UMcNHJ07nt7HMoDQQwpst+LelGT7+/NZAqIiIyCBlj+D//43OEXF+7fieTNdjO+6oAsMB/vbmT92Ix9re2ov6S4lIgJSIiMkidNu44Hr7oS8w7bgpuW8+SYzsO6X1YIptlX7yVbXv28m40RjKT7Y/mjkgKpERERAaxj44ey/0XfJEHLriUz086mcpA6IjnRHx+In4/nrXsjjWzu7mZeEr5pophSGU2FxERGYkcx+H04ydTXVJGyDzPA2+83mVZA5xz4lQCPhcnY0hmMjRGY+xrbuGkqjGUhYL4lCLhmNFPUkREZIiYWDmKb559Dmccf3ynxw1QEQzxxVNmEEskeb+5hZZ0Fp/j0JxM8ef9B9kTbSad1VDfsaJASkREZAgJ+f3830u+yN98/BMEPpQS6JQxVdxx7vmEfX4OtCbwOYaxJREigQAB1yWdzfJOU5Q39+6nNZXWRPRjQOkPikzpD0REpFiak0n+682dHEwkGB8ppbq0jEjAT2OsmYOtCcaXllARDtOaSrM7FqcsEKIk6COZyTCurJRxZWWMLg3jaqivg55+f2uOlIiIyBBVFgpxySkziCaT7I42s7s5RjqXJZHJMLakhMpwmFQ2x/6WBH7XYXRpiKDPR86zZHIeuw40kcykmVBZjr+ThNdyZAqkREREhjDXcRgdiRD2+wn4XFrSaRwMQZ8PYwxNrUlS2SzjSksIOA7NiRT7m1vJZHL4HId4MkU661FVXkJ5+MhvBEp76ssTEREZBsJ+P1NGVTKutBS/63Iw3srB1iSJdJoxJREqwiESmSwH4q34XIdRJWHGlEXwuw6x1hRvv3+QvbG45k31knqkREREhomAz0d1eRlhv5/3fDEy2RwuUBIMQlvvVCaXY2xZCX7XoSWZZl+shYqQh8XSkkyTyeQYVRYhHPAP9O0MCQqkREREhhFjDKMiYfyuy754nHg6zcHWBI4xtCTTjC4JUx4K0prK8H40juu4VJSGcI0hnsiwLxanJZmmenQ5peHgQN/OoKdASkREZBgqDQYI+SsJ+vy8dzCK6zgYaxkVCQPQ1JIg5+WoKivDZxxaU2kOxlvxwkH2NSdoSWaYNLaC8tKQJqJ3Q4GUiIjIMOVzHCaUlxLyuextbiXWkiSWSJP1POLJNGPLw5SFA7Qk07wfbcF1DWWRICUeJNMZGptitKTTVI8qx+9TMNUZBVIiIiLDmDGGykiYcCCAawx7Y3FKgn4oKWF0JEzOs8RaE1jrMaa0DNdxaE2mibUmsdayr6mVZCpL9egySsNBjOl+weSRRoGUiIjIMGeMIeT3MWXsKEqDAfY2t9Bi0iQyWVqTGZoTGcZVRCgJ+oknUuyNtuB3XSKBAK6TJdaSAGsZVRZhbEWJgqnDKJASEREZIXyuw9jyEiJBP57nEU+kCfodRpeFqSyJkPM8oi1JMDCmLILPdUikLYlkBpuz7G9qIZnKUFVZSiiot/pAeaRERERGFNdxKAuHOPm4cUwaU0nY78c1DtlcPohqTqQZXRoiHPQRb02x90AL6UyWYMAPBvY1tdB4oJl4IjXQtzIoqEdKRERkBAr6fYwfVUrI7yOR2k8inQFgbGmEikiYTMajqSWB4xrGVpYS8LskkoZEMk0u63GgKc6Jx42hrCSEbwRPRFePlIiIyAgV8Pmoqixl+qTxjC0voTwUxO/zYYFYa4pEIk1lSZCgz6W5Jcneg3ESqQx+n0trIsOuxibePxAnlc4O9K0MGPVIiYiIjHAVJSH8rsP7vjjNjU20JDJkczlGV5RSURomk8kRbU7gug5jKkrw+xxe376b51+qZ+/+ZkpLQlwwfwYXL6ihvGxkrddnrBbVKapYLEZFRQXRaJTy8vKBbo6IiEiXPM9jz4E4TS0JWhNpwFISDrJ3fzNN8QRjRpVQ6g/w7/c/w2vbd2MMHIoiDDC6soR7vnUpUyaOGcjbOCZ6+v2toT0REREBwHEcxo8uo3pMOaXhAM2taQ7GE6TSWUaVlVAZCbP5yT/y2vbdwAdBFIAFDkRb+YdvbyLWnBiYGxgACqRERESkwHEMlaVhTjhuDCdMGE15KEhlaZiySJBszuPpLX/q8lxrLe82NrHpsZd4r7EJz/P6seUDQ3OkREREpIOg38fx4yo4GGulMRfjYHOC/fviJJOZbs8zBra+8S5jKsv46InjOPHEKgL+4ftWnwIpERER6ZTf51I1qpRQwEdgX4x4c7JH5xmg6WCcVxNp4vEkH/3oeMrLwsVt7ABRICUiIiJdMsZQXhrG7/MRdF3CIT+JbnqlrIUJY8sx1pJMZtj+p93Eoq2cfPIEqqtHDbvlZTRHSkRERI4oHPIz+fjRfO7cj9NVKGQMjKksobI0hOM62GwOrMfuPU289PKf+dO23aSHWc4pBVIiIiLSI67j8LUvfZozak4A8oHT4UrCQc4546P4HUMykcFa8HLgeBA70Mqrf2zgla0NxJqGz1t9yiNVZMojJSIiw43nWZ74f2/ws0de5L09UYIBH9NOqOKkSWMJuQYwZHIWP4ZMNkckEsDL5XBsPvI67rhKpp0ykXETKnCcwTnUN+B5pG6//Xbmzp1LJBKhsrKy0zLXXnsts2bNIhgMUlNT0+H4qlWrMMZ02EpKSrq99n//939zzjnnUFlZyahRo1iwYAFbt25tV+aVV15h3rx5hEIhJk2axD//8z93qOfBBx9k+vTphEIhPv7xj/PYY4/1+P5FRESGK8cxXDB/Bv/+T3/N3au+yJc/P5tPTDuOkrAfx+cjk7EEfC5ZL4fP55BszeDl8rmmjDG8s+sALz//Jm9ubxzyQ31FC6TS6TSLFi3ia1/7Wrflli1bxuLFizs9dtNNN7F79+5224wZM1i0aFGX9cXjcS644AImT57M888/z29/+1vKyspYsGABmUx+clwsFuP8889nypQpvPjii9xxxx2sWrWKH/zgB4V6nnvuOS677DKWL1/Oyy+/zMUXX8zFF1/Mq6++2oefhoiIyPATDPo56YTxzP7LqUyorsTzDMZAKOgnl83hug5gsXgYLMlEGmstjmOIxhK8uOVPPPeb19j3fnSgb6XPij60t379eq6//nqampq6LLNq1Soeeugh6urquq1r69at1NTU8OyzzzJv3rxOy7zwwgucfvrpNDQ0MGnSJAD++Mc/MnPmTHbs2MFJJ53E9773PW655RYaGxsJBAIArFixgoceeoht27YBsHjxYlpaWnjkkUcKdZ955pnU1NRw77339vj+NbQnIiIjQaw5wc439/DWm++TTmZwMFg8kon8IsfZdBbXdbA5S9O+GFt//xbv1u8DwO93mb9wJsuuP5/RVYPju3LAh/aKYc2aNUybNq3LIArg5JNPZsyYMaxdu5Z0Ok0ikWDt2rWccsopnHDCCQBs2bKFT3/604UgCmDBggVs376dgwcPFsqce+657epesGABW7Zs6baNqVSKWCzWbhMRERnuysvCnHrKRE47bTLlFRFy1mJwCAX9eDmLz/WBB+/9eT+/euC/ee/tfYVzM5kcTz1cxzWX/juN7x4YwLvovSETSCWTSTZs2MDy5cu7LVdWVsbTTz/Nj370I8LhMKWlpTz++OP86le/wufLp81qbGxk/Pjx7c479LmxsbHbMoeOd2X16tVUVFQUtkO9YiIiIsNdMOjjxBPGc/rpJzLx+FEYLNbme5z8AYd0OsOWJ1/FepYPj4dZa2na38Ld33yI2MGWgbmBPuhVILVixYpOJ38fvh0aGjvWNm3aRHNzM0uWLOm2XCKRYPny5Xzyk5/k97//Pb/73e/42Mc+xmc/+1kSieK/brly5Uqi0Whh27VrV9GvKSIiMlg4jmH8+EpOm3UC06ZXE/Q7GJMf0tv3bhOpRHfJPC2v/Pdb/O7Xr9K4a/+QWKuvV5nNb7zxRpYuXdptmalTpx5Ne7q0Zs0aLrzwwg69RB/24x//mLfffpstW7bgOE5h36hRo/jlL3/JX//1XzNhwgT27NnT7rxDnydMmFD4387KHDrelWAwSDAY7NW9iYiIDDfl5RE+/okplJWF2fHGe0QPttC0vwXjGKzX9fRsz7O8+cZ72EyWj35sMpM/Op5gONBl+YHWq0CqqqqKqqqqYrWlS/X19dTW1rJ58+Yjlm1tbcVxnHYp6A99PhTZzpkzh1tuuYVMJoPf7wfgySef5OSTT2bUqFGFMk899RTXX399oZ4nn3ySOXPmHMM7ExERGb58PpePfHQCJaVBdrz+Hq8F3qbDmF4nAj6XeCzBq394k9jBOB+dOYnyUaXFb3AfFG2OVENDA3V1dTQ0NJDL5airq6Ouro54PF4os3PnTurq6mhsbCSRSBTKpNPpdnWtW7eO6upqFi5c2OE6mzZtYvr06YXP5513HgcPHuSaa67hjTfe4LXXXuPLX/4yPp+Ps88+G4DLL7+cQCDA8uXLee2113jggQf47ne/yw033FCo57rrruPxxx/nzjvvZNu2baxatYoXXniBr3/968f6RyUiIjJsGWOoPm40NadPZf6Cjx8xjiqvDBOOBLCeJZvJ8Nbr7/LSs9t55633B+VQX9EWLb711lu5//77C59PO+00AGpraznrrLMAuOqqq3jmmWc6lKmvry+8Yed5HuvXr2fp0qW4rtvhOtFolO3btxc+T58+nYcffpjbbruNOXPm4DgOp512Go8//jjV1dUAVFRU8MQTT3DNNdcwa9Ysxo4dy6233spXv/rVQj1z587lxz/+Md/4xje4+eab+ehHP8pDDz3Exz72sWPzAxIRERlByisizD//4zz96FZe/N2OLgOqmTWTcZz8KFIuZ3Fdw55d+2iNttLa1MKUGRMJhvz92/huaImYIlMeKRERkQ8kE2n+6e9/zAu/25GfhmPyk8yNMZx+5kmcNH08rnHIZTOEwn7SySz+gAsWAgEfE08az0drplA5trjfqT39/lYgVWQKpERERDp69cV6Hvnp8+xtjBKOBJg6bQLBgIvf55JKJDEAOTCug+OAweIYg8Uy7vhRnFRzAsedOK7T0apjoaff30Ub2hMRERHpysdmncjJH5/In3fs4c3X36VpfzNYyOWyuI6LsR4Zm8PxLOlUjlBJkGwmg+M4NP55Py3NCVoOtnDixycTDA3cW30KpERERGRA+AN+PjLjeMorS9jxxwbefXtvvtfJtWSTHoGQj0wig+tzSLemcF2DtTkc47Dnrfd59enXGFM9inlfmM300z86IPegQEpEREQGjDGGccePIlIWorQiQv0b75BuzeAPueTSOVyfA9Yjay3Ws7TEErz4+Mu8VVdfmLC+8Z9/ySfOOpWVG65lTPXofm3/kFkiRkRERIav0vIwM2adwMzZH6ViTAnGGly/H8fnks16BAI+MokUtT96lrfq3u7w1t8fn32dv5/3/9ESa+3XdiuQEhERkUHBH/BzwinH8RfzZ3DcR8ZjjIdjIBDwk81keGdHI3t37aOz9+Q8z9JYv5dHvv9kv7ZZgZSIiIgMGo7jUHXcKGo+OY1T/vIj+P0u1svhD/p586W3OGzhkg4slv+6r7b/GosCKRERERmESitLmD7rRP7i7BmMGl+Bl/VojbV2nxndQtP70X5rI2iyuYiIiAxS/oCfyScfT6Qswp9efJOSyhIONkY7HdoDwMDo6lH92kb1SImIiMigZYyh6vjRzPz0KfyPy+d1HUS1uWDZ2f3UsjwFUiIiIjLolVWWsujGz3HyGSfll5b5EMd1mDz9eD7zlXP7tV0KpERERGRICIYC3PHUNzn3S5/G9X2wNIzjGOZ87i+585lvESkN92ubtNZekWmtPRERkWOvaW+UPz77Bp7nccqZH2XcpKpjWr/W2hMREZFhq7KqgnlfOHOgm6GhPREREZG+UiAlIiIi0kcKpERERET6SIGUiIiISB8pkBIRERHpIwVSIiIiIn2kQEpERESkjxRIiYiIiPSRAikRERGRPlJm8yI7tAJPLBYb4JaIiIhITx363j7SSnoKpIqsubkZgEmTJg1wS0RERKS3mpubqaio6PK4Fi0uMs/zeO+99ygrK8MYM9DNGRFisRiTJk1i165dWih6ENLzGdz0fAY3PZ/+Y62lubmZ4447DsfpeiaUeqSKzHEcJk6cONDNGJHKy8v1h2YQ0/MZ3PR8Bjc9n/7RXU/UIZpsLiIiItJHCqRERERE+kiBlAw7wWCQb37zmwSDwYFuinRCz2dw0/MZ3PR8Bh9NNhcRERHpI/VIiYiIiPSRAikRERGRPlIgJSIiItJHCqRERERE+kiBlAwa3/ve95g5c2Yh0dycOXP41a9+VTh+1llnYYxpt1199dXd1hmPx/n617/OxIkTCYfDzJgxg3vvvbddmWQyyTXXXMOYMWMoLS3lC1/4Anv27CnKPQ5lA/V8+lLvSFWMZ7Rnzx6WLl3KcccdRyQS4YILLmDHjh3tyuh3qGcG6vnod6jIrMggsXnzZvvoo4/aP/3pT3b79u325ptvtn6/37766qvWWmvnz59vv/KVr9jdu3cXtmg02m2dX/nKV+xHPvIRW1tba+vr6+33v/9967qu/eUvf1koc/XVV9tJkybZp556yr7wwgv2zDPPtHPnzi3qvQ5FA/V8+lLvSHWsn5HnefbMM8+08+bNs3/4wx/stm3b7Fe/+lU7efJkG4/HC+X0O9QzA/V89DtUXAqkZFAbNWqUXbNmjbU2/8fguuuu69X5p556qv3Wt77Vbt9f/MVf2FtuucVaa21TU5P1+/32wQcfLBx/4403LGC3bNlydI0fAYr9fPpar3zgaJ7R9u3bLVD4orfW2lwuZ6uqqux//ud/Wmv1O3S0iv18+lKv9I6G9mRQyuVy/PSnP6WlpYU5c+YU9m/YsIGxY8fysY99jJUrV9La2tptPXPnzmXz5s28++67WGupra3lT3/6E+effz4AL774IplMhnPPPbdwzvTp05k8eTJbtmwpzs0NA/31fPparxybZ5RKpQAIhUKFfY7jEAwG+e1vfwvod6iv+uv59KVe6aWBjuREDvfKK6/YkpIS67quraiosI8++mjh2Pe//337+OOP21deecX+6Ec/sscff7y95JJLuq0vmUzaK6+80gLW5/PZQCBg77///sLxDRs22EAg0OG8008/3f7DP/zDsbuxYaK/n09f6x3JjuUzSqfTdvLkyXbRokX2wIEDNpVK2e985zsWsOeff761Vr9DvdXfz6cv9UrvKJCSQSWVStkdO3bYF154wa5YscKOHTvWvvbaa52Wfeqppyxgd+7c2WV9d9xxh502bZrdvHmz3bp1q/23f/s3W1paap988klrrb4Eequ/n09f6x3JjvUzeuGFF+wnPvEJC1jXde2CBQvswoUL7QUXXGCt1e9Qb/X38+lrvdJzCqRkUDvnnHPsV7/61U6PxeNxC9jHH3+80+Otra3W7/fbRx55pN3+5cuX2wULFlhrP/iDcvDgwXZlJk+ebO+6666jv4FhrtjPpy/1SntH84wO19TUZN9//31rrbVnnHGG/du//VtrrX6Hjlaxn8/R1itHpjlSMqh5nleYB/BhdXV1AFRXV3d6PJPJkMlkcJz2/5m7rovneQDMmjULv9/PU089VTi+fft2Ghoa2s1bkM4V+/n0pV5p72ie0eEqKiqoqqpix44dvPDCC1x00UWAfoeOVrGfz9HWKz0w0JGcyCErVqywzzzzjK2vr7evvPKKXbFihTXG2CeeeMLu3LnTfutb37IvvPCCra+vt7/85S/t1KlT7ac//el2dZx88sn2F7/4ReHz/Pnz7amnnmpra2vtW2+9Ze+77z4bCoXsf/zHfxTKXH311Xby5Mn2N7/5jX3hhRfsnDlz7Jw5c/rtvoeKgXg+Pa1X8orxjDZu3Ghra2vtm2++aR966CE7ZcoU+/nPf77dOfod6pmBeD76HSo+BVIyaCxbtsxOmTLFBgIBW1VVZc855xz7xBNPWGutbWhosJ/+9Kft6NGjbTAYtCeddJL9X//rf3XIhQLY++67r/B59+7ddunSpfa4446zoVDInnzyyfbOO++0nucVyiQSCfu3f/u3dtSoUTYSidhLLrnE7t69u1/ueSgZiOfT03olrxjP6Lvf/a6dOHGi9fv9dvLkyfYb3/iGTaVS7c7R71DPDMTz0e9Q8RlrrR3YPjERERGRoUlzpERERET6SIGUiIiISB8pkBIRERHpIwVSIiIiIn2kQEpERESkjxRIiYiIiPSRAikRERGRPlIgJSIiItJHCqRERERE+kiBlIiIiEgfKZASERER6SMFUiIiIiJ99P8DqN96ZMq/od8AAAAASUVORK5CYII=", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "from scipy.interpolate import interp1d\n", "\n", "taup_model = TauPyModel(model='1066a')\n", "\n", "ray_path = taup_model.get_ray_paths_geo(source_depth_in_km=10,\n", " source_latitude_in_deg=35.766,\n", " source_longitude_in_deg=-117.605,\n", " receiver_latitude_in_deg=35.98249,\n", " receiver_longitude_in_deg=-117.80885,\n", " phase_list=['P', 'S'])\n", "\n", "ray_path[0].path\n", "\n", "# Define the number of points N\n", "N = 100\n", "\n", "# Create the interpolator objects for latitude, longitude, and depth\n", "interp_latitude = interp1d(np.linspace(0, ray_path[0].path['lat'].max(), len(ray_path[0].path['lat'])), ray_path[0].path['lat'])\n", "interp_longitude = interp1d(np.linspace(0, ray_path[0].path['lon'].max(), len(ray_path[0].path['lon'])), ray_path[0].path['lon'])\n", "interp_depth = interp1d(np.linspace(0, ray_path[0].path['depth'].max(), len(ray_path[0].path['depth'])), ray_path[0].path['depth'])\n", "\n", "# Resample the ray path to N points\n", "resampled_latitude = interp_latitude(np.linspace(0, ray_path[0].path['lat'].max(), N))\n", "resampled_longitude = interp_longitude(np.linspace(0, ray_path[0].path['lon'].max(), N))\n", "resampled_depth = interp_depth(np.linspace(0, ray_path[0].path['depth'].max(), N))\n", "\n", "plt.scatter(resampled_latitude, resampled_longitude, c=resampled_depth, cmap='viridis', alpha=0.1)\n", "plt.scatter(ray_path[0].path['lat'], ray_path[0].path['lon'], c=ray_path[0].path['depth'], cmap='viridis')\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": 101, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([(825.43624587, 0. , 0. , 10. , 35.766 , -117.605 ),\n", " (825.43624587, 0.26845098, 0.00012109, 11. , 35.7715199 , -117.61017979),\n", " (825.43624587, 1.63173458, 0.00177268, 11.01091139, 35.84678737, -117.68090241),\n", " (825.43624587, 2.99501819, 0.00342427, 11. , 35.92201332, -117.75175935),\n", " (825.43624587, 3.26346917, 0.00354535, 10. , 35.92752691, -117.75695956),\n", " (825.43624587, 4.33607208, 0.00402859, 6.00385879, 35.94952856, -117.77772004),\n", " (825.43624587, 4.47129397, 0.00408945, 5.5 , 35.9522991 , -117.78033535),\n", " (825.43624587, 5.20940225, 0.00442139, 2.74941477, 35.96740958, -117.79460335),\n", " (825.43624587, 5.57829224, 0.00458712, 1.37456112, 35.97495353, -117.80172933),\n", " (825.43624587, 5.94707298, 0.0047527 , 0. , 35.98248996, -117.80884997)],\n", " dtype=[('p', '" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "def compute_velocity_model(azimuth, elevation, interpolate):\n", " filename = list(output_csv.temp_files)[0]\n", " \n", " df = pd.read_csv(filename)\n", " filename = filename.split('/')[-1]\n", " \n", " # Current EQ location\n", " eq_lat = float(filename.split(\"_\")[0])\n", " eq_lon = float(filename.split(\"_\")[1])\n", " eq_depth = float(filename.split(\"_\")[2])\n", "\n", " # Define the region of interest (latitude, longitude, and depth ranges)\n", " lat_range = (np.min([df.st_lat.min(), eq_lat]), np.max([df.st_lat.max(), eq_lat]))\n", " lon_range = (np.min([df.st_lon.min(), eq_lon]), np.max([df.st_lon.max(), eq_lon]))\n", " depth_range = (0, 50)\n", "\n", " # Define the number of nodes in each dimension\n", " n_lat = 10\n", " n_lon = 10\n", " n_depth = 10\n", " num_points = 100\n", "\n", " taup_model = TauPyModel(model='1066a')\n", "\n", " # Create the grid\n", " lat_values = np.linspace(lat_range[0], lat_range[1], n_lat)\n", " lon_values = np.linspace(lon_range[0], lon_range[1], n_lon)\n", " depth_values = np.linspace(depth_range[0], depth_range[1], n_depth)\n", "\n", " # Initialize the velocity model with constant values\n", " initial_velocity = 0 # km/s, this can be P-wave or S-wave velocity\n", " velocity_model = np.full((n_lat, n_lon, n_depth), initial_velocity, dtype=float)\n", "\n", " # Loop through the stations and update the velocity model\n", " for i in range(len(df)):\n", " if ~np.isnan(df['velocity_p, km/s'].iloc[i]):\n", "\n", " ray_path = taup_model.get_ray_paths_geo(source_depth_in_km=eq_depth,\n", " source_latitude_in_deg=eq_lat,\n", " source_longitude_in_deg=eq_lon,\n", " receiver_latitude_in_deg=df.st_lat.iloc[i],\n", " receiver_longitude_in_deg=df.st_lon.iloc[i],\n", " phase_list=['P', 'S'])\n", "\n", " # Create the interpolator objects for latitude, longitude, and depth\n", " interp_latitude = interp1d(np.linspace(0, ray_path[0].path['lat'].max(), len(ray_path[0].path['lat'])), ray_path[0].path['lat'])\n", " interp_longitude = interp1d(np.linspace(0, ray_path[0].path['lon'].max(), len(ray_path[0].path['lon'])), ray_path[0].path['lon'])\n", " interp_depth = interp1d(np.linspace(0, ray_path[0].path['depth'].max(), len(ray_path[0].path['depth'])), ray_path[0].path['depth'])\n", "\n", " # Resample the ray path to N points\n", " lat_values_interp = interp_latitude(np.linspace(0, ray_path[0].path['lat'].max(), num_points))\n", " lon_values_interp = interp_longitude(np.linspace(0, ray_path[0].path['lon'].max(), num_points))\n", " depth_values_interp = interp_depth(np.linspace(0, ray_path[0].path['depth'].max(), num_points))\n", "\n", " # Loop through the interpolated coordinates and update the grid cells with the average P-wave velocity\n", " for lat, lon, depth in zip(lat_values_interp, lon_values_interp, depth_values_interp):\n", " lat_index = find_closest_index(lat_values, lat)\n", " lon_index = find_closest_index(lon_values, lon)\n", " depth_index = find_closest_index(depth_values, depth)\n", " \n", " if velocity_model[lat_index, lon_index, depth_index] == initial_velocity:\n", " velocity_model[lat_index, lon_index, depth_index] = df['velocity_p, km/s'].iloc[i]\n", " else:\n", " velocity_model[lat_index, lon_index, depth_index] = (velocity_model[lat_index, lon_index, depth_index] +\n", " df['velocity_p, km/s'].iloc[i]) / 2\n", "\n", " # Create the figure and axis\n", " fig = plt.figure(figsize=(8, 8))\n", " ax = fig.add_subplot(111, projection='3d')\n", "\n", " # Set the plot limits\n", " ax.set_xlim3d(lat_range[0], lat_range[1])\n", " ax.set_ylim3d(lon_range[0], lon_range[1])\n", " ax.set_zlim3d(depth_range[1], depth_range[0])\n", "\n", " ax.set_xlabel('Latitude')\n", " ax.set_ylabel('Longitude')\n", " ax.set_zlabel('Depth (km)')\n", " ax.set_title('Velocity Model')\n", " \n", " # Create the meshgrid\n", " x, y, z = np.meshgrid(\n", " np.linspace(lat_range[0], lat_range[1], velocity_model.shape[0]+1),\n", " np.linspace(lon_range[0], lon_range[1], velocity_model.shape[1]+1),\n", " np.linspace(depth_range[0], depth_range[1], velocity_model.shape[2]+1),\n", " indexing='ij'\n", " )\n", "\n", " # Create the color array\n", " norm = plt.Normalize(vmin=2, vmax=8)\n", " colors_vel = plt.cm.plasma(norm(velocity_model)) \n", " \n", " # Plot the voxels\n", " if interpolate:\n", " interpolated_velocity_model = interpolate_vel_model(velocity_model, initial_velocity, lat_values, lon_values, depth_values, n_lat, n_lon, n_depth)\n", " colors_interp = plt.cm.plasma(norm(interpolated_velocity_model))\n", " ax.voxels(x, y, z, interpolated_velocity_model > 0, facecolors=colors_interp, alpha=0.5, edgecolor='k')\n", " \n", " ax.voxels(x, y, z, velocity_model > 0, facecolors=colors_vel, alpha=1, edgecolor='black')\n", "\n", " # Set the view angle\n", " ax.view_init(elev=elevation, azim=azimuth)\n", "\n", "\n", "compute_velocity_model(0, 40, False)" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "35.766_-117.605_10.0_2019-07-04 17:33:49_3.csv\n", "35.766_-117.605_10.0_2019-07-04T17:33:49-00_3.csv\n", "35.766_-117.605_2019-07-04 17:33:49_3.csv\n", "35.766_-117.605_2019-07-04 17:33:49_9.csv\n", "current_vel_model.csv\n", "testt\n" ] } ], "source": [ "!ls data/velocity" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "35.766_-117.605_10.0_2019-07-04 17:33:49_3.csv\n", "35.766_-117.605_10.0_2019-07-04T17:33:49-00_3.csv\n", "35.766_-117.605_2019-07-04 17:33:49_3.csv\n", "35.766_-117.605_2019-07-04 17:33:49_9.csv\n", "current_vel_model.csv\n", "testt\n" ] } ], "source": [ "!ls '/Users/anovosel/Documents/phase-hunter/data/velocity'" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "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.2" }, "orig_nbformat": 4 }, "nbformat": 4, "nbformat_minor": 2 }