1.GeoPandas.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Hello \"GeoPandas\" World!\n",
    "\n",
    "Let us start by importing required modules"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import geopandas as gpd\n",
    "import matplotlib.pyplot as plt\n",
    "import pandas as pd\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Open shapefile\n",
    "plot_locations = gpd.read_file('data/sites.shp')\n",
    "type(plot_locations)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Like PANDAS (DataFrame), the plot_locations variable is a **\"GeoDataFrame\"** and comes with several in-built methods / functions. \n",
    "\n",
    "Let us print the first (head) and last (tail) 6 points data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#change head to tail and see the difference\n",
    "print(plot_locations.head(6))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(plot_locations.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#to see the bounds\n",
    "print(plot_locations.total_bounds)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#to view CRS of object\n",
    "print(plot_locations.crs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#to view the geometry type\n",
    "print(plot_locations.geom_type)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Plot a Shapefile**\n",
    "\n",
    "Next, you can visualize the data in your Python geodata.frame object using the .plot() method. \n",
    "\n",
    "Notice that you can create a plot using the geopandas base plotting using the syntax:\n",
    "\n",
    "```python\n",
    "dataframe_name.plot()\n",
    "```\n",
    "\n",
    "You can call .plot() without setting up a figure and axis object like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_locations.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In general it is good practice to setup an axis object so you can plot different layers together. \n",
    "\n",
    "When you do that you need to provide the plot function with the axis object that you want it to plot on. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax1 = plt.subplots(figsize=(10, 10))\n",
    "\n",
    "plot_locations.plot(ax=ax1)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Categorical Symbology and more** ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(10, 10))\n",
    "\n",
    "# Plot the data and add a legend\n",
    "plot_locations.plot(    column='COVER',\n",
    "                        categorical=True,\n",
    "                        legend=True,\n",
    "                        markersize=50,\n",
    "                        cmap=\"Oranges\",\n",
    "                        ax=ax)\n",
    "# Add a title\n",
    "ax.set_title('Plot Locations\\nBear River, WY')\n",
    "# plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can use the cmap argument to adjust the colors of our plot. \n",
    "\n",
    "Below you used a colormap that is a part of the matplotlib colormap library.\n",
    "\n",
    "![Colormaps](https://matplotlib.org/_images/lightness_01.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The keyword [marker](https://matplotlib.org/3.1.1/api/markers_api.html) can be used to change the symbol."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(10, 10))\n",
    "plot_locations.plot(column='COVER',\n",
    "                         categorical=True,\n",
    "                         legend=True,\n",
    "                         marker='*',\n",
    "                         markersize=65,\n",
    "                         cmap='OrRd', \n",
    "                         ax=ax)\n",
    "ax.set_title('Plot Locations\\nBear River, WY')\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you get an error running the cell below, please see this issue [(#134)](https://github.com/pyproj4/pyproj/issues/134) with pyproj.\n",
    "\n",
    "**Workaround**: open anaconda prompt, activate your environment and then start the jupyter notebook."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import geopandas\n",
    "world = geopandas.read_file(geopandas.datasets.get_path('naturalearth_lowres'))\n",
    "print(world.crs)\n",
    "world = world.to_crs({'init': 'epsg:3395'})\n",
    "print(world.crs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Classwork:** Import Point & Line Shapefiles\n",
    "\n",
    "Using the steps above, import the data/sites and data/roads shapefiles into Python. \n",
    "\n",
    "Answer the following questions:\n",
    "\n",
    "    What type of Python spatial object is created when you import each layer?\n",
    "    What is the CRS and extent for each object?\n",
    "    Do the files contain, points, lines or polygons?\n",
    "    How many spatial objects are in each file?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import crop boundary\n",
    "roads = gpd.read_file(\"data/road.shp\")\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(10, 10))\n",
    "\n",
    "# First setup the plot using the crop_extent layer as the base layer\n",
    "roads.plot(           color='lightgrey',\n",
    "                      ax=ax)\n",
    "\n",
    "# Add another layer using geopandas syntax .plot, and calling the ax variable as the axis argument\n",
    "plot_locations.plot(     column='COVER', \n",
    "                         categorical=True,\n",
    "                         marker='D',\n",
    "                         legend=True,\n",
    "                         markersize=50,\n",
    "                         cmap='Set2', \n",
    "                         ax=ax)\n",
    "# Clean up axes\n",
    "ax.set_title('Plot Locations along North Bear Lake Boulevard\\nBear River, WY')\n",
    "\n",
    "# ax.set_axis_off()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Projection"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Geopandas has strong **CRS transformation** and mapping options.\n",
    "\n",
    "The example below shows the land boundaries in 2 projection systems."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from matplotlib.ticker import ScalarFormatter\n",
    "import numpy as np\n",
    "from shapely.geometry import Point"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Import graticule & world bounding box shapefile data\n",
    "worldBound = gpd.read_file(\"data/land.shp\")\n",
    "bbox = gpd.read_file(\"data/bbox.shp\")\n",
    "graticule = gpd.read_file(\"data/graticule.shp\")\n",
    "\n",
    "# Create numpy array of x,y point locations\n",
    "add_points = np.array([[-77.08, 38.89],            #Washington DC\n",
    "                       [ 77.14, 28.53],            #New Delhi\n",
    "                       [-0.24, 51.53]])            #London\n",
    "\n",
    "# Turn points into list of x,y shapely points \n",
    "city_locations = [Point(xy) for xy in add_points]\n",
    "\n",
    "# Create geodataframe using the points\n",
    "city_locations = gpd.GeoDataFrame(city_locations, \n",
    "                                  columns=['geometry'],\n",
    "                                  crs=worldBound.crs)\n",
    "\n",
    "# Reproject graticules and bounding box  to robinson\n",
    "worldBound_robin = worldBound.to_crs('+proj=robin')\n",
    "bbox_robinson = bbox.to_crs('+proj=robin')\n",
    "graticule_robinson = graticule.to_crs('+proj=robin')\n",
    "city_locations_robin = city_locations.to_crs('+proj=robin')\n",
    "\n",
    "# # Setup plot with 2 \"rows\" one for each map and one column\n",
    "fig, (ax0, ax1) = plt.subplots(2, 1, figsize=(13, 12))\n",
    "\n",
    "# First plot\n",
    "bbox.plot(ax=ax0,\n",
    "          alpha=.1,\n",
    "          color='grey')\n",
    "graticule.plot(ax=ax0,\n",
    "               color='lightgrey')\n",
    "worldBound.plot(ax=ax0,\n",
    "                color='black')\n",
    "city_locations.plot(ax=ax0,\n",
    "                    markersize=100,\n",
    "                    color='springgreen')\n",
    "ax0.set(title=\"World Map - Geographic (long/lat degrees)\")\n",
    "\n",
    "# Second plot\n",
    "bbox_robinson.plot(ax=ax1,\n",
    "                   alpha=.1,\n",
    "                   color='grey')\n",
    "graticule_robinson.plot(ax=ax1,\n",
    "                        color='lightgrey')\n",
    "worldBound_robin.plot(ax=ax1,\n",
    "                      color='black')\n",
    "city_locations_robin.plot(ax=ax1,\n",
    "                          markersize=100,\n",
    "                          color='springgreen')\n",
    "ax1.set(title=\"World Map Projected - Robinson (Meters)\")\n",
    "\n",
    "for axis in [ax1.xaxis, ax1.yaxis]:\n",
    "    formatter = ScalarFormatter()\n",
    "    formatter.set_scientific(False)\n",
    "    axis.set_major_formatter(formatter)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Dissolve Polygons with Geopandas"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us try to dissolve all the polygons in the \"ganga-states.shp\" file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import states shapefile data\n",
    "states = gpd.read_file(\"data/ganga-states.shp\")\n",
    "print(states)\n",
    "states.plot(cmap=\"Greys\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# There is a column called LSAD containing zeroes. Let us merge all the features based on this column.\n",
    "# Select the columns that you wish to use for the dissolve and keep.\n",
    "states1 = states[['LSAD', 'geometry']]\n",
    "basin  = states1.dissolve(by='LSAD')\n",
    "\n",
    "# View the result\n",
    "print(basin)\n",
    "basin.reset_index().plot() "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#to save the output to a shapefile is just 1 line of code!\n",
    "basin.to_file(\"data/basin.shp\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Spatial Analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us perform an intersection query.\n",
    "\n",
    "Question: To find all the features in \"populated_places.shp\" which lie inside the newly created \"basin.shp\"."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "places = gpd.read_file('data/populated_places.shp')\n",
    "print('There are totally', places.shape[0], ' features in the file')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot the data\n",
    "fig, ax = plt.subplots(figsize=(12, 8))\n",
    "places.plot(markersize=0.5, \n",
    "            ax=ax)\n",
    "basin.plot(alpha=0.7,\n",
    "           ax=ax)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![clip](https://www.earthdatascience.org/images/courses/earth-analytics/spatial-data/vector-clip.png)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Let us define a function to clip points with polygon \n",
    "#It involves 2 steps\n",
    "# - Create a mask where every point that overlaps the polygon that you wish to clip to is set to true\n",
    "# - Apply that mask to filter the geopandas dataframe.\n",
    "def clip(point, polygon):\n",
    "    poly = polygon.geometry.unary_union\n",
    "    points_clip = point[point.geometry.intersects(poly)]\n",
    "    return points_clip"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Clip the data using the clip_data module\n",
    "clipped_points = clip(places, basin)\n",
    "clipped_points.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('There are ', clipped_points.shape[0], ' number of points inside ganga basin')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## More complex analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From the highway line data (highway.shp), let us find the cumulative length in each state (ganga-states.shp).\n",
    "\n",
    "Further let us find state-wise, road density (i.e., length of road in km per unit area in sq. km) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "highway = gpd.read_file('data/highway.shp')\n",
    "states = gpd.read_file('data/ganga-states.shp')\n",
    "\n",
    "# Plot the data\n",
    "fig, ax = plt.subplots(figsize=(12, 8))\n",
    "states.plot(alpha=1,\n",
    "            facecolor=\"none\",\n",
    "            ax=ax)\n",
    "highway.plot(ax=ax)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The sjoin function performs spatial join in geopandas.\n",
    "\n",
    "The \"how\" keyword has 3 possible values (left, right, inner) [[Details]](https://github.com/geopandas/geopandas/blob/master/geopandas/tools/sjoin.py#L18).\n",
    "\n",
    "The \"op\" keyword has 10 possible values (\\_\\_eq\\_\\_, equals, almost_equals, contains, crosses, disjoint, intersects, overlaps, touches, within) [[Details]](https://shapely.readthedocs.io/en/stable/manual.html?highlight=binary%20predicates#binary-predicates)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Roads within region\n",
    "highway_ganga = gpd.sjoin(highway, \n",
    "                          states,\n",
    "                          how='inner',\n",
    "                          op='intersects')\n",
    "\n",
    "# Notice once you have joins the data - you have attributes \n",
    "# from the regions_object (i.e. NAME_1) attached to each road feature\n",
    "highway_ganga[[\"FID\", \"index_right\", \"NAME_1\"]].head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot the results\n",
    "fig, ax = plt.subplots(figsize=(12, 8))\n",
    "states.plot(alpha=1,\n",
    "            facecolor=\"none\",\n",
    "            ax=ax)\n",
    "highway_ganga.plot(column='NAME_1',\n",
    "                   ax=ax,\n",
    "                   legend=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To do distance and area measurements, it is recommended that the data should be reprojected to projected coordinate system. But, for a large geographical area like India, which projection system is ideal? It spreads across 6 UTM zones (42-47). Thankfully, there is a standardized LCC based projection system [recommended by NRSC/ISRO](https://swat.tamu.edu/docs/swat/india-dataset/2012/Data_sources.txt). The EPSG code for it is [7755](http://epsg.io/7755)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "highway_ganga_lcc = highway_ganga.to_crs({'init': 'epsg:7755'}) \n",
    "states_lcc = states.to_crs({'init': 'epsg:7755'}) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Turn off scientific notation\n",
    "pd.options.display.float_format = '{:.4f}'.format\n",
    "\n",
    "# Calculate the total length of road \n",
    "highway_ganga_lcc['rdlength'] = highway_ganga_lcc.length\n",
    "\n",
    "# Sum existing columns\n",
    "sub = highway_ganga_lcc[['rdlength', 'NAME_1']].groupby('NAME_1').sum()\n",
    "sub"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "states_lcc['area'] = states_lcc.area/1e6\n",
    "states_lcc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sub = sub.merge(states_lcc, on='NAME_1')\n",
    "sub[['NAME_1', 'rdlength', 'area']]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sub['highway_length_by_area'] = sub['rdlength']/sub['area']\n",
    "sub[['NAME_1', 'highway_length_by_area']].sort_values(by=['highway_length_by_area'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Homework**: Try out [handling missing spatial attribute data](https://www.earthdatascience.org/courses/earth-analytics-python/spatial-data-vector-shapefiles/missing-data-vector-data-in-python/)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:conda-geo_env] *",
   "language": "python",
   "name": "conda-env-conda-geo_env-py"
  },
  "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.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}