{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import holoviews as hv\n",
    "import geoviews as gv\n",
    "import geoviews.feature as gf\n",
    "import cartopy\n",
    "import cartopy.feature as cf\n",
    "\n",
    "from geoviews import opts\n",
    "from cartopy import crs as ccrs\n",
    "\n",
    "gv.extension('matplotlib', 'bokeh')\n",
    "\n",
    "gv.output(dpi=120, fig='svg')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Cartopy and shapely make working with geometries and shapes very simple, and GeoViews provides convenient wrappers for the various geometry types they provide. In addition to Path and Polygons types, which draw geometries from lists of arrays or a geopandas DataFrame, GeoViews also provides the ``Feature`` and ``Shape`` types, which wrap cartopy Features and shapely geometries respectively.\n",
    "\n",
    "### Feature\n",
    "\n",
    "The Feature Element provides a very convenient means of overlaying a set of basic geographic features on top of or behind a plot. The ``cartopy.feature`` module provides various ways of loading custom features, however geoviews provides a number of default features which we have imported as ``gf``, amongst others this includes coastlines, country borders, and land masses. Here we demonstrate how we can plot these very easily, either in isolation or overlaid:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "(gf.ocean + gf.land + gf.ocean * gf.land * gf.coastline * gf.borders).cols(3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These default features simply wrap around cartopy Features, therefore we can easily load a custom ``NaturalEarthFeature`` such as graticules at 30 degree intervals: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "graticules = cf.NaturalEarthFeature(\n",
    "    category='physical',\n",
    "    name='graticules_30',\n",
    "    scale='110m')\n",
    "\n",
    "(gf.ocean() * gf.land() * gv.Feature(graticules, group='Lines') * gf.borders * gf.coastline).opts(\n",
    "    opts.Feature('Lines', projection=ccrs.Robinson(), facecolor='none', edgecolor='gray'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The scale of features may be controlled using the ``scale`` plot option, the most common options being `'10m'`, `'50m'` and `'110m'`. Cartopy will downloaded the requested resolution as needed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gv.output(backend='bokeh')\n",
    "\n",
    "(gf.ocean * gf.land.options(scale='110m', global_extent=True) * gv.Feature(graticules, group='Lines') + \n",
    " gf.ocean * gf.land.options(scale='50m', global_extent=True) * gv.Feature(graticules, group='Lines'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Zoom in using the bokeh zoom widget and you should see that the right hand panel is using a higher resolution dataset for the land feature.\n",
    "\n",
    "Instead of displaying a ``Feature`` directly it is also possible to request the geometries inside a ``Feature`` using the ``Feature.geoms`` method, which also allows specifying a ``scale`` and a ``bounds`` to select a subregion:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gf.land.geoms('50m', bounds=(-10, 40, 10, 60))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When working interactively with higher resolution datasets it is sometimes necessary to dynamically update the geometries based on the current viewport. The ``resample_geometry`` operation is an efficient way to display only polygons that intersect with the current viewport and downsample polygons on-the-fly."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gv.operation.resample_geometry(gf.coastline.geoms('10m')).opts(width=400, height=400, color='black')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Try zooming into the plot above and you will see the coastline geometry resolve to a higher resolution dynamically (this requires a live Python kernel)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Shape\n",
    "\n",
    "The ``gv.Shape`` object wraps around any shapely geometry, allowing finer grained control over each polygon. We can, for example, access the geometries on the ``LAND`` feature and display them individually. Here we will get the geometry corresponding to the Australian continent and display it using shapely's inbuilt SVG repr (not yet a HoloViews plot, just a bare SVG displayed by Jupyter directly):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "land_geoms = gf.land.geoms(as_element=False)\n",
    "land_geoms[21]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Instead of letting shapely render it as an SVG, we can now wrap it in the ``gv.Shape`` object and let matplotlib or bokeh render it, alone or with other GeoViews or HoloViews objects:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "australia = gv.Shape(land_geoms[21])\n",
    "alice_springs = gv.Text(133.870,-21.5, 'Alice Springs')\n",
    "\n",
    "australia * gv.Points([(133.870,-23.700)]).opts(color='black', width=400) * alice_springs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also supply a list of geometries directly to a Polygons or Path element:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gv.Polygons(land_geoms) + gv.Path(land_geoms)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This makes it possible to create choropleth maps, where each part of the geometry is assigned a value that will be used to color it. However, constructing a choropleth by combining a bunch of shapes can be a lot of effort and is error prone. For that reason, the Shape Element provides convenience methods to load geometries from a shapefile. Here we load the boundaries of UK electoral districts directly from an existing shapefile:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hv.output(backend='matplotlib')\n",
    "\n",
    "shapefile = '../assets/boundaries/boundaries.shp'\n",
    "gv.Shape.from_shapefile(shapefile, crs=ccrs.PlateCarree())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To combine these shapes with some actual data, we have to be able to merge them with a dataset. To do so we can inspect the records the cartopy shapereader loads:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "shapes = cartopy.io.shapereader.Reader(shapefile)\n",
    "list(shapes.records())[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As we can see, the record contains a ``MultiPolygon`` together with a standard geographic ``code``, which we can use to match up the geometries with a dataset. To continue we will require a dataset that is also indexed by these codes. For this purpose we load a dataset of the 2016 EU Referendum result in the UK:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "referendum = pd.read_csv('../assets/referendum.csv')\n",
    "referendum = hv.Dataset(referendum)\n",
    "referendum.data.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The ``from_records`` function optionally also supports merging the records and dataset directly. To merge them, supply the name of the shared attribute on which the merge is based via the ``on`` argument. If the name of attribute in the records and the dimension in the dataset match exactly, you can simply supply it as a string, otherwise supply a dictionary mapping between the attribute and column name. In this case we want to color the choropleth by the `'leaveVoteshare'`, which we define via the `value` argument.\n",
    "\n",
    "Additionally we can request one or more indexes using the ``index`` argument. Finally we will declare the coordinate reference system in which this data is stored, which will in most cases be the simple Plate Carree projection.  We can then view the choropleth, with each shape colored by the specified value (the percentage who voted to leave the EU):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hv.output(backend='bokeh')\n",
    "\n",
    "gv.Shape.from_records(shapes.records(), referendum, on='code', value='leaveVoteshare',\n",
    "                      index=['name', 'regionName']).opts(tools=['hover'], width=350, height=500)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## GeoPandas"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "GeoPandas extends the datatypes used by pandas to allow spatial operations on geometric types, which makes it a very convenient way of working with geometries with associated variables. GeoViews ``Path``, ``Contours`` and ``Polygons`` Elements natively support projecting and plotting of\n",
    "geopandas DataFrames using both ``matplotlib`` and ``bokeh`` plotting extensions. We will load an example dataset of Airbnb rentals in Chicago, which also includes some additional data about the city's communities:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import geodatasets as gds\n",
    "import geopandas as gpd\n",
    "\n",
    "data = gpd.read_file(gds.get_path('geoda airbnb'))\n",
    "data[[\"community\", \"population\", \"num_spots\", \"geometry\"]].head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can simply pass the GeoPandas DataFrame to a Polygons, Path or Contours element and it will plot the data for us. The ``Contours`` and ``Polygons`` will automatically color the data by the first specified value dimension defined by the ``vdims`` keyword (the geometries may be colored by any dimension using the ``color`` plot option):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "poly_plot = gv.Polygons(data, vdims=[\"population\", \"community\", \"num_spots\"]).opts(width=600, height=600)\n",
    "gv.tile_sources.OSM * poly_plot"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we will switch the color by number of spots (``num_spots``) and activating the hover tool to reveal information about the plot. The switch will work in both Matplotlib and Nokeh, but the bokeh version will be more interactive:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gv.tile_sources.OSM * poly_plot.opts(color=\"num_spots\", cmap='tab20', tools=['hover'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The \"Working with Bokeh\" GeoViews notebook shows how to enable hover data that displays information about each of these shapes interactively."
   ]
  }
 ],
 "metadata": {
  "language_info": {
   "name": "python",
   "pygments_lexer": "ipython3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}