{ "cells": [ { "cell_type": "markdown", "metadata": { "id": "iSo6_GrBhoj4" }, "source": [ "# 気象庁が公開している地震データを解析する\n", "\n", "例えば、過去の地震データは \n", "https://www.data.jma.go.jp/svd/eqev/data/bulletin/shindo.html \n", "に公開されている。\n", "\n", "ここでは、それをダウンロードし、ファイル内容を解析して、震源地のマップやマグニチュードごとの頻度をプロットすることを目指す。" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "id": "otiZtC2LiMGd" }, "outputs": [], "source": [ "import urllib # ファイルのダウンロードに用いる\n", "import zipfile # Zipファイルの解凍に用いる\n", "import datetime # 日時データを取り扱うために用いる\n", "\n", "import numpy as np\n", "import matplotlib.pyplot as plt" ] }, { "cell_type": "markdown", "metadata": { "id": "9FBz4iuah-LK" }, "source": [ "## データの例\n", "\n", "2019年の地震データは以下のURLに公開されている。ダウンロードして、どのようなファイル構成になっているか見てみよう\n", "\n", "https://www.data.jma.go.jp/svd/eqev/data/bulletin/data/shindo/i2019.zip" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "id": "oIg13A5pac7J" }, "outputs": [], "source": [ "url = 'https://www.data.jma.go.jp/svd/eqev/data/bulletin/data/shindo/i2019.zip'" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "colab": { "base_uri": "https://localhost:8080/" }, "id": "MK_llfs0ao6p", "outputId": "edc8c7f1-0854-4170-db9e-0048fafc8164" }, "outputs": [ { "data": { "text/plain": [ "('i2019.zip', )" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# ファイルをダウンロードするコマンド。ファイルをダウンロードして、i2019.zipの名前で保存する\n", "urllib.request.urlretrieve(url, filename='i2019.zip')" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "id": "ajBY_GJ8azJ6" }, "outputs": [], "source": [ "# ダウンロードしたZipファイルを展開する\n", "with zipfile.ZipFile('i2019.zip', 'r') as zip_ref:\n", " zip_ref.extract('i2019.dat')" ] }, { "cell_type": "markdown", "metadata": { "id": "HF64dQWicgi_" }, "source": [ "## ファイルフォーマット\n", "\n", "ファイルのフォーマット（構成）は以下のページに公開されている。 \n", "https://www.data.jma.go.jp/svd/eqev/data/bulletin/data/shindo/format_j.txt \n", "ここでは、地震発生位置・震源地・マグニチュードのみを扱うことにする" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "id": "hNiLDPq3bXqs" }, "outputs": [], "source": [ "def read_earthquake_data(file):\n", " '''\n", " 地震データを読み込むプログラム\n", "\n", " 日時、緯度、経度、マグニチュードをリストとして返す\n", " '''\n", " # 展開したファイルを開く。ここで公開されているファイルには`shift_jist`というエンコーディングが使われているので注意する。\n", " # エンコーディングについては各自調べること\n", " with open(file, 'r', encoding='shift_jis') as f:\n", " lines = f.readlines()\n", "\n", " dates = []\n", " latitudes = []\n", " longtitudes = []\n", " magnitudes = []\n", " for line in lines:\n", " if line[0] in 'ABCD': # 震源レコード\n", " # date\n", " year = int(line[1:5])\n", " month = int(line[5:7])\n", " day = int(line[7:9])\n", " hour = int(line[9:11])\n", " min = int(line[11:13])\n", " try:\n", " sec = int(line[13:15])\n", " except ValueError:\n", " sec = 0\n", " # location\n", " latitude_deg = int(line[21:24])\n", " latitude_min = float(line[24:28]) * 0.01\n", " longtitude_deg = int(line[33:36])\n", " longtitude_min = float(line[36:40]) * 0.01\n", " # magnitude\n", " magnitude = line[52:54]\n", " try: \n", " magnitude = float(magnitude) * 0.1\n", " magnitudes.append(float(magnitude))\n", " dates.append(datetime.datetime(year, month, day, hour, min, sec))\n", " latitudes.append(latitude_deg + latitude_min / 60.0)\n", " longtitudes.append(longtitude_deg + longtitude_min / 60.0)\n", " except ValueError:\n", " pass\n", " return dates, latitudes, longtitudes, magnitudes" ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "id": "YETub7nibwwe" }, "outputs": [], "source": [ "dates, latitudes, longtitudes, magnitudes = read_earthquake_data('i2019.dat')" ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "colab": { "base_uri": "https://localhost:8080/", "height": 351 }, "id": "B33iE8jxjWaI", "outputId": "a90617af-762a-437e-8a3a-27dfb59bc4b6" }, "outputs": [ { "data": { "text/plain": [ "Text(0, 0.5, 'magnitude')" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAA24AAAE9CAYAAABz1DEXAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAABZcElEQVR4nO29f7ReR3nf+x3pyDYG2QghLDmAVXEpBBliI4HlWwpZaXNjp9ALgXuB/CDOAgz9CSs3TSm9K3FZDSFZadL0QuM6QIEWAwXCLeViJ4Q6wW4sEQkMlnCNQbGNQY7lH9gCjH2ONPeP877ye7b2j/nxzMwze38/a50l6X33nnlmnpnZe46e7zzGWgtCCCGEEEIIIXpZV9oAQgghhBBCCCH9cONGCCGEEEIIIcrhxo0QQgghhBBClMONGyGEEEIIIYQohxs3QgghhBBCCFEON26EEEIIIYQQopyl0gYs8uQnP9lu3769tBmEEEIIIYQQUoQDBw7ca63d0vxc1cZt+/bt2L9/f2kzCCGEEEIIIaQIxpg72j5nqCQhhBBCCCGEKIcbN0IIIYQQQghRDjduhBBCCCGEEKIcbtwIIYQQQgghRDncuBFCCCGEEEKIcrhxI4QQQgghhBDlcONGCCGEEEIIIcrhxo0QQgghhBBClMONGyGEEJUcuOMBvOe6b+DAHQ+oLpMQQgjJwVJpAwghhJAmB+54AD/33r14dOUETltahw+/YQ92nbdJXZmEEEJILvg/boQQQtSx9/B9eHTlBE5YYHnlBPYevk9lmYQQQkguuHEjhBCijj07NuO0pXVYb4ANS+uwZ8dmlWUSQgghuTDW2tI2nGT37t12//79pc0ghBCigAN3PIC9h+/Dnh2bxUIaU5RJCCGESGKMOWCt3d38nBo3QgghKtl13ibxzVWKMrXDzSohhIwDbtwIIYSQkcIDWQghZDxQ40YIIYSMFB7IQggh44EbN0IIIWSk8EAWQggZDwyVJMWh/oJ0wbFBSBy7ztuED79hD+cRIYSMAG7cSFGovyBdcGwQIsMUD2QhhJAxwlBJUhTqL0gXHBuEEEIIIY/BjRspCvUXpAuODUIIIYSQx2ACblIc6phIFxwbhBBCCJkaTMBN1EL9BemCY4OQ8vT9AoW/XCGEkHxw40YIIYSQVvoOCeIBQoQQkhdq3AghhBDSSt8hQTxAiBBC8sKNGyGEEEJa6TskiAcIEUJIXng4CSGECELNDxkb1LgRQkheeDgJIYQkhpofMkb6DgniAUKEEJIPhkoSQogQ1PwQQgghJBXcuBFCiBDU/BBCCCEkFQyVJKQF6jZkmUp/7jpvEz78hj2ja+tU/EcIIYRohhs3QhpQpyTL1PpzbJqfqfmPEEII0QpDJQlpQJ2SLOzPuqH/CCGEEB1w40ZIA+qUZGF/1g39RwghhOiAedwIaYGaHlnYn3VD/xFCCCH56Mrjxo0bIYQQQiYBfwlBCKkBJuAmhBBCyGThQTuEkNpJqnEzxjzRGPMJY8z/NMbcYoy5OGV9hBBCCCFt8KAdQkjtpP4ft98HcK219lXGmNMAnJm4PkIIIYSQU5gftLO8coIH7RBCqiTZxs0YcxaAFwO4DACstY8CeDRVfVphPP20of8JIQDXAg3sOm8TPvyGPfQDIaRaUv6P2w4ARwH8R2PMjwE4AOAt1trvJ6xTFYynnzb0PyEE4FqgiV3nbWLfE0KqJaXGbQnA8wH8gbX2QgDfB/C25kXGmMuNMfuNMfuPHj2a0Jz8MJ5+2tD/hBCAawEhhBAZUm7c7gJwl7V23+zfn8DqRm4N1tqrrLW7rbW7t2zZktCc/DBx7bSh/wkhANcCQgghMiTN42aMuR7AG6y1txpjrgDweGvtP+u6fox53KhrmDb0PyEE4FpACCHEnSIJuI0xFwB4L4DTABwG8EvW2ge6rq9548aHMpkKHOuE1A3nMCGE6KZIAm5r7U0ATql0bFB4TqYCxzohdcM5TAgh9ZI0AfdUoPCcTAWOdULqhnOYEELqhRs3ASg8J1OBY52QuuEcJoSQekmqcfOFGjcSC/2Qnpg+pn8IKQ/nISHEFa4XZShyOIkvNW/cSHmo3dAN/UMIIYTUA5/b5ejauDFUkowGajd0Q/8QQggh9cDntj64cSOjgdoN3dA/hBBCSD3wua0PhkqSUcFYbN3QP2vR2B8abSKEEFIGPhPKQI0bIYQoQqN2QKNNhBBCyNSgxo0QQhShUTug0SZCCCGErMKNGyGEFECjdkCjTYQQQghZhaGShBBSCI3aAY02EUIIIVOiK1RyqYQxpCxDL2Z8cSMhlBo3PvXmsNGnjl3nbVI3xzTaVAtcO8vBvidjRnJ8c67UDTduE2Po8AEeTkBCKDVufOrNYSPnz3Sh78vBvidjRnJ8c67UDzVuE2Po8AEeTkBCKDVufOrNYSPnz3Sh78vBvidjRnJ8c67UDzduE2Po8AEeTkBCKDVufOrNYSPnz3Sh78vBvidjRnJ8c67UDw8nmSC1ady02UPaocYtXx1EJ22+53joR6p/2M9kzExR41aLnalgAm5SJYzHJoTUCtevftg/hJA2uDYwATepFMZjE0JqhetXP+wfQkgbXBu64caNqIbx2ISQWuH61Q/7hxDSBteGbhgqSdQz9ThnQki9cP3qh/1DCGlj6msDNW6EQN9BFrm4et+duObgEVx6/jb87EVPL20OIYSQBcb0vGnic2jPmPshFeyzcdK1cWMCbjIZtCVrzsXV++7E2z91MwDg+tvuBQBu3gghRAljet40aWsbgNb2jrkfUsE+mx7UuJHJoC1Zcy6uOXik99+EEELKMabnTZO2tnW1d8z9kAr22fTgxo1MBm3JmnNx6fnbev9NCCGkHGN63jRpa1tXe8fcD6lgn00PatwqhPHM4WjQuJXwHzVucXDOEUJS0lxjxrTmUOOWFo19ptGm2uDhJCOB8cx1Q//VB31GCMkJ1xxSMxy/MjAB90hgPHPd0H/1QZ8RQnLCNYfUDMdvWrhxqwzGM9cN/Vcf9BkhJCdcc0jNcPymhaGSFcLY4bqh/+qDPiOE5IRrDqkZjt94qHEjg3CiEUJqJnQN49pXP/QhmRoc8+OGCbhJLxSTEkJqJnQN49pXP/QhmRoc89OFGjcCgGJSQkjdhK5hXPvqhz4kU4Njfrpw40YAUExKCKmb0DWMa1/90IdkanDMTxdq3AIYa1yx1nZptYvIQP/KUnt/xtg/dY3bWNoRwpTbngv2sS7G6I9cbaqh73g4iRCMK84L+3vc0L+y1N6ftdtfEvYdSQnHF0lNrjFWy1hmAm4hGFecF/b3uKF/Zam9P2u3vyTsO5ISji+SmlxjrPaxnPRUSWPM7QCOATgOYKVt51gb87ji5ZUTjCvOAPt73NC/stTen7XbXxL2HUkJxxdJTa4xVvtYThoqOdu47bbW3utyfQ2hkkAdsbFj4up9d+Kag0dw6fnb8KytG6vp+wN3PIBPfukuGAA/8/ynqrRXQzz5mOfTUNtStL32/qzd/pw0+2rsfaexfRptSsWU2gpMr70a0PBOooUiGrexbtxIPhZjkZfWGcAYrBzXHZcMrNr92qtuxKPHV+fXaUvr8JE36rKX8eRpGWr3VPuFyDC18aOxvRptIjLQt6Q0pTRuFsCfGGMOGGMu7zDscmPMfmPM/qNHjyY2h9TGmljk4xbLlcQl7z18H5aPP/ZLEY32Mp48LUPtnmq/EBmmNn40tlejTUQG+pZoJfXG7W9Za58P4FIA/8gY8+LmBdbaq6y1u621u7ds2ZLYHFIba3KVrDfYUEnekj07NmPDenPy3xrtzZUHZqr5ZobaPdV+ITJMbfxobK9Gm4gM9C3RSrZ0AMaYKwB8z1r7O13XMFSStLEYiwxAfVzyHGrc8tejjRIaNzIdpjZ+NLZXo01EBvqWlCS7xs0Y83gA66y1x2Z//xyAd1hrr+26Z8obt/kCsenM0/DADx7tXShqSDK7eKDIz170dG8bat2sNeHCP37afEy/EyIL5xQZM77jO8d84JwrS9fGLWU6gHMAfMoYM6/n6r5N25SZi2AfWT4BC2CdQacYNlQwm1Noe/W+O/H2T90MALj+ttVzaX72oqc721DrgSRNKG4eP20+BkC/EyII11IyZnzHd475wDmnl2QaN2vtYWvtj81+dlprfyNVXbUzF8HO/++zTwwbKpjNKbS95uCR1n+72lDrgSRNKG4eP20+pt8JkYVziowZ3/GdYz5wzukl9eEkxIG5CHbujHU9YthQwWxOoe2l529r/berDbUeSNKE4ubx0+Zj+p0QWTinyJjxHd855gPnnF6yHU7iAjVueTRuLnWEMq/j2MPLOHTkoZMaN9+6+zRuGmPBNdZN+pHyTZuek36fFvR3etjH+qGPwtH4XlPDAWtzxqg1L5KA25cpb9xykTJuuatsyTo1xoKT+pAaFxxfhGOAEM6DMVKLT8eqNS+VgJsoI2XcclfZknVqjAUn9SE1Lji+CMcAIZwHY6QWn05Na57yVEmikHnc8vLKCfG45a6yJev0LStle0m9SI0Lji/CMUAI58EYqcWnXXbWYHsIDJWcICnjfrvKlqxTYyw4qQ+pccHxRTgGCOE8GCO1+JQat0Jw41YHWiZDjsNWXOov3Q8kDxK/MNAs9s41npv11DKPmPC2n5ptJ3VS6y/f+uor8Yt1l+tqm9+12dtGiQTcZIRoEav6JC1PWX/pfiB5kDgUBwBee9WNePT46i/LPn7gLnzkjTrGTa7x3Kzn1166E+/4zCH184gJb/up2XZSJ7UeMNVXX4nD41yuA+o66GPs6xEPJyFeaBF8+iQtT1l/6X4geZA4FGfv4fuwfPyxCAdN4ybXeG7Wc83BI1XMIya87adm20md1HrAVF99JQ6Pc7mutvldm72+cONGvNCSlNEnaXnK+kv3A8mDRILUPTs2Y8N6c/IaTeMm13hu1nPp+duqmEdMeNtPzbaTOpEac7nHbl99KW1xLbvr2VXT/K7NXl+ocZsQzZjftsTBIeWUIlbjFtsO1/4reTBLSlumRi6NWykdXIqx4aKVqGVMjk3jJq1jCb03dZu5Zo6XsWnchtb+2GdQjMat+b6jfV4s2gdAta1d8HCSidOM+b3s4u248guHT37/zlc812vzVjuxMdAx8eK5ko+ntIWk4cAdD6zRwZ22tE6NDs4Xjje9aNGxpB4jXDNJLQyNNQmdtdS7Ry3aZKDuOcwE3BOnGfN77aG713x/zcEjhSwrQ2wMdEy8uHabSTk06+B84XjTixYdS+o6uWaSWhgaaxI6aynbatEmA+Ocw9y4TYRmzO8lO7eu+f7S87cVsqwMsTHQMfHi2m0m5dCsg/OF400vWnQsqevkmklqYWisSeispWyrRZsMjHMOM1RSkJIxvy55QTadeRoOfufBk/HTt959LEjj5lN3jvu7yhkqN5fdktdJ2Vwq951WSuXV6bJFQuOWS7PWR6iOVoocc6pWtORqSj1GNDyfXJ9F0uuy9rHt+t6i8VmVYv74vLMAw7otSd1pKW1yjvmnFWrcElMyjtY1L8jSOgMYg5XjcjZq0RDUHIMN5B0/Ncd8p6BUXp2UpLC7pMYilb2lbZw6U+h/V+2SdE5S7X3rsu6WytM6RGmNqIZ8bznQZEsJqHFLTMk4Wue8IMctloVt1KIhqDkGG8g7fsYY8x2D8/ypqK9S2F1SYxGCS/2lbZw6U+h/V+2SdE5S7X3rsu6WytM6RGmNaMq6NI0bTbZoghs3IUrG0TrnBVlvsEHYRi0agppjsIG842eMMd8xOM+fivoqhd0lNRYhuNRf2sapM4X+d9UuSeck1d63LutuqTytQ5TWiKasS9O40WSLJpxCJY0xBsDPAdhhrX2HMebpALZaa78oaUzNoZKAfo2ba1y0ZN057u8qp7a45tI5nKaMJo2bFBo0bqX7jho3/Uyh/6lxa4caN9n6ayi7ZltyE6VxM8b8AYATAH7CWvujxphNAP7EWvsCSSNr37iVZsoDvI3cCWFr6X8tG+3c9aTyq2SS05TUMj4Bv5eimtolwdjGDf3ajcZnUaqyfcdB6l9IS9jue02KelPcu3j/fGMdusEO3QBPYX3o2rgtOd5/kbX2+caYLwOAtfYBY8xpohaSKKYu4mwS2h+578uNlsNkcteTyq+hh+LkHi+1jE/AT/hfU7skGNu46Sp/an5tQ+OzKFXZvuMg9aFrErb7XpOi3hT3Lt4/PzzGAEGHyIQe8jL19cFV47ZsjFmPVd/AGLMFq/8DR5RAEedaQvsj93250XKYTO56Uvk19FCc3OOllvEJ+An/a2qXBGMbN/RrNxqfRanK9h0Haz5PcOiahO2+16SoN8W9i/fP4/VCD5HxWesl7a8d143bvwPwKQBPMcb8BoAbALwzmVXEG4o41xLaH7nvy42Ww2Ry15PKr6GH4uQeL7WMT8BP+F9TuyQY27ihX7vR+CxKVbbvOFjzeYJD1yRs970mRb0p7l28f76BMLM/fQ+R8VnrJe2vHec8bsaYZwP4O1j10eettbdIGzNljZuvjiakDN/raqcmjVtOn8wT3u7cdhY2Pm6D95jKdZCHtBYvp8atmUC7ZBx/SpF/Tr1LCi1U6LrrOrZ8bEuxfjevDUl2nVqbA7Trk6aifUuxnrrOn5Dycx+k0vX54lh+1taNomMixSFM1Li119m2JlHjtkrQ4STGmCf1FWqtvV/AtpNMdePmq6OJieedemywRnL6pBmb3heT7muXprFV0pZQ3dtQOSFtqFHvkovQddfVvz79k0MDc9nF23HlFw6f/P6dr3iu8+ZNkpLaHG2U1AwD4QmjS/uB61p9SI/BKRCagPsAgP2zP48C+DqA22Z/PyBt5FTx1dHExPNOPTZYIzl90oxNl4wj1zS2StoSqnsbKiekDTXqXXIRuu66+tenf3JoYK49dPea7685eCS6Dgm7cmpztJGrPaFaotJ2l6i/dNvGivQYnDK9Gzdr7d+w1u4A8McAXmatfbK1djOAlwL4oxwGTgFfHU1MPO/UY4M1ktMnzdj0vph0X7s0ja2StoTq3obKCWlDjXqXXISuu67+9emfHBqYS3ZuXfP9pedvi65Dwq6c2hxt5GpPqJaotN0l6i/dtrEiPQanjGsetwPW2l2Nz/a3/RdeDFMNlQRkNG5SdZH85PSJj0YhRax/Lkrakktfl6uMEmXnoEaNW0z7QjRuKahFD5yDXO2R0rj1lZcTrmv1IT0Gx05sAu4/BnA9gP+M1ZM/fx7Ai621PyVppNaNW8mHXU2DOuUhCEN1lnghDqEmf7rSdtCAZDLOrjpy9mHIAQJj9HUMMQchsA9XCTkQoAuJ55oW/3AOlqdUX0u+n/kcipHzPSc3Us/bnIdXjZHYBNyvBfDrWE0JAABfmH02eq7edyfe/qmbAQDX33YvAGTbvNUkkvU59EK6Tk0JSjXWm5LFNs0ToS6vhI0Bl8MgSiRb7fOb6wEWY/B1DCH9wT5cS3ONNQhLegvIPNe0+IdzsDyl+lry/cylDSXec3Ij9bxNNSY4rx3zuFlr77fWvsVae+Hs5y3SJ0pqpSngzinorkm46XPohXSdmhKUaqw3JWvaNEuEKpmMs6uOnH3Y5zfXAyzG4OsYQvqDfbiW5hobs9ZKPNe0+IdzsDyl+lry/cylDSXec3Ij9bxNNSY4rx03bsaY64wx/735k9o4DTQF3DkF3TUJN30OvZCuU1OCUo31pmRNm2aJUEPHgNNhEAWSrfb5zfUAizH4OoaQ/mAfrqW5xprZnyFrrcRzTYt/OAfLU6qvJd/PXNpQ4j0nN1LP21RjgvPaXeO2eDDJGQBeCWDFWvurksaMSeMmdahAM5FvaDmp44EXbd157tlRsd+u9sb0j1aNm6SWa6gsqT7w1bjF6sXa2uJqo7T2ztXmXKG4UuNFKqRl6ACPoTqv3ncnPvaXd+Kcs87Am17yjKgDPly/A2Q0lCkPG6HGrdsOAN6+dB2bUocOdZUfc50UoeuZ1PuKj01AWY3bwe88mKS9pWnOJ4n3LACnzJ2+cofu1bDupCbqcJKOAv/cWvuSaMsW0Lpx80UiBlcqjjdHPLBkHa5ljTHOWVLLNVQWUCbxZQm/TU2XEDteJPrGRy/SlwD7tVfdiEePrz6jTltah4+8MSyJtet3UhrKMa5PNSClpQROXR+bn3UlXpe0M/c4itELTk3TpNUuaVK8z/7aS3fiik8f7FzbS+vatRCagHt+85MWfp5sjPkpAFsHb5woEjG4UnG8OeKBJetwLWuMcc6SWq6hsqak85ucLiF2vAj0jY9epC8B9vLxx36xGJPE2vk7IQ3lGNenGpDSUrp81pV4XdLO3OMoRi84NU2TVrukSfE+e83BI71re4o1eUy4nip5ADh5kNUKgL8C8HqXG40x6wHsB/Bta+1LQ4ysjXkM7vLKiejEuTFlSJaTqw7XsnK0KzeLbVo/+y3T8eNh7XMpq0T/lfDbvM5Hl0/gBMatSwjtV2m/uJQ3dM2eHZuxYb05+VtZ1yTWbeW5fhc773zaT+QJ6feue4Y+u/T8bfjL2+8P8rHW51xXfRLzWdqm0mi1S5oU77OXnr8N+w7f17m2p1iTx4Srxu0Ma+0PG5+dbq19xOHeXwawG8BZQxu3WkMlu2LhY+O9Q2J82+qV1m+E6gFiym/7vvY8KkN9BpTXuIXaE6o1CsVVf3HvsUewZePp0TpMSVz6KkWy9JD7fdY6lzVw6Jo2jVtof7neJ6VbSampe9dnb8G1h+7GJTu34m0//aPO9iz2dUi9KZDSg8foYVz1bG3P1tDnvKs2K4UWua8futojoVlNaXvfvSW0dyXLkmTondJXZx6jcYvR9krqU3MRm4D7S9ba5w991nLfUwF8EMBvAPjlMW7cYmLCQ+voiq0f0oTk1N5NLd7dlxraERpnrkWXIdGW1Lhor7Ro87q0QF1rTqwuxlV7JKW1yTVuY8fiuz57C678wuGT/37zi3cMbt6az4al9QbrFMyB1M+k1Fpx7WPNpz4fm2p4fgH+OtlSaO3P0DU4pS0hz0PXd2htBGncjDFbZydKPs4Yc6Ex5vmznx8HcKZDvf8WwK8CONFTx+XGmP3GmP1Hjx51KFIXMTHhoXV0xdYPaUJSxCrnjs0fS1x5De0IjTPXosvovEZRzLyL9kqLNq9L99O15sTqYly1Rz5l+rYvBbFj8dpDd/f+u6vORT+tKJkDqZ9JUj6Vfs6rXyOVrfGhDL0TaUFrf4auwSltCXkeur5D18LQ4SQ/BeB3ADwVwO8C+Dezn18G8Pa+G40xLwVwj7X2QN911tqrrLW7rbW7t2zZ4my4FuaxuClzyDTLuvT8bZ11blhvTt7nmnNKor2h16WqXzs1tGONjR75XHK3zaW+0Lakps/2+Xdacga12dq35nj7xWG9ii3Tt30piB2Ll+zc2vvvrjoX/bSkZA6kfiZJ+VT6Oa9+jVS2xocy9E6kBa39GboGp7Ql5Hno+g5dC66hkq+01n7Sq2BjfhPAL2D1MJMzAJwF4I+stT/fdU+NoZJAeN6TmDp8NG4utsbYEntdqvq1U0M7UmjcUtsprddLjZTGLQe+et5QfdFQfTFl+rYvBdS4rbUr5TMptUZM+1jzqc/HphqeX0BajZskWvszp5bR1RZq3Ho2bsaYn7fW/mdjzP8F4JQLrbW/61j5jwP4ldo1brkf6qFJjEPKB+Qf4iUnRg5Bt29y4JL4vNjk3NznKi/lWJTqv9CXT00vjpKb49S/+NGyGY5FwwtIjN8lNgsShyX4oKHPpW3xWcdi1zzptS5kUzp0jwYf5x7XTRYP0HnW1o1F+0Pje0lqujZuQ+kAHj/78wnyJtVFbuF6nwAzxpa28lMc1lBSbJtakA6sLmhv/9TNsysfxHW33oOPXn6xysXAR7wvIdqV9n1seSnHolT/hR6woOlwhNhDN1Lbt1iulgNfYtFwqEGM3yUOxAjpA4lnqIaDDVI96/rWsdg1T3qtCzl4ZchmDT7OPa6bLL7jXH/bvdiw3uD4CVukPyTapcGnUvRq3Ky1/2H2579q+3GtxFr7Z7XncMstXO8TYMbY0lp+AqF6SbFtakE6AFxz8Miaa5ePW7UC1752pBDtSvs+tryUY1Gq/3x8JLEWpOiTNWVKJpAX9JnL+loTJdfZVhs8/e5jf9e1IX0g8QzVMHZSPev61rHYNU96rXPpA1+bNfg497hu0vaOU6o/JNqlwadSDB1OAgAwxmwxxrzdGHOVMeb985/Uxmkit3C9T4AZY0tr+QmE6iXFtqkF6QBw6fnb1ly7Yb1RK3Dta0cK0a6072PLSzkWpfrPx0cSa0GKPllTZuSakspnLutrTZRcZ1tt8PS7j/1d14b0gcQzVMPBBqmedX3rWOyaJ73WufSBr80afJx7XDdpe8cp1R8S7dLgUylcDyf5CwDXAzgA4Pj8c98DS4bQrnGbC8IveNoT8cxzNiaJ4V+8PkTj5hLHPaRxW/x7Ks1OTL902RWrr/EV4c41bqcvrcMzz9k4KHx28VmI3S519R1m4NOvUvqBvjKGkieH9FEqzV1bwubcGre5DmHntrOw8XEbvOZTqGi/y1chfuqbF++65hZ86/4f4OUX/MiacSshUpdMau57nY/dQ37Krb0eGp8AvNrYlojaZXz5aJx813aXvih50MvQPHNpb8y633Wtq8/77HdZ65pJy12SmTev6VvDF230nauh88BnjLr4ZqjeofF/9b478f4bDgPG4O8++yl46JGVLAe8SOlY28ZoDYfULBKbgPsma+0FKQxbRPPGba2mCTAATt8gG8Mfi4TeJoXmyaUe335s2hVrZ2qtkEt8vpTWsFnXZRdvX5Ow952veG7nwy2kfZLx50PaIw1x6pp0UjG2xI55ifb32dCVaDpX/8fO7a7rfOw+cIeOBMKx67VP3wEQ9a92PWdM3a5rJIDkz/XUmixJjdvivX3zK9S++ZhYXnEbwxL94HPN0PfN991c+rYSmk3NBCXgXuAzxpifFrapKprxvhbyMfyxxMaeS5URYqtvPzbtirUz1E+u9znF5wtpDZt1NRP0NsdybPsk48+HtEca4tRdbdVuS+yYl2h/nw1diaZz9X/s3O66zsfuvYd1JBCOXa99+k7avynWjBTrdlDdjmtkjud6SD/73OP0DPV8Rg3Nr2D7ZmPCdQxL9IPPNUPfl9K3Sc3VXO+xpXDduL0Fq5u3h40xDxljjhljHkppmDaa8b4G8jH8scTGnkuVEWKrbz827Yq1M9RPrvc5xecLaQ2bdTUT9DbHcmz7JOPPh7RHGuLUXW3VbkvsmJdof58NXYmmc/V/7Nzuus7H7j07dCQQjl2vffpO2r8p1owU63ZQ3Y5rZI7nekg/+9zj9Az1fEYNza9g+2ZjwnUMS/SDzzVD35fSt0nN1VzvsaVwCpXMheZQSSBOS5Ir/l1SaxRTRkg9sXbF2plKyzd0XQ6Nm4sGwNdu32tc64mNzc9BjMZKky2xY16i/X02dGkzc/V/DRq3XMSu1z7XSfs3xZqRYt0OqTtU05e6T1KsQTH6p757pTSk0hq3mGtDxsgipXK4SY3LXO+xKYnVuD2/5eMHAdxhrV0RsA+A/o1bG7GTw/WamLJiH4LS9g29hEiIqH1s9l20cy0AMQt5179DhdYSGyif8dtlZ+gYS/UwSGF/Srt8/J7ypdm1jqFfOGh8GEvb5LNxK7XJK7mJkSD25TzkmerjJ23j3LUPJH+hV/qXM1L1dB0K1/ackDxAKfT62PJCf2mceuxom1NtxG7c9gJ4PoC5WvG5AL4CYDOAN1tr/0TCyNo2brECUNdrYsoC4oTe0vYNCe2HRKUuh5LM2+xis68w2bXsWGLEyl195HsoQuzY9m3P/JouO0PH2BUvkxEmu/ZJjP0hSPpqqA0StgFua1JTIN88VCdVf8YgbZPLmhlyrSSLbc59UIcEoWttaFt9/aRtnLv2gYTdqdueq2/b1tWuvvNZg1M8h0Pa5XrIievBaKnHjrY51UXs4SS3A7jQWrvLWrsLwAUADgL4uwB+W8rI2ogVgLpeE1PW/LNQobe0fUNC+2Y5Q4eQtIlOfWz2FSa7lh2LTz2ufeZ7KELs2Pa9fmisho4xKWGya5/E2J/SLh+/pzwYwrWOpkC++e9cc9EHaZtc1syQayVZ0+bMB3VIELzWBrbV10/axrlrH0jYnbrtuZ/na9a8jr7zWYNTPIdD2hW6hqe0s68MbXPKF9eN27OttYfm/7DWfg2rG7nDPfeMnlgBqOs1MWXNPwsVekvbNyS0b5YzdAhJm+jUx2ZfYbJr2bH41OPaZ76HIsSObd/rh8Zq6BiTEia79kmM/Snt8hlLKQ+GcK2jKZBv/jvXXPRB2iaXNTPkWknWtDnzQR0SBK+1gW319ZO2ce7aBxJ2p2577uf5mjWvo+981uAUz+GQdrkecuJ6MFrqsaNtTvniGir5MQD3A/jo7KNXA3gygF8AcIO19gUSxmgNlZzHo9977BFs2Xg6dp57dm+McjP+20XT5aLJidX3pNK4LZa/mMiyq2+G2hrz73mdPoLgsWvcAODKP/8m7nnoh7h4x+bWg3WG2gm4JXmNbU/XWPWN+S+tcTv28DIOHXnolLh+LRq3tn/P+2znuWfj4HcexL3HHgEAbNl4epReqm9sLfbTXAC/6OPf/9Ov44u3348Xbn8SPvT6i04pO+bQHUm6+q+v73zW866y+vzokhg6Vvs4L99Ht5N63WyzeajOq/fdiY/95Z0456wz8KaXPMP5mbH4LuCzBh57eBk3Hr7vZH1AdxJqYPgdwqW9Q2uADy4atxC9ZYnnba7nBdA9V/qed23JwbvqiH0Oh5Q19P40f/d49QuePrhGS+hHXdszBY3b4wD8QwAvwupJ+DcA+PcAfgjgTGvt9ySM1Lhxa8ajz3FNEOyS+C/knj57c8budsW7u8b/S9pbS9xyTkJ0hdqSXpeu3xft9ratN1d8+uDJMbK03mCdR/JYCTvaEtYaAIvL7jwBd1c7SvVzc44tzf43ZSVizrlck/r7NltcNDp99+f0lWuduTRnXdqlPt22VB9LvmO4ENJHpeZzqr7xGX8atFjSernS+rsxEKVxs9Y+bK39N9baV1hrX26t/R1r7Q+stSekNm1aacajz3HVk7joa0Lu6bM3Z+xuZ7y7Y/y/pL21xy2nwEVPMdRvpfu1dP2+aLe3bb1ZHCMrnsljRexoqbO57DYTcmvp5+YcWzluT27agHT659Tft13rotHpuz+nr1zrdFkjQ8rtuq85r3p120J9LPmO4dPWVMm4JUnVNz7jT4MWy7cu6fcGLet5DTht3IwxzzTGfMIY8zVjzOH5T2rjNNCMR5/jmiDYRV8Tck+fvTljdzvj3R3j/yXtrT1uOQUueoqhfivdr6Xr90W7vW3rzeIYWfJMHitiR0udzWW3mZBbSz8359jSenPyf92AdPrn1N+3Xeui0em7P6evXOt0WSNDyu26rzmvenXbQn0s+Y7h09ZUybglSdU3PuNPgxbLty7p9wYt63kNuIZK3gDg1wH8HoCXAfil2b2/LmmMxlBJoF/jFqIn6aojRfx5rlCDtrjkxb+HxuXH2DL1/2af46I18NVI5aZ0/b5ot9dVG5UzwXVbnZ87dHdrAu6udpSirf9i55zLNam/b7vWRc8W0pYUuNbpq8cKbcuQhrdL4ybRx5LvGD5t9Sm/1HxO1Tc+40+DFktaty5d3tSI1bgdsNbuMsbcbK197uyz6621f1vSSI0bt9AJlWIiuk6Srgfr4t81b5C6ypZaXH2Evy72xQiqSyxssf3o6p/YdoSI28dC15hKJaJ3qb/v87bvgfBf3HS94Ia2pevAjj4bfedw6IYmBC2HsgB+G9DUvwjIRa0vmKnePWLLl7DD5Tlcyj6fdbHEXOlaT2LebXypaaOfi66N25Lj/T80xqwDcJsx5h8D+DaAp0gaqJFQ0WgKsamrELRLPC6ZGDWliLSrbCkBcbOf5nz8wF1OyWpd7ADcE4DnFu/G9qOrf2IPKmgeGODqnzHQ1j8ATukPqcTiLvW7+Hjx+5jDiQC3BN2ubVnst4/t/xbWOayJvgcLhB7aEcJiUtvrb7sXAIpt3lznsmRS99LUeohCqneP2PIl7O9bR0ofyOKzLpaYK13riUufpn7mSN8zFlzzuL0VwJkA/imAXQB+HsDrEtmkhlDRaAqxqasQtFM8LpgYNaWItKvs5uehAuJmP80JFVDHJgDPLd6N7UdX//SV43Lt3sNlkglroK1/2voj1QEDoT5e833E4USta1nEeG87eEUqefDguptg3IYmtU2B61yW8qcGUj7/UpLq3SO2fFd836uk3hli7fNZF0vMla71xKVPUz9zpO8ZC64bNwvgPwH4NIDdAP4mgD9MZZQWQkWjKcSmrkLQTvG4YGLUlCLSrrKbn4cKiJv9NCdUQB2bADy3eDe2H13901eOy7V7dpRJJqyBtv5p649UBwyE+njN9xGHE7WuZRHjve3gFankwYPrboJxG5rUNgWuc1nKnxpI+fxLSap3j9jyXfF9r5J6Z4i1z2ddLDFXutYTlz5N/cyRvmcsuGrcbgXwzwDcDODE/HNr7R2SxlDj1n2NS4LQ+T23/fUx3PSt7+KSnVvxkzu3emtOXPogVH/k0i+uGr13ffYWXHvoblzwtCfimedsdO5PACcPmwHWJhh2iUV3EZLP2zAUp+47FrpsGiq3rwzfOPYQ/VPTty794pt4uESsewrNYd9nvgdghNrfNV4W/y6lcWvTVzR1qH3JwIfasOnM03Ddrffgr45+Dzu2PAE//qyndGrRQvSvbWPVte2hzNe++aEtJbUeuTRuEnPe1Q4fLaekVqyEBsulTtcE5bVo3Pr8H/O+U4vGrWlPV5JsaY1b7HtxW3ldfgl5z9ZG7OEkN1hrX5TEsgU0btxK0Izdvezi7bjyC49lX3jnK57bqWlYjFceujbElhBN01BZXd8B3fHUzXYaAKdviEtwHqLRCalHov9T2iEdO+4bt++jpcqhXRhqV27/5qy7pG5l/v2iRg1YmyDZVVdyAv0JsRdt8dWX5NZaaBn/OZFos6t/U/nTd13TkqS8OQe75k8thLyHjLX9qfW4XfVK1RXzflqDRi4qATeAXzfGvNcY81pjzM/Mf4RtJDOasbvNxLN9mgZp/UOo3sWnrK7v+q5vtssiPsG5Tyy6S/skYuol+llKxxaKb9z+UP25tQtdSPRTaBm56k6tI3Dx9XIjE7fLOtScy0MJsRfv8dWX5NZaaBn/OZFos6t/U/nTd13L4UfXNWBMeuOQ95DRtj+xHrezXqG6Yt5PS8w3KVw3br8E4AIAl2A1j9vLALw0kU2Tpxm720w826dpkNY/hOpdfMrq+q7v+ma7DOITnPvEoru0TyKmXqKfpXRsofjG7Q/Vn1u70IVEP4WWkavu1DoCF18vatQAt3WoOZeHEmIv3uOrL8mttdAy/nMi0WZX/6byp++6lsOPrmvAmPTGIe8ho21/Yj1uZ71CdcW8n5aYb1K4hkqezN+WEoZKdmt7fPL2uMajD9kxFJMtob1a/KwtprsvZnwel33xjs3Y+LgNrdqHuV1DCU9D2z0Upy2RbymkjFg9iFTsd4i2LUSDEqK7kyCFxk1Kh+NTV6q57WNTW3tf9759+B/fuBfr1xlsO/tx2LDeYMeWJ6xZ1/r0eHP9w85zzz6pWdt57tmd60zfOG2zc17fohavTYe32NY+TYarZlGLxrPLnjZi18LUGrcQHZLr/fPvjj28jENHHursg5R+dNGC9um+JHNqhmhJJXF9hgMQt9NH57hYv2vC9rb1qml7zFh3HaM++jhfPdqizV32t41ZbevmELEatz8E8HvW2q+lMG7O1DduUrHksbG7vvfH1JdSK5ciftvH3pIaKA3kHEc5ystBqM7Kp+yh/pCek6G2/f6ffh1fmOUVatK2NrrooIBTdbNtn7no35p5Mod0dH3re806IhffS2uvpYld1wF0Pnfm47BkHjtN+rrmXJqjbczn1GS1fT8fR8sra8dNc13rW+di+zhUWyr9nhTyXjcGLXCsxu1FAG4yxtxqjPmqMeZmY8xXZU0kUrHUsbG7vvfH1Bdyr3OMeoL4bR97JWKoa47DzjmOcpSXA1cdTkzZKfIxpVhzvnj7/Z3Xh2ofuzQsrnqIznxtDjq6vvVdau0vgUv/aco910bsut733JmPw5y5uVxs9vk+hS3N/zLQNuZT9ImXH2bjqDlumuta3zoX28cua2pIO32vDXmvC7W9Blw3bpcAeCaA/w2P6dtelsqoqSIVSx0bu+t7f0x9Ifc6x6gniN/2sVcihrrmOOyc4yhHeTlw1eHElJ0iH1OKNeeF25/UeX2o9rFLw+Kqh+jLkzmko+tb36XW/hK49J+m3HNtxK7rfc+d+TgsmcduqH0518rmXJqjbcyn6BMvP8zGUXPcNNe1vnUuto9d1tSQdvpeG/JeF2p7DTiFSuZi6qGSgFwseWzsbk5NlIRWp+s7QD6fko+9KTRQNVFKW5eqvBy4atxiyg7RJ0iV7XP/6963D3/xzXvxhDM24DW7n4aHHlnpXRt9tB8u2ou2svvyTA6t3VIaN23k0LilJnZdd9G4pczNNYSPtiq1baU1bq6k6BMNGrcYe1M8Q2I0bjVpgX2J0rjlQuPGrc3RzQeQz6LcTJwaa0ssqTc5XXWlEvfmfvlpWxik6g/trxoWJ9d+qmkj5yP+v/XuY8leYkNE3D6fD33X1WafTdIc35f9ruvnn+/cdtaag4xC2iS9VsYI9n3sz4Gv331f2Eqvg6FtkKgPiB9zKTaQKZ7rQPvGV/oXw10v9l2/nHHZJHT1bYn5ObR+hm7CQjZEOa8ZC10bt6USxtRCm2jy1ruPnRRZX3/bvbjzvu/jAzfe7iQ8ftdnbzmZSHv+p+vmLbVQNnUiRin7+8ppCvw/fuCupGLnNvHrFZ8+KFJ/aH/lFJmH4uon6bak7JuustsOtVhngJUTq/ddPzt8Q2rz1jWnu+b3kN1DBwD1Cezb1kRg+CCQOYsHWrj0U9f1zc8NgNM3hLVJeq0M6X+fcnLi6/e2z6TWgNTPy5DxHFOfxJjr84/0+I0pp+twF+nDz9qe381DZHwPwujq2xLzc2j9dLWpq598Dv1wqUvqmingqnGbJG2iyaao+tpDdzsLj5uJtJv/9rUllhDBp0hdEeX3lbP3cF6Bf5v4Var+0P5KMU6kcfWTdFtS9k1X2fPPF9eH+aZtjuRBDZ1zumN+D9ntmqS2zYa2NdHHB74HWnRd3/zcIrxN0mtlSP/7lJMTX7+72KxpHQxtg0h9AmOuzz8iNko97zoOd5E+/Kzt+d17AFFP3UN9W2J+Dq2frjZ19ZOPf6TmuoZ1TgPcuPXQJppsiqov2bnVWXjcTKTd/LevLbGECD5F6ooov6+cPTvyCvzbxK9S9Yf2V4pxIo2rn6TbkrJvusqef764Piw1Vl3Jgxo653TH/B6y2zVJbZsNbWuijw98D7Tour75uUF4m6TXypD+9yknJ75+d7FZ0zoY2gaR+gTGXJ9/RGyUet51HO4iffhZ2/O79wCinrqH+rbE/BxaP11t6uonH/9IzXUN65wGqHEboC1GuE/jNiQC9dG4NWN5Y3QaXbH5iwlqr7v1Htzz0A/x6hc8vTckybeOZtx6bHy9hBA3ld6u2aehCZSbGjDAXd9Q6hAAn7j3eZL4echD15gL8VOsvsrnGhctRCmNW3MuLNo3/3vffMmlcevTO77uffvwxdvvxwu3Pwkfev1FgzZ0jf35unvB056IZ56zsdVHPgeLLCb0dj0ope+ztrpTaNx8tSq+Nsz7f/PjT8N93390jR9C9WGh61mofse3zBS6sbb6gGG9qmt5Q4ddhNoYo0GsSePWt/Z0jf0UCcyH+iNE4+bymes64avRpA5uLTycJICSOgOX+GtX7clQbL5L8tihNrrG/7d9lkJ7UFoP4qOZcY29d7G3VAy4j83za3+4vDZmUCIpb6yOwyfOvmRC3SFcxpcG+5t6x8W1Z1ETDABvfvGOk7/s8hnnMXOxq5wh7YuPHirXnPXVqvhqizRrqFKXn3Pdla6rxDOjZq1S39yO1fWG1p9CV6l97mqrNwWxCbgnSUmdgVP8taNOYzA23yF5bHQdCTUBJf0UWkdI7H2upN8h+Ng8v7aJhNZrSGvg246+OHtJrYg0LuNLg/17D3frHfs0wT7jPGYudl43oH3xWftyzdmhORqiXWkrX6OGKnX5Oddd6bpKPDNKPack6Jvbsbre0PqlSPkOW+pdpKax5UqyjZsx5gxjzBeNMV8xxhwyxvyrVHWloqTOwCn+2lGnMRib75A8NrqOhJqAkn4KrSMk9j5X0u8QfGyeX2saZUhovYa0Br7t6IuzL5lQdwiX8aXB/j07uvWOfZpgn3EeMxc7rxvQvvisfbnm7NAcDdGutJWvUUOVuvyc6650XSWeGaWeUxL0ze1YXW9o/VKkfIct9S5S09hyJVmopDHGAHi8tfZ7xpgNAG4A8BZr7d6ue7SFSgJpdAahdbvE0sfoz1xjsEPrCNVuhPSV9PUpbAqJvXext1R8t4/N8++OPbyMQ0ceEtV6xWpOQjVu2sIxXPpfg/19+o8+TbDknAjRQwH9mg2ftS/XnE2tcUsxrlL3jVT5Oddd6bpKPDNq1iH5ah1r8lfKd9hc1Dy2FimqcTPGnInVjds/sNbu67pO48bNheYgWRSEPmvrxsENU+kBFvuwHio35iE+ZEuf+FaqXS7C4tjNbAiaXr5Tkkv8n6NMjb9g8K0vdsPkc7/LYU4udUgc2CPp1xL0rYdA3MEBQ3M01TNmyDaN82foPsD9l6g1EJsE2uX5mwsJ/8ZoxrSvMVJI9lfN70dFEnAbY9YDOADgfwHwnr5NW600hZCXXbz9pKD++tvuxYb1BsdP2M5DQVImvQ6xPybZZVu5MUL1IZFpX4JJqXY1D09oSxbdZafPAQW+aDpgIiUp25lCxOxyKIimQ3R86/OxKXb8Lx5OMv+z7X/dhuporsmAf8JzSb+WoG89dHkOubTfNfmw1DNmyDaN82fovuZBYW3Pm5qITQLt8vzNhYR/Q8dhDWuMFJL9Ndb3o6SHk1hrj1trLwDwVAAvNMac37zGGHO5MWa/MWb/0aNHU5qThKYQsimoXz5uew8FSZn0OsT+mGSXbeXGCNWHRKZ9CSal2rX38HCy6C47+0TMGvq3BlK2M4WIua9M3/pyi6xd6vOxKXb89x1O4lNH876QQ3Ak/VqC3vXQ4Tnk0n7X5MNSz5gh2zTOn8H7PA4Kq4HYJNAuz99ciPg30P4a1hgpJPtrrO9HWU6VtNZ+F8CfAbik5burrLW7rbW7t2zZksMcUZpCyKagfsN603soSMqk1yH2xyS7bCs3Rqg+JDLtSzAp1a49O4aTRXfZ2Sdi1tC/NZCynSlEzH1l+taXW2TtUp+PTbHjv+9wEp86mveFHIIj6dcS9K6HDs8hl/a7Jh+WesYM2aZx/gze53FQWA3EJoF2ef7mQsS/gfbXsMZIIdlfY30/Snk4yRYAy9ba7xpjHgfgTwD8lrX2M1331KRxW4yfbSZiXdRl/OTOrYMaN5f47RhdyVCMeDMWOOTwk7bvAbeE26HakQN3PIAr//ybnUnDJXQu8777xl8fwyMrJ3oTRbsm0c2hcfPxqaRNrjb71BMbq55bm9SnHfIdDzH2DR3s4TJe2+aMj00+97d9/taPfhmfvfkINqxfh1/Yc16rzm0qGrfYOpo+B/wSOvc9R+bf3XvsEWzZePpgMveQ55nUoVySuLTLVcM1f8ZcvGMzHnpkRfSdIOR6KcaocfN9Fkn0/dDZCVJ1Sr2z+Opeu9oau16PUeOWcuP2PAAfBLAeq/+z91+ste/ou6eWjVtf/Czgp2FKrSv5tZfuxBWfPtia4FbajpDY5Fxx4zG2ucRJa4pBb9ptgF77c9leoo9K+kV6bvswlLy6K+n1Iov6FMA/QbqEHu7/vPIvsBAptaYdU0KDTqaEzk+zVkXymfnD5bV5LYfmmna97Jgp/RyL1aS61hOjy/fVvdamGc5F9gTc1tqvWmsvtNY+z1p7/tCmrSb64mdT6FhidCXXHDziFCMuYUdIbHKuuPEY21zipDXFoDftHrI/l+0l+qikX6Tntg9Dyatd1oQhfcoQMevWfC1d3LQ12zElNOhk+u5PNY41a1Ukn5lNhuZajucfaaf4cyxSk+pcT+A77WI5rrrXHGvJmMiicRsbffGzKXQsMbqSS8/f5hQjLmFHSGxyrrjxGNtc4qQ1xaA37Z57v8v+XLaX6KOSfpGe2z4MJa92WROG9ClDxKxb87V0wcxT2jElNOhk+u5PNY41a1Ukn5mNYT4413I8/0g7xZ9jkZpU53oC32kXy3HVveZYS8ZEljxurtQSKgmsjZ+97tZ71uisXOLtAX99ga8mYDHu3ldH12XfkK0xsdC+9yxqKnaee/bJOOYuG1Prq/r0Dl02DdUbo2dZtDtW4yYRL/6uz96C//emb+PpTzoT//zSHx0M7Qkd7zFlDeGrt3CpW6KtV++7Ex/7yztxzlln4E0vecYpWltXjVsTnzyJbd+1aYB99EhX77sTv/u5W/H9R47jp3aeg3/7mguzanVc53BJjZuPHkR6XfHVwYTqjGPXnlRaVgl9+vz7Yw8v49CRhwb7pqTOShqNNrURo9GMravtcyDsndH3eTSvx0Wn1rQpRuMm9R5Zy/jqomgCbldq2rjNcdWB+MYol0S7rX0aLi32xvShthhvCZ1Jn+aqq75QbVTKF2cXbVgqutraXIM2rDf46OUXJ7XLVaPQHPs+Obza+vuKl8nmAHNtY98cLjlfY3WIIYS2t4StQBr/lPK5tmdDDLW0JfczJlVd0rpITe+JUto8bWTXuE0FVx2Ib4xySbTb2qvhUmJvTB9qi/GW0Jm45ORq1heqjUqFqzYsZf1tbW2uOcvHbXK7nDUKjbHvk8Orrb+lc4A5t7FnDpecr7E6xBBC21vCViB/vsaUaHs2xFBLW3I/Y1LV5Vv20PWa3hOltHm1wI1bJK46EN8Y5ZJot7VXw6XE3pg+1BbjLaEzccnJ1awvVBuVCldtWMr629raXHM2rDfJ7XLWKDTGvk8Or7b+ls4B5tzGnjlccr7G6hBDCG1vCVuB/PkaU6Lt2RBDLW3J/YxJVZe0LlLTe6KUNq8WGCopgGvcfozmKTfabe3TcAE67C2pcZNGSuPWpbnqqk9C4yaJr8YtRf2uGrdStjS/A4Zz2fXV0ezvqWrcupDIeeRLaHtL2Aqk07iV8Lm2Z0MMtbSl1JqTSz8Xer2m90Rq3ApR68YtN6VeaqXFpLltjL3f51CFWHuk+sS1HI0LnOvhF9JJuSWuL430oSgx9QP9Avc+XDf7qdamoXqkRPRTQeIXQFI2uB6iUNK2UuWXPGxH8t6U7wMhpLIn1ftCynW1xDpY29rbtXFbKmEMCafUwQ2uhxH41pPLxlRtTNEvUn3iWo5GkbhLv4YeliIt0taGi725RPBz0fryir+vFg+0mf/ZtnlLtTYN1eNz0IqLrWNH4pAjKRuaPtSQ4Dv12ND0PMjxzpDyfSCEVPakel/oW9+09a3WOlNBjVtllDq4wfkwAs96ctkYe3/Xdyn6RapPXMvRKOJ16dfQw1KkRdracLE3mwh+JloP8ZXrgTap1qahekIOSqltLEkSO28lbWj6UEOC79RjQ9PzIMc7Q8r3gRBS2ZPqfaFvfdPWt1rrTAU3bpVR6uAG58MIPOvJZWPs/V3fpegXqT5xLUejiNelX0MPS5EWaWvDxd5sIviZaD3EV64H2qRam4bqCTkopbaxJEnsvJW0oelDDQm+U48NTc+DHO8MKd8HQkhlT6r3hb71TVvfaq0zFdS4RVBCmOl7CIGklieVlmt++MDG05ecEpCG2u9yX5/Ooeua1Bo3IDzpuWv9zUMDNOg+JDRuUpqA0LkuqQOQrldCn+XSv3O7Q8bS6963D1+8/X68cPuT8KHXX+Td3rbDTXzb4tr+eTtL6m762huDhPZ30bad5559ymFSIYl6Q+yf19f0YVuy+Ob9qZOAa9G45ThEpq+OofpzadzmduzcdhYeemTl5Ni49e5jJ+171taNzmMj1u7cGrc+jXHIWrPo13m/udjsMx5d5nqJ95kQeDiJMG06jtTJB1MnMM0dA7xY3zoDrJx47LtcyVmbduTUDUrYJmGPZt1HDCl8lXOslFhjXO0J1Xj5ELve+WjfYtuSa10YalOKBPES2t+u/h3SvaVY33w1URLavJr0NaU1bqUStXfZOPf7nPXrgOML7yob1husHLdiz0otY2XID752hvrV5z5tz8xYmIBbmDYdR+rY2dQJTHPHAC/Wt7hpA/IlZ23aoSnhs0t9EvZo1n3EkMJXOcdKiTXG2Z5AjZcPseudj/Ytti251oWhNqVIEO/Tj126m67+nV/XtdakWN98NVFDNkraoIEctvbVUSpRe5Om3+ccb7yrLM82bYDMs1LLWBnyg6+doX71uU/bMzMV3LgF0qbjSB07mzqBae4Y4MX6lhojMVdy1qYdmhI+u9QnYY9m3UcMKXyVc6yUWGOc7QnUePkQu975aN9i25JrXRhqU4oE8T792KW76erf+XVda02K9c1XEzVko6QNGshha18dpRK1N2n6fc76xgcb1hvRZ6WWsTLkB187Q/3qc5+2Z2YqGCoZQSmNW8rY89wxwIv1LcaN5w6N8Gl3yT7y1biF1qFB4yZBCl/lHCsl1hhXe3JoBmLXOx9tVkrtmSQ1a9z6rsulcYvRFqXWuGkhh619dZRK1N6kS/8YqnHzrbf0WJHSGrqWJ3GftmdmDNS4FUJSNJrigdxWT85Bn/vFuvQDKTcSB1CktCFVX5V6CR/TQ6MNiUMBUr/At10n7RdNc5yMB58xfuWffxP3PPRDvPoFTy+mARvDL/dcSDnfY3+BVPp9R4qabM0FE3AXQEpk2hSdf/zAXSKi8z57cwg7cx8eUVp0nZsQW6TtL+EPyUNBfO4fmzC6SewhD33fSR1S0XYdAFG/aJrjZDz4jPFX/4e/OKkL/8pdqwc35Ny8SRzYUgsp53vsIUml33ekqMlWDVDjlhApkWkq0XlbPTmFnbkPjygtus5NiC3S9pfwR2y5ofePXRgde8hD33dDB0DE1C3tF01znIwHnzG+UvAwr7kNYzjAyoWU871Ztu8hSaXfd6SoyVYNcOOWECmRaSrReVs9OYWduQ+PKC26zk2ILdL2l/BHbLmh949dGB17yEPfd0MHQMTULe0XTXOcjAefMb5U8DCvuQ1jOMDKhZTzvVm27yFJpd93pKjJVg1Q4xaAT0xyTNxuU5cRo3EL0YcA6ZMXpjjsIOZAAlctTIkDArpoE+7G6A+monHr66MxadwWx+Ni4mMJ7Z/rwRSLn837uyuRffPzxSS4Gx+3wbvu5nyUOARJQjta+uCcoes1juWx4zOm+zRuEmtu18EcXfM5x/PCx35JUj7Tm+3pe567vm+keLa52D5kVx8xB9KMVR/Hw0mEiI1JDq2ndELHMerD+nzZp4UZsjNnO9qSUz5r60bGiw8wFY1GUx8LyLXVZ5yHrptdyVdDNW+pk4Z3kcMO33XHVR8zRr1mTaTyq8u8aSafPm1JJql7iO0+9ue2K2U9ALz6pG9tyfVe2XzuDI2bGLvGrI9jAm4hYmOSQ+uJKVeirDHqw3p92aOFGbIzZzvaklMyXnyYqWg09h5eq48F5NrqM85C182u5KuhmrfUScO7yGGH77x31seMUK9ZE8n86jBvmr/WT+1/Tc/WEvUOanQd+qRvbcn1Xtl87vg+H3zsmuL7DjdunsTGJIfWE1OuRFlj1If1+rJHCzNkZ852tCWnZLz4MFPRaOzZsVYfC8i11Wecha6bXclXQzVvqZOGd5HDDt9576yPGaFesyaS+dVh3jRfEFP7X9OztUS9gxpdhz7pW1tyvVc2nzu+zwcfu6b4vsNQyQAktApdZSzGJANyuoLU2rGUdaSMX27GVS/+e55YM0QnIKGRcr3vXZ+9BdceuhuX7NyKt/30j54s55Nfugv3HnsEWzaeHhyXn7Lvc48XlzmXK+wmp4ZDUuMWY2voutmlffD1c6y2LHa85NS4udo4Jo3bWHUuQJh2sUtf5Ko9mtd57OFl3Hj4Ppxz1hl400uekXzN0qRxKzEHQt6BfNaW1L5aXIOuu/Ue51yDMXaNde5T46aIrpjkmvU2qeKMU8Yva9G/dNkTk3fNN8Zcyp6cZfvoAbTEwWvVcIyF1OtQDWv0FMfQFNvch4TOzfXasfb9WNuVEmpjZaHGTRFdMck1621SxRmnjF/Won/pssel/q57fGPMpezJWbaPHkBLHLxWDcdYSL0O1bBGT3EMTbHNfUjo3FyvHWvfj7VdKVnTZ9TGJoMbtwJ0xSTXrLdJFWecMn5Zi/6lyx6X+rvu8Y0xl7InZ9k+egAtcfBaNRxjIfU6VMMaPcUxNMU29yGhc3O9dqx9P9Z2pWRNn1EbmwyGShZCi95GklRxxjl1VqVjpSX0koufx+aeGbPGrRSaNBxjJPU6VMMaPcUxNMU29yGlc3K5dqx9P9Z2paQ2baxmqHGbMLmEqrFoEKemOBgiVd83N2aA3yKpyfeaqKlfusaAhs2Faz9q7W/fX+r4vuACaV9qNP7SRbr9uX9ZmGOsaps3WudnSbT5SALfdxWpOaJxndJC18ZtqYQxJB9SIuWSdqa8N0U5LuXF1tU8fORj+7+FdR5CYE2+10RN/dI1BpZXyh+g4dqPWvvb9+Ai30McUgv3NR4sJN3+3Adi5Rir2uaN1vlZEm0+ksD3XQWQOQRM4zpVA9S4jRwpkXJJO1Pem6Icl/Ji62oePrLiKQTW5HtN1NQvXWNAwwEarv2otb+bdg0dXOR9iENi4b7Gg4Wk25/6IJoSBxxpmzda52dJtPlIAt93Fak5onGdqgFu3EaOlEi5pJ0p701Rjkt5sXU1Dx9Z8hQCa/K9Jmrql64xoOEADdd+1NrfTbuGDi7yPsQhsXBf48FC0u1PfRBNiQOOtM0brfOzJNp8JIHvu4rUHNG4TtUANW4ToCuBLZBX06X13lQ2LOqP2hIfS2rcdp57Ng5+50GvhNs+ceupkD6MRdsYkKTLX1o0bqEH/Wg9TMZXQ+qjcdt05mk4+J0How4OcrFfm3aEGrfwupvf+8xxrWuaRiQPYyl9vkDfu1+sraU1bkPPm+b6XePY5eEkEyVnnO+YY4qH6NPE5NKzxCYHzu2/kPpK6k+GmHK8fk59aC5ya6gIGUJyzSSnMibt4NX77sTbP3XzyX+/8xXPHdy81cKQBnlIk1wLTMA9UXLG+Y45pniIXk1MJj1LrLYpt/9C6iupPwm1TXvZEuTUh+Yit4aKkCEk10xyKmPSDl5z8Ejvv2tmSIM8pEmunWQbN2PM04wx1xljbjHGHDLGvCVVXaSbnHG+Y44pHqJXE5NJzxKrbcrtv5D6SupPQm3TXrYEOfWhucitoSJkCMk1k5zKmLSDl56/rfffNTOkQR7SJNdOslBJY8w2ANustV8yxmwEcADAy621X+u6h6GSacgZoz61ePhF+mKugTw5m2K1TdS4xaFRV5QLTXpZKXJrqAgZQnLNJKcypvx4Phq32nDRuNU+5otr3Iwx/xXAu621n+u6ZqobtzEMsDZ8Rb5A2s1NH1Ibn9j6Uwt6SRok/TG1FzMtm6Pa+rA2e30ZU/s0t8XHNi0bgRIHUWlpOxn/2j6naAJuY8x2ABcC2JejvpoYq3BYW0JaF1tjD/eIrV9D0krij6Q/pnb4gJYDQGrrw9rs9WVM7dPcFh/bFg+7uP62ewGgyAZGqj9rbDsZ/9ruQvLDSYwxTwDwSQBvtdY+1PL95caY/caY/UePHk1tjjrGKhx2adeaaxIf4OFia6nExb5jYKxjplYk/TG1wwe0HABSWx/WZq8vY2qf5rb42KblsAup/qyx7WT8a7sLSTduxpgNWN20fdha+0dt11hrr7LW7rbW7t6yZUtKc1QyVuGwS7vWXJP4AA8XW0slLvYdA2MdM7Ui6Y+pHT6g5QCQ2vqwNnt9GVP7NLfFxzYth11I9WeNbSfjX9tdSHk4iQHwQQD3W2vf6nLP2DRurocoSCZ81AQ1bv71jz1me6xQ4xZGyiSpGubUlA+riWVM7dPcljbtVpe9uXVeqQ+icnn/mF9z7OFlHDrykEqNW8kDVUqMbd93aM3zr4/sh5MYY14E4HoANwM4Mfv47dbaz3bdM6aNm2Si4DHG6BJCps3Y17Wxt4/UT9sYBaBi3OacP5Lva7kpmTQcKD9WavBRKNkTcFtrb7DWGmvt86y1F8x+OjdtY6MrrnZqGhZCCGlj7Ova2NtH6qdtjGoZtzntkHxfy00uG7WOFQ025CbLqZJTZB5Xu7xyojVRcPPzkLIIIaRWxr6ujb19pH66xqiGcZtz/ki+r+Uml41ax0oNPpImWx43F8YUKgnIxmfXGqNLCCFdjH1dG3v7SP1o0S252pa7Li190cfUNG4udo2B4gm4Xah545ZSZJ+TsU4AScbWR2NrTxdTaWdt8AAPQsLIMb5LHdrDuasbbf7RZo8ERRNwj50DdzyA1151Ix49vroJ/viBu/CRN9YnkByzyFOKsfXR2NrTxVTaWRsp/UKfkzGTY3ynqMOlTM5d3WjzjzZ7UpM8AfcU2Hv4Piwff+x/LmsVSE5R5OnL2PpobO3pYirtrI2UfqHPyZjJMb5T1OFSJueubrT5R5s9qeHGTYA9OzZjw3pz8t+1CiTHmKhQmrH10dja08VU2lkbKf1Cn5Mxk2N8p6jDpUzOXd1o8482e1JDjZsQrho37XG4EvZpb2MsY2tf7qSqpRib38ZCn19ifUafkzFDjRvJiUvC8lKEjhfN44yHkyhgCnG4U2jjmKC/iFY4NgkhRAdjXI+1tyl7Am5yKlOIw51CG8cE/UW0wrFJCCE6GON6XGubuHHLyBTicKfQxjFBfxGtcGwSQogOxrge19omhkpmphlPqzm+NpQxtmnM0F9EKxybhJAYpryGSLd9jH2puU3UuClEe3wtIYQQQkiNTPkda8ptHwvUuCmk1vhaQgghhBDNTPkda8ptHzvcuBWk1vhaQgghhBDNTPkda8ptHzsMlSyM5vhaQgghhJBamfI71pTbPgaocSOEkAnBhzYhhBBSJ10bt6USxhBCCEkHhemEEELI+KDGjRBCRgaF6YQQQsj44MaNEEJGBoXphBBCyPhgqCQhhIyMXedtwoffsIcaN1DrR9zgOKkX+o5MCW7cCCFkhOw6b9PkX2Ko9SMucJzUC31HpgZDJQkhhIwSav2ICxwn9ULfkanBjRshhJBRQq0fcYHjpF7oOzI1mMeNEELIaKH+hbjAcVIv9B0ZI0zATQghhBBCCCHK6dq4MVSSEEIIIYQQQpTDjRshhBBCCCGEKIcbN0IIIYQQQghRDjduhBBCCCGEEKIcbtwIIYQQQgghRDncuBFCCCGEEEKIcrhxI4QQQgghhBDlcOM2Ag7c8QDec903cOCOB0qbQgghhBBCCEnAUmkDSBwH7ngAP/fevXh05QROW1qHD79hD3adt6m0WYQQQgghhBBB+D9ulbP38H14dOUETlhgeeUE9h6+r7RJhBBCCCGEEGGSbdyMMe83xtxjjDmYqg4C7NmxGactrcN6A2xYWoc9OzaXNokQQgghhBAiTMpQyQ8AeDeADyWsY/LsOm8TPvyGPdh7+D7s2bGZYZKEEEIIIYSMkGQbN2vtF4wx21OVTx5j13mbuGEjhBBCCCFkxFDjRgghhBBCCCHKKb5xM8ZcbozZb4zZf/To0dLmEEIIIYQQQog6im/crLVXWWt3W2t3b9mypbQ5hBBCCCGEEKKO4hs3QgghhBBCCCH9pEwH8BEANwJ4ljHmLmPM61PVRQghhBBCCCFjJuWpkq9NVTYhhBBCCCGETAmGShJCCCGEEEKIcrhxI4QQQgghhBDlcONGCCGEEEIIIcox1trSNpzEGHMUwB2l7WjwZAD3ljaCZIP+ni70/XSh7wnAcTBl6Ptpo9H/51lrT8mTpmrjphFjzH5r7e7SdpA80N/Thb6fLvQ9ATgOpgx9P21q8j9DJQkhhBBCCCFEOdy4EUIIIYQQQohyuHEb5qrSBpCs0N/Thb6fLvQ9ATgOpgx9P22q8T81boQQQgghhBCiHP6PGyGEEEIIIYQoZ3QbN2PM04wx1xljbjHGHDLGvGX2+ZOMMZ8zxtw2+3PT7PPNs+u/Z4x5d6OsVxtjvjor57d76txljLnZGPMNY8y/M8aY2ecvNsZ8yRizYox5Vcp2TxFlvn7z7PObjDE3GGOek7LtRJ3/LzPGHJ35/yZjzBtStn3qKPP97y34/evGmO8mbDpZQNk4OM8Y8/lZGX9mjHlqyrZPnUK+/w1jzLeMMd9rfM53vYwE+P4njTEHZvP2gDHmJxbKap3PLXXqec+31o7qB8A2AM+f/X0jgK8DeA6A3wbwttnnbwPwW7O/Px7AiwC8GcC7F8rZDOBOAFtm//4ggL/TUecXAVwMwAC4BsCls8+3A3gegA8BeFXpvhnbjzJfn7Vwzd8HcG3p/hn7jzL/X7ZYJn+m4/vGNf8EwPtL989UfjSNAwAfB/CLs7//BID/VLp/xvxTyPd7ZvV+r/H5dvBdT7PvLwRw7uzv5wP49kJZg+t633UlfD+6/3Gz1h6x1n5p9vdjAG4B8CMA/nesTkjM/nz57JrvW2tvAPDDRlE7AHzdWnt09u8/BfDKZn3GmG1YfWm/0a568UMLZd9urf0qgBNiDSQnUebrhxYufTwAikcTo8n/JC+Kff9aAB8JbxnxQdk4eA6Az8/+ft3MBpKI3L6flbHXWnuk5XO+62UkwPdfttZ+Z/b5IQBnGGNOd13Xtb3nj27jtogxZjtWd9r7AJwzn3CzP58ycPs3ADzbGLPdGLOEVSc9reW6HwFw18K/75p9RjKiwdfGmH9kjPkmVn/r80/DWkJC0OB/AK+chdt8whjTdj9JgBLfwxhzHoC/AeC/+7eCxKJgHHwFj73wvwLARmPMZv+WEF8y+Z4oJMD3rwTwZWvtI3B/f1f1nj/ajZsx5gkAPgngrY3/DXHCWvsAgH8A4GMArgdwO4CVtqrabvetj4SjxdfW2vdYa58B4J8D+L997SBhKPH/fwOw3Vr7PKz+xvaDLdcSYZT4fs5rAHzCWnvc1w4Sh5Jx8CsAXmKM+TKAlwD4dkcZRJCMvifK8PW9MWYngN8C8Kb5Ry2Xtb2/q3rPH+XGzRizAavO/LC19o9mH//17L875//tec9QOdba/2atvchaezGAWwHcZoxZbx4Tor8DqzvvRRHyUwF8p608Io9SX38UDKHLghb/W2vvm/0GDwD+EMAuifaRbrT4foHXgGGS2dEyDqy137HW/oy19kIA/3L22YNCzSQtZPY9UYSv783qYUGfAvA6a+03Zx+3zmft7/mj27jNTnp5H4BbrLW/u/DVpwH84uzvvwjgvzqU9ZTZn5sA/EMA77XWHrfWXjD7+bXZf8ceM8bsmdX9OpeySTyafG2MeeZCcX8PwG2RzSMDKPP/toXi/j5WY+5JIjT5fnbvswBsAnCjQPOII5rGgTHmycaY+TvVvwDwfoEmkg5y+17WehKDr++NMU8E8P8B+BfW2v8xv7hrPqt/z7cKToiR/MHqqUEWwFcB3DT7+Wmsnhz0eay+UH8ewJMW7rkdwP0AvofVnfVzZp9/BMDXZj+v6alzN4CDAL4J4N14LLH5C2blfR/AfQAOle6fMf0o8/XvY1X0ehNWhek7S/fP2H+U+f83Z/7/ysz/zy7dP2P+0eT72XdXAHhX6X6Z2o+mcQDgVbP6vg7gvQBOL90/Y/4p5Pvfnt13YvbnFbPP+a6n2PdYla58f+HamwA8ZfZd57reqFPNe/68YkIIIYQQQgghShldqCQhhBBCCCGEjA1u3AghhBBCCCFEOdy4EUIIIYQQQohyuHEjhBBCCCGEEOVw40YIIYQQQgghyuHGjRBCyKQwxlxhjPmVnu9fbox5Tk6bCCGEkCG4cSOEEELW8nIA3LgRQghRBfO4EUIIGT3GmH8J4HUAvgXgKIADAB4EcDmA0wB8A8AvALgAwGdm3z0I4JWzIt4DYAuAHwB4o7X2f2Y0nxBCCOHGjRBCyLgxxuwC8AEAFwFYAvAlAFcC+I/W2vtm1/xrAH9trf1/jDEfAPAZa+0nZt99HsCbrbW3GWMuAvCb1tqfyN8SQgghU2aptAGEEEJIYv42gE9Za38AAMaYT88+P3+2YXsigCcA+OPmjcaYJwD4XwF83Bgz//j01AYTQgghTbhxI4QQMgXawks+AODl1tqvGGMuA/DjLdesA/Bda+0FySwjhBBCHODhJIQQQsbOFwC8whjzOGPMRgAvm32+EcARY8wGAD+3cP2x2Xew1j4E4K+MMf8HAJhVfiyf6YQQQsgq1LgRQggZPQuHk9wB4C4AXwPwfQC/OvvsZgAbrbWXGWP+FoA/BPAIgFcBOAHgDwBsA7ABwEette/I3ghCCCGThhs3QgghhBBCCFEOQyUJIYQQQgghRDncuBFCCCGEEEKIcrhxI4QQQgghhBDlcONGCCGEEEIIIcrhxo0QQgghhBBClMONGyGEEEIIIYQohxs3QgghhBBCCFEON26EEEIIIYQQopz/H5qXWY+Mr8OdAAAAAElFTkSuQmCC", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# マグニチュードを日時の関数としてプロットしてみる\n", "plt.figure(figsize=(15, 5))\n", "plt.plot(dates, magnitudes, '.')\n", "plt.xlabel('date')\n", "plt.ylabel('magnitude')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 地震の大きさと頻度" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "Text(0, 0.5, 'frequency')" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEGCAYAAACKB4k+AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAARM0lEQVR4nO3de6xlZX3G8e8joALKoDI1CowDHkSJVdRTmhZLmkIVCqMWaRWlDZYyxVs0jWlBjZbSC7XVqo2KIyC0okgQLSMgGpXgXWYochGNBFFGbBHBAawVwV//2GsWx2Eua2b2Omv2Pt9PcsJea99+i+TMc97Let9UFZIkATxs6AIkSdsPQ0GS1DIUJEktQ0GS1DIUJEmtHYcuYFvssccetXTp0qHLkKSJsnr16juqavGGnpvIUEiyDFg2MzPDqlWrhi5HkiZKku9t7LmJ7D6qqpVVtXzRokVDlyJJU2UiQ0GS1A9DQZLUMhQkSS1DQZLUMhQkSS1DQZLUMhQkSa2JvHlNmi9LT76k0+tuOf3IniuR5octBUlSy5aCJoJ/sUvzw5aCJKllS0FTxRaFtG0MBS1IXcNDWmjsPpIktQwFSVLLUJAktQwFSVLLUJAktQwFSVLLKanaYt4LIE0vQ0G96RIeBoe0fdmuuo+SvCjJB5L8Z5LnDV2PJC00vYdCkrOT3J7k+vXOH57k20luSnIyQFV9oqpOBI4HXtJ3bZKkXzUfLYVzgMPnnkiyA/Ae4AjgAODYJAfMecmbm+clSfOo91CoqiuBO9c7fRBwU1XdXFX3AecDL8zIPwGXVdXVfdcmSfpVQw007wncOud4DfCbwGuBw4BFSWaq6oz135hkObAcYMmSJfNQqjQeztrSJBgqFLKBc1VV7wbevak3VtUKYAXA7Oxs9VCbJC1YQ80+WgPsPed4L+C2gWqRJDWGailcBeyXZB/gB8BLgZcNVIsa7jEgaT6mpH4E+Aqwf5I1SU6oqvuB1wCXAzcCF1TVDVvwmcuSrFi7dm0/RUvSAtV7S6Gqjt3I+UuBS7fyM1cCK2dnZ0/cltokSb9qu7qjWZI0LENBktSayFBwTEGS+jGRoVBVK6tq+aJFi4YuRZKmykSGgiSpH4aCJKllKEiSWu68tgB4p7KkriYyFJIsA5bNzMwMXYo0dq6mqiFNZCh4R7NkeKgfExkKmh52bUnbFweaJUktQ0GS1DIUJEmtiQwF1z6SpH5MZCi49pEk9cPZR9IYOItK02IiWwqSpH4YCpKklqEgSWoZCpKk1kSGglNSJakfExkKTkmVpH5MZChIkvphKEiSWoaCJKnlHc3SlHMzHm0JWwqSpJahIElqGQqSpNZEjikkWQYsm5mZGbqU3nTpB7YPWNK4TWRLwZvXJKkfExkKkqR+GAqSpJahIElqGQqSpJahIElqGQqSpJahIElqGQqSpJahIElqTWQouEezJPVjIkPBZS4kqR8TuSCepGG4Yc/0m8iWgiSpH4aCJKll99EE69qUl6SubClIklqGgiSpZShIklqGgiSpZShIklqGgiSptdlQSLIqyauTPGY+CpIkDadLS+GlwBOBq5Kcn+T5SdJzXZKkAWw2FKrqpqp6E/AU4MPA2cD3k5ya5LF9FyhJmj+dxhSSPAN4O/DPwMeAY4C7gc/1V5okab5tdpmLJKuBnwBnASdX1c+bp76W5OAea9tUTcuAZTMzM0N8/TZxaQpJ27MuLYU/qqpDq+rDcwIBgKo6uqe6Nsn9FCSpH11C4c+T7L7uIMljkvxdfyVJkobSJRSOqKqfrDuoqruAP+itIknSYLqEwg5JHrHuIMnOwCM28XpJ0oTqsp/Ch4DPJvkgUMCfAef2WpUkaRCbDYWqeluS64BDgQCnVdXlvVcmSZp3nXZeq6rLgMt6rkWSNLAuax8dneQ7SdYmuTvJPUnuno/iJEnzq0tL4W3Asqq6se9iJEnD6jL76H8MBElaGLq0FFYl+SjwCaC9o7mqLuqrKEnSMLqEwm7A/wLPm3OuAENBkqZMlympr5iPQiRJw+uySupTgPcBj6+qpzfLaL+gqlz/SNIGDbUa8C2nHznI906TLgPNHwBOAX4BUFXXMtqNTZI0ZbqEwi5V9fX1zt3fRzGSpGF1CYU7kjyZ0eAySY4BfthrVZKkQXSZffRqYAXw1CQ/AL4LHNdrVZKkQXSZfXQzcFiSXYGHVdU9/ZclSRpCl9lHb1nvGICq+tueapI0APcPF3TrPvrpnMePBI4CXPZCkqZQl+6jt889TvIvwMW9VSRJGkyX2Ufr2wXYd9yFSJKG12VM4Tqa6ajADsBiYOzjCUn2Bd4ELKqqY8b9+ZKkzesypnDUnMf3M1pKu9PNa0nObt5/e1U9fc75w4F3MQqZM6vq9GaW0wlJLuxcvSRprLp0H90z5+dnwG5JHrvuZzPvPQc4fO6JJDsA7wGOAA4Ajk1ywJYWLkkavy4thauBvYG7gAC7A99vnis2Mb5QVVcmWbre6YOAm5qWAUnOB14IfLNLwUmWA8sBlixZ0uUtkqSOurQUPsVoO849qupxjLqDLqqqfapqawac9wRunXO8BtgzyeOSnAE8K8kpG3tzVa2oqtmqml28ePFWfL0kaWO6tBR+o6pOWndQVZclOW0bvjMbOFdV9WPgpA08J0maJ11C4Y4kbwY+xKi76Djgx9vwnWsYdUetsxdw2zZ8niRpTLp0Hx3LaBrqx5ufxc25rXUVsF+SfZI8nNHeDFt0M1ySZUlWrF27dhvKkCStb7OhUFV3VtXrgN+pqmdX1eur6s4uH57kI8BXgP2TrElyQjOd9TXA5YyWy7igqm7YkqKramVVLV+0aNGWvE2StBldbl77beBM4FHAkiTPBP6iql61ufdW1QZbFFV1KXDpFtYqSepZl+6jfwWeTzOOUFXfAA7psyhJ0jA6rX1UVbeud+qBHmqRJA2sSyjc2nQhVZKHJ3kDAy+d7UCzJPWjSyicxGhLzj0ZTSc9sDkejAPNktSPTQ40N+sUvbOqXj5P9UiSBrTJlkJVPQAsbu4nkCRNuS53NN8CfCnJxczZmrOq3tFXUZKkYWy0pZDkP5qHLwE+2bz20XN+BuNAsyT1Y1MtheckeRKjZbL/bZ7q6aSqVgIrZ2dnTxy6FkmaJpsKhTMYLZu9D7BqzvmwmX0UJEmTaaPdR1X17qp6GvDBqtp3zs/W7qMgSdrOdVkQ75XzUYgkaXidlrmQJC0MXaakbneSLAOWzczMDF2KpAm09ORLNvuaW04/ch4q2f5MZEvBZS4kqR8TGQqSpH4YCpKklqEgSWoZCpKklqEgSWpNZCi4IJ4k9WMiQ8EpqZLUj4kMBUlSPwwFSVLLUJAktQwFSVLLUJAktQwFSVLLUJAktdxPYUy6rM8uaXIM9Ts99D4OE9lS8OY1SerHRIaCJKkfhoIkqWUoSJJahoIkqWUoSJJahoIkqWUoSJJahoIkqWUoSJJaLnOxGS5fIWk+df03p6/lMCaypeAyF5LUj4kMBUlSPwwFSVLLUJAktQwFSVLLUJAktQwFSVLLUJAktQwFSVLLUJAktQwFSVLLUJAktQwFSVLLUJAktQwFSVLL/RQkTQ33P9l2E9lScD8FSerHRIaCJKkfhoIkqWUoSJJahoIkqWUoSJJahoIkqWUoSJJahoIkqWUoSJJahoIkqWUoSJJahoIkqWUoSJJahoIkqWUoSJJahoIkqWUoSJJahoIkqWUoSJJahoIkqWUoSJJahoIkqbXj0AWsk2RX4L3AfcAVVXXewCVJ0oLTa0shydlJbk9y/XrnD0/y7SQ3JTm5OX00cGFVnQi8oM+6JEkb1nf30TnA4XNPJNkBeA9wBHAAcGySA4C9gFublz3Qc12SpA3otfuoqq5MsnS90wcBN1XVzQBJzgdeCKxhFAzXsImwSrIcWA6wZMmSra5t6cmXbPV7JWlaDTHQvCcPtghgFAZ7AhcBL07yPmDlxt5cVSuqaraqZhcvXtxvpZK0wAwx0JwNnKuq+inwivkuRpL0oCFaCmuAvecc7wXcNkAdkqT1DBEKVwH7JdknycOBlwIXb8kHJFmWZMXatWt7KVCSFqq+p6R+BPgKsH+SNUlOqKr7gdcAlwM3AhdU1Q1b8rlVtbKqli9atGj8RUvSAtb37KNjN3L+UuDSPr9bkrTlXOZCktQyFCRJrYkMBQeaJakfqaqha9hqSX4EfK/DS/cA7ui5nCFM63XB9F6b1zV5pvHanlRVG7z7d6JDoaskq6pqdug6xm1arwum99q8rskzzde2IRPZfSRJ6oehIElqLZRQWDF0AT2Z1uuC6b02r2vyTPO1PcSCGFOQJHWzUFoKkqQODAVJUmuqQ2Fje0RPuiR7J/l8khuT3JDkdUPXNA5JHpnk60m+0VzXqUPXNE5JdkjyX0k+OXQt45TkliTXJbkmyaqh6xmXJLsnuTDJt5rftd8auqb5MNVjCkkOAe4F/r2qnj50PeOS5AnAE6rq6iSPBlYDL6qqbw5c2jZJEmDXqro3yU7AF4HXVdVXBy5tLJL8JTAL7FZVRw1dz7gkuQWYraqpusErybnAF6rqzGaZ/12q6icDl9W7qW4pVNWVwJ1D1zFuVfXDqrq6eXwPoyXI9xy2qm1XI/c2hzs1P1PxV0uSvYAjgTOHrkWbl2Q34BDgLICqum8hBAJMeSgsBEmWAs8CvjZwKWPRdLFcA9wOfKaqpuK6gHcCfwX8cuA6+lDAp5OsTrJ86GLGZF/gR8AHmy6/M5PsOnRR88FQmGBJHgV8DHh9Vd09dD3jUFUPVNWBjLZpPSjJxHf7JTkKuL2qVg9dS08OrqpnA0cAr266bSfdjsCzgfdV1bOAnwInD1vS/DAUJlTT5/4x4LyqumjoesataapfARw+bCVjcTDwgqbv/Xzg95J8aNiSxqeqbmv+ezvwceCgYSsaizXAmjkt1QsZhcTUMxQmUDMgexZwY1W9Y+h6xiXJ4iS7N493Bg4DvjVoUWNQVadU1V5VtZTRnuSfq6rjBi5rLJLs2kx2oOleeR4w8bP9quq/gVuT7N+cOhSY6IkcXfW6HefQmj2ifxfYI8ka4K1VddawVY3FwcCfANc1/e8Ab2y2OZ1kTwDOTbIDoz9YLqiqqZq+OYUeD3x89HcKOwIfrqpPDVvS2LwWOK+ZeXQz8IqB65kXUz0lVZK0Zew+kiS1DAVJUstQkCS1DAVJUstQkCS1DAWpR0lOSvKnzePjkzxxKz7jliR7jL866aGm+j4FaWhVdcacw+MZ3dh12zDVSJtnS0ELSpKlzfr4Zya5Psl5SQ5L8qUk30lyUPPz5WYhtC+vu6s1yS5JLkhybZKPJvlaktnmuXuT/H2zF8RXkzy+Of83Sd6Q5BhGy2af1+w7sPPcFkCS2SRXNI8fl+TTzfe/H8ic+o9r9py4Jsn7mxv9pLExFLQQzQDvAp4BPBV4GfBc4A3AGxktrXFIsxDaW4B/aN73KuCuqnoGcBrwnDmfuSvw1ap6JnAlcOLcL6yqC4FVwMur6sCq+tkm6nsr8MXm+y8GlgAkeRrwEkYL0B0IPAC8fGv+B0gbY/eRFqLvVtV1AEluAD5bVZXkOmApsIjRchv7MVoWeqfmfc9lFCZU1fVJrp3zmfcB65bkWA38/jbUdwhwdPM9lyS5qzl/KKMguqpZVmJnRkuMS2NjKGgh+vmcx7+cc/xLRr8TpwGfr6o/bParuKJ5PmzcL+rBNWMeoNvv1v082Fp/5HrPbWj9mQDnVtUpHT5b2ip2H0kPtQj4QfP4+Dnnvwj8MUCSA4Bf38LPvQd49JzjW3iwC+rFc85fSdMtlOQI4DHN+c8CxyT5tea5xyZ50hbWIG2SoSA91NuAf0zyJWDuQO57gcVNt9FfA9cCa7fgc88Bzlg30AycCrwryRcYtS7WORU4JMnVjJai/j5Aswf3mxntcnYt8BlGK8tKY+MqqVJHzUyfnarq/5I8mdFf7k+pqvsGLk0aG8cUpO52AT7f7HoX4JUGgqaNLQVJUssxBUlSy1CQJLUMBUlSy1CQJLUMBUlS6/8B4fSz9IrNQlEAAAAASUVORK5CYII=", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "plt.hist(magnitudes, bins=31)\n", "plt.yscale('log')\n", "plt.xlabel('magnitude')\n", "plt.ylabel('frequency')" ] }, { "cell_type": "markdown", "metadata": { "id": "uRS_KHxCjvLp" }, "source": [ "## 2005年からのファイルを全てダウンロード" ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "id": "Dp_sO5N-jx1R" }, "outputs": [], "source": [ "dates_all, latitudes_all, longtitudes_all, magnitudes_all = [], [], [], []\n", "for year in range(2005, 2020):\n", " # 文字列に対する.format()関数を行うことで、文字列中の{}の部分に値を代入することができる\n", " url = 'https://www.data.jma.go.jp/svd/eqev/data/bulletin/data/shindo/i{}.zip'.format(year)\n", " # ファイルをダウンロードし、ダウンロードしたZipファイルを展開する\n", " '''\n", " ここに続きのプログラムを入力する\n", " '''" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "colab": { "base_uri": "https://localhost:8080/", "height": 351 }, "id": "0m9F_4tHkXnJ", "outputId": "9ec317fb-f060-4b85-86b3-4fc930123043" }, "outputs": [], "source": [ "# マグニチュードを日時の関数としてプロットしてみる\n", "plt.figure(figsize=(15, 5))\n", "plt.plot(dates_all, magnitudes_all, '.', markersize=2)\n", "plt.xlabel('date')\n", "plt.ylabel('magnitude')" ] }, { "cell_type": "markdown", "metadata": { "id": "ZvgijwtJkiBt" }, "source": [ "## マグニチュードと頻度の関係をプロットする" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "colab": { "base_uri": "https://localhost:8080/", "height": 295 }, "id": "jbwf3lNekhFD", "outputId": "988e2179-23bc-49c4-eb7e-8a19edc84dd4" }, "outputs": [], "source": [ "plt.hist(magnitudes, bins=31)\n", "plt.yscale('log')\n", "plt.xlabel('magnitude')\n", "plt.ylabel('frequency')" ] }, { "cell_type": "markdown", "metadata": { "id": "RHqgIbEPk5_1" }, "source": [ "## 地震発生位置をプロットする\n", "\n", "地図をプロットできるパッケージ folium を使う\n", "\n", "参考 \n", "https://newtechnologylifestyle.net/python_gpsprot/\n", "\n" ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "colab": { "base_uri": "https://localhost:8080/" }, "id": "h85EX-tygosn", "outputId": "061b0764-b3c5-4e0c-cb8e-862d143665d3" }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Requirement already satisfied: folium in /home/keisukefujii/miniconda3/envs/project/lib/python3.9/site-packages (0.12.1)\r\n", "Requirement already satisfied: jinja2>=2.9 in /home/keisukefujii/miniconda3/envs/project/lib/python3.9/site-packages (from folium) (3.0.0)\r\n", "Requirement already satisfied: numpy in /home/keisukefujii/miniconda3/envs/project/lib/python3.9/site-packages (from folium) (1.20.2)\r\n", "Requirement already satisfied: branca>=0.3.0 in /home/keisukefujii/miniconda3/envs/project/lib/python3.9/site-packages (from folium) (0.4.2)\r\n", "Requirement already satisfied: requests in /home/keisukefujii/miniconda3/envs/project/lib/python3.9/site-packages (from folium) (2.25.1)\r\n", "Requirement already satisfied: MarkupSafe>=2.0.0rc2 in /home/keisukefujii/miniconda3/envs/project/lib/python3.9/site-packages (from jinja2>=2.9->folium) (2.0.1)\r\n", "Requirement already satisfied: urllib3<1.27,>=1.21.1 in /home/keisukefujii/miniconda3/envs/project/lib/python3.9/site-packages (from requests->folium) (1.26.4)\r\n", "Requirement already satisfied: idna<3,>=2.5 in /home/keisukefujii/miniconda3/envs/project/lib/python3.9/site-packages (from requests->folium) (2.10)\r\n", "Requirement already satisfied: chardet<5,>=3.0.2 in /home/keisukefujii/miniconda3/envs/project/lib/python3.9/site-packages (from requests->folium) (4.0.0)\r\n", "Requirement already satisfied: certifi>=2017.4.17 in /home/keisukefujii/miniconda3/envs/project/lib/python3.9/site-packages (from requests->folium) (2021.5.30)\r\n" ] } ], "source": [ "!pip install folium" ] }, { "cell_type": "code", "execution_count": 11, "metadata": { "id": "utOlXckbg8s7" }, "outputs": [], "source": [ "import folium" ] }, { "cell_type": "code", "execution_count": 12, "metadata": { "id": "4rMbLfWCgvHf" }, "outputs": [], "source": [ "map = folium.Map(location=[latitudes[0], longtitudes[0]], zoom_start=4)\n", "# 最初の100個の地震をプロットする\n", "for i in range(100):\n", " folium.Marker(location=[latitudes[i], longtitudes[i]]).add_to(map)" ] }, { "cell_type": "code", "execution_count": 13, "metadata": { "colab": { "base_uri": "https://localhost:8080/", "height": 1000 }, "id": "zgtjRuDihKBE", "outputId": "3e99801c-6293-4d27-e5c7-6e6246fd0f49" }, "outputs": [ { "data": { "text/html": [ "
Make this Notebook Trusted to load map: File -> Trust Notebook
" ], "text/plain": [ "" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "map" ] }, { "cell_type": "markdown", "metadata": { "id": "a5qFs3IIlPj5" }, "source": [ "## 課題（提出不要）\n", "\n", "+ 2019年以外のデータもダウンロードし、2011年の東北地方太平洋沖地震の大きさやその影響を調べよ\n", "+ 各地方ごとの地震の頻度を比べる \n", " 九州・四国・中国・近畿・中部・関東・東北・北海道など\n", "\n", "+ 地震の頻度は、規模 M に対して $\\rho(M) \\propto M^{-\\alpha}$ というふうにべき乗分布で表されることが知られている。ここでマグニチュード$m$は $m = \\log_{10} M + C$ ($C$は定数)として表される。べき指数 $α$ を上記のデータから求めてみよ。" ] }, { "cell_type": "markdown", "metadata": { "id": "0yxQpq01frUz" }, "source": [ "## その他の応用\n", "\n", "+ 重力波データの解析 \n", "https://www.gw-openscience.org/GW150914data/GW150914_tutorial.html\n", "\n", "+ 太陽フレアの観測データ \n", "https://www.ngdc.noaa.gov/stp/solar/solarflares.html \n", "https://www.swpc.noaa.gov/products/alerts-watches-and-warnings \n", "ftp://ftp.swpc.noaa.gov/pub/indices/events/ \n", "\n", "+ neuralnetwork playground \n", "https://playground.tensorflow.org/" ] } ], "metadata": { "colab": { "collapsed_sections": [], "name": "地震データ解析", "provenance": [] }, "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.7.7" } }, "nbformat": 4, "nbformat_minor": 1 }