diff --git a/notebooks/tutorial/HadISD/1_HadISD_Download.ipynb b/notebooks/tutorial/HadISD/1_HadISD_Download.ipynb index 31a7c932..ad495a08 100644 --- a/notebooks/tutorial/HadISD/1_HadISD_Download.ipynb +++ b/notebooks/tutorial/HadISD/1_HadISD_Download.ipynb @@ -30,7 +30,8 @@ "from tqdm.auto import tqdm\n", "import tarfile\n", "import gzip\n", - "import shutil" + "import shutil\n", + "from pathlib import Path" ] }, { @@ -45,13 +46,13 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 2, "id": "2eadaf27", "metadata": {}, "outputs": [], "source": [ - "%run Data_Config.ipynb\n", - "print(f\"Data will be downloaded to: {download_dir}\") " + "# ruff: noqa: F821\n", + "%run Data_Config.ipynb" ] }, { @@ -60,95 +61,93 @@ "metadata": {}, "source": [ "### Download HadISD Data\n", - "The following code will download the HadISD data files. Some files take longer to download than others depending on time of day. To download different WMO datasets, you can change `wmo_id_range` in the `Data_Config.ipynb` notebook .\n", + "The following code will download the HadISD data files. Some files take longer to download than others depending on time of day. To download different WMO datasets, you can change `wmo_id_ranges` in the `Data_Config.ipynb` notebook.\n", "\n", "The full list of available data can be found here:\n", - "https://www.metoffice.gov.uk/hadobs/hadisd/v340_2023f/download.html" + "https://www.metoffice.gov.uk/hadobs/hadisd/v340_2023f/download.html\n", + "\n", + "Station data has been split up into ranges to make downloads more managable. You may download as much or as little as you like. To get started we reccomend just downloading a few station ranges to get an idea of how to use HadISD data with PyEarthTools. " ] }, { "cell_type": "code", - "execution_count": null, - "id": "feb8d671", + "execution_count": 3, + "id": "8ddbebda", "metadata": {}, "outputs": [], "source": [ - "# Explain why stations are split into ranges, file size, and how it's not neccesssary to download all stations. " + "wmo_id_ranges = wmo_id_ranges # This has been defined in HadISD_data_config.ipynb" ] }, { "cell_type": "code", - "execution_count": null, + "execution_count": 4, "id": "11a188d4", "metadata": {}, - "outputs": [], + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Downloading HadISD data for WMO range: ['500000-549999', '722000-722999', '800000-849999']\n" + ] + } + ], "source": [ - "print(f\"Downloading HadISD data for WMO range: {wmo_id_range}\")" + "print(f\"Downloading HadISD data for WMO range: {wmo_id_ranges}\")" ] }, { "cell_type": "code", - "execution_count": null, - "id": "8ddbebda", - "metadata": {}, - "outputs": [], - "source": [ - "wmo_id_range = wmo_id_range # This has been defined in HadISD_data_config.ipynb\n", - "\n", - "wmo_str = f\"WMO_{wmo_id_range}\"\n", - "url = f\"https://www.metoffice.gov.uk/hadobs/hadisd/v340_2023f/data/{wmo_str}.tar.gz\"\n", - "tar_name = f\"{wmo_str}.tar\"\n", - "filename = download_dir / tar_name" - ] - }, - { - "cell_type": "code", - "execution_count": null, + "execution_count": 5, "id": "08ac36fd", "metadata": {}, "outputs": [], "source": [ - "# Get remote file size using HTTP HEAD\n", - "head = requests.head(url, allow_redirects=True)\n", - "remote_size = int(head.headers.get('content-length', 0))\n", - "\n", - "local_size = filename.stat().st_size if filename.exists() else 0\n", - "\n", - "if filename.exists() and local_size == remote_size:\n", - " print(f\"File already fully downloaded: {filename} ({local_size/1024**2:.2f} MB)\")\n", - "else:\n", - " headers = {}\n", - " mode = 'wb'\n", - " initial_pos = 0\n", - " if filename.exists() and local_size < remote_size:\n", - " headers['Range'] = f'bytes={local_size}-'\n", - " mode = 'ab'\n", - " initial_pos = local_size\n", - " print(f\"Resuming download for {filename.name} at {local_size/1024**2:.2f} MB...\")\n", - " else:\n", - " print(f\"Starting download for {filename.name}...\")\n", - "\n", - " response = requests.get(url, stream=True, headers=headers)\n", - " total = remote_size\n", - "\n", - " with open(filename, mode) as f, tqdm(\n", - " desc=f\"Downloading {filename.name}\",\n", - " total=total,\n", - " initial=initial_pos,\n", - " unit='B', unit_scale=True, unit_divisor=1024\n", - " ) as bar:\n", - " for chunk in response.iter_content(chunk_size=8192):\n", - " if chunk:\n", - " f.write(chunk)\n", - " bar.update(len(chunk))\n", - "\n", - " final_size = filename.stat().st_size\n", - " if final_size == remote_size:\n", - " print(f\"Download complete: {filename} ({final_size/1024**2:.2f} MB)\")\n", + "def download_wmo_range(wmo_id_range, download_dir):\n", + " wmo_str = f\"WMO_{wmo_id_range}\"\n", + " url = f\"https://www.metoffice.gov.uk/hadobs/hadisd/v340_2023f/data/{wmo_str}.tar.gz\"\n", + " tar_name = f\"{wmo_str}.tar\"\n", + " filename = Path(download_dir) / tar_name\n", + "\n", + " head = requests.head(url, allow_redirects=True)\n", + " remote_size = int(head.headers.get('content-length', 0))\n", + " local_size = filename.stat().st_size if filename.exists() else 0\n", + "\n", + " if filename.exists() and local_size == remote_size:\n", + " print(f\"File already fully downloaded: {filename} ({local_size/1024**2:.2f} MB)\")\n", " else:\n", - " print(f\"Warning: Download incomplete. Local size: {final_size}, Remote size: {remote_size}\")\n", + " headers = {}\n", + " mode = 'wb'\n", + " initial_pos = 0\n", + " if filename.exists() and local_size < remote_size:\n", + " headers['Range'] = f'bytes={local_size}-'\n", + " mode = 'ab'\n", + " initial_pos = local_size\n", + " print(f\"Resuming download for {filename.name} at {local_size/1024**2:.2f} MB...\")\n", + " else:\n", + " print(f\"Starting download for {filename.name}...\")\n", + "\n", + " response = requests.get(url, stream=True, headers=headers)\n", + " total = remote_size\n", + " with open(filename, mode) as f, tqdm(\n", + " desc=f\"Downloading {filename.name}\",\n", + " total=total,\n", + " initial=initial_pos,\n", + " unit='B', unit_scale=True, unit_divisor=1024\n", + " ) as bar:\n", + " for chunk in response.iter_content(chunk_size=8192):\n", + " if chunk:\n", + " f.write(chunk)\n", + " bar.update(len(chunk))\n", + "\n", + " final_size = filename.stat().st_size\n", + " if final_size == remote_size:\n", + " print(f\"Download complete: {filename} ({final_size/1024**2:.2f} MB)\")\n", + " else:\n", + " print(f\"Warning: Download incomplete. Local size: {final_size}, Remote size: {remote_size}\")\n", "\n", - "# Possibly also add check to see if netcdf files esist for the downloaded tar file, if so then don't download again" + " return filename, tar_name\n" ] }, { @@ -161,77 +160,239 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 6, "id": "fb79a81c", "metadata": {}, "outputs": [], "source": [ - "extract_dir = download_dir / tar_name.replace('.tar', '')\n", - "extract_dir.mkdir(exist_ok=True)\n", - "\n", - "extracted_files = list(extract_dir.glob('*'))\n", - "if extracted_files:\n", - " print(f\"Extraction directory '{extract_dir}' already contains {len(extracted_files)} files. Skipping extraction.\")\n", - "elif filename.exists():\n", - " with tarfile.open(filename, \"r:gz\") as tar:\n", - " tar.extractall(path=extract_dir)\n", + "def extract_wmo_tar(filename, tar_name, download_dir):\n", + " extract_dir = Path(download_dir) / tar_name.replace('.tar', '')\n", + " extract_dir.mkdir(exist_ok=True)\n", " extracted_files = list(extract_dir.glob('*'))\n", " if extracted_files:\n", - " print(f\"Extraction successful. {len(extracted_files)} files found in {extract_dir}.\")\n", - " # Delete the tar file after extraction\n", - " filename.unlink()\n", - " print(f\"Deleted tar file: {filename}\")\n", + " print(f\"Extraction directory '{extract_dir}' already contains {len(extracted_files)} files. Skipping extraction.\")\n", + " elif filename.exists():\n", + " with tarfile.open(filename, \"r:gz\") as tar:\n", + " tar.extractall(path=extract_dir)\n", + " extracted_files = list(extract_dir.glob('*'))\n", + " if extracted_files:\n", + " print(f\"Extraction successful. {len(extracted_files)} files found in {extract_dir}.\")\n", + " filename.unlink()\n", + " print(f\"Deleted tar file: {filename}\")\n", + " else:\n", + " print(f\"Warning: No files extracted to {extract_dir}. Tar file will not be deleted.\")\n", + " raise RuntimeError(\"Extraction failed, tar file not deleted.\")\n", " else:\n", - " print(f\"Warning: No files extracted to {extract_dir}. Tar file will not be deleted.\")\n", - " raise RuntimeError(\"Extraction failed, tar file not deleted.\")\n", - "else:\n", - " print(f\"No tar file found and extraction directory is empty. Nothing to extract.\")\n", - " raise FileNotFoundError(f\"Missing tar file: {filename}\")\n" + " print(\"No tar file found and extraction directory is empty. Nothing to extract.\")\n", + " raise FileNotFoundError(f\"Missing tar file: {filename}\")\n", + " return extract_dir" ] }, { "cell_type": "code", - "execution_count": null, - "id": "53161550", + "execution_count": 7, + "id": "4e43dcc4", "metadata": {}, "outputs": [], "source": [ - "# Create subfolder for netcdf\n", - "netcdf_dir = download_dir / \"netcdf\"\n", - "netcdf_dir.mkdir(parents=True, exist_ok=True)" + "# Move extracted .nc files into netcdf_dir after extraction\n", + "def move_netcdf_files(extract_dir, download_dir):\n", + " netcdf_dir = Path(download_dir) / \"netcdf\"\n", + " netcdf_dir.mkdir(parents=True, exist_ok=True)\n", + " num_files = 0\n", + " for gz_path in extract_dir.glob('*.nc.gz'):\n", + " nc_path = gz_path.with_suffix('') # Remove .gz extension\n", + " with gzip.open(gz_path, 'rb') as f_in, open(nc_path, 'wb') as f_out:\n", + " f_out.write(f_in.read())\n", + " gz_path.unlink()\n", + " shutil.move(str(nc_path), netcdf_dir / nc_path.name)\n", + " num_files += 1\n", + " print(f\"{num_files} .nc files have been extracted, cleaned up, and moved to the netcdf directory: {netcdf_dir}\")\n", + "\n", + " # Delete extraction directory\n", + " try:\n", + " shutil.rmtree(extract_dir)\n", + " print(f\"Deleted extraction directory: {extract_dir}\")\n", + " except Exception as e:\n", + " print(f\"Could not delete extraction directory {extract_dir}: {e}\")" + ] + }, + { + "cell_type": "markdown", + "id": "932e8906", + "metadata": {}, + "source": [ + "### Idempotent Checks" ] }, { "cell_type": "code", - "execution_count": null, - "id": "4e43dcc4", + "execution_count": 8, + "id": "dcb9b902", "metadata": {}, "outputs": [], "source": [ - "# Move extracted .nc files into netcdf_dir after extraction\n", - "num_files = 0\n", - "for gz_path in extract_dir.glob('*.nc.gz'):\n", - " nc_path = gz_path.with_suffix('') # Remove .gz extension\n", - " with gzip.open(gz_path, 'rb') as f_in, open(nc_path, 'wb') as f_out:\n", - " f_out.write(f_in.read())\n", - " gz_path.unlink() # Delete the .gz file after extraction\n", - " shutil.move(str(nc_path), netcdf_dir / nc_path.name)\n", - " num_files += 1\n", - "\n", - "print(f\"{num_files} .nc files have been extracted, cleaned up, and moved to the netcdf directory: {netcdf_dir}\")\n", - "\n", - "# Delete the extraction directory after processing\n", - "try:\n", - " shutil.rmtree(extract_dir)\n", - " print(f\"Deleted extraction directory: {extract_dir}\")\n", - "except Exception as e:\n", - " print(f\"Could not delete extraction directory {extract_dir}: {e}\")" + "def netcdf_files_exist_for_range(wmo_id_range, netcdf_dir):\n", + " \"\"\"Check if any .nc files for the given WMO range exist in the netcdf directory.\"\"\"\n", + " start, end = map(int, wmo_id_range.split('-'))\n", + " nc_files = list(Path(netcdf_dir).glob(\"*.nc\"))\n", + " for nc_file in nc_files:\n", + " try:\n", + " # Extract the first 6 digits from the station part of the filename\n", + " station_part = nc_file.stem.split('_')[-1]\n", + " wmo_number = int(station_part.split('-')[0])\n", + " if start <= wmo_number <= end:\n", + " return True\n", + " except Exception as e:\n", + " print(f\"Skipping file {nc_file.name}: {e}\")\n", + " continue\n", + " print(f\"No NetCDF files found for WMO range {wmo_id_range}.\")\n", + " return False" ] + }, + { + "cell_type": "code", + "execution_count": 9, + "id": "36ad9920", + "metadata": {}, + "outputs": [], + "source": [ + "def is_tar_fully_downloaded(wmo_id_range, download_dir):\n", + " \"\"\"Check if the tar file exists and is fully downloaded (size matches remote).\"\"\"\n", + " wmo_str = f\"WMO_{wmo_id_range}\"\n", + " tar_name = f\"{wmo_str}.tar\"\n", + " tar_path = Path(download_dir) / tar_name\n", + " url = f\"https://www.metoffice.gov.uk/hadobs/hadisd/v340_2023f/data/{wmo_str}.tar.gz\"\n", + "\n", + " if not tar_path.exists():\n", + " return False\n", + "\n", + " # Get remote file size\n", + " head = requests.head(url, allow_redirects=True)\n", + " remote_size = int(head.headers.get('content-length', 0))\n", + " local_size = tar_path.stat().st_size\n", + "\n", + " return local_size == remote_size" + ] + }, + { + "cell_type": "markdown", + "id": "77f044b9", + "metadata": {}, + "source": [ + "### Loop through each WMO range, download if necessary, extract, and move files" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "id": "ffcc5730", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "No NetCDF files found for WMO range 500000-549999.\n", + "Starting download for WMO_500000-549999.tar...\n" + ] + }, + { + "data": { + "application/vnd.jupyter.widget-view+json": { + "model_id": "274ec4d5674b4781a8657b2c98db1d66", + "version_major": 2, + "version_minor": 0 + }, + "text/plain": [ + "Downloading WMO_500000-549999.tar: 0%| | 0.00/411M [00:00\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "
<xarray.Dataset> Size: 327GB\n",
+       "Dimensions:                (station: 770, time: 473352, flagged: 19, test: 71,\n",
+       "                            reporting_v: 19, reporting_t: 1116, reporting_2: 2)\n",
+       "Coordinates:\n",
+       "    elevation              (station) float64 6kB 1.676e+03 401.0 ... 117.7 36.0\n",
+       "    latitude               (station) float64 6kB 38.47 42.33 ... 45.93 30.48\n",
+       "    longitude              (station) float64 6kB 102.2 120.7 ... 126.6 -87.19\n",
+       "  * station                (station) object 6kB '526740-99999' ... '722223-13...\n",
+       "  * time                   (time) datetime64[ns] 4MB 1970-01-01 ... 2023-12-3...\n",
+       "Dimensions without coordinates: flagged, test, reporting_v, reporting_t,\n",
+       "                                reporting_2\n",
+       "Data variables: (12/25)\n",
+       "    cloud_base             (station, time) float64 3GB dask.array<chunksize=(1, 59169), meta=np.ndarray>\n",
+       "    dewpoints              (station, time) float64 3GB dask.array<chunksize=(1, 59169), meta=np.ndarray>\n",
+       "    flagged_obs            (station, time, flagged) float64 55GB dask.array<chunksize=(1, 29585, 3), meta=np.ndarray>\n",
+       "    high_cloud_cover       (station, time) float64 3GB dask.array<chunksize=(1, 59169), meta=np.ndarray>\n",
+       "    low_cloud_cover        (station, time) float64 3GB dask.array<chunksize=(1, 59169), meta=np.ndarray>\n",
+       "    mid_cloud_cover        (station, time) float64 3GB dask.array<chunksize=(1, 59169), meta=np.ndarray>\n",
+       "    ...                     ...\n",
+       "    stnlp                  (station, time) float64 3GB dask.array<chunksize=(1, 59169), meta=np.ndarray>\n",
+       "    temperatures           (station, time) float64 3GB dask.array<chunksize=(1, 59169), meta=np.ndarray>\n",
+       "    total_cloud_cover      (station, time) float64 3GB dask.array<chunksize=(1, 59169), meta=np.ndarray>\n",
+       "    wind_gust              (station, time) float64 3GB dask.array<chunksize=(1, 59169), meta=np.ndarray>\n",
+       "    winddirs               (station, time) float64 3GB dask.array<chunksize=(1, 59169), meta=np.ndarray>\n",
+       "    windspeeds             (station, time) float64 3GB dask.array<chunksize=(1, 59169), meta=np.ndarray>\n",
+       "Attributes: (12/39)\n",
+       "    Conventions:                 CF-1.6\n",
+       "    Metadata_Conventions:        Unidata Dataset Discovery v1.0, CF Discrete ...\n",
+       "    acknowledgement:             RJHD was supported by the Joint BEIS/Defra M...\n",
+       "    cdm_data_type:               station\n",
+       "    creator_email:               robert.dunn@metoffice.gov.uk\n",
+       "    creator_name:                Robert Dunn\n",
+       "    ...                          ...\n",
+       "    station_id:                  526740-99999\n",
+       "    station_information:         Where station is a composite the station id ...\n",
+       "    summary:                     Quality-controlled, sub-daily, station datas...\n",
+       "    time_coverage_end:           2002-05-13T09:00Z\n",
+       "    time_coverage_start:         1964-01-01T00:00Z\n",
+       "    title:                       HadISD
" + ], + "text/plain": [ + " Size: 327GB\n", + "Dimensions: (station: 770, time: 473352, flagged: 19, test: 71,\n", + " reporting_v: 19, reporting_t: 1116, reporting_2: 2)\n", + "Coordinates:\n", + " elevation (station) float64 6kB 1.676e+03 401.0 ... 117.7 36.0\n", + " latitude (station) float64 6kB 38.47 42.33 ... 45.93 30.48\n", + " longitude (station) float64 6kB 102.2 120.7 ... 126.6 -87.19\n", + " * station (station) object 6kB '526740-99999' ... '722223-13...\n", + " * time (time) datetime64[ns] 4MB 1970-01-01 ... 2023-12-3...\n", + "Dimensions without coordinates: flagged, test, reporting_v, reporting_t,\n", + " reporting_2\n", + "Data variables: (12/25)\n", + " cloud_base (station, time) float64 3GB dask.array\n", + " dewpoints (station, time) float64 3GB dask.array\n", + " flagged_obs (station, time, flagged) float64 55GB dask.array\n", + " high_cloud_cover (station, time) float64 3GB dask.array\n", + " low_cloud_cover (station, time) float64 3GB dask.array\n", + " mid_cloud_cover (station, time) float64 3GB dask.array\n", + " ... ...\n", + " stnlp (station, time) float64 3GB dask.array\n", + " temperatures (station, time) float64 3GB dask.array\n", + " total_cloud_cover (station, time) float64 3GB dask.array\n", + " wind_gust (station, time) float64 3GB dask.array\n", + " winddirs (station, time) float64 3GB dask.array\n", + " windspeeds (station, time) float64 3GB dask.array\n", + "Attributes: (12/39)\n", + " Conventions: CF-1.6\n", + " Metadata_Conventions: Unidata Dataset Discovery v1.0, CF Discrete ...\n", + " acknowledgement: RJHD was supported by the Joint BEIS/Defra M...\n", + " cdm_data_type: station\n", + " creator_email: robert.dunn@metoffice.gov.uk\n", + " creator_name: Robert Dunn\n", + " ... ...\n", + " station_id: 526740-99999\n", + " station_information: Where station is a composite the station id ...\n", + " summary: Quality-controlled, sub-daily, station datas...\n", + " time_coverage_end: 2002-05-13T09:00Z\n", + " time_coverage_start: 1964-01-01T00:00Z\n", + " title: HadISD" + ] + }, + "execution_count": 7, + "metadata": {}, + "output_type": "execute_result" + } + ], "source": [ "ds_combined = load_combined_dataset(zarr_output_dir)\n", "ds_combined" ] - }, - { - "cell_type": "markdown", - "id": "d448a05f", - "metadata": {}, - "source": [ - "# Data Organization: NetCDF and Zarr stores\n", - "\n", - "To keep your workflow clear and reproducible, we recommend storing both the raw NetCDF files and the processed Zarr data in separate subfolders inside your main WMO directory. For example:\n", - "\n", - "- `HadISD_data/WMO_080000-099999/netcdf/` (raw NetCDF files)\n", - "- `HadISD_data/WMO_080000-099999/zarr/` (processed Zarr stores with harmonized time coordinates)\n", - "\n", - "This makes it obvious which data is raw and which is ready for fast, parallel analysis." - ] } ], "metadata": { "kernelspec": { - "display_name": "Python 3 (ipykernel)", + "display_name": "pyearthtools", "language": "python", "name": "python3" }, @@ -312,7 +3812,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.13.2" + "version": "3.13.5" } }, "nbformat": 4, diff --git a/notebooks/tutorial/HadISD/3_HadISD_XGBoost_Pipeline.ipynb b/notebooks/tutorial/HadISD/3_HadISD_XGBoost_Pipeline.ipynb index 8d2cf388..fc7dc079 100644 --- a/notebooks/tutorial/HadISD/3_HadISD_XGBoost_Pipeline.ipynb +++ b/notebooks/tutorial/HadISD/3_HadISD_XGBoost_Pipeline.ipynb @@ -19,8 +19,9 @@ "outputs": [], "source": [ "import numpy as np\n", - "import pandas as pd\n", - "from pathlib import Path\n", + "import matplotlib.pyplot as plt\n", + "from xgboost import XGBClassifier\n", + "from sklearn.metrics import classification_report, confusion_matrix\n", "\n", "import pyearthtools.pipeline as petpipe\n", "import pyearthtools.data as petdata\n", @@ -36,10 +37,11 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 2, "metadata": {}, "outputs": [], "source": [ + "# ruff: noqa: F821\n", "%run Pipeline_Config.ipynb" ] }, @@ -53,9 +55,17 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 3, "metadata": {}, - "outputs": [], + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Total number of stations: 770\n" + ] + } + ], "source": [ "hadisd = HadISDIndex()\n", "all_stations = hadisd.get_all_station_ids()\n", @@ -65,13 +75,21 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 4, "metadata": {}, - "outputs": [], + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "List of first ten stations: ['501360-99999', '502460-99999', '503530-99999', '504340-99999', '504420-99999', '504680-99999', '505270-99999', '505480-99999', '505570-99999', '505640-99999']\n" + ] + } + ], "source": [ "# Select first n stations\n", "first_ten_stations = all_stations_ordered[:10]\n", - "print(f\"List of first ten stations:\", first_ten_stations)" + "print(\"List of first ten stations:\", first_ten_stations)" ] }, { @@ -90,16 +108,721 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 5, "metadata": {}, - "outputs": [], + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
Pipeline\n",
+       "\tDescription                    `pyearthtools.pipeline` Data Pipeline\n",
+       "\n",
+       "\n",
+       "\tInitialisation                 \n",
+       "\t\t exceptions_to_ignore           None\n",
+       "\t\t iterator                       None\n",
+       "\t\t sampler                        None\n",
+       "\tSteps                          \n",
+       "\t\t HadisdDataClass.HadISDIndex    {'HadISDIndex': {'station': "['501360-99999', '502460-99999', '503530-99999', '504340-99999', '504420-99999', '504680-99999', '505270-99999', '505480-99999', '505570-99999', '505640-99999']", 'variables': 'None'}}\n",
+       "\t\t values.AddFlaggedObs           {'AddFlaggedObs': {'flagged_labels': "['temperatures', 'dewpoints', 'slp', 'stnlp' ... 'precip12_depth', 'precip15_depth', 'precip18_depth', 'precip24_depth']"}}\n",
+       "\t\t values.SetMissingToNaN         {'SetMissingToNaN': {'varname_val_map': {'total_cloud_cover': '-999.0', 'low_cloud_cover': '-999.0', 'mid_cloud_cover': '-999.0', 'high_cloud_cover': '-999.0', 'winddirs': '-999.0'}}}\n",
+       "\t\t variables.Drop                 {'Drop': {'__args': '()', 'variables': "'flagged_obs'"}}\n",
+       "\t\t coordinates.Drop               {'Drop': {'__args': '()', 'coordinates': "['latitude', 'longitude', 'elevation']", 'ignore_missing': 'False'}}\n",
+       "\t\t __main__.TrainTestSplit        {'TrainTestSplit': {}}\n",
+       "\t\t __main__.FeatureTargetSplit    {'FeatureTargetSplit': {}}\n",
+       "\t\t conversion.ToNumpy             {'ToNumpy': {'reference_dataset': 'None', 'run_parallel': 'False', 'saved_records': 'None', 'warn': 'True'}}\n",
+       "\t\t __main__.MedianImputePerStation {'MedianImputePerStation': {}}\n",
+       "\t\t __main__.TransposeAndFlattenX  {'TransposeAndFlattenX': {}}\n",
+       "\t\t conversion.ToNumpy[1]          {'ToNumpy': {'reference_dataset': 'None', 'run_parallel': 'False', 'saved_records': 'None', 'warn': 'True'}}\n",
+       "\t\t __main__.SelectAndFlattenY     {'SelectAndFlattenY': {}}\n",
+       "\t\t __main__.FilterNaNSamples      {'FilterNaNSamples': {}}\n",
+       "\t\t __main__.FeatureTargetSplit[1] {'FeatureTargetSplit': {}}\n",
+       "\t\t conversion.ToNumpy[2]          {'ToNumpy': {'reference_dataset': 'None', 'run_parallel': 'False', 'saved_records': 'None', 'warn': 'True'}}\n",
+       "\t\t __main__.MedianImputePerStation[1] {'MedianImputePerStation': {}}\n",
+       "\t\t __main__.TransposeAndFlattenX[1] {'TransposeAndFlattenX': {}}\n",
+       "\t\t conversion.ToNumpy[3]          {'ToNumpy': {'reference_dataset': 'None', 'run_parallel': 'False', 'saved_records': 'None', 'warn': 'True'}}\n",
+       "\t\t __main__.SelectAndFlattenY[1]  {'SelectAndFlattenY': {}}\n",
+       "\t\t __main__.FilterNaNSamples[1]   {'FilterNaNSamples': {}}
" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/html": [ + "

Graph

" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "image/svg+xml": [ + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "HadISDIndex_8fcd7115-0d04-46ad-930b-ca03e14913a3\n", + "\n", + "HadisdDataClass.HadISDIndex\n", + "\n", + "\n", + "\n", + "AddFlaggedObs_9b13e396-f451-47e2-bf2a-cd1fde1498a6\n", + "\n", + "values.AddFlaggedObs\n", + "\n", + "\n", + "\n", + "HadISDIndex_8fcd7115-0d04-46ad-930b-ca03e14913a3->AddFlaggedObs_9b13e396-f451-47e2-bf2a-cd1fde1498a6\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "SetMissingToNaN_24f2948c-1801-416a-9800-63a82673caf7\n", + "\n", + "values.SetMissingToNaN\n", + "\n", + "\n", + "\n", + "AddFlaggedObs_9b13e396-f451-47e2-bf2a-cd1fde1498a6->SetMissingToNaN_24f2948c-1801-416a-9800-63a82673caf7\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "Drop_3fcf5737-43ba-4619-b446-dd0ca4cd8f4a\n", + "\n", + "variables.Drop\n", + "\n", + "\n", + "\n", + "SetMissingToNaN_24f2948c-1801-416a-9800-63a82673caf7->Drop_3fcf5737-43ba-4619-b446-dd0ca4cd8f4a\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "Drop_3bb14fd5-7ee9-494f-966c-352afdf69348\n", + "\n", + "coordinates.Drop\n", + "\n", + "\n", + "\n", + "Drop_3fcf5737-43ba-4619-b446-dd0ca4cd8f4a->Drop_3bb14fd5-7ee9-494f-966c-352afdf69348\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "TrainTestSplit_efda0281-ace4-45f9-97b6-c5a90ab7cef1\n", + "\n", + "__main__.TrainTestSplit\n", + "\n", + "\n", + "\n", + "Drop_3bb14fd5-7ee9-494f-966c-352afdf69348->TrainTestSplit_efda0281-ace4-45f9-97b6-c5a90ab7cef1\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "FeatureTargetSplit_18774330-4f23-463d-801a-a66de3b568a1\n", + "\n", + "__main__.FeatureTargetSplit\n", + "\n", + "\n", + "\n", + "TrainTestSplit_efda0281-ace4-45f9-97b6-c5a90ab7cef1->FeatureTargetSplit_18774330-4f23-463d-801a-a66de3b568a1\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "FeatureTargetSplit_891d48df-c84a-4c0f-b14f-1a5f1707d584\n", + "\n", + "__main__.FeatureTargetSplit\n", + "\n", + "\n", + "\n", + "TrainTestSplit_efda0281-ace4-45f9-97b6-c5a90ab7cef1->FeatureTargetSplit_891d48df-c84a-4c0f-b14f-1a5f1707d584\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "ToNumpy_607494ac-d92e-463b-936a-ce83b20b6187\n", + "\n", + "conversion.ToNumpy\n", + "\n", + "\n", + "\n", + "FeatureTargetSplit_18774330-4f23-463d-801a-a66de3b568a1->ToNumpy_607494ac-d92e-463b-936a-ce83b20b6187\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "ToNumpy_ca59c6d8-f1c6-4b20-9ef2-acc334da1ad5\n", + "\n", + "conversion.ToNumpy\n", + "\n", + "\n", + "\n", + "FeatureTargetSplit_18774330-4f23-463d-801a-a66de3b568a1->ToNumpy_ca59c6d8-f1c6-4b20-9ef2-acc334da1ad5\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "MedianImputePerStation_bd3a0a4a-f8eb-4f00-9a2c-314c013ac8e0\n", + "\n", + "__main__.MedianImputePerStation\n", + "\n", + "\n", + "\n", + "ToNumpy_607494ac-d92e-463b-936a-ce83b20b6187->MedianImputePerStation_bd3a0a4a-f8eb-4f00-9a2c-314c013ac8e0\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "TransposeAndFlattenX_fd90c437-479d-49fe-a223-54429c389000\n", + "\n", + "__main__.TransposeAndFlattenX\n", + "\n", + "\n", + "\n", + "MedianImputePerStation_bd3a0a4a-f8eb-4f00-9a2c-314c013ac8e0->TransposeAndFlattenX_fd90c437-479d-49fe-a223-54429c389000\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "FilterNaNSamples_82e45d93-8672-4efd-bd48-7e2aaf5fb7b4\n", + "\n", + "__main__.FilterNaNSamples\n", + "\n", + "\n", + "\n", + "TransposeAndFlattenX_fd90c437-479d-49fe-a223-54429c389000->FilterNaNSamples_82e45d93-8672-4efd-bd48-7e2aaf5fb7b4\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "SelectAndFlattenY_2ad4a151-252a-47cf-8d3d-e09727602866\n", + "\n", + "__main__.SelectAndFlattenY\n", + "\n", + "\n", + "\n", + "ToNumpy_ca59c6d8-f1c6-4b20-9ef2-acc334da1ad5->SelectAndFlattenY_2ad4a151-252a-47cf-8d3d-e09727602866\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "SelectAndFlattenY_2ad4a151-252a-47cf-8d3d-e09727602866->FilterNaNSamples_82e45d93-8672-4efd-bd48-7e2aaf5fb7b4\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "ToNumpy_036f6dbc-cd85-4f8c-bf7d-23e4a3f4d971\n", + "\n", + "conversion.ToNumpy\n", + "\n", + "\n", + "\n", + "FeatureTargetSplit_891d48df-c84a-4c0f-b14f-1a5f1707d584->ToNumpy_036f6dbc-cd85-4f8c-bf7d-23e4a3f4d971\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "ToNumpy_e4b5574b-2dff-4fe9-8edc-9dc90d46fa50\n", + "\n", + "conversion.ToNumpy\n", + "\n", + "\n", + "\n", + "FeatureTargetSplit_891d48df-c84a-4c0f-b14f-1a5f1707d584->ToNumpy_e4b5574b-2dff-4fe9-8edc-9dc90d46fa50\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "MedianImputePerStation_876cb65b-1c97-4b0a-855d-607bf960730c\n", + "\n", + "__main__.MedianImputePerStation\n", + "\n", + "\n", + "\n", + "ToNumpy_036f6dbc-cd85-4f8c-bf7d-23e4a3f4d971->MedianImputePerStation_876cb65b-1c97-4b0a-855d-607bf960730c\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "TransposeAndFlattenX_6bcf84c0-b50d-4f66-b725-a8aa32ce6fd0\n", + "\n", + "__main__.TransposeAndFlattenX\n", + "\n", + "\n", + "\n", + "MedianImputePerStation_876cb65b-1c97-4b0a-855d-607bf960730c->TransposeAndFlattenX_6bcf84c0-b50d-4f66-b725-a8aa32ce6fd0\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "FilterNaNSamples_c58d7fad-bd15-437d-a663-b7209d458754\n", + "\n", + "__main__.FilterNaNSamples\n", + "\n", + "\n", + "\n", + "TransposeAndFlattenX_6bcf84c0-b50d-4f66-b725-a8aa32ce6fd0->FilterNaNSamples_c58d7fad-bd15-437d-a663-b7209d458754\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "SelectAndFlattenY_69d7fd56-8e65-4603-9b90-314f0b76a9f5\n", + "\n", + "__main__.SelectAndFlattenY\n", + "\n", + "\n", + "\n", + "ToNumpy_e4b5574b-2dff-4fe9-8edc-9dc90d46fa50->SelectAndFlattenY_69d7fd56-8e65-4603-9b90-314f0b76a9f5\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "SelectAndFlattenY_69d7fd56-8e65-4603-9b90-314f0b76a9f5->FilterNaNSamples_c58d7fad-bd15-437d-a663-b7209d458754\n", + "\n", + "\n", + "\n", + "\n", + "\n" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], "source": [ "data_prep_pipe = petpipe.Pipeline(\n", " petdata.archive.hadisd(station = first_ten_stations), # Use all stations\n", - " SqueezeStationCoordinates(),\n", " petdata.transforms.values.AddFlaggedObs(flagged_labels),\n", " petdata.transforms.values.SetMissingToNaN(varname_val_map),\n", " petdata.transforms.variables.Drop(\"flagged_obs\"),\n", + " petdata.transforms.coordinates.Drop(['latitude', 'longitude', 'elevation']),\n", " TrainTestSplit(test_size=0.2, random_state=42, dim='time'), # returns (ds_train, ds_test)\n", " (\n", " # branch 1: train\n", @@ -114,7 +837,7 @@ " ),\n", " # branch y_train\n", " (\n", - " petpipe.operations.xarray.conversion.ToNumpy(),\n", + " petpipe.operations.xarray.conversion.ToNumpy(),\n", " SelectAndFlattenY(test_number=33),\n", " ),\n", " 'map'\n", @@ -148,9 +871,18 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 6, "metadata": {}, - "outputs": [], + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/Users/joelmiller/Projects/Python/PyEarthTools_Fork/PyEarthTools/packages/data/src/pyearthtools/data/indexes/_indexes.py:487: IndexWarning: Could not find time in dataset to select on. Petdt('1969-01-01T00')\n", + " warnings.warn(\n" + ] + } + ], "source": [ "# Select data based on date\n", "(X_train, y_train),(X_test, y_test) = data_prep_pipe[\"1969-01-01T00\"] # Curretnly we have to select a date outside of the date range to select the entire dataset (fix coming soon)" @@ -158,9 +890,20 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 7, "metadata": {}, - "outputs": [], + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "X_train shape: (920121, 22)\n", + "y_train shape: (920121,)\n", + "X_test shape: (232093, 22)\n", + "y_test shape: (232093,)\n" + ] + } + ], "source": [ "# Print the shapes of the resulting arrays\n", "print(f\"X_train shape: {X_train.shape}\")\n", @@ -178,12 +921,21 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 8, "metadata": {}, - "outputs": [], + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/Users/joelmiller/miniconda3/envs/pyearthtools/lib/python3.13/site-packages/xgboost/training.py:183: UserWarning: [17:41:03] WARNING: /Users/runner/work/xgboost/xgboost/src/learner.cc:738: \n", + "Parameters: { \"use_label_encoder\" } are not used.\n", + "\n", + " bst.update(dtrain, iteration=i, fobj=obj)\n" + ] + } + ], "source": [ - "from xgboost import XGBClassifier\n", - "\n", "# Calculate scale_pos_weight for class imbalance\n", "scale_pos_weight = (len(y_train) - np.sum(y_train)) / np.sum(y_train)\n", "#scale_pos_weight = num_zeros / num_ones \n", @@ -213,25 +965,71 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 9, "metadata": {}, - "outputs": [], + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + " precision recall f1-score support\n", + "\n", + " 0.0 0.99 0.97 0.98 223592\n", + " 1.0 0.53 0.76 0.63 8501\n", + "\n", + " accuracy 0.97 232093\n", + " macro avg 0.76 0.87 0.80 232093\n", + "weighted avg 0.97 0.97 0.97 232093\n", + "\n", + "[[217897 5695]\n", + " [ 2013 6488]]\n" + ] + } + ], "source": [ "# compare predictions with true labels\n", - "from sklearn.metrics import classification_report, confusion_matrix\n", "print(classification_report(y_test, y_pred))\n", - "print(confusion_matrix(y_test, y_pred)) \n", - "\n", - "# Plot confusion matrix\n", - "import matplotlib.pyplot as plt\n", - "import seaborn as sns\n", + "print(confusion_matrix(y_test, y_pred)) " + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAArAAAAJOCAYAAABP3DTfAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjMsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvZiW1igAAAAlwSFlzAAAPYQAAD2EBqD+naQAAZllJREFUeJzt3Qd4FOXWwPGTBAi9k4TepXekqDRpKiIIKE0NXbiAFKUp0q94QWmCcBURVFBABaWIIFWlSS8KAlKlqhTpJfs95/WbvbtJSDIhy2aS/89nDTvz7sxsP3vmzJkAl8vlEgAAAMAhAv29AQAAAIAdBLAAAABwFAJYAAAAOAoBLAAAAByFABYAAACOQgALAAAARyGABQAAgKMQwAIAAMBRCGABAADgKASwydyBAwekQYMGkilTJgkICJCFCxcm6PKPHDliljtz5swEXa6T1a5d21wSyuXLl6VTp04SFhZmHuvevXtLUqKvHb1f+lpCVO3atZMCBQp4TdPHa9iwYX7bpqRAHz99HBMCn4NAwiOATQQOHTokL774ohQqVEhSp04tGTNmlIcfflgmTpwo165d8+m6w8PDZffu3fLvf/9bPv74Y6lcubIkpS92/dLQxzO6x1GDd52vl7feesv28k+ePGm+5Hbs2CH+9MYbb5gvxm7dupnn8Pnnn/fp+jRYevLJJ6Odt2bNGvN4fv755z7dBv0BULp0aa9pN2/eNO+ZChUqmOc8c+bMUqpUKenSpYvs27cvSkBsXfQ9lytXLmnYsKFMmjRJ/v77b1vbcuzYMenatat5XIKDgyUkJESefvppWb9+vSQWui36Wr1w4YJPlr9gwQJ5/PHHJXv27JIqVSrzeD777LOyatUq8aU5c+bIhAkTJDF8zkR3WbZsmV+3DUjKUvh7A5K7JUuWyDPPPGO++F544QXzpaxfxD/88IP069dP9u7dK++9955P1q1B3YYNG+S1116THj16+GQd+fPnN+tJmTKl+EOKFCnk6tWrsmjRIvOF6mn27NkmeLl+/Xq8lq0B7PDhw03gUr58+Tjfbvny5ZKQNEioVq2aDB06VJKz5s2byzfffCOtW7eWzp07y61bt0zgunjxYnnooYekePHiXuNHjBghBQsWNONOnz5tgm/NXo8bN06+/vprKVu2bKzr/PHHH+WJJ54w/9YseMmSJc2yNEh+5JFHZMqUKeaHxf2m7zl97XsGsPpa1WBLA/uE4nK5pEOHDub+6g+Hvn37mj0Bp06dMkFt3bp1zWOkj7+vAtg9e/b4fa+Dfn5Pnz49yvRy5cr5ZXuA5IAA1o8OHz4srVq1MkGeBiE5c+Z0z+vevbscPHjQBLi+cu7cOfM3Ib/QIrMyXP78YtFs9qeffholgNUvv0aNGskXX3xxX7ZFA+m0adOaDFVCOnv2rAmcEsrt27clIiIiwbfTl3766ScTqOqehFdffdVr3uTJk6PNPGrG0HOPw6BBg8z7ULPLTz31lPzyyy+SJk2au67z/Pnz0qJFCzNGg7TChQu752kgpxndnj17msBOf2DcT/frPff222+b4NUK/D13uesPY90j4BlI+5P+UNXXdGBgwu941Pv43HPPJfhyAdwdJQR+NGbMGFO/+MEHH3gFr5YiRYpIr169vAKLkSNHmi9KDcw086df1jdu3Ih2F69mcatUqWK+zLQ84aOPPnKP0d2JGjgrzfTqF49VRxddTd3dasJWrFhhMk0aBKdPn16KFSvmFUDcrfZLA4UaNWpIunTpzG2bNGliAobo1qeBvJU50lrd9u3bm2Awrtq0aWMyc55BjAY8WkKg8yL766+/5JVXXpEyZcqY+6S7ozXY2blzp3uMZusefPBB82/dHmuXoXU/rV3cW7dulZo1a5rA1XpcItfAahmHPkeR778GQFmyZDGZ3ph21+sPIf2hY22DVSuqgW3Hjh0lNDTULF+zQbNmzfJahvX8aAmF7oq1Xls///yzJJSjR4/Kv/71L/Pa0GAvW7ZsZq9DdDWtusfh0UcfNePy5Mkjo0aNMsF0XMpwlP5YiSwoKMisMy503a+//rrZ5k8++STGsf/9739NtnXs2LFewavS7bcea830xlZXGV2d71dffWV+YOnueH1OdB36/r9z506s98OzBlb/6ntcacbZ83VSq1atu2YJ9fnS12BMWd7Ro0ebzLa+fqK7X1rOop9Blt9++80891mzZjXvCQ3sI/9It17X8+bNMz9I9HWgr1/N5upngUXfQ3pbfa6s+2R9blnL+Oyzz2Tw4MGSO3dus75Lly6Z+fPnz5dKlSqZ50nLHjT4/P333+V+2rVrl/lcs0rHNHOt2ew///wzyli9P/pjS8fp60Bfe/H5PAaSksTx0ziZ0t3a+uEV191ruotSvxQ16/Pyyy/Lpk2bzBeIBj66u86TftDrOA1gNECaMWOG+bDUD22tC2zWrJn5kOvTp4/Z5aq7QfUDzw4NNjRQ1l2t+iWtX7K6Xs1GxeS7774zAaHed/0Q1i/Cd955xwQf27ZtixI8a+ZUv3j1vup83VWndYb/+c9/4rSdel+1RvHLL780XxBW9lW/eCtWrBhlvH7J6sFs+kWr6z1z5oz5wtAvew3sNKAoUaKEuc9DhgwxNZYajCvP51K/iPR+apZdvyA1kIyO1m1qQK/Pk5Z0aMCl69NSA81g6fqio9ug8/U51C95fU2oHDlymMdUv+D1+dDyEL0f+qWtrwEN5D1/GKkPP/zQZKj0vujzqAFGTHS3+x9//BFl+sWLF6NM0x8LugtbHwfdTg2cpk6darZPH08NLJQGg3Xq1DE/1AYOHGh+3Gj5TExZUIv1Y0zLQvR1dC9ZPw269EtfH38tRYjp/asBReTMvkUfcw0m9PWuj63drKgGtfqe1Gyu/tXXiL7eNAjToDmu9PX/66+/mr0Q48ePNwGb9TrR+6r3UXfDe9YU63Omt9Hg7270B7L+2NPsq75mY6PvI31/6I/Pl156yfyo0M8zzXZrzbTWDXt68803TbZUf0zq60p/8Ldt29Z87lkZXp1+4sQJc79U5M8wDfg166rL0B/6+m99XPVHp/4A1c8U3S59D+rn1vbt2+O9Ryry+0HLpvQH991osKmfNbotGrxa5WL6d+PGje7gVLfpscceM0kOLQPRHzD62aPPX0J8HgOO5YJfXLx40aUPf5MmTeI0fseOHWZ8p06dvKa/8sorZvqqVavc0/Lnz2+mrVu3zj3t7NmzruDgYNfLL7/snnb48GEzbuzYsV7LDA8PN8uIbOjQoWa8Zfz48eb6uXPn7rrd1jo+/PBD97Ty5cu7QkJCXH/++ad72s6dO12BgYGuF154Icr6OnTo4LXMp59+2pUtW7a7rtPzfqRLl878u0WLFq66deuaf9+5c8cVFhbmGj58eLSPwfXr182YyPdDH78RI0a4p/30009R7pulVq1aZt60adOinacXT99++60ZP2rUKNdvv/3mSp8+vatp06auuNDnqlGjRl7TJkyYYJb3ySefuKfdvHnTVb16dbPsS5cuue+XjsuYMaN5jcR1fXqbmC7z5893j7969WqUZWzYsMGM++ijj9zTevfubaZt2rTJPU23KVOmTGa6bqtFH79SpUq5r0dERLgf89DQUFfr1q1dU6ZMcR09ejTKuvX50nH6/N2NrrNChQoxPg6ZM2d2lStXLsYxL730klnXrl27on0PRd4mz/sY3eP24osvutKmTWteozG9X3VZui6Lvr4jL19duHDBlTp1ateAAQOibLe+dy5fvnzX+zZx4kSzzAULFrjiwnp+v//+e/e0v//+21WwYEFXgQIF3O+51atXm3ElSpRw3bhxI8r6du/e7Z6mr/voPqusZRQqVMjrcdT3gH72lC5d2nXt2jX39MWLF5vxQ4YMcU+723MVmT7+0b0HPN/j0X0ORvf8fvrpp1E+uxs3bmye899//9097cCBA64UKVLY/jwGkhJKCPzE2pWVIUOGOI1funSp+avZGE9W1i3ybjitibSygkp/revuJP3Fn1CsTIXu6ozLbl6lB3foUfuaCfTM8mnWoH79+u776Umzp570fml203oM40JLBXQ3nGb5NJOlf6MrH1CaubDq5DTboeuydsdpBjiudDmaXYkLbWWmnSg0c6IZM83WaRY2vvRx1KyOZtc9M0Ka+dKylbVr10Y5ACpyRicmVatWNRmkyJfoujl4ZlA1c6uPp5bH6OvH8/HUbdZdyp67nHWbNOsWG81Wffvtt6bkQMsuNNuodeSamW3ZsqXto+/1+Y6tG4HOj+39a82329kg8uOmt9cMn772NYPp2VXhXmiGUMt39PH6J+795zU/d+5cadq0qcmCJ+RnmD63mpX2fJw1669Z+chlK/re8azDtj7P7HyG6V4Nz8dxy5YtprRGS1o8M+JaqqF7ZOJ7zIEuK/J7QeuDY+K5XZqh1+fXqpW23hf6XGgGX58Lzz0x+v7RvTv3+nkMOBkBrJ9oXaWdLzat89KgSj+4PGmQoh9cOt9Tvnz5oixDv9j1wJOEooGB7q7V0gbdPa67iLVuLaYPT2s7NRiMbpe4fohfuXIlxvui90PZuS9aIqFftPrFrLuZdfdh5MfSotuvuySLFi1qglDd5aqBlNasRbeL/G607s7OgVAa/GlQrwG+tnPSMon40sdZtz/yASv6GFvzI+/utkMfk3r16kW5aIlKZFrOoLu+8+bN6/V4alDp+Xha2xxZdK+V6OiydbeyltRo3bAGZRoQ6GvSbpcNDfLjEpzGJchV8XkudZew7lbXIFM/L/Qxsw4UsvM6jI12P9FWYN9//725rgGT7laPrR1bfD7D7va+t+Yn9Ps+8us6ps8fDWAjb0NcaQlFXN4LnrT8Qkt59LNTg1l9fq3ttZ5fDbb1/RPdZ1XkafH5PAacjADWT/TDX39Ra+2ZHXFtrH23mjQryxKfdUQ+eEQ/dNetW2e+8PTLTgM8/RDVTGpcDjSJq3u5L57BjWY2teZO64Xvln21+qpqplsPvtIDeTSzpxkVrR2282UQl9pNT1rrpl9YSnvz3k92t9UOPRJfD8bRWlH9QtXaUn08tQbSV1+uWi+oX+D6+tSgWNertbVxoTWVGkDc7QeO516O/fv3RzmI0pO+J/RHjP6YsfPe0uBea671wEHNymu9rT5mVt13Qj5ueqCWBjzWQWv6V38YaxAWE6stma9eqwnxvvfl6/pe6fvh/fffd9fn6/vC6hsbn+f3fn0eA4kFAawfacG9Hj2tB+7ERneF6oeaHjnvSTMl+mVnHcSSEDTTEd0u1+iyE5rh06ODtYWO7gLUQEV30a9evfqu90PpF39kultUs3Mx7ba8Fxq0apCoGSMNbu5GDyjRg4m0O4SO0937+mUe+TFJqLP0KM066y5TDYp0l6oesKIH0sSXPs76Won8RWjtek7I10ts9PHUXbm6S1UPLNQvVN2NHPnxtLY5suheK3GlZRNannK3g86iowfGqZiOwFeNGzc2u3714Ljo6G5xzWrq+9wKpKwsYuT7Hvm9peUuWmqhBxxplk6Xoa9B6/Z2xfRa1UBR3xv6PGl2Uw9g1NKT2A7M0ufQKteIS4Ckz+/d3vfWfLvsvgdj+vzRaffrfaGP88qVK83Binpglmba9X2hB7Z60sy9lid4dl+wRDfN7ucx4GQEsH7Uv39/E6zpLh8NRCPT4FaPjlVWs/TIZ53RDyqrhiuhaJsWzUDpL3iL1Zg88i6wyKyG/nfLSmlmTMdoJtTzS1wz0ZqBsO6nL2hQqkcla19QzTDdjX5xR87yaJASuc2OFWgnxNmNBgwYYHbj6uOiz6l2YtCgL6bsXkz0cdQ6Xy2ZsGgGUrs9aN2hZvful+geT92OyEGPbrMefb1582avXsVa8hEbDXz18YtMnxv9gaiBVlxqfPXLXl8juis3ttpbrVnW15G2qIpcl6mBrdVeTd/nFqvdlmbKPH+8RG5vZgWPno+bnuDk3XfflfiI7bWqGTsNqvQ+aflEXHqaavcIfd1qyYb+jS4zqtlc6/nU51f/7fmDXe+7Hnmvr/f49DLW+2WnnEJbUWlQOG3aNK/3lrbZ0/uRkJ+jMYnu+Y3u890qTdAfFZ7t9DR41W2+189jwMloo+VH+mWm7Zx0N4/WgXmeiUvbDlltj5T2atSARj/srd2L+mWgX3xa4K/BWULRrKN+IWlWQA/60YNGtO3RAw884HXQje7a1C9i/dDXzIXu/tYvWG2V5HmgRmTaAkgPQKhevbpp82W10dJaP1+ev12zEzG1BbJotkvvmwYg2vZHd5FqEBU5O6LPn9Yf65eh1kPql6ke3GS3nlSDJn3c9ExaVlsvbWulbaa0J6lmY+3SLK4eBKavH+1FqwGCZti0pY5+Scb1wJuEoI+nZjX1+dUgRQMY3c0ZuTerBno6TlsGadbRaqOlry3PH1PR0V3tmkXU15Ue7KO1xPqDQ98f+sWv9zlyRlEDAM3+aWCvPyD1edDd9Lo+PRNXbG2vNCjWx1QDM33eIp+JS4Na/bGkrwmLZvO1tlNf9xr46jZpizsNrj0DcH3d6fL1Pa/vQQ2E9bGxs/vck1WPqTXC+v7WzLRmkK3AVk+2oJ89+pmjn0XRtZeLjnW2QM2ua5ZPM+wa1OtjoEGXfkZZp9TVbKNma/U50vukz5E+P9rHWE8mEp8TDOj90h9pWvKjde3640zv193o/dYyDH1v62eoZpqtNlr6HtGWdPerhExLlPS9rXsHtMREf8DrYxGZfibqPK1v1bO66Q8/fV3p8+V5Guv4fh4DjuXvNghwuX799VdX586dTSuZVKlSuTJkyOB6+OGHXe+8845Xu5xbt26Z1k/adiZlypSuvHnzugYNGuQ15m5tlaJr33S3Nlpq+fLlptWMbk+xYsVMO6bIbWVWrlxp2oDlypXLjNO/2r5I70/kdURuNfXdd9+Z+5gmTRrTwklbxfz8889eY6z1RW4LE13LodjaaN3N3dpoabuxnDlzmu3T7dS2T9G1v/rqq69cJUuWdLe0se5n5DZPnjyXo+2s9PmqWLGieX499enTx7QW03XH5G7P95kzZ1zt27d3Zc+e3Tw/ZcqUifI8xPQasLs+z/ZFnm20zp8/794ObeHVsGFD1759+8xy9DnypO2m9LHR1k65c+d2jRw50vXBBx/E2kZL7+ubb75ppuvzps9HlixZXI8++qjr888/j/b1Y130sdG2avXr1zetmqwWY3F15MgRV5cuXVz58uVzvw70oq/x6GzdutVVtWpVs169zbhx46J9Tf/444+uatWqmdegvrf69+/vbremj7OdNlpKH0t9TPU1Fd37Z8yYMWb6G2+84bJLH+MGDRq4smbNah4DfQ5atmzpWrNmjde4Q4cOmZZ22oJMn+MqVaqYFlaxvYbu9lmibb7atGljlqfzrMfhbsuwzJ0717RJ09Z4us1t27Z1nThxwmuMnTZacf2c8dx2XZ+2BNRt17ZtzzzzjOvkyZPRPnf6Wavbq6+ZwoULu6ZPn24+o/QxtPN5DCQlAfo/fwfRAJBUaG2jZmU166VZXqecklezkJqB1Nrd6LqYIHHRPW+a/Y6ubhxIDqiBBYAEpAfR6K5x3aWuu6qdkCPQbdSDFnW3OsFr4qNlVp40aNW+up6npAaSGzKwAJBM6UFUWu+rwba2dNIm+HpqVyQuevCr1rNrHb52rNBjEvTALO2qEl3vZCA5IIAFgGRKywX0oEM9GFHPTqVtl5D4aCZff2TowXHa01oPgNV+1XE92A5IighgAQAA4CjUwAIAAMBRCGABAADgKI4+kYGeJlOblGtT9oQ8rScAAPAvrXDUU3/nypUrXie68CU9256edMiXtAVfbCdUSc4cHcBq8Jo3b15/bwYAAPCR48ePmzOKJabgNU2GbCK3r/p0PXpWOz07G0FsEgxgrdNhpioZLgFBzmgWDiD+fltp/7S6AJzp778vSfHC+e/rqa/jwmReb1+V4JLhIr6KPe7clNM/zzLrIoBNggGsVTagwSsBLJD06TnkASQvibZEMEVqn8UeroDEVTKRGPEIAQAAwFEcnYEFAADwC00M+yo7nEiTzokJGVgAAAA4ChlYAAAAu7RO1Ve1qtTAxopHCAAAAI5CBhYAAMAurX/1WQ0sRbCxIQMLAAAARyEDCwAAYBc1sH7FIwQAAABHIQMLAABgFzWwfkUGFgAAAI5CBhYAAMA2H9bAkl+MFY8QAAAAHIUMLAAAgF3UwPoVGVgAAAA4ChlYAAAAu+gD61c8QgAAAHAUMrAAAAB2UQPrV2RgAQAA4ChkYAEAAOyiBtaveIQAAADgKGRgAQAA7KIG1q/IwAIAAMBRyMACAADYRQ2sX/EIAQAAwFHIwAIAAMSrBtZXGVhqYGNDBhYAAACOQgYWAADArsCAfy6+WjZiRAYWAAAAjkIGFgAAwC66EPgVjxAAAAAchQwsAACAXZyJy6/IwAIAAMBRyMACAADYRQ2sX/EIAQAAwFEIYAEAAOJbA+uriw2jR4+WBx98UDJkyCAhISHStGlT2b9/v9eY69evS/fu3SVbtmySPn16ad68uZw5c8ZrzLFjx6RRo0aSNm1as5x+/frJ7du3vcasWbNGKlasKMHBwVKkSBGZOXNmlO2ZMmWKFChQQFKnTi1Vq1aVzZs3296W2BDAAgAAONjatWtNQLhx40ZZsWKF3Lp1Sxo0aCBXrlxxj+nTp48sWrRI5s+fb8afPHlSmjVr5p5/584dE7zevHlT1q9fL7NmzTLB6ZAhQ9xjDh8+bMbUqVNHduzYIb1795ZOnTrJt99+6x4zd+5c6du3rwwdOlS2bdsm5cqVk4YNG8rZs2fjvC1xEeByuVziUJcuXZJMmTJJcJnOEhCUyt+bA8DHzm2c5O9NAHAfv+Nzh2SRixcvSsaMGSXRxR51RkhAitQ+WYfr9nW5sXpIvO/7uXPnTAZVg8OaNWua5eTIkUPmzJkjLVq0MGP27dsnJUqUkA0bNki1atXkm2++kSeffNIEk6GhoWbMtGnTZMCAAWZ5qVKlMv9esmSJ7Nmzx72uVq1ayYULF2TZsmXmumZcNRs8efJkcz0iIkLy5s0rPXv2lIEDB8ZpW+KCDCwAAEAipMGy5+XGjRtxut3FixfN36xZs5q/W7duNVnZevXquccUL15c8uXLZ4JGpX/LlCnjDl6VZk51vXv37nWP8VyGNcZahmZvdV2eYwIDA811a0xctiUuCGABAAASYQ2sZi4122tdtNY1NhEREWbX/sMPPyylS5c2006fPm0yqJkzZ/Yaq8GqzrPGeAav1nxrXkxjNMi9du2a/PHHH6YUIboxnsuIbVvigjZaAAAAidDx48e9Sgj0wKnYdO/e3ezi/+GHHyQpI4AFAABIhH1gNXi1UwPbo0cPWbx4saxbt07y5Mnjnh4WFmZ272utqmfmU4/813nWmMjdAqzOAJ5jIncL0Ou6jWnSpJGgoCBziW6M5zJi25a4oIQAAADAwVwulwleFyxYIKtWrZKCBQt6za9UqZKkTJlSVq5c6Z6mbba0bVb16tXNdf27e/dur24B2tFAg9OSJUu6x3guwxpjLUNLA3RdnmO0pEGvW2Pisi1xQQYWAADArnj0a7W1bBu6d+9ujur/6quvTC9Yq5ZU62Y1M6p/O3bsaNpb6YFdGpRqVwANGK2j/rXtlgaqzz//vIwZM8YsY/DgwWbZVulC165dTXeB/v37S4cOHUywPG/ePNOZwKLrCA8Pl8qVK0uVKlVkwoQJpp1X+/bt3dsU27bEBQEsAACAg02dOtX8rV27ttf0Dz/8UNq1a2f+PX78eNMRQE8aoN0MtHvAu+++6x6ru/61/KBbt24mmEyXLp0JREeMGOEeo5ldDVa1j+vEiRNNmcL06dPNsiwtW7Y0bbe0f6wGweXLlzcttjwP7IptW+KCPrAAHIM+sEDykej7wNZ7UwJS+qgP7K3rcuO7gYnuvicm1MACAADAUSghAAAAcHANbHJEBhYAAACOQgYWAAAgXhlYX/WBJQMbGzKwAAAAcBQysAAAAInwTFy4Ox4hAAAAOAoZWAAAALvoQuBXBLAAAAB2UULgVzxCAAAAcBQysAAAAHZRQuBXZGABAADgKGRgAQAA7KIG1q94hAAAAOAoZGABAADsogbWr8jAAgAAwFHIwAIAANgUEBBgLj5auG+Wm4SQgQUAAICjkIEFAACwiQysf5GBBQAAgKOQgQUAALBLk6S+SpSSgI0VGVgAAAA4ChlYAAAAm6iB9S8ysAAAAHAUMrAAAAA2kYH1LzKwAAAAcBQysAAAADaRgfUvMrAAAABwFDKwAAAANpGB9S8ysAAAAHAUMrAAAAB2cSYuvyIDCwAAAEchAwsAAGATNbD+RQYWAAAAjkIGFgAAIB5JUt9lYH2z2KSEDCwAAAAchQwsAACATQH6n89qVUnBxoYMLAAAAByFDCwAAIBNdCHwLzKwAAAAcBQysAAAAHZxJi6/IgMLAAAARyEDCwAAYJcPa2Bd1MDGigwsAAAAHIUMLAAAQCLqQuC7/rJJBxlYAAAAOAoZWAAAAJvIwPoXGVgAAAAHW7dunTRu3Fhy5cplgt+FCxdGG2xHvowdO9Y9pkCBAlHmv/nmm17L2bVrl9SoUUNSp04tefPmlTFjxkTZlvnz50vx4sXNmDJlysjSpUu95rtcLhkyZIjkzJlT0qRJI/Xq1ZMDBw7Yvs8EsAAAAPHtA+uriw1XrlyRcuXKyZQpU6Kdf+rUKa/LjBkzTIDavHlzr3EjRozwGtezZ0/3vEuXLkmDBg0kf/78snXrVhP8Dhs2TN577z33mPXr10vr1q2lY8eOsn37dmnatKm57Nmzxz1Gg95JkybJtGnTZNOmTZIuXTpp2LChXL9+3dZ9poQAAADAwR5//HFzuZuwsDCv61999ZXUqVNHChUq5DU9Q4YMUcZaZs+eLTdv3jTBb6pUqaRUqVKyY8cOGTdunHTp0sWMmThxojz22GPSr18/c33kyJGyYsUKmTx5sglYNfs6YcIEGTx4sDRp0sSM+eijjyQ0NNRkjVu1ahXn+0wGFgAAwKa77ZZPqIuvnDlzRpYsWWKypJFpyUC2bNmkQoUKJsN6+/Zt97wNGzZIzZo1TfBq0czp/v375fz58+4xWhLgScfodHX48GE5ffq015hMmTJJ1apV3WPiigwsAABAIqS77T0FBweby72YNWuWybQ2a9bMa/pLL70kFStWlKxZs5pSgEGDBpkyAs2wKg08CxYs6HUbzZxa87JkyWL+WtM8x+h0a5zn7aIbE1cEsAAAAImwC4EeKOVp6NChpu70XmgJQNu2bc1BVp769u3r/nfZsmVNpvXFF1+U0aNH33PQ7AsEsAAAAInQ8ePHJWPGjO7r9xpIfv/992aX/9y5c2Mdq7v1tYTgyJEjUqxYMVMbq+UHnqzrVt3s3cZ4zremaRcCzzHly5e3dV+ogQUAAEiENbAavHpe7jWA/eCDD6RSpUqmY0Fs9ACtwMBACQkJMderV69u2nXdunXLPUYP0NLgVssHrDErV670Wo6O0elKSxA0iPUco2US2o3AGhNXZGABAAAc7PLly3Lw4EH3dT1YSgNQrWfNly+fO1DUHq1vv/12lNvrAVQaRGpnAq2P1et9+vSR5557zh2ctmnTRoYPH24O/howYIBpjaVdB8aPH+9eTq9evaRWrVpmHY0aNZLPPvtMtmzZ4m61pYF57969ZdSoUVK0aFET0L7++uumf62227KDABYAAMDBZ+LasmWLCT4j17OGh4fLzJkzzb81mNQ2VtqnNTLN7Op8ra+9ceOGCSw1gPWsi9VuAcuXL5fu3bubLG727NnNCQmsFlrqoYcekjlz5pg2Wa+++qoJUrU9VunSpd1j+vfvb/rW6u0uXLggjzzyiCxbtixKTW6sj5FL741D6a8JfUCDy3SWgKD/tXUAkDSd2zjJ35sA4D5+x+cOySIXL170qgNNLLFHaLuPJTBVWp+sI+LmVTkz8/lEd98TEzKwAAAAdsXjjFm2lo0YcRAXAAAAHIUMLAAAgINrYJMjMrAAAABwFDKwAAAANpGB9S8ysAAAAHAUMrAAAAA2kYH1LzKwAAAAcBQysAAAAHbRB9avyMACAADAUcjAAgAA2EQNrH+RgQUAAICjkIEFAACwiQysf5GBRYJ7pUMD+eGTfnL2h7fk6MrRMm9cZymaP8RrTIdmD8u37/eSM9+PlWvbJ0um9Gm85teoVNRMj+5SqWQ+97h61UvI2lkvm3UdWzVaPn2rk+TLmdVrWS8+W1O2fzFY/towTnYueF3aPFnFa75uR3Tr+XJSV588PkBy88bI4ZIhdZDXpWLZkl5jNm3cII0a1pPQrBkkV47M0rBubbl27Zp7/o7t2+SpJxpIntCski9XDun5rxfl8uXLXsuIvA69fD7vs/t2PwHcP2RgkeBqVCwi0+auk617j0qKFEEyvEdjWTy1h1RoNkquXr9pxqRNnVJWrP/ZXEa+1CTKMjbu/E0K1BvkNW3Iv56UOlWKydafj5nr+XNlk/nju8ikT1ZJu9dmSab0qWXMK83ls7c7y0Nt/mPGdH7mERnRs7F0H/mpbNl7VB4sXUCmvN5aLly6KkvX7TFjWr38vqRKGeReT9ZM6WTz3EHy5YrtPn2cgOSkRMlSsmjpcvf1oBQpvILXZk89IX37DZS3xk808/bs2imBgf/kWE6dPGmC12YtnpW3Jrwjf1+6JANe6StdO7eXTz6d77Weqe99IPUbPOa+nilz5vty/5D8BIgPM7C0IXBGADtlyhQZO3asnD59WsqVKyfvvPOOVKninSWDczTp8a7X9S5DP5Hjq96UCiXzyo/bDplpk+escWdao3Pr9h058+ff7uspUgTKk7XLytTP1rqnVSyZV4ICA2XYlMXicrnMtAkfrTRBrY6/fTtC2jSqIh988aN8vnybmX/k9z+lUql88nK7+u4A9vylq17rfqZhJRNoE8ACCSdFihQSGhYW7byB/V+Wrv/qKS/3G+Ce9sADxdz//mbpYkmRMqWMmzjZHdROnPyuVKtcXg4dOiiFCxfxCljvth4gIVFCkMxLCObOnSt9+/aVoUOHyrZt20wA27BhQzl79qy/Nw0JJGP61Obv+YvegaIdT9YqK9kypZOPv9ronrbt5+MS4YqQF5pUk8DAALMeDVhXbdpvgleVKmUKuX7zlteyrl2/JZVL5zdBbnTCmz4k87/d5s4WA7h3hw4ekKIF80iZ4kWkY/hzcvzYP3tSzp09K1s2b5IcISFSt/YjUihfTnmsXh1Z/+MP7tvevHlTUqVM5Q5eVeo0/5QdbfAYp17u3VPy5w6R2o9Uk49mznD/uAWQtPg9gB03bpx07txZ2rdvLyVLlpRp06ZJ2rRpZcaMGf7eNCTQr8ixr7SQ9dsPyc+HTsV7OeFNq8uKDb/I72cvuKcdPfmnPPmvKaZE4eKmCXLm+7ckd2hmea7//1473234Rdo1fUgqlMhrrlcsmU/aPf2QCWyzZ04fZT2VS+WX0kVzycwF6+O9rQC8Va5SRaa9P0MWfL1Uxr8zRY4cPSIN69aSv//+Ww4f/s2MeWPUcGnXvqMZU75CBWn8eH05ePCAmVerdh05c+a0TBj3lglmz58/L0MH/1NipHvuLIOHDJdZn3wmXy/5Vp5q2kz69uoh096d7Kd7jWRzIgNfXZB4Swj0g2jr1q0yaND/ah31F3a9evVkw4YNUcbfuHHDXCyXLl26b9uK+Jkw6FkpVSSn1G0/Pt7LyB2SWepXLyHPDfD+UROaLYO8+3obmb1ok8xbtlXSpwuWId2elDlvdZRGXf/50hr9/jIJzZZR1s56RXSPzNm//jbjX25fXyIiXNEGyrt//d3UywJIGA0aPu7+d+kyZaXyg1Wl1AMF5cvP50mx4iXM9A4du8jz4e3Nv8uVryBrVq+Sj2d+KMNHvWHqZ/87/UMZNOAVGfb6qxIUFCRdu/eUkNBQr6zsgFcHu/+ty7h65YpMHPeWdOve877eXwBJPAP7xx9/yJ07dyQ0NNRrul73/FVtGT16tGTKlMl9yZv3n6waEqfxA56RJ2qUloadJ3llTu16vkk1+fPiFVm8dpfX9Bdb1pRLl6/JaxO/kp37T5j62g6vzZJHqxaXKmUKmDHXb9ySrsNnS9aH+kjxRkOl6OOvy9FTf5rbnTvvfQRz2tSpTP3rrIVRfzwBSDiZM2eWIkUfkN8OHZKwsJxmWvES/wSylmLFi8uJ4/+UGahnW7WRQ0dPyq+/HZejJ8/Jq4OHyh/nzknBggVjzPz+/vsJr8QHkNA1sL66IJGXENihmdqLFy+6L8ePH/f3JiGG4PWpR8vJYy9OMrv678ULT1WTOYs3u+taPQPOyFnUOxH/jNGaWE96Ww2idbwGqd98vzdKbVyz+hUkOFUK+XTpT/e0vQBipu2vDv92SMJy5pT8BQpIzly55MCvv3qNOXjggOTNlz/KbTXrmj59evli/lxJnTq11Klb/67r2b1zp2TJkkWCg4N9cj8AJNMSguzZs5tdQWfOnPGartfDojmKVD+E+CByRtlAy8cryzN93pPLV66bXf3q4uXrJiOqdJru2i+cL7u5rnWnf1+5LsdPn/fqClC7ygNSME92+TCamlQNQnu2rSODujxmSggypA2W4T2eMgHzjn0nzJgi+ULMAVs/7TkiWTKklZeef1RKFs4lnV7/OMry2jWtLovW7JK/Ll7x2WMDJEevDuwnTzzxpAlIT506KW+MHCaBQUHS4tlWJtPUq88rZlqZsmWlTLnyMufjj+TX/fvk4znz3Mv479QpUrVadUmXPr2sXvmdDB7UX4aPGm2yuWrpkkVy9swZqVK1mgSnTi2rV66Qt8aMlpd6v+zHe46kjC4EyTiATZUqlVSqVElWrlwpTZs2NdMiIiLM9R49evhz03AP9MQBasX03l7TOw/5WD5ZtMn8u1OLGjK46xPued/N6BNljNIDsDbsOCS/HvH+kaPW/vSrtHt1lvQJryd9w+ubrgGbdh2Wp7q/6w6Ug4ICpNfzj8oD+UNNa651W36VOu3elmOn/vJalp5o4eGKRdy1swASzsnfT0j78Lby159/SvYcOaT6Qw/LqrXrJUeOHGZ+95695Pr16zKw38ty/vxfUrpsOflqybdSqHBh9zK2/rRZ/j1ymFy5fFkeKFZcJk6eKq3bPu+enzJlSnn/v1NlUP+Xzd6VQoWLyOj/vCXtOnb2y30G4FsBLj/3GNE2WuHh4fLf//7X9H6dMGGCzJs3T/bt2xelNjYyPYhLa2GDy3SWgKBU922bAfjHuY2T/L0JAO4T/Y7PHZLFlAxmzJhREgsr9ijY43MJDE7rk3VE3Lgqhye3SHT3PTHx+4kMWrZsKefOnZMhQ4aYA7fKly8vy5YtizV4BQAAQPLk9wBWabkAJQMAAMAptEzVdzWwPllskuKoLgQAAABAosjAAgAAOIrJwPpu2YgZGVgAAAA4ChlYAAAAm+gD619kYAEAAOAoZGABAADi1YXAd8tGzMjAAgAAwFHIwAIAANgUGBhgLr7g8tFykxIysAAAAHAUMrAAAAA2UQPrX2RgAQAA4ChkYAEAAGyiD6x/kYEFAACAo5CBBQAAsIkaWP8iAwsAAABHIQMLAABgEzWw/kUGFgAAAI5CBhYAAMAmMrD+RQYWAAAAjkIGFgAAwCa6EPgXGVgAAAA4CgEsAACATQH6X4CPLmIvBbtu3Tpp3Lix5MqVy9x+4cKFXvPbtWsXZR2PPfaY15i//vpL2rZtKxkzZpTMmTNLx44d5fLly15jdu3aJTVq1JDUqVNL3rx5ZcyYMVG2Zf78+VK8eHEzpkyZMrJ06VKv+S6XS4YMGSI5c+aUNGnSSL169eTAgQNiFwEsAACAg125ckXKlSsnU6ZMuesYDVhPnTrlvnz66ade8zV43bt3r6xYsUIWL15sguIuXbq451+6dEkaNGgg+fPnl61bt8rYsWNl2LBh8t5777nHrF+/Xlq3bm2C3+3bt0vTpk3NZc+ePe4xGvROmjRJpk2bJps2bZJ06dJJw4YN5fr167buMzWwAAAADq6Bffzxx80lJsHBwRIWFhbtvF9++UWWLVsmP/30k1SuXNlMe+edd+SJJ56Qt956y2R2Z8+eLTdv3pQZM2ZIqlSppFSpUrJjxw4ZN26cO9CdOHGiCZT79etnro8cOdIExJMnTzYBq2ZfJ0yYIIMHD5YmTZqYMR999JGEhoaarHGrVq3ifJ/JwAIAACRxa9askZCQEClWrJh069ZN/vzzT/e8DRs2mLIBK3hVums/MDDQZEmtMTVr1jTBq0Uzp/v375fz58+7x+jtPOkYna4OHz4sp0+f9hqTKVMmqVq1qntMXJGBBQAASIR9YHW3feQsql7s0qxos2bNpGDBgnLo0CF59dVXTcZWg8agoCATVGpw6ylFihSSNWtWM0/pX729J82cWvOyZMli/lrTPMd4LsPzdtGNiSsCWAAAgERID5TyNHToUFN3apfnrnk9sKps2bJSuHBhk5WtW7euOBEBLAAAQCKsgT1+/LjpCmCJT/Y1OoUKFZLs2bPLwYMHTQCrtbFnz571GnP79m3TmcCqm9W/Z86c8RpjXY9tjOd8a5p2IfAcU758ebGDGlgAAIBESINXz0tCBbAnTpwwNbBWEFm9enW5cOGC6S5gWbVqlURERJj6VGuMdia4deuWe4weoKU1tVo+YI1ZuXKl17p0jE5XWoKgQaznGC2T0Dpba0xcEcACAADY5LMesPGorb18+bLpCKAX62Ap/fexY8fMPO0KsHHjRjly5IgJHrUDQJEiRcwBVqpEiRKmTrZz586yefNm+fHHH6VHjx6m9EA7EKg2bdqYA7i0RZa225o7d67pOtC3b1/3dvTq1ct0M3j77bdl3759ptxhy5YtZlnWY9a7d28ZNWqUfP3117J792554YUXzDq03ZYdlBAAAAA42JYtW6ROnTru61ZQGR4eLlOnTjUnIJg1a5bJsmqwqP1ctcWVZ0ZX22RpoKklBdp9oHnz5qZfq2e3gOXLl0v37t2lUqVKpgRBT0jg2Sv2oYcekjlz5pg2WXqgWNGiRU17rNKlS7vH9O/f3/St1dvp9jzyyCMm6NUTH9gR4NKmXA6laWd9QIPLdJaAoP+1dQCQNJ3b+L8PUwBJm37H5w7JIhcvXvSqA00ssUfF1xdLUOp0PlnHnetXZNvIJxPdfU9MKCEAAACAo1BCAAAAkAj7wOLuyMACAADAUcjAAgAA2OXDPrC6bMSMDCwAAAAchQwsAACATdTA+hcZWAAAADgKGVgAAACbNEnqq0QpCdjYkYEFAACAo5CBBQAAsIkaWP8iAwsAAABHIQMLAABgEzWw/kUGFgAAAI5CBhYAAMAmamD9iwwsAAAAHIUMLAAAgE1kYP2LDCwAAAAchQwsAACATXQh8C8ysAAAAHAUMrAAAAA2UQPrX2RgAQAA4ChkYAEAAGyiBta/yMACAADAUcjAAgAA2EQNrH8RwAIAANikIabPSgh8s9gkhRICAAAAOAoZWAAAAJsCAwLMxVfLRszIwAIAAMBRyMACAADYRBst/yIDCwAAAEchAwsAAGATbbT8iwwsAAAAHIUMLAAAgE2BAf9cfLVsxIwMLAAAAByFDCwAAIBdpgsBp+LyFzKwAAAAcBQysAAAADbRB9a/yMACAADAUcjAAgAA2BTw///5atmIGRlYAAAAOAoZWAAAAJvoA+tfZGABAADgKGRgAQAAbNIesL7qA+uz/rJJCBlYAAAAOAoZWAAAAJvoA+tfZGABAADgKGRgAQAAbAoMCDAXXy0bMSMDCwAAAEchgAUAAIhnDayvLnasW7dOGjduLLly5TIdDBYuXOied+vWLRkwYICUKVNG0qVLZ8a88MILcvLkSa9lFChQwN1Zwbq8+eabXmN27dolNWrUkNSpU0vevHllzJgxUbZl/vz5Urx4cTNG17l06VKv+S6XS4YMGSI5c+aUNGnSSL169eTAgQP27jABLAAAgLNduXJFypUrJ1OmTIky7+rVq7Jt2zZ5/fXXzd8vv/xS9u/fL0899VSUsSNGjJBTp065Lz179nTPu3TpkjRo0EDy588vW7dulbFjx8qwYcPkvffec49Zv369tG7dWjp27Cjbt2+Xpk2bmsuePXvcYzTonTRpkkybNk02bdpkguqGDRvK9evXbd1namABAAAc3Af28ccfN5foZMqUSVasWOE1bfLkyVKlShU5duyY5MuXzz09Q4YMEhYWFu1yZs+eLTdv3pQZM2ZIqlSppFSpUrJjxw4ZN26cdOnSxYyZOHGiPPbYY9KvXz9zfeTIkWbduj4NWDX7OmHCBBk8eLA0adLEjPnoo48kNDTUZI1btWoV5/tMBhYAACAR0qyn5+XGjRsJstyLFy+aIDlz5sxe07VkIFu2bFKhQgWTYb19+7Z73oYNG6RmzZomeLVo5lSzuefPn3eP0ZIATzpGp6vDhw/L6dOnvcZogF21alX3mLgiAwsAAGDT/egDq3WmnoYOHWp2298L3VWvNbG6qz9jxozu6S+99JJUrFhRsmbNakoBBg0aZMoINMOqNPAsWLCg17I0c2rNy5Ili/lrTfMco9OtcZ63i25MXBHAAgAAJELHjx/3CjKDg4PvaXm3bt2SZ5991uzKnzp1qte8vn37uv9dtmxZk2l98cUXZfTo0fe8Xl+ghAAAACCefWB9dVEavHpe7iWQvPX/wevRo0dNXapnYBwd3a2vJQRHjhwx17U29syZM15jrOtW3ezdxnjO97xddGMSNAP79ddfx3mB0R3VBgAAAP+49f/Bq7arWr16talzjY0eoBUYGCghISHmevXq1eW1114zy0qZMqWZpoFwsWLFTPmANWblypXSu3dv93J0jE5XWoKggaqOKV++vJmmtb3ajaBbt24JH8BqC4S40ILgO3fu2NoAAAAAp9Ecqa/Ol2V3uZcvX5aDBw+6r+vBUhqAaj2r9ltt0aKFaaG1ePFiE6dZ9aY6X0sF9AAqDSLr1KljOhHo9T59+shzzz3nDk7btGkjw4cPNy2ytIZWW2Np14Hx48e719urVy+pVauWvP3229KoUSP57LPPZMuWLe5WWxonanA7atQoKVq0qAlotb2X9qaNa6xpK4CNiIiwtVAAAADcH1u2bDHBZ+R61vDwcHPQl7Un3cp6WjQbW7t2bVOaoMGmjtVOBxpYagDrWRer3QKWL18u3bt3l0qVKkn27NnNCQmsFlrqoYcekjlz5pg2Wa+++qoJUrU9VunSpd1j+vfvb/rW6u0uXLggjzzyiCxbtsyc+MCOAJdW8t7DkWx2V5iQNO2sD2hwmc4SEPS/tg4AkqZzGyf5exMA3Mfv+NwhWUzLp9jqNf0RezSf9r2kTJPeJ+u4de2yfNG1RqK7744+iEtTz9qYNnfu3JI+fXr57bffzHRNAX/wwQe+2EYAAAAg/gHsv//9b5k5c6Y5FZhnM1tND0+fPt3u4gAAABwnMMC3FyRwAKun/NJi3LZt20pQUJB7up6Dd9++fXYXBwAAAPj2RAa///67FClSJNoDvbS1AgAAQFKnR9TrxVfLRgJnYEuWLCnff/99lOmff/65OXcuAAAAkKgysNoyQdsyaCZWs65ffvml7N+/35QWaH8xAACA5IBEqYMysE2aNJFFixbJd999J+nSpTMB7S+//GKm1a9f3zdbCQAAAMQ3A6tq1KhhTg0GAACQHFED68AA1jrrg2ZerbpYPSsDAAAAkOgC2BMnTkjr1q3lxx9/lMyZM5tpeiowPX2YnoYsT548vthOAACARMOX/VrpA+uDGthOnTqZdlmaff3rr7/MRf+tB3TpPAAAACBRZWDXrl0r69evl2LFirmn6b/feecdUxsLAACQ1FED67AMbN68eaM9YcGdO3ckV65cCbVdAAAAQMIEsGPHjpWePXuag7gs+u9evXrJW2+9ZXdxAAAAjhPg4wsSoIQgS5YsXunsK1euSNWqVSVFin9ufvv2bfPvDh06SNOmTeOySAAAAMB3AeyECRPit3QAAIAkKDAgwFx8tWwkQACrp44FAAAAHH0iA3X9+nW5efOm17SMGTPe6zYBAAAkapok9VWilASsDw7i0vrXHj16SEhIiKRLl87Ux3peAAAAgEQVwPbv319WrVolU6dOleDgYJk+fboMHz7ctND66KOPfLOVAAAAibAPrK8uSOASgkWLFplAtXbt2tK+fXtz8oIiRYpI/vz5Zfbs2dK2bVu7iwQAAAB8l4HVU8cWKlTIXe+q19Ujjzwi69ats7s4AAAAx9bA+uqCBA5gNXg9fPiw+Xfx4sVl3rx57sxs5syZ7S4OAAAA8G0JgZYN7Ny5U2rVqiUDBw6Uxo0by+TJk83pZceNG2d3cQAAAI5DH1iHBbB9+vRx/7tevXqyb98+2bp1q6mDLVu2bEJvHwAAAJBwfWCVHrylFwAAgOSCPrAOCGAnTZoU5wW+9NJL97I9AAAAwL0HsOPHj4/LMNO3jAAWAAAkdb7s10of2AQKYK2uA4nVsTVvcQpbIBm4fSfC35sAAEgKNbAAAADJsQ9poA+XjZgRwAIAANhECYF/EeQDAADAUcjAAgAA2KRJ0kDaaPkNGVgAAAAk/QD2+++/l+eee06qV68uv//+u5n28ccfyw8//JDQ2wcAAJDoaPbVlxckcAD7xRdfSMOGDSVNmjSyfft2uXHjhpl+8eJFeeONN+wuDgAAAPBtADtq1CiZNm2avP/++5IyZUr39Icffli2bdtmd3EAAACO7ULgqwsSOIDdv3+/1KxZM8r0TJkyyYULF+wuDgAAAPBtABsWFiYHDx6MMl3rXwsVKmR3cQAAAI5DDazDAtjOnTtLr169ZNOmTSbFffLkSZk9e7a88sor0q1bN99sJQAAABDfPrADBw6UiIgIqVu3rly9etWUEwQHB5sAtmfPnnYXBwAA4DhapuqrUlVKYH0QwGrW9bXXXpN+/fqZUoLLly9LyZIlJX369HYXBQAAANy/M3GlSpXKBK4AAADJTWBAgLn4atlI4AC2Tp06MbZ3WLVqld1FAgAAAL4LYMuXL+91/datW7Jjxw7Zs2ePhIeH210cAACAI4+CD/ThspHAAez48eOjnT5s2DBTDwsAAAD4UoIF+c8995zMmDEjoRYHAACQ6LsQ+OqC+xTAbtiwQVKnTp1QiwMAAEAcrFu3Tho3biy5cuUyxyktXLjQa77L5ZIhQ4ZIzpw5JU2aNFKvXj05cOCA15i//vpL2rZtKxkzZpTMmTNLx44do+xZ37Vrl9SoUcPEe3nz5pUxY8ZE2Zb58+dL8eLFzZgyZcrI0qVLbW+LTwLYZs2aeV2efvppqVatmrRv315efPFF2xsAAADgNIHyTxcCn1zEXgr2ypUrUq5cOZkyZUq08zXQnDRpkkybNs2ciCpdunTSsGFDuX79unuMBq979+6VFStWyOLFi01Q3KVLF/f8S5cuSYMGDSR//vyydetWGTt2rCkffe+999xj1q9fL61btzbB7/bt26Vp06bmosdJ2dmWuAhwaShsgwaqngIDAyVHjhzy6KOPmjt2P+mDmSlTJjnz50XziwFA0nb7ToS/NwHAffyOzx2SRS5eTFzf8Vbs0e/zbRKczjc98G9cuSxjW1SM130PCAiQBQsWmMBRaZinmdmXX37ZnHRK6XJDQ0Nl5syZ0qpVK/nll19Ma9SffvpJKleubMYsW7ZMnnjiCTlx4oS5/dSpU815AE6fPm1aqVont9Js7759+8z1li1bmmBaA2CLJjm1AYAGrHHZFp8cxHXnzh0TwGpKOEuWLHZuCgAAkGQ45Uxchw8fNkGn7qq3aABetWpVU/6pQaP+1bIBK3hVOl6TlJol1b3tOkbPvmoFr0ozp//5z3/k/PnzJi7UMX379vVav46xShrisi0+KSEICgoyWdYLFy7YuRkAAADike31vNy4ccP2Mk6fPm3+apbTk1635unfkJAQr/kpUqSQrFmzeo2Jbhme67jbGM/5sW2Lz2pgS5cuLb/99pvdmwEAACQZgQG+vSg9UEozlNZl9OjR/r7bzu0DO2rUKFO3MHLkSKlUqZIpvvWUmOpUAAAAnOr48eNecVVwcLDtZYSFhZm/Z86cMUf+W/S6dXIqHXP27Fmv292+fdt0JrBur3/1Np6s67GN8Zwf27YkeAZ2xIgRpjBXC3p37twpTz31lOTJk8fUPOhFayeoiwUAAMmB1qn6qguBVQOrwavnJT4BbMGCBU3guHLlSvc0LUfQ2tbq1aub6/pXy0O1u4Bl1apVEhERYepTrTHamUDPwGrRjgXFihVzx386xnM91hhrPXHZlgTPwA4fPly6du0qq1evtrUCAAAA+M7ly5fl4MGD7ut6sNSOHTtMDWu+fPmkd+/eZg960aJFTRD5+uuvm24AVqeCEiVKyGOPPSadO3c23QI0SO3Ro4c5qErHqTZt2phYUFtkDRgwwLTGmjhxotcZWnv16iW1atWSt99+Wxo1aiSfffaZbNmyxd1qSzskxLYtCR7AWt22dMMAAACSs8TUhWDLli1Sp04d93WrE0B4eLhpT9W/f3+zF137umqm9ZFHHjFtsjxPQDV79mwTtNatW9d0H2jevLnp12rRGtzly5dL9+7dTQlp9uzZzQkJPHvFPvTQQzJnzhwZPHiwvPrqqyZI1Q4EevyUJS7bkqB9YPXOaI2C9nxNLOgDCyQv9IEFko/E3gf21YXbJHW6DD5Zx/Urf8sbTePXBza5sHUQ1wMPPGDSvzHRgl8AAICkzLNbgC+WjQQMYLX2QX91AAAAAI4IYLWYN3KjWwAAgOQm4P//89WykUBttGIrHQAAAADuB9tdCAAAAJI7amAdEsBqM1sAAADAcaeSBQAASO7IwDqkBhYAAABIDMjAAgAA2KQHt/vqAHcOnI8dGVgAAAA4ChlYAAAAm6iB9S8ysAAAAHAUMrAAAAA2aZmqr0pVKYGNHRlYAAAAOAoZWAAAAJsCAwLMxVfLRszIwAIAAMBRyMACAADYRBcC/yIDCwAAAEchAwsAAGCXD7sQ6LIRMzKwAAAAcBQysAAAADYFSoC5+GrZiBkZWAAAADgKGVgAAACbOBOXf5GBBQAAgKOQgQUAALCJPrD+RQYWAAAAjkIGFgAAwKbAgABz8dWyETMysAAAAHAUMrAAAAA20YXAv8jAAgAAwFHIwAIAAMTnTFy+qoHlTFyxIgMLAAAARyEDCwAAYBM1sP5FBhYAAACOQgYWAAAgHhlAX2UByS7GjscIAAAAjkIGFgAAwKaAgABz8dWyETMysAAAAHAUMrAAAAA2aY7UV3lS8q+xI4AFAACwSU9i4LMTGVBCECtKCAAAAOAoZGABAADigTyp/5CBBQAAgKOQgQUAALCJU8n6FxlYAAAAOAoZWAAAAJs4kYF/kYEFAACAo5CBBQAAiEcG0FdZQLKLseMxAgAAcLACBQq4Sxo8L927dzfza9euHWVe165dvZZx7NgxadSokaRNm1ZCQkKkX79+cvv2ba8xa9askYoVK0pwcLAUKVJEZs6cGWVbpkyZYrYnderUUrVqVdm8ebNP7jMZWAAAAAfXwP70009y584d9/U9e/ZI/fr15ZlnnnFP69y5s4wYMcJ9XQNVi95Wg9ewsDBZv369nDp1Sl544QVJmTKlvPHGG2bM4cOHzRgNfGfPni0rV66UTp06Sc6cOaVhw4ZmzNy5c6Vv374ybdo0E7xOmDDBzNu/f78JihMSGVgAAAAHy5Ejhwk+rcvixYulcOHCUqtWLa+A1XNMxowZ3fOWL18uP//8s3zyySdSvnx5efzxx2XkyJEmm3rz5k0zRoPSggULyttvvy0lSpSQHj16SIsWLWT8+PHu5YwbN84Eyu3bt5eSJUua2+h6Z8yYkeD3mQAWAADApgAfX9SlS5e8Ljdu3Ih1u27evGkC0Q4dOnhlcjVrmj17dildurQMGjRIrl696p63YcMGKVOmjISGhrqnaeZU17l37173mHr16nmtS8fodGu9W7du9RoTGBhorltjEhIlBAAAAIlQ3rx5va4PHTpUhg0bFuNtFi5cKBcuXJB27dq5p7Vp00by588vuXLlkl27dsmAAQPMbv0vv/zSzD99+rRX8Kqs6zovpjEa5F67dk3Onz9vShGiG7Nv3z5JaASwAAAAibAG9vjx4167+vXgqdh88MEHpgRAg1VLly5d3P/WTKvWrdatW1cOHTpkSg2ciAAWAAAgEdLg1TOAjc3Ro0flu+++c2dW70YPsFIHDx40AazWxEbuFnDmzBnzV+dZf61pnmN0+9KkSSNBQUHmEt0YaxkJiRpYAACAePaB9dUlPj788ENztL92C4jJjh07zF/NxKrq1avL7t275ezZs+4xK1asMMGpHoxljdHOA550jE5XqVKlkkqVKnmNiYiIMNetMQmJABYAAMDhIiIiTAAbHh4uKVL8bwe7lgloRwE9wOrIkSPy9ddfmxZZNWvWlLJly5oxDRo0MIHq888/Lzt37pRvv/1WBg8ebPrIWmUL2j7rt99+k/79+5ua1nfffVfmzZsnffr0ca9LW2i9//77MmvWLPnll1+kW7ducuXKFdOVIKFRQgAAAODgPrBKSwf0ZATafcCTZkZ1nvZk1WBSDwxr3ry5CVAtuutfW29pwKnZ0nTp0plA2LNvrLbQWrJkiQlYJ06cKHny5JHp06e7e8Cqli1byrlz52TIkCHmoC9tybVs2bIoB3YlhACXy+USh9Ij3zJlyiRn/rxoq0YEgDPdvhPh700AcB+/43OHZJGLFxPXd7wVe3zy46+SNn0Gn6zj6uW/5bmHH0h09z0xIQMLAABgk2e/Vl8sGzGjBhYAAACOQgYWAADAJi1T9VEJrM+Wm5SQgQUAAICjkIEFAACwKVACzMVXy0bMyMACAADAUcjAAgAA2EQNrH+RgQUAAICjkIEFAACwKeD///PVshEzMrAAAABwFDKwAAAANlED619kYAEAAOAoZGABAADiUafqq36t1MDGjgwsAAAAHIUMLAAAgE3UwPoXGVgAAAA4ChlYAAAAm8jA+hcZWAAAADgKGVgAAACbOBOXf5GBBQAAgKOQgQUAALApMOCfi6+WjZiRgQUAAICjkIEFAACwiRpY/yIDCwAAAEchAwsAAGATfWD9iwwsAAAAHIUMLAAAgE2aJPVdDSxiQwYWAAAAjkIGFgAAwCb6wPoXGVgAAAA4ChlYAAAAm+gD619kYAEAAOAoBLDwi7H/GS0PV3tQcmTJIPlyhcgzzZvKr/v3e425fv269O7ZXXKHZpPsmdNLq2eby5kzZ7zG9O39kjxUpZJkShcsVSuVj7IeXWbDenUkf+5QyZw+tZR4oJAMGzJYbt265fP7COB/Tv7+u3Rq97zky5VDcmROJ1UrlZNtW7dEO7ZXj26SIXWQTHlnotf0Awd+lZYtmkr+3CGSK0dmqV+npqxbs9przNYtP8mTj9WXPKFZJW9YNmn65GOye9dOn943JO8+sL66IBEHsOvWrZPGjRtLrly5JCAgQBYuXOjPzcF99P26tdK1W3dZ+8NGWfzNCrl965Y8+UQDuXLlintM/5f7yJIli2T2Z/Nl+cq1curkSWn1TLMoy3qhXQdp8UzLaNeTMmVKafvcC7Jo6XLZuXe/jH17gnz4wfsycvhQn94/AP9z/vx5qV+nhqRImVK+/GqJ/LR9j7zx5ljJnDlLlLFff7VAftq8SXLmyhVl3jNPPyW3b9+WJcu+k3UbfpIyZcvKM82ekjOnT5v5ly9flqefekLy5M0rq77fIMtXrZP0GTJI08aP86MVSGL8WgOrwUq5cuWkQ4cO0qxZ1MAESdfXS5Z5XX/vg5kmE7t921Z5pEZNuXjxosz88AOZ+fEcqV3n0X/GTP9QypcpIZs2bpSq1aqZaeMmTDJ///jjnOzZvSvKegoWKmQulvz588u6tWvkxx++9/E9BGAZ//YYyZ0nr0x7f4Z7WoGCBaPN0vbr20sWLvpGWjRt7DXvjz/+kEMHD8iUae9L6TJlzbTho0bL+/+dKj/v3SOhYWHy6/59cv6vv2TwkOEmiFWDXhsi1SqXl2PHjkrhwkV8fl+R3PrA+m7ZSMQZ2Mcff1xGjRolTz/9tD83A4nApYsXzd8sWbKavxrIasbk0br13GOKFS8uefPlk00bN8R7PYcOHpQVy5dJjZq1EmCrAcTF0sWLpGKlSvJ8m2elYN4webhqJbMnxFNERIR07hAuvfq8IiVKloqyjGzZsknRB4rJp7M/NskPzcTOmP6e5AgJkfIVK5kxOj9rtmzy0cwZcvPmTbl27Zr5d7HiJSR//gL37f4C8D1qYOF3+sXV7+XeUv2hh6VU6dJm2unTpyVVqlSSOXNmr7EhIaFy5sw/uwvtqF3jIVMDW7pEUXn44RoyZNiIBNt+ADE7cvg3mf7eNClcuKjJrnbs/KL0f7m3zP54lnvMuLfGSIoUQdKte89ol6FlZloKtGvnDsmZPZNkz5RWJk8aLwu+XipZsvxTipAhQwb5ZvkqmfvpbFNnG5Yto6xY/q0pW0iRgqY7SFiBEiCBAT66kINNWgHsjRs35NKlS14XOJ8eqLV37x75aPZnPlvHx3PmyobN20xJwjffLJHx497y2boARP2RWq5CRRk28t9SrnwF6dCpi7Tr0Ek+mP6ee4/L1CmTZNr7H5pANToul0v69u4hOXLkkG9XrpU1P2yUJxs3kWebN5HTp06ZMZpx7d61s1St/pCsWrdeVqz+XkqWKiUtnm5s5gFIOhz1k3T06NEyfPhwf28GElDvl3rI0qWL5btV6yRPnjzu6WFhYWYX4IULF7yysGfPnpHQ0DDb68n7//VwJUqWlIg7d6R7ty7Su8/LEhQUlED3BMDdhIXllOLFS3hN05KgrxZ+af69/scf5NzZs1Ki6P9289+5c0deHfCKvPvORNn762+ydvUqWbZ0iRw//adkzJjRjClfoaKsXvmdzP7kI3m53wCZ99kcOXr0iKxc+6MEBv6Tn5kxa7bpRrBk0VfS4tlW9/V+I2mjBta/HJWBHTRokDm4x7ocP37c35uEeNJsigavesTxsuWrohzQUaFiJdNBYPWqlV4tsY4fOyZVq1W/52yQ1tfqXwC+V636Q3Lg11+9ph08cEDy5stv/t2qzXOyccsOWb95m/uiXQh69X1FFiz+xoy5eu2q+WsFppaAwED3e1mzrDrfM4trXef9DiQtjsrABgcHmwuSRtnA3M/myPwvvzJtbrTmVWXKlEnSpElj/rZr31EG9OsrWbNmlQwZMkrf3j1N8Gp1ILAOytLWOdpG59r1a7Jzxw53plVraD+dM9sEwqVLlzGvna1bt8jrgweZtls6HYDvdX+pt9Sr/Yjp/9ysxTOy9afN5iCuSVOmuQ/Q0ounlClSmr0tDzxQzFyvUrW6ZM6SRV7s1E4Gvvq6pE6TRmbOmC5HjxyWxx5/woypU7eeDB7UX/r26iEv/quHuCIiZNzY/5j615q16vjhniNJIwWbfANYDTwOHjzovn748GHZsWOHCVjy5cvnz02Dj73336nmb4O6tb2nT/9Qng9vZ/495u3xJnvS+tnmpv65XoOGMvGdd73Gd3uxk+kpa6n2YAXzd9+Bw5K/QAHzxaVfYNoAXbO++fLnl27/6iE9e/W5D/cSgKpU+UGZM+8LGfb6a/KfN0ZK/gIF5c2x46Rl67ZxXkb27NnNAVsjhg6WRo/VM72ji5csJZ99vkDKlC1nxhQrVlzmffGVjP73SKlX62Hz+VG2XAX58uulEpYzpw/vIYD7LcCl3+p+smbNGqlTJ+qv4vDwcJk5c2ast9eDuDRTd+bPi+6aKABJ1+077AYGkgv9js8dksWUDCam73gr9li5/Ziky+Cb7bry9yWpWyFforvviYlfM7C1a9c2WTEAAAAgSdbAAgAAJAoB2p/Yd8tGEupCAAAAAJCBBQAAsIkmBP5FBhYAACC+EayvLjYMGzbM9Dv2vBQvXtw9//r169K9e3fTri59+vTSvHlzOXPmjNcyjh07Jo0aNZK0adNKSEiI9OvXT27fvh3l4PuKFSuatpRFihSJ9oD7KVOmSIECBSR16tRStWpV2bx5s/gCASwAAIDDlSpVSk6dOuW+/PDDD+55ffr0kUWLFsn8+fNl7dq1cvLkSWnWrJnXme80eNUzYK5fv15mzZplgtMhQ4Z4tTrVMdo9Slue9u7dWzp16iTffvute8zcuXOlb9++MnToUNm2bZuUK1dOGjZsKGfPnk1abbTuFW20gOSFNlpA8pHY22it3nlc0vuojdblvy9JnXJ543zfhw0bJgsXLjSBZWS6jBw5csicOXOkRYsWZtq+ffukRIkSsmHDBqlWrZp888038uSTT5rANjQ01IyZNm2aDBgwQM6dO2dODKT/XrJkiezZs8e97FatWplTvi9btsxc14zrgw8+KJMnTzbX9Qx4eir3nj17ysCBAyUhkYEFAABwuAMHDkiuXLmkUKFC0rZtW1MSoLZu3WpOn16vXj33WC0v0BNGaQCr9G+ZMmXcwavSzKkG63v37nWP8VyGNcZahmZvdV2eY/RkInrdGpOQOIgLAADApgAfttGylqsBpCetPdVLZJr51F3+xYoVM+UDw4cPlxo1aphsqZ6qXTOomTNn9rqNBqvWadz1r2fwas235sU0Rrfx2rVrcv78eVOKEN0YzfgmNAJYAACAREh3v3vS2lItF4js8ccfd/+7bNmyJqDNnz+/zJs3T9KkSSNJEQEsAABAImyjdfz4ca8a2Oiyr9HRbOsDDzwgBw8elPr165vd+1qr6pmF1S4EYWFh5t/6N3K3AKtLgeeYyJ0L9LpunwbJQUFB5hLdGGsZCYkaWAAAgERIg0PPS1wD2MuXL8uhQ4ckZ86cUqlSJUmZMqWsXLnSPX///v2mRrZ69ermuv7dvXu3V7eAFStWmHWWLFnSPcZzGdYYaxlapqDr8hyjB3HpdWtMQiIDCwAA4OAzGbzyyivSuHFjUzagnQS01ECzoa1btzYdEzp27GjaW2XNmtUEpdoVQINK7UCgGjRoYALV559/XsaMGWPqXQcPHmx6x1pBc9euXU13gf79+0uHDh1k1apVpkRBOxNYdB3h4eFSuXJlqVKlikyYMEGuXLki7du3T9jHhwAWAADA2U6cOGGC1T///NO0zHrkkUdk48aN5t9q/PjxpiOAnsDgxo0bpnvAu+++6769BruLFy+Wbt26mcA2Xbp0JhAdMWKEe0zBggVNsKo9ZSdOnCh58uSR6dOnm2VZWrZsadpuaf9YDYLLly9vWmxFPrArIdAHFoBj0AcWSD4Sex/YdbtP+LQPbM0yeRLdfU9MqIEFAACAo1BCAAAAkAj7wOLuyMACAADAUcjAAgAAOLcJQbJEBhYAAACOQgYWAADALlKwfkUGFgAAAI5CBhYAAMCmgP//z1fLRszIwAIAAMBRyMACAADYRB9Y/yIDCwAAAEchAwsAAGATTQj8iwwsAAAAHIUMLAAAgF2kYP2KDCwAAAAchQwsAACATfSB9S8ysAAAAHAUMrAAAAA20QfWv8jAAgAAwFHIwAIAANhEEwL/IgMLAAAARyEDCwAAYBcpWL8iAwsAAABHIQMLAABgE31g/YsMLAAAAByFDCwAAIBN9IH1LzKwAAAAcBQysAAAADbRhMC/yMACAADAUcjAAgAA2EUK1q/IwAIAAMBRyMACAADYRB9Y/yIDCwAAAEchAwsAAGCXD/vAkoCNHRlYAAAAOAoZWAAAAJtoQuBfZGABAADgKGRgAQAA7CIF61dkYAEAAOAoZGABAABsog+sf5GBBQAAgKOQgQUAALApwId9YH3WXzYJIQMLAAAARyEDCwAAYBNNCPyLDCwAAAAchQwsAACAXaRg/YoMLAAAAByFDCwAAIBN9IH1LzKwAAAADjZ69Gh58MEHJUOGDBISEiJNmzaV/fv3e42pXbu2BAQEeF26du3qNebYsWPSqFEjSZs2rVlOv3795Pbt215j1qxZIxUrVpTg4GApUqSIzJw5M8r2TJkyRQoUKCCpU6eWqlWryubNmxP8PhPAAgAAxKcENsBHF5vbsnbtWunevbts3LhRVqxYIbdu3ZIGDRrIlStXvMZ17txZTp065b6MGTPGPe/OnTsmeL1586asX79eZs2aZYLTIUOGuMccPnzYjKlTp47s2LFDevfuLZ06dZJvv/3WPWbu3LnSt29fGTp0qGzbtk3KlSsnDRs2lLNnz0pCCnC5XC5xqEuXLkmmTJnkzJ8XJWPGjP7eHAA+dvtOhL83AcB9/I7PHZJFLl5MXN/xVuyx5/BZyeCj7fr70iUpXTAk3vf93LlzJoOqgW3NmjXdGdjy5cvLhAkTor3NN998I08++aScPHlSQkNDzbRp06bJgAEDzPJSpUpl/r1kyRLZs2eP+3atWrWSCxcuyLJly8x1zbhqNnjy5MnmekREhOTNm1d69uwpAwcOlIRCBhYAACCeTQh8dbGCZc/LjRs34rRtFy9eNH+zZs3qNX327NmSPXt2KV26tAwaNEiuXr3qnrdhwwYpU6aMO3hVmjnV9e7du9c9pl69el7L1DE6XWn2duvWrV5jAgMDzXVrTELhIC4AAIBESDOXnnS3/LBhw2K8jWY8ddf+ww8/bAJVS5s2bSR//vySK1cu2bVrl8mmap3sl19+aeafPn3aK3hV1nWdF9MYDXKvXbsm58+fN6UI0Y3Zt2+fJCQCWAAAAJuselVfLVsdP37cq4RAD5yKTffu3c0u/h9++MFrepcuXdz/1kxrzpw5pW7dunLo0CEpXLiwOA0lBAAAAImwiECDV89LbAFsjx49ZPHixbJ69WrJkydPjGO1VlUdPHjQ/A0LC5MzZ854jbGu67yYxui2pUmTxpQnBAUFRTvGWkZCIYAFAABwMJfLZYLXBQsWyKpVq6RgwYKx3ka7CCjNxKrq1avL7t27vboFaEcDDU5LlizpHrNy5Uqv5egYna70QK9KlSp5jdGSBr1ujUkolBAAAAAkwhKCuOrevbvMmTNHvvrqK9ML1qpZ1W4JmhnVMgGd/8QTT0i2bNlMDWyfPn1Mh4KyZcuasdp2SwPV559/3rTX0mUMHjzYLNvK/GrfWO0u0L9/f+nQoYMJlufNm2c6E1i0hVZ4eLhUrlxZqlSpYroeaDuv9u3bJ+RDRAALAADgZFOnTnW3yvL04YcfSrt27Uxm9LvvvnMHk3pwWPPmzU2AatFd/1p+0K1bN5MtTZcunQlER4wY4R6jmV0NVjX4nThxoilTmD59uulEYGnZsqVpu6X9YzUI1tZd2mIr8oFd94o+sAAcgz6wQPKR2PvA7jt6zqd9YIvnz5Ho7ntiQg0sAAAAHIUSAgAAAAfXwCZHZGABAADgKGRgAQAAbAr4//98tWzEjAwsAAAAHIUMLAAAgF3/O2GWb5aNGJGBBQAAgKOQgQUAALCJBKx/kYEFAACAo5CBBQAAsIk+sP5FBhYAAACOQgYWAADAJvrA+hcZWAAAADgKGVgAAAC7aEPgV2RgAQAA4ChkYAEAAGwiAetfZGABAADgKGRgAQAAbKIPrH+RgQUAAICjkIEFAACwzXd9YKmCjR0ZWAAAADgKGVgAAACbqIH1LzKwAAAAcBQCWAAAADgKASwAAAAchRpYAAAAm6iB9S8ysAAAAHAUMrAAAADx6gLrm1Sp7/rLJh1kYAEAAOAoZGABAABsogbWv8jAAgAAwFHIwAIAANikSVJfJUpJwMaODCwAAAAchQwsAACAXaRg/YoMLAAAAByFDCwAAIBN9IH1LzKwAAAAcBQysAAAADbRB9a/yMACAADAUcjAAgAA2EQTAv8iAwsAAABHIQMLAABgFylYvyIDCwAAAEchAwsAAGATfWD9iwwsAAAAHIUMLAAAgE30gfUvRwewLpfL/P370iV/bwqA++D2nQh/bwKA++Tvvy95fdcnNpd8GHv4ctlJhaMD2L///tv8LVIwr783BQAA+Oi7PlOmTJJYpEqVSsLCwqSoj2MPXYeuC9ELcCXWnzZxEBERISdPnpQMGTJIAPn2ZEN/mebNm1eOHz8uGTNm9PfmAPAh3u/Jl4YnGrzmypVLAgMT1yE7169fl5s3b/p0HRq8pk6d2qfrcDJHZ2D1BZ0nTx5/bwb8RL/M+EIDkgfe78lTYsq8etLAkuDSvxLXTxoAAAAgFgSwAAAAcBQCWDhOcHCwDB061PwFkLTxfgeQ5A7iAgAAQPJDBhYAAACOQgALAAAARyGABQAAgKMQwMJxpkyZIgUKFDA9+KpWrSqbN2/29yYB8IF169ZJ48aNTSN7PVnNwoUL/b1JABIJAlg4yty5c6Vv377mqORt27ZJuXLlpGHDhnL27Fl/bxqABHblyhXzHtcfrQDgiS4EcBTNuD744IMyefJk9+mE9TSTPXv2lIEDB/p78wD4iGZgFyxYIE2bNvX3pgBIBMjAwjH0vNNbt26VevXqeZ1OWK9v2LDBr9sGAADuHwJYOMYff/whd+7ckdDQUK/pev306dN+2y4AAHB/EcACAADAUQhg4RjZs2eXoKAgOXPmjNd0vR4WFua37QIAAPcXASwcI1WqVFKpUiVZuXKle5oexKXXq1ev7tdtAwAA90+K+7gu4J5pC63w8HCpXLmyVKlSRSZMmGBa7bRv397fmwYggV2+fFkOHjzovn748GHZsWOHZM2aVfLly+fXbQPgX7TRguNoC62xY8eaA7fKly8vkyZNMu21ACQta9askTp16kSZrj9iZ86c6ZdtApA4EMACAADAUaiBBQAAgKMQwAIAAMBRCGABAADgKASwAAAAcBQCWAAAADgKASwAAAAchQAWAAAAjkIACwAAAEchgAWQoNq1aydNmzZ1X69du7b07t3bL2dxCggIkAsXLtx1jM5fuHBhnJc5bNgwc/a3e3HkyBGzXj0lKgAgfghggWQSVGrQpJdUqVJJkSJFZMSIEXL79m2fr/vLL7+UkSNHJljQCQBACn9vAID747HHHpMPP/xQbty4IUuXLpXu3btLypQpZdCgQVHG3rx50wS6CSFr1qwJshwAACxkYIFkIjg4WMLCwiR//vzSrVs3qVevnnz99ddeu/3//e9/S65cuaRYsWJm+vHjx+XZZ5+VzJkzm0C0SZMmZhe45c6dO9K3b18zP1u2bNK/f39xuVxe641cQqAB9IABAyRv3rxmmzQb/MEHH5jl1qlTx4zJkiWLycTqdqmIiAgZPXq0FCxYUNKkSSPlypWTzz//3Gs9GpQ/8MADZr4ux3M740q3S5eRNm1aKVSokLz++uty69atKOP++9//mu3Xcfr4XLx40Wv+9OnTpUSJEpI6dWopXry4vPvuu7a3BQBwdwSwQDKlgZ5mWi0rV66U/fv3y4oVK2Tx4sUmcGvYsKFkyJBBvv/+e/nxxx8lffr0JpNr3e7tt9+WmTNnyowZM+SHH36Qv/76SxYsWBDjel944QX59NNPZdKkSfLLL7+YYFCXqwHhF198Ycbodpw6dUomTpxormvw+tFHH8m0adNk79690qdPH3nuuedk7dq17kC7WbNm0rhxY1Nb2qlTJxk4cKDtx0Tvq96fn3/+2az7/fffl/Hjx3uNOXjwoMybN08WLVoky5Ytk+3bt8u//vUv9/zZs2fLkCFDzI8BvX9vvPGGCYRnzZple3sAAHfhApDkhYeHu5o0aWL+HRER4VqxYoUrODjY9corr7jnh4aGum7cuOG+zccff+wqVqyYGW/R+WnSpHF9++235nrOnDldY8aMcc+/deuWK0+ePO51qVq1arl69epl/r1//35Nz5r1R2f16tVm/vnz593Trl+/7kqbNq1r/fr1XmM7duzoat26tfn3oEGDXCVLlvSaP2DAgCjLikznL1iw4K7zx44d66pUqZL7+tChQ11BQUGuEydOuKd98803rsDAQNepU6fM9cKFC7vmzJnjtZyRI0e6qlevbv59+PBhs97t27ffdb0AgJhRAwskE5pV1UynZlZ1l3ybNm3MUfWWMmXKeNW97ty502QbNSvp6fr163Lo0CGz21yzpFWrVnXPS5EihVSuXDlKGYFFs6NBQUFSq1atOG+3bsPVq1elfv36XtM1C1yhQgXzb810em6Hql69utg1d+5ckxnW+3f58mVzkFvGjBm9xuTLl09y587ttR59PDVrrI+V3rZjx47SuXNn9xhdTqZMmWxvDwAgegSwQDKhdaFTp041QarWuWqw6SldunRe1zWAq1SpktklHlmOHDniXbZgl26HWrJkiVfgqLSGNqFs2LBB2rZtK8OHDzelExpwfvbZZ6ZMwu62aulB5IBaA3cAQMIggAWSCQ1Q9YCpuKpYsaLJSIaEhETJQlpy5swpmzZtkpo1a7ozjVu3bjW3jY5meTVbqbWrehBZZFYGWA8Os5QsWdIEqseOHbtr5lYPmLIOSLNs3LhR7Fi/fr05wO21115zTzt69GiUcbodJ0+eND8CrPUEBgaaA99CQ0PN9N9++80EwwAA3+AgLgDR0gAse/bspvOAHsR1+PBh06f1pZdekhMnTpgxvXr1kjfffNOcDGDfvn3mYKaYergWKFBAwsPDpUOHDuY21jL1oCilAaR2H9Byh3PnzpmMpu6Wf+WVV8yBW3oglO6i37Ztm7zzzjvuA6O6du0qBw4ckH79+pld+XPmzDEHY9lRtGhRE5xq1lXXoaUE0R2Qpp0F9D5oiYU+Lvp4aCcC7fCgNIOrB53p7X/99VfZvXu3aV82btw4W9sDALg7AlgA0dIWUevWrTM1n3qEv2Y5tbZTa2CtjOzLL78szz//vAnotBZUg82nn346xuVqGUOLFi1MsKstprRW9MqVK2aelghoAKgdBDSb2aNHDzNdT4SgR/JrYKjboZ0QtKRA22op3UbtYKBBsbbY0m4FevS/HU899ZQJknWderYtzcjqOiPTLLY+Hk888YQ0aNBAypYt69UmSzsgaBstDVo146xZYw2mrW0FANy7AD2SKwGWAwAAANwXZGABAADgKASwAAAAcBQCWAAAADgKASwAAAAchQAWAAAAjkIACwAAAEchgAUAAICjEMACAADAUQhgAQAA4CgEsAAAAHAUAlgAAAA4CgEsAAAAxEn+D5mqnRFXIN6KAAAAAElFTkSuQmCC", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "# Plot confusion matrix using only matplotlib\n", "def plot_confusion_matrix(cm, classes, title='Confusion Matrix', cmap=plt.cm.Blues):\n", " plt.figure(figsize=(8, 6))\n", - " sns.heatmap(cm, annot=True, fmt='d', cmap=cmap,\n", - " xticklabels=classes, yticklabels=classes)\n", + " plt.imshow(cm, interpolation='nearest', cmap=cmap)\n", " plt.title(title)\n", + " plt.colorbar()\n", + " tick_marks = np.arange(len(classes))\n", + " plt.xticks(tick_marks, classes)\n", + " plt.yticks(tick_marks, classes)\n", + "\n", + " # Annotate cells with counts\n", + " thresh = cm.max() / 2.\n", + " for i in range(cm.shape[0]):\n", + " for j in range(cm.shape[1]):\n", + " plt.text(j, i, format(cm[i, j], 'd'),\n", + " ha=\"center\", va=\"center\",\n", + " color=\"white\" if cm[i, j] > thresh else \"black\")\n", + "\n", " plt.ylabel('True label')\n", " plt.xlabel('Predicted label')\n", + " plt.tight_layout()\n", " plt.show()\n", "\n", "cm = confusion_matrix(y_test, y_pred)\n", @@ -241,7 +1039,7 @@ ], "metadata": { "kernelspec": { - "display_name": "Python 3 (ipykernel)", + "display_name": "pyearthtools", "language": "python", "name": "python3" }, @@ -255,7 +1053,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.13.2" + "version": "3.13.5" } }, "nbformat": 4, diff --git a/notebooks/tutorial/HadISD/Data_Config.ipynb b/notebooks/tutorial/HadISD/Data_Config.ipynb index db28f31b..ad00c678 100644 --- a/notebooks/tutorial/HadISD/Data_Config.ipynb +++ b/notebooks/tutorial/HadISD/Data_Config.ipynb @@ -57,27 +57,42 @@ "metadata": {}, "outputs": [], "source": [ - "# A sample list of WMO number ranges. Users can find more at the official HadISD download page.\n", + "# For any station ranges you don't want to download, you can comment them out here\n", "wmo_id_ranges = [\n", - " \"000000-029999\",\n", - " \"080000-099999\",\n", - " \"200000-249999\",\n", - " \"720000-721999\",\n", + " #\"000000-029999\",\n", + " #\"030000-049999\",\n", + " #\"050000-079999\",\n", + " #\"080000-099999\",\n", + " #\"100000-149999\",\n", + " #\"150000-199999\",\n", + " #\"200000-249999\",\n", + " #\"250000-299999\",\n", + " #\"300000-349999\",\n", + " #\"350000-399999\",\n", + " #\"400000-449999\",\n", + " #\"450000-499999\",\n", + " \"500000-549999\",\n", + " #\"550000-599999\",\n", + " #\"600000-649999\",\n", + " #\"650000-699999\",\n", + " #\"700000-709999\",\n", + " #\"710000-714999\",\n", + " #\"715000-719999\",\n", + " #\"720000-721999\",\n", + " \"722000-722999\",\n", + " #\"723000-723999\",\n", + " #\"724000-724999\",\n", + " #\"725000-725999\",\n", + " #\"726000-726999\",\n", + " #\"727000-729999\",\n", + " #\"730000-799999\",\n", + " \"800000-849999\",\n", + " #\"850000-899999\",\n", + " #\"900000-949999\",\n", + " #\"950000-999999\",\n", "]" ] }, - { - "cell_type": "code", - "execution_count": null, - "id": "35617ad2", - "metadata": {}, - "outputs": [], - "source": [ - "\n", - "# User sets the WMO number range to download\n", - "wmo_id_range = \"080000-099999\" # Change this to the desired WMO range, either from the sample list or from the HadISD page." - ] - }, { "cell_type": "markdown", "id": "7aec321a", @@ -97,15 +112,15 @@ "# Set the date range to reindex the time coordinate\n", "DATE_RANGE = (\"1970-01-01T00\", \"2023-12-31T23\")\n", "# Set the input directory to the folder with raw NetCDFs\n", - "input_dir = download_dir / f\"WMO_{wmo_id_range}\" / \"netcdf\"\n", + "input_dir = download_dir / \"netcdf\"\n", "# Set the Zarr output directory to a sibling folder under the same WMO directory\n", - "zarr_output_dir = download_dir / f\"WMO_{wmo_id_range}\" / \"zarr\"" + "zarr_output_dir = download_dir / \"zarr\"" ] } ], "metadata": { "kernelspec": { - "display_name": "Python 3 (ipykernel)", + "display_name": "pet_tutorial", "language": "python", "name": "python3" }, @@ -119,7 +134,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.13.2" + "version": "3.11.11" } }, "nbformat": 4, diff --git a/notebooks/tutorial/HadISD/HadISD_QC_Exploration.ipynb b/notebooks/tutorial/HadISD/HadISD_QC_Exploration.ipynb index c26523fa..0b3694cf 100644 --- a/notebooks/tutorial/HadISD/HadISD_QC_Exploration.ipynb +++ b/notebooks/tutorial/HadISD/HadISD_QC_Exploration.ipynb @@ -15,10 +15,7 @@ "metadata": {}, "outputs": [], "source": [ - "import datetime\n", "import numpy as np\n", - "import pandas as pd\n", - "from pathlib import Path\n", "\n", "import pyearthtools.pipeline as petpipe\n", "import pyearthtools.data as petdata\n", @@ -27,12 +24,13 @@ }, { "cell_type": "code", - "execution_count": 16, + "execution_count": null, "id": "4a7da841", "metadata": {}, "outputs": [], "source": [ - "# %run HadISD_config.ipynb" + "# ruff: noqa: F821\n", + "%run Pipeline_Config.ipynb" ] }, { @@ -89,17 +87,7 @@ "metadata": {}, "outputs": [], "source": [ - "y = data_prep_pipe[\"1969-01-01T07\"]\n", - "y" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "3d3819b9", - "metadata": {}, - "outputs": [], - "source": [ + "x = data_prep_pipe[\"1969-01-01T07\"]\n", "x" ] }, @@ -110,7 +98,7 @@ "metadata": {}, "outputs": [], "source": [ - "qc = y[\"quality_control_flags\"].values\n" + "qc = x[\"quality_control_flags\"].values\n" ] }, { @@ -214,21 +202,13 @@ "# for qc, test 12, time 826, station 0, print the value of the test\n", "print(\"QC value for test 12, time 826, station 0:\", qc[0, 826, 12])" ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "8c8f1bc9", - "metadata": {}, - "outputs": [], - "source": [] } ], "metadata": { "kernelspec": { - "display_name": "Python 3 (ipykernel)", + "display_name": "pyearthtools", "language": "python", - "name": "python3" + "name": "pyearthtools" }, "language_info": { "codemirror_mode": { @@ -240,7 +220,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.13.2" + "version": "3.13.5" } }, "nbformat": 4, diff --git a/notebooks/tutorial/HadISD/Pipeline_Config.ipynb b/notebooks/tutorial/HadISD/Pipeline_Config.ipynb index e0dcc486..4b5298b3 100644 --- a/notebooks/tutorial/HadISD/Pipeline_Config.ipynb +++ b/notebooks/tutorial/HadISD/Pipeline_Config.ipynb @@ -2,12 +2,12 @@ "cells": [ { "cell_type": "code", - "execution_count": 1, + "execution_count": null, "id": "1abab3c3", "metadata": {}, "outputs": [], "source": [ - "# import pyearthtools.pipeline as petpipe" + "import pyearthtools.pipeline as petpipe" ] }, { @@ -15,7 +15,7 @@ "id": "2307aa00", "metadata": {}, "source": [ - "# Lists and Dictionaries Required For Sustom Pipeline Steps" + "# Lists and Dictionaries Required For Custom Pipeline Steps" ] }, { @@ -62,6 +62,16 @@ "source": [ "# Custom operation to remove redundent coordinates\n", "class SqueezeStationCoordinates(petpipe.Operation):\n", + " \"\"\"\n", + " Squeeze singleton dimensions from specified station-based coordinates in an xarray.Dataset.\n", + "\n", + " This operation is useful for removing unnecessary singleton dimensions (e.g., shape (n, 1))\n", + " from coordinates like latitude, longitude, and elevation, ensuring they are 1D and indexed\n", + " by 'station'.\n", + "\n", + " Args:\n", + " coords (tuple of str): Names of coordinates to squeeze. Defaults to (\"latitude\", \"longitude\", \"elevation\").\n", + " \"\"\"\n", " def __init__(self, coords=(\"latitude\", \"longitude\", \"elevation\")):\n", " super().__init__()\n", " self.coords = coords\n", @@ -74,7 +84,7 @@ " # Undo function added otherwise pyearthtools will complain\n", " def undo_func(self, ds):\n", " # No undo operation needed for this operation\n", - " return ds" + " return ds\n" ] }, { @@ -314,8 +324,6 @@ "metadata": {}, "outputs": [], "source": [ - "import xarray as xr\n", - "import numpy as np\n", "from pyearthtools.data.transforms.values import AddFlaggedObs\n", "\n", "def test_add_flagged_obs():\n", diff --git a/packages/data/src/pyearthtools/data/indexes/_indexes.py b/packages/data/src/pyearthtools/data/indexes/_indexes.py index d08abd03..b4c358ff 100644 --- a/packages/data/src/pyearthtools/data/indexes/_indexes.py +++ b/packages/data/src/pyearthtools/data/indexes/_indexes.py @@ -66,7 +66,7 @@ LOG = logging.getLogger("pyearthtools.data") - + class Index(CallRedirectMixin, CatalogMixin, metaclass=ABCMeta): """ Base Level Index to define the structure diff --git a/packages/data/src/pyearthtools/data/indexes/utilities/fileload.py b/packages/data/src/pyearthtools/data/indexes/utilities/fileload.py index 2ed6f9b9..eee5cf2f 100644 --- a/packages/data/src/pyearthtools/data/indexes/utilities/fileload.py +++ b/packages/data/src/pyearthtools/data/indexes/utilities/fileload.py @@ -83,7 +83,7 @@ def get_config(mf: bool = False): return xr.open_mfdataset( filter_files(location), decode_timedelta=True, # TODO: should we raise a warning? It seems to be required for almost all our data. - compat='override', + compat="override", **get_config(True), ) diff --git a/packages/data/src/pyearthtools/data/transforms/variables.py b/packages/data/src/pyearthtools/data/transforms/variables.py index 93d0117e..fec0253f 100644 --- a/packages/data/src/pyearthtools/data/transforms/variables.py +++ b/packages/data/src/pyearthtools/data/transforms/variables.py @@ -87,11 +87,20 @@ def apply(self, dataset: xr.Dataset) -> xr.Dataset: if self._variables is None: return dataset - var_included = set(dataset.data_vars).difference(set(self._variables)) + # 3/9/2025 - old logic was replaced with a simple drop of the variables + # A new issue will be raised to review how coordinate protection should + # work because people need a way to drop coords when needed. - if not var_included: - return dataset - return dataset[var_included] + # Calculate the difference between the data variables on the dataset + # and the variables requested for drop. This leaves coordinate variables + # unaffected + # var_included = set(dataset.data_vars).difference(set(self._variables)) + + # if not var_included: + # return dataset + # return dataset[var_included] + + return dataset.drop_vars(self._variables) class Select(Transform): diff --git a/packages/data/tests/indexes/test_indexes.py b/packages/data/tests/indexes/test_indexes.py index b8554205..14509f93 100644 --- a/packages/data/tests/indexes/test_indexes.py +++ b/packages/data/tests/indexes/test_indexes.py @@ -7,6 +7,7 @@ from pyearthtools.data.exceptions import DataNotFoundError + def test_Index(monkeypatch): monkeypatch.setattr("pyearthtools.data.indexes.Index.__abstractmethods__", set()) diff --git a/packages/nci_site_archive/tests/test_radar_proj.py b/packages/nci_site_archive/tests/test_radar_proj.py index 6a75f080..545323f2 100644 --- a/packages/nci_site_archive/tests/test_radar_proj.py +++ b/packages/nci_site_archive/tests/test_radar_proj.py @@ -24,13 +24,18 @@ import platform import xarray as xr -from site_archive_nci._Rainfields3 import ( - ErrorRadarProj, - ProjErrorStatus, - ProjKind, - RadarProj, - WarnRadarProj, -) +try: + from site_archive_nci._Rainfields3 import ( + ErrorRadarProj, + ProjErrorStatus, + ProjKind, + RadarProj, + WarnRadarProj, + ) + +except ImportError: + pytest.skip(allow_module_level=True) + PYPROJ_SAMPLE = pyproj.Proj("+proj=aea +lat_1=-36 +lat_2=-18 +lon_0=132 +units=m") EXPECTED_KEYS = [ diff --git a/packages/pipeline/src/pyearthtools/pipeline/operations/xarray/join.py b/packages/pipeline/src/pyearthtools/pipeline/operations/xarray/join.py index c4146509..b0e83742 100644 --- a/packages/pipeline/src/pyearthtools/pipeline/operations/xarray/join.py +++ b/packages/pipeline/src/pyearthtools/pipeline/operations/xarray/join.py @@ -46,15 +46,15 @@ def unjoin(self, sample: Any) -> tuple: class LatLonInterpolate(Joiner): - ''' + """ Makes additional assumptions about how interpolation should work and how the data is structured. In this case, interpolation is primarily - expected to occur according to latitude and longitude, performing + expected to occur according to latitude and longitude, performing no broadcasting, and iterating over other dimensions instead. It assumed the dimensions 'latitude', 'longitude', 'time', and 'level' will be present. 'lat' or 'lon' may also be used for convenience. - ''' + """ _override_interface = "Serial" @@ -78,15 +78,15 @@ def __init__( self._merge_kwargs = merge_kwargs def raise_if_dimensions_wrong(self, dataset): - ''' + """ Raise exceptions if the supplied dataset does not meet requirements - ''' + """ - if not hasattr(self, 'required_dims'): - if 'lat' in dataset.coords: - self.required_dims = ['lat', 'lon'] + if not hasattr(self, "required_dims"): + if "lat" in dataset.coords: + self.required_dims = ["lat", "lon"] else: - self.required_dims = ['latitude', 'longitude'] + self.required_dims = ["latitude", "longitude"] present_in_coords = [d in dataset.coords for d in self.required_dims] if not all(present_in_coords): @@ -100,12 +100,12 @@ def raise_if_dimensions_wrong(self, dataset): # raise ValueError(f"Cannot perform a GeoMergePancake on the data variables {data_var} without {self.required_dims} as a dimension") def maybe_interp(self, ds): - ''' + """ This method will only interpolate the datasets if the latitudes and longitudes don't already match. This means, for example, you can't use it to interpolate between time steps or vertical levels alone. The primary purpose here is lat/lon interpolation, not general model interpolation or arbitrarily-dimensioned data interpolation. - ''' + """ ds_coords_ok = [ds[coord].equals(self.reference_dataset[coord]) for coord in self.required_dims] @@ -115,7 +115,6 @@ def maybe_interp(self, ds): return ds - def _join_two_datasets(self, sample_a: xr.Dataset, sample_b: xr.Dataset) -> xr.Dataset: """ Used to reduce a sequence of joinable items. Only called by the public interface join method. @@ -144,7 +143,7 @@ def join(self, sample: tuple[Union[xr.Dataset, xr.DataArray], ...]) -> xr.Datase return merged def unjoin(self, sample: Any) -> tuple: - raise NotImplementedError("Not Implemented") + raise NotImplementedError("Not Implemented") class GeospatialTimeSeriesMerge(Joiner): diff --git a/packages/pipeline/src/pyearthtools/pipeline/operations/xarray/normalisation.py b/packages/pipeline/src/pyearthtools/pipeline/operations/xarray/normalisation.py index 3496dd79..68fe514f 100644 --- a/packages/pipeline/src/pyearthtools/pipeline/operations/xarray/normalisation.py +++ b/packages/pipeline/src/pyearthtools/pipeline/operations/xarray/normalisation.py @@ -157,10 +157,10 @@ class Deviation(xarrayNormalisation): """Deviation Normalisation""" def __init__( - self, - mean: FILE | xr.Dataset | xr.DataArray | float, + self, + mean: FILE | xr.Dataset | xr.DataArray | float, deviation: FILE | xr.Dataset | xr.DataArray | float, - debug=False + debug=False, ): """ Each argument take take a Dataset, DataArray, float or file object. @@ -173,7 +173,9 @@ def __init__( self.record_initialisation() if debug: - import pdb; pdb.set_trace() + import pdb + + pdb.set_trace() if isinstance(mean, xr.Dataset): self.mean = mean @@ -187,7 +189,7 @@ def __init__( if isinstance(deviation, xr.Dataset): self.deviation = deviation elif isinstance(deviation, xr.DataArray): - self.deviation = deviation + self.deviation = deviation elif isinstance(deviation, float): self.deviation = deviation else: diff --git a/packages/tutorial/src/pyearthtools/tutorial/HadisdDataClass.py b/packages/tutorial/src/pyearthtools/tutorial/HadisdDataClass.py index 1fba6cf9..0a5921ba 100644 --- a/packages/tutorial/src/pyearthtools/tutorial/HadisdDataClass.py +++ b/packages/tutorial/src/pyearthtools/tutorial/HadisdDataClass.py @@ -136,62 +136,32 @@ def __init__( self.record_initialisation() - # def get_all_station_ids(self, root_directory: Path | str) -> list[str]: - # """ - # Retrieve all station IDs by scanning the dataset directory. - - # Args: - # root_directory (Path | str): The root directory containing station data. - - # Returns: - # list[str]: A list of all station IDs. - # """ - # root_directory = Path(root_directory) - # station_ids = [] - # for folder in cached_iterdir(root_directory): - # if folder.is_dir(): - # for file in cached_iterdir(folder): - # if file.suffix == ".nc": # Check for NetCDF files - # # Extract the station ID from the filename - # station_id = file.stem.split("_")[-1] # Assuming station ID is the last part of the filename - # station_ids.append(station_id) - # return station_ids - def get_all_station_ids(self, root_directory: Path | str = None) -> list[str]: """ - Retrieve all station IDs by scanning the dataset directory. + Retrieve all station IDs by scanning the Zarr directory. Args: - root_directory (Path | str, optional): The root directory containing station data. - Defaults to HADISD_HOME/netcdf. + root_directory (Path | str, optional): The directory containing Zarr files. + Defaults to HADISD_HOME/zarr. Returns: list[str]: A list of all station IDs. """ - HADISD_HOME = self.ROOT_DIRECTORIES["hadisd"] if root_directory is None: - # Search all WMO folders for netcdf subfolders - wmo_folders = [f for f in Path(HADISD_HOME).iterdir() if f.is_dir() and f.name.startswith("WMO_")] - station_ids = [] - for wmo_folder in wmo_folders: - netcdf_dir = wmo_folder / "netcdf" - if cached_exists(netcdf_dir): - for file in cached_iterdir(netcdf_dir): - if file.suffix == ".nc": - station_id = file.stem.split("_")[-1] - station_ids.append(station_id) - return station_ids + zarr_dir = Path(HADISD_HOME) / "zarr" else: - root_directory = Path(root_directory) - if not cached_exists(root_directory): - raise DataNotFoundError(f"Root directory does not exist: {root_directory}") - station_ids = [] - for file in cached_iterdir(root_directory): - if file.suffix == ".nc": - station_id = file.stem.split("_")[-1] - station_ids.append(station_id) - return station_ids + zarr_dir = Path(root_directory) + + if not cached_exists(zarr_dir): + raise DataNotFoundError(f"Zarr directory does not exist: {zarr_dir}") + + station_ids = [] + for file in cached_iterdir(zarr_dir): + if file.suffix == ".zarr": + station_id = file.stem.split("_")[-1] + station_ids.append(station_id) + return station_ids def filesystem(self, *args, date_range=("1970-01-01T00", "2023-12-31T23"), **kwargs) -> dict[str, Path]: """ @@ -222,66 +192,18 @@ def filesystem(self, *args, date_range=("1970-01-01T00", "2023-12-31T23"), **kwa if not isinstance(station_ids, list) or not all(isinstance(sid, str) for sid in station_ids): raise TypeError(f"Expected station_ids to be a str or list[str], but got: {type(station_ids)}") - # Define the station ranges and corresponding folders - STATION_RANGES = [ - (0, 29999, "WMO_000000-029999"), - (30000, 49999, "WMO_030000-049999"), - (50000, 79999, "WMO_050000-079999"), - (80000, 99999, "WMO_080000-099999"), - (100000, 149999, "WMO_100000-149999"), - (150000, 199999, "WMO_150000-199999"), - (200000, 249999, "WMO_200000-249999"), - (250000, 299999, "WMO_250000-299999"), - (300000, 349999, "WMO_300000-349999"), - (350000, 399999, "WMO_350000-399999"), - (400000, 449999, "WMO_400000-449999"), - (450000, 499999, "WMO_450000-499999"), - (500000, 549999, "WMO_500000-549999"), - (550000, 599999, "WMO_550000-599999"), - (600000, 649999, "WMO_600000-649999"), - (650000, 699999, "WMO_650000-699999"), - (700000, 709999, "WMO_700000-709999"), - (710000, 714999, "WMO_710000-714999"), - (715000, 719999, "WMO_715000-719999"), - (720000, 721999, "WMO_720000-721999"), - (722000, 722999, "WMO_722000-722999"), - (723000, 723999, "WMO_723000-723999"), - (724000, 724999, "WMO_724000-724999"), - (725000, 725999, "WMO_725000-725999"), - (726000, 726999, "WMO_726000-726999"), - (727000, 729999, "WMO_727000-729999"), - (730000, 799999, "WMO_730000-799999"), - (800000, 849999, "WMO_800000-849999"), - (850000, 899999, "WMO_850000-899999"), - (900000, 949999, "WMO_900000-949999"), - (950000, 999999, "WMO_950000-999999"), - ] - # Map station IDs to their file paths paths = {} for station_id in station_ids: - wmo_number = station_id[:6] # Extract the first 6 digits of the station ID - station_numeric = int(wmo_number) # Convert the WMO number to an integer - - # Find the parent folder dynamically - parent_folder = None - for start, end, folder in STATION_RANGES: - if start <= station_numeric <= end: - parent_folder = folder - break - - if parent_folder is None: - raise ValueError(f"Station ID {station_id} does not fall within any defined range.") - - # Construct the expected filename - date_range = "19310101-20240101" # Hardcoded for now; adjust if dataset is updated + date_range_str = "19310101-20240101" # Hardcoded for now; adjust if dataset is updated version = "hadisd.3.4.0.2023f" - filename_nc = f"{version}_{date_range}_{station_id}.nc" - filename_zarr = f"{version}_{date_range}_{station_id}.zarr" + filename_zarr = f"{version}_{date_range_str}_{station_id}.zarr" + + # filename_nc = f"{version}_{date_range_str}_{station_id}.nc" # Uncomment to test with netcdf + # file_path_nc = Path(HADISD_HOME) / "netcdf" / filename_nc # Uncomment to test with netcdf # Construct the full path - _file_path_nc = Path(HADISD_HOME) / parent_folder / "netcdf" / filename_nc - file_path_zarr = Path(HADISD_HOME) / parent_folder / "zarr" / filename_zarr + file_path_zarr = Path(HADISD_HOME) / "zarr" / filename_zarr # Check if the file exists (comment out if testing with single netcdf) if not file_path_zarr.exists(): @@ -290,6 +212,7 @@ def filesystem(self, *args, date_range=("1970-01-01T00", "2023-12-31T23"), **kwa # Add the file path to the dictionary paths[station_id] = ( file_path_zarr # Change to file_path_zarr to test with zarr files or remove "_zarr" to test with netcdf files + # file_path_nc # Uncomment to test with netcdf files ) return paths diff --git a/packages/utils/src/pyearthtools/utils/data/converter.py b/packages/utils/src/pyearthtools/utils/data/converter.py index a8f2e0c2..e8bac483 100644 --- a/packages/utils/src/pyearthtools/utils/data/converter.py +++ b/packages/utils/src/pyearthtools/utils/data/converter.py @@ -158,9 +158,16 @@ def _distill_dataset(self, dataset: XR_OBJECT) -> dict[DISTILL_KEYS, Any]: try: dims[use_shape.index(size)] = coord except ValueError as e: - raise RuntimeError( - "Cannot record coordinate, currently converter can only handle datasets with variables of the same dimensions." - ) from e + + msg = ( + f"Cannot record coordinate '{coord}', currently the conversion can only handle data variables with " + f"the same dimensionality as the dataset coords {list(dataset.coords)}. " + f"Data variable {variables[0]} with dimensions of {dataset[variables[0]].dims} was used to estimate the shape required. " + f"You may need to drop unused coordinates, drop mismatching data variables, or broadcast your data variables onto the " + "coordinates of the dataset yourself as the proper approach is user-defined." + ) + + raise RuntimeError(msg) from e use_shape[use_shape.index(size)] = 1e10 while None in dims: