{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "3d8b50c7",
   "metadata": {},
   "source": [
    "# Searching for Stations\n",
    "## Where in the world are these stations?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c1ff413c",
   "metadata": {},
   "source": [
    "One common frustration is finding relevant stations for analysis and validation\n",
    "\n",
    "**What are some methods you have used to identify stations that fit your research?**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9231cb26",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Imports\n",
    "from datetime import datetime\n",
    "import pandas as pd\n",
    "import geopandas as gpd\n",
    "from pathlib import Path\n",
    "\n",
    "from metloom.pointdata import SnotelPointData, CDECPointData, MesowestPointData"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fb184ca5",
   "metadata": {},
   "outputs": [],
   "source": [
    "try:\n",
    "    import folium\n",
    "except Exception as e:\n",
    "    !pip install folium\n",
    "    \n",
    "try:\n",
    "    import mapclassify\n",
    "except Exception as e:\n",
    "    !pip install mapclassify"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ffc9a610",
   "metadata": {},
   "source": [
    "## Search for points in an area"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "69501ab6",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Find your area\n",
    "sf_path = Path(\"./data/outline.shp\").expanduser()\n",
    "sf = gpd.read_file(str(sf_path))\n",
    "sf[\"name\"] = [\"Our outline\"]\n",
    "variables = [SnotelPointData.ALLOWED_VARIABLES.SNOWDEPTH]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "91fbe619",
   "metadata": {},
   "outputs": [],
   "source": [
    "# What does the area look like\n",
    "sf.explore()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0493e54e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Find all the points in the area for our variables\n",
    "points = SnotelPointData.points_from_geometry(sf, variables)\n",
    "print(len(points))\n",
    "# This is an iterator\n",
    "print(type(points))\n",
    "# It contains the points in a `points` attribute\n",
    "print(points.points)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "65e808fc",
   "metadata": {},
   "source": [
    "### What if I want a station nearby?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "17b66820",
   "metadata": {},
   "outputs": [],
   "source": [
    "# buffer - add buffer (in degrees) to extents\n",
    "print(SnotelPointData.points_from_geometry(sf, variables, buffer=0.5).points)\n",
    "\n",
    "# within_geometry - doesn't do anything since our geometry is already a square\n",
    "print(SnotelPointData.points_from_geometry(sf, variables, within_geometry=False).points)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "af4d1bd0-5024-4415-875e-ac0d3ea30eab",
   "metadata": {},
   "outputs": [],
   "source": [
    "# imports for a static map\n",
    "try:\n",
    "    import contextily as ctx\n",
    "except Exception:\n",
    "    !pip install contextily\n",
    "    import contextily as ctx\n",
    "import matplotlib.pyplot as plt\n",
    "from shapely.geometry import Point"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8644756d-d3b9-4a23-9af2-a5bc0e1dd8c8",
   "metadata": {
    "tags": [
     "nbsphinx-thumbnail",
     "nbsphinx-gallery"
    ]
   },
   "outputs": [],
   "source": [
    "# Make a static plot of the points\n",
    "\n",
    "# turn that iterator into a dataframe\n",
    "df = points.to_dataframe()\n",
    "# look at what is in the dataframe\n",
    "df.head(10)\n",
    "\n",
    "def remove_z(geom):\n",
    "    if geom.is_empty:\n",
    "        return geom\n",
    "    else:\n",
    "        return Point(geom.x, geom.y)\n",
    "\n",
    "df.geometry = df.geometry.apply(remove_z)\n",
    "df = df.set_crs(\"EPSG:4326\")\n",
    "\n",
    "\n",
    "print(df)\n",
    "ax = sf.to_crs(\"EPSG:4326\").plot(alpha=0.1)\n",
    "ax = df.plot(ax=ax, color=\"red\")\n",
    "ctx.add_basemap(\n",
    "    ax,\n",
    "    # source=ctx.providers.OpenStreetMap.Mapnik,\n",
    "    source=ctx.providers.USGS.USTopo,\n",
    "    crs=df.crs,\n",
    "    # zoom=5\n",
    ")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d8495b3d",
   "metadata": {},
   "source": [
    "## Our plot could be more interesting\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a3215268",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "# plot the shapefile\n",
    "m = sf.explore(\n",
    "    tooltip=False, color=\"grey\", highlight=False, style_kwds={\"opacity\": 0.2}, popup=[\"name\"]\n",
    ")\n",
    "# plot the points on top of the shapefile\n",
    "df.explore(m=m, tooltip=[\"name\", \"id\"], color=\"red\", marker_kwds={\"radius\":4})"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9067a52d",
   "metadata": {},
   "source": [
    "## Extra\n",
    "**How would we get the daily data for the first point from points?**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dce90f1d",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.18"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}