{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Contextily 1: Intro"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Overview:\n",
    "1. Read in the NYS Mesonet sites table with GeoPandas\n",
    "1. Plot the NYS Mesonet sites on a map\n",
    "1. Use the [contextily](https://contextily.readthedocs.io/en/latest/) package to overlay an image as a basemap, served via a remote tile server"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Prerequisites\n",
    "\n",
    "| Concepts | Importance | Notes |\n",
    "| --- | --- | --- |\n",
    "| (Geo)Pandas| Necessary | |\n",
    "| Cartopy| Necessary | |\n",
    "\n",
    "* **Time to learn**: 10 minutes\n",
    "***"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Imports"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import cartopy.crs as ccrs\n",
    "import cartopy.feature as cfeature\n",
    "import geopandas as gp\n",
    "import pandas as pd\n",
    "from metpy.plots import USCOUNTIES\n",
    "import contextily as ctx\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Read in the NYS Mesonet sites table with Geopandas"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, open the NYSM sites table using Pandas, as we've done before."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nysm_sites = pd.read_csv('/spare11/atm533/data/nysm_sites.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nysm_sites.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the Shapefiles Intro notebook, we used Geopandas to read in a shapefile as a Dataframe. Geopandas is not limited to shapefiles; it can also georeference existing Pandas dataframes that contain relevant columns or rows. The NYSM dataframe has latitude and longitude columns; we will pass these in to the `GeoDataFrame` method which creates the geo-referenced dataframe."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-info\">\n",
    "    <b>Note:</b> A <code>GeoDataFrame</code> requires a column named <i>geometry</i>. Besides polygons and multipolygons, another valid <i>geometry</i> is <b>Points</b>. We thus create a <i>Geometry</i> series using Geopandas' <code>points_from_xy</code> method. This method accepts two series... in this case, the longitude and latitude columns from the NYSM sites dataframe.\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gdf = gp.GeoDataFrame(nysm_sites,geometry=gp.points_from_xy(nysm_sites.lon,nysm_sites.lat))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Examine the `GeoDataFrame` and its *geometry* column"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gdf"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As we did in the first Shapefiles notebook, let's take a look at the geometry attribute of one of the sites."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gdf.loc[gdf.stid == 'VOOR'].geometry"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It's a `POINT` ... whose x- and y- coordinates are the corresponding longitude and latitude in degrees."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Also as in the previous notebook, we can use Geopandas `plot` method to get a quick visualization."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gdf.plot(figsize=(12,9));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Make a plot using Matplotlib and Cartopy. \n",
    "To plot the points corresponding to NYS Mesonet sites, we use the same Geopandas `plot` method, but need to specify the axes. In this case, we first define a set of `GeoAxes`, as we always do when using Cartopy, and then specify the name of that `GeoAxes` object as the value for the `ax` argument."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(15,10))\n",
    "ax = plt.axes(projection=ccrs.PlateCarree())\n",
    "\n",
    "res='10m'\n",
    "plt.title('NYS Mesonet Sites')\n",
    "gdf.plot(ax=ax,color='tab:red')\n",
    "\n",
    "ax.coastlines(resolution=res)\n",
    "ax.add_feature (cfeature.LAKES.with_scale(res), alpha = 0.5)\n",
    "ax.add_feature (cfeature.STATES.with_scale(res))\n",
    "ax.add_feature(USCOUNTIES,edgecolor='grey');\n",
    "ax.set_extent([-80, -72,40.5,45.2], ccrs.PlateCarree())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Add a background GeoTIFF raster image. \n",
    "We'll use the *contextily* package's `add_basemap` method for that. It will request GeoTIFFs from a remote server that not only provides a source of background maps, but also makes them available at various zoom levels. Depending on how we define the extent of our map, one or more *tiles* will be downloaded from the remote image server. The default basemap comes from [Openstreetmap](https://www.openstreetmap.org/)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-info\">\n",
    "    <b>Note:</b> The default coordinate reference system (<i>crs</i>) in Contextily is <b>WebMercator</b> (<a href=\"https://epsg.io/3857\">EPSG 3857</a>). We can either transform the remote images to our desired projection, or set up our map using WebMercator.\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We specify a zoom level here, although if you do not specify one, Contextily will determine an appropriate default value. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-info\">\n",
    "    <b>Tip:</b> The larger the zoom value, the greater detail in the background image (<i>but larger values will take <b>longer</b> to download and render!</i>).\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(15,10))\n",
    "ax = plt.axes(projection=ccrs.epsg('3857'))\n",
    "\n",
    "res='10m'\n",
    "plt.title('NYS Mesonet Sites')\n",
    "gdf.plot(ax=ax,color='tab:red',transform=ccrs.PlateCarree())\n",
    "ax.coastlines(resolution=res)\n",
    "ax.add_feature (cfeature.LAKES.with_scale(res), alpha = 0.5)\n",
    "ax.add_feature (cfeature.STATES.with_scale(res))\n",
    "ax.add_feature(USCOUNTIES,edgecolor='grey', linewidth=1 );\n",
    "ax.set_extent([-80, -71.4,40.5,45.2], ccrs.PlateCarree())\n",
    "ctx.add_basemap(ax,zoom=9,transform=ccrs.epsg('3857'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Try a lower zoom value."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(15,10))\n",
    "ax = plt.axes(projection=ccrs.epsg('3857'))\n",
    "\n",
    "res='10m'\n",
    "plt.title('NYS Mesonet Sites')\n",
    "gdf.plot(ax=ax,color='tab:red',transform=ccrs.PlateCarree())\n",
    "ax.coastlines(resolution=res)\n",
    "ax.add_feature (cfeature.LAKES.with_scale(res), alpha = 0.5)\n",
    "ax.add_feature (cfeature.STATES.with_scale(res))\n",
    "ax.add_feature(USCOUNTIES,edgecolor='grey', linewidth=1 );\n",
    "\n",
    "ax.set_extent([-80, -71.4,40.5,45.2], ccrs.PlateCarree())\n",
    "ctx.add_basemap(ax,zoom=7,transform=ccrs.epsg('3857'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plot a map whose extent covers just the Greater Capital District.\n",
    "Do not specify a zoom level. Let contextily's `add_ basemap` figure determine the optimum one."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(15,10))\n",
    "ax = plt.axes(projection=ccrs.epsg('3857'))\n",
    "\n",
    "res='10m'\n",
    "plt.title('NYS Mesonet Sites')\n",
    "gdf.plot(ax=ax,color='tab:red',transform=ccrs.PlateCarree())\n",
    "ax.coastlines(resolution=res)\n",
    "\n",
    "ax.add_feature (cfeature.LAKES.with_scale(res), alpha = 0.5)\n",
    "ax.add_feature (cfeature.STATES.with_scale(res))\n",
    "ax.add_feature(USCOUNTIES,edgecolor='grey', linewidth=1 );\n",
    "\n",
    "ax.set_extent([-74.4, -73.4,42.3,43.1], ccrs.PlateCarree())\n",
    "ctx.add_basemap(ax)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-warning\">\n",
    "    <b>Question:</b> How might you add the NYSM site ID's to the plot?\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "## Summary\n",
    "* Besides shapefiles, Geopandas can geo-reference Pandas dataframes.\n",
    "* The contextily package allows for the plotting of remotely-served basemaps.\n",
    "* These basemap servers support tiling; thus, zooming in or out will reveal greater or less features.\n",
    "\n",
    "### What's Next?\n",
    "We'll explore additional basemaps besides Openstreetmaps ... such as high-resolution MODIS satellite composites from NASA."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Reference:\n",
    "1. [Contextily documentation](https://contextily.readthedocs.io/en/latest/)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 August 2023 Environment",
   "language": "python",
   "name": "aug23"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}