11_OPTIONAL_Basemap_with_Contextily.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 11. Adding Basemaps with Contextily\n",
    "\n",
    "If you work with geospatial data in Python, you most likely are familiar with the fantastic [GeoPandas](https://geopandas.org/) library.  GeoPandas leverages the power of [Matplotlib](https://matplotlib.org/) to enable users to make maps of their data.  However, until recently, it has not been easy to add basemaps to these maps. Basemaps are the contextual map data, like Google Maps, on top of which geospatial data are often displayed.\n",
    "\n",
    "\n",
    "The new Python library [contextily](https://github.com/geopandas/contextily), which stands for *context map tiles*, now makes it possible and relatively straight forward to add basemaps to Geopandas maps. Below we walk through a few common workflows for doing this.\n",
    "\n",
    "First, let's load are libraries. This assumes you have the following Python libraries installed in your environment:\n",
    "\n",
    "- pandas\n",
    "- matplotlib\n",
    "- geopandas (and all dependancies)\n",
    "- contextily\n",
    "- descartes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "import pandas as pd\n",
    "import geopandas as gpd\n",
    "import contextily as cx\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Read data into a Geopandas GeoDataFrame\n",
    "\n",
    "Fetch the census places data to map. Census places includes cities and other populated places. Here we fetch the 2019 cartographic boundary (`cb_`) file of California (`06`) places."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ca_places = \"https://www2.census.gov/geo/tiger/GENZ2019/shp/cb_2019_06_place_500k.zip\"\n",
    "places = gpd.read_file(ca_places)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Use the geodatarame `plot` method to make a quick map."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "places.plot();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we can see those cities, let's take a look at the data in the geodataframe."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "places.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can subset the data by selecting a row or rows by place name. Let's select the city of Berkeley, CA."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "berkeley = places[places['NAME']=='Berkeley']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "berkeley.plot();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Use Contextily to add a basemap\n",
    "\n",
    "Above we can see the map of the boundary of the city of Berkeley, CA.  The axis labels display the longitude and latitude coordinates for the bounding extent of the city.\n",
    "\n",
    "Let's use `contextily` in it's most simple form to add a basemap to provide the geographic context for Berkeley. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ax = berkeley.to_crs('EPSG:3857').plot(figsize=(9, 9))\n",
    "cx.add_basemap(ax)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are a few important things to note about the above code.\n",
    "\n",
    "- We use `matplotlib` to define the plot canvas as `ax`.\n",
    "- We then add the contextily basemap to the map with the code `cx.add_basemap(ax)`\n",
    "\n",
    "Additionally, we **dynamically transform the coordinate reference system**, or CRS, of the Berkeley geodataframe from geographic lat/lon coordinates to `web mercator` using the method **to_crs('EPSG:3857')**. [Web mercator](https://en.wikipedia.org/wiki/Web_Mercator_projection) is the default CRS used by all web map tilesets. It is referenced by a the code `EPSG:3857` where [EPSG](https://en.wikipedia.org/wiki/EPSG_Geodetic_Parameter_Dataset) stands for the the initials of the organization that created these codes  (the European Petroleum Survey Group).\n",
    "\n",
    "Let's clean up the map by adding some code to change the symbology of the Berkeley city boundary. This will highlight the value of adding a basemap.\n",
    "\n",
    "First, let's map the boundary with out a fill color."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "berkeley.plot(edgecolor=\"red\", facecolor=\"none\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, let's build on those symbology options and add the contextily basemap."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ax = berkeley.to_crs('EPSG:3857').plot(edgecolor=\"red\", \n",
    "                                       facecolor=\"none\", # or a color \n",
    "                                       alpha=0.95,       # opacity value for colors, 0-1\n",
    "                                       linewidth=2,      # line, or stroke, thickness\n",
    "                                       figsize=(9, 9)\n",
    "                                      )\n",
    "cx.add_basemap(ax)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Mapping Point Data\n",
    "\n",
    "Let's expand on this example by mapping a point dataset of BART station locations.\n",
    "\n",
    "First we fetch these data from a D-Lab web mapping tutorial."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bart_url = 'https://raw.githubusercontent.com/dlab-berkeley/Geospatial-Fundamentals-in-Python/master/notebook_data/transportation/bart.csv'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bart = pd.read_csv(bart_url)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bart.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Converting Point Data in a dataframe to Geospatial Data in a geodataframe\n",
    "\n",
    "Because these data are in a CSV file we read them into a Pandas DataFrame.\n",
    "\n",
    "In order to map these data we need to convert these data to a GeoPandas GeoDataFame. To do this, we need to specify:\n",
    "\n",
    "- the data, here the geodataframe `bart`\n",
    "- the coordinate data, here `bart['X']` and `bart['Y']`\n",
    "- the CRS of the bart coordinate data, here `EPSG:4326`\n",
    "\n",
    "The CRS code 'EPSG:4326' stands for the World Geodectic System of 1984, or WGS84. This is the most commonly used CRS for geographic (lat/lon) coordinate data.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Convert the DataFrame to a GeoDataFrame. \n",
    "bart_gdf = gpd.GeoDataFrame(bart, geometry=gpd.points_from_xy(bart['lon'], \n",
    "                                                              bart['lat']), \n",
    "                            crs='EPSG:4326') \n",
    "\n",
    "# and take a look\n",
    "bart_gdf.plot();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we have the BART data in a geodataframe we can use the same commands as we did above to map it with a contextily basemap."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ax = bart_gdf.to_crs('EPSG:3857').plot(figsize=(9, 9))\n",
    "cx.add_basemap(ax)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We have the full range of `matplotlib` style options to enhance the map, a few of which are shown in the example below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ax = bart_gdf.to_crs('EPSG:3857').plot(\n",
    "                                    color=\"red\",\n",
    "                                    edgecolor=\"black\",\n",
    "                                    markersize=50, \n",
    "                                    figsize=(9, 9))\n",
    "\n",
    "ax.set_title('Bay Area Bart Stations')\n",
    "cx.add_basemap(ax)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Changing the Basemap\n",
    "\n",
    "By default `contextiley` returns maptiles from the OpenStreetmap Mapnik basemap. However, ther are other available tilesets from different providers.  These tilesets are stored in the contextily `cx.providers` dictionary.\n",
    "\n",
    "That's a large dictionary and you can view it. Alternatively, and more simply, you can access the list of the providers in this dictionary using the command `cs.providers.keys`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# change basemap - can be one of these\n",
    "# first see available provider names\n",
    "cx.providers.keys()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "Once you have the list of providers, you can find the names of their specific tilesets. \n",
    "\n",
    "Below, we retrieve the list of the tilesets available from the provider `CartoDB`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Then find the names of the tile sets for a specific provider\n",
    "cx.providers.CartoDB.keys()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can specify a different tileset using the **source** argument to the `add_basemap` method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cx.providers.Esri.keys()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Change the basemap provider and tileset\n",
    "ax = bart_gdf.to_crs('EPSG:3857').plot(figsize=(9, 9))\n",
    "cx.add_basemap(ax,  source=cx.providers.NASAGIBS.ModisTerraTrueColorCR)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Learning More\n",
    "\n",
    "Above, we prove a very short introduction to the excellent `contextily` library. You can find more detailed information on the `contextily` homepage, available at: [https://github.com/geopandas/contextily](https://github.com/geopandas/contextily). We especially encourage you to check out the notebook examples provided in that github repo.\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "geo_env2",
   "language": "python",
   "name": "geo_env2"
  },
  "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.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}