{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import xarray as xr\n",
    "import hvplot.pandas  # noqa\n",
    "import hvplot.xarray  # noqa\n",
    "import cartopy.crs as ccrs\n",
    "\n",
    "from bokeh.sampledata.airport_routes import airports"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Installation\n",
    "\n",
    "The plot API also has support for geographic data built on top of Cartopy and GeoViews. Both can be installed using conda with:\n",
    "\n",
    "    conda install -c pyviz geoviews\n",
    "    \n",
    "or if the cartopy dependency has been satisfied in some other way, GeoViews may also be installed using pip:\n",
    "\n",
    "    pip install geoviews\n",
    "\n",
    "## Usage\n",
    "\n",
    "Only certain hvPlot types support geographic coordinates, currently including: 'points', 'polygons', 'paths', 'image', 'quadmesh', 'contour', and 'contourf'. As an initial example, consider a dataframe of all US airports (including military bases overseas):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "airports.head(3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plotting points\n",
    "\n",
    "If we want to overlay our data on geographic maps or reproject it into a geographic plot, we can set ``geo=True``, which declares that the data will be plotted in a geographic coordinate system.  The default coordinate system is the ``PlateCarree`` projection, i.e., raw longitudes and latitudes. If the data is in another coordinate system, you will need to [declare an explicit ``crs``](#Declaring-a-CRS) as an argument, in which case `geo=True` is assumed.  Once hvPlot knows that your data is in geo coordinates, you can use the ``tiles`` option to overlay a the plot on top of map tiles."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "airports.hvplot.points('Longitude', 'Latitude', geo=True, color='red', alpha=0.2,\n",
    "                       xlim=(-180, -30), ylim=(0, 72), tiles='ESRI')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Declaring a CRS\n",
    "\n",
    "To declare a geographic plot we have to supply a ``cartopy.crs.CRS`` (or coordinate reference system).  Coordinate reference systems are described in the [GeoViews documentation](http://geoviews.org/user_guide/Projections.html) and the full list of available CRSs is in the [cartopy documentation](https://scitools.org.uk/cartopy/docs/v0.15/crs/projections.html). "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Geopandas\n",
    "\n",
    "Since a GeoPandas ``DataFrame`` is just a Pandas DataFrames with additional geographic information, it inherits the ``.hvplot`` method. We can thus easily load shapefiles and plot them on a map:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import geopandas as gpd\n",
    "\n",
    "cities = gpd.read_file(gpd.datasets.get_path('naturalearth_cities'))\n",
    "\n",
    "cities.hvplot(global_extent=True, frame_height=450, tiles=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The GeoPandas support allows plotting ``GeoDataFrames`` containing ``'Point'``, ``'Polygon'``, ``'LineString'`` and ``'LineRing'`` geometries, but not ones containing a mixture of different geometry types. Calling ``.hvplot`` will automatically figure out the geometry type to plot, but it also possible to call ``.hvplot.points``, ``.hvplot.polygons``, and ``.hvplot.paths`` explicitly.\n",
    "\n",
    "To draw multiple GeoDataFrames onto the same plot, use the ``*`` operator:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "world = gpd.read_file(gpd.datasets.get_path('naturalearth_lowres'))\n",
    "\n",
    "world.hvplot(geo=True) * cities.hvplot(geo=True, color='orange')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is possible to declare a specific column to use as color with the ``c`` keyword:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "world.hvplot(geo=True) + world.hvplot(c='continent', geo=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Spatialpandas\n",
    "\n",
    "Spatialpandas is another powerful library for working with geometries and is optimized for rendering with datashader, making it possible to plot millions of individual geometries very quickly:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import spatialpandas as spd\n",
    "\n",
    "spd_world = spd.GeoDataFrame(world)\n",
    "\n",
    "spd_world.hvplot(datashade=True, project=True, aggregator='count_cat', c='continent', color_key='Category10')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Declaring an output projection"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The ``crs=`` argument specifies the *input* projection, i.e. it declares how to interpret the incoming data values. You can independently choose any *output* projection, i.e. how you want to map the data points onto the screen for display, using the ``projection=`` argument. After loading the same temperature dataset explored in the [Gridded Data](Gridded_Data.ipynb) section, the data can be displayed on an Orthographic projection:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "air_ds = xr.tutorial.open_dataset('air_temperature').load()\n",
    "\n",
    "air_ds.hvplot.quadmesh(\n",
    "    'lon', 'lat', 'air', projection=ccrs.Orthographic(-90, 30),\n",
    "    global_extent=True, frame_height=540, cmap='viridis',\n",
    "    coastline=True\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you don't need to pass any keyword arguments to a given projection and you don't have cartopy.crs (ccrs) imported, you can use the string representation: e.g. ``'LambertConformal'`` instead of ``ccrs.LambertConformal()``. Note that it is case sensitive!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "air_ds.hvplot.quadmesh(\n",
    "    'lon', 'lat', 'air', projection='LambertConformal',\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that when displaying raster data in a projection other than the one in which the data is stored, it is more accurate to render it as a ``quadmesh`` rather than an ``image``. As you can see above, a QuadMesh will project each original bin or pixel into the correct non-rectangular shape determined by the projection, accurately showing the geographic extent covered by each sample. An Image, on the other hand, will always be rectangularly aligned in the 2D plane, which requires warping and resampling the data in a way that allows efficient display but loses accuracy at the pixel level. Unfortunately, rendering a large QuadMesh using Bokeh can be very slow, but there are two useful alternatives for datasets too large to be practical as native QuadMeshes.\n",
    "\n",
    "The first is using the ``rasterize`` or ``datashade`` options to regrid the data before rendering it, i.e., rendering the data on the backend and then sending a more efficient image-based representation to the browser. One thing to note when using these operations is that it may be necessary to project the data **before** rasterizing it, e.g. to address wrapping issues. To do this provide ``project=True``, which will project the data before it is rasterized (this also works for other types and even when not using these operations). Another reason why this is important when rasterizing the data is that if the CRS of the data does not match the displayed projection, all the data will be projected every time you zoom or pan, which can be very slow. Deciding whether to ``project`` is therefore a tradeoff between projecting the raw data ahead of time or accepting the overhead on dynamic zoom and pan actions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rasm = xr.tutorial.open_dataset('rasm').load()\n",
    "\n",
    "\n",
    "\n",
    "rasm.hvplot.quadmesh(\n",
    "    'xc', 'yc', crs=ccrs.PlateCarree(), projection=ccrs.PlateCarree(),\n",
    "    ylim=(0, 90), cmap='viridis', project=True, geo=True,\n",
    "    rasterize=True, coastline=True, frame_width=800, dynamic=False,\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Another option that's still relatively slow for larger data but avoids sending large data into your browser is to plot the data using ``contour`` and ``contourf`` visualizations, generating a line or filled contour with a discrete number of levels:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rasm.hvplot.contourf(\n",
    "    'xc', 'yc', crs=ccrs.PlateCarree(), projection=ccrs.PlateCarree(),\n",
    "    ylim=(0, 90), frame_width=800, cmap='viridis', levels=10,\n",
    "    coastline=True\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As you can see, hvPlot makes it simple to work with geographic data visually.  For more complex plot types and additional details, see the [GeoViews](http://geoviews.org) documentation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Geographic options\n",
    "\n",
    "The API provides various geo-specific options:\n",
    "\n",
    "- ``coastline`` (default=False): Whether to display a coastline on top of the plot, setting ``coastline='10m'/'50m'/'110m'`` specifies a specific scale\n",
    "- ``crs`` (default=None): Coordinate reference system of the data specified as Cartopy CRS object, proj.4 string or EPSG code\n",
    "- ``features`` features (default=None): A list of features or a dictionary of features and the scale at which to render it. Available features include 'borders', 'coastline', 'lakes', 'land', 'ocean', 'rivers' and 'states'. Available scales include '10m'/'50m'/'110m'.\n",
    "- ``geo`` (default=False): Whether the plot should be treated as geographic (and assume PlateCarree, i.e. lat/lon coordinates)\n",
    "- ``global_extent`` (default=False): Whether to expand the plot extent to span the whole globe\n",
    "- ``project`` (default=False): Whether to project the data before plotting (adds initial overhead but avoids projecting data when plot is dynamically updated)\n",
    "- ``tiles`` (default=False): Whether to overlay the plot on a tile source. Tiles sources can be selected by name, the default is 'Wikipedia'.\n",
    "Other options are: 'CartoDark', 'CartoEco', 'CartoLight', 'CartoMidnight', 'EsriImagery', 'EsriNatGeo', 'EsriReference''EsriTerrain', 'EsriUSATopo', 'OSM', 'StamenLabels', 'StamenTerrain', 'StamenTerrainRetina', 'StamenToner', 'StamenTonerBackground', 'StamenWatercolor'"
   ]
  }
 ],
 "metadata": {
  "language_info": {
   "name": "python",
   "pygments_lexer": "ipython3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}