03a_Intro_to_GIS_Vector_Data.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "EfKSXVSMZgh3"
   },
   "source": [
    "# Intro to GIS with Python\n",
    "## What is GIS?\n",
    "GIS stands for _geographic information system_. Colloquially, it's the process of presenting and analyzing data on maps. GIS allows us to visualize and characterize the nature of spatially distributed data, including weather, infrastructure, and populations. As you can imagine, this is key for disaster response scenarios for both diagnosing the situation, as well as planning and monitoring the response.\n",
    "\n",
    "There are dozens of different GIS software options, both free and commercial. In this course, we will focus on free, python-based tools and packages. The principles taught in this course should carry over to most common GIS implementations. \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "U1VkoUkxZVeg"
   },
   "outputs": [],
   "source": [
    "import geopandas as gpd\n",
    "import contextily as ctx # for basemaps\n",
    "from shapely.geometry import Point, LineString, Polygon\n",
    "from matplotlib import pyplot as plt\n",
    "import folium"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "I4FdxWjGoa3Q"
   },
   "source": [
    "## Reading in Vector data\n",
    "For this lesson we are using data in the ESRI [Shapefile](https://doc.arcgis.com/en/arcgis-online/reference/shapefiles.htm) format. \n",
    "\n",
    "Geopandas supports reading a number of different GIS vector file formats: https://geopandas.org/en/stable/docs/user_guide/io.html\n",
    "\n",
    "Geopandas uses [fiona](https://fiona.readthedocs.io/en/stable/fiona.html) to handle reading and writing vector file types."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# print out the supported file types\n",
    "import fiona; fiona.supported_drivers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will first look at some flood risk assessment data from the Philippines. This data is originally from the Humanitarian Data Exchange: https://data.humdata.org/dataset/wfp-geonode-ica-philippines-flood-risk\n",
    "\n",
    "It is currently stored in this folder as a .zip, which we will unzip using the command-line `unzip` command with the argument `-d` to provide a destination folder name for the unzipped files."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!unzip phl_ica_floodrisk_geonode_mar2014.zip -d philippines_flood_risk"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "SIt2VIrloies"
   },
   "outputs": [],
   "source": [
    "# path to shapefile\n",
    "filepath = \"philippines_flood_risk/phl_ica_floodrisk_geonode_mar2014.shp\"\n",
    "\n",
    "# Read file using gpd.read_file()\n",
    "data = gpd.read_file(filepath)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 487
    },
    "colab_type": "code",
    "id": "uSY2yeJupiBZ",
    "outputId": "fb033a3b-fba1-4fcb-d483-3b91079c5e20"
   },
   "outputs": [],
   "source": [
    "data.head() #look at top entries - looks like a pandas dataframe"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data.columns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 34
    },
    "colab_type": "code",
    "id": "nVgTgQ1fpjue",
    "outputId": "4f9ef7d2-f5e6-428f-9b83-03d580c2fdc2"
   },
   "outputs": [],
   "source": [
    "# Note the column 'geometry' is full of shapely Polygon objects\n",
    "type(data['geometry'].iloc[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that the data are in (lon, lat) ordering --- this is because the convention is (x, y) for computers, but (lat, lon) for coordinates. This is a frequent cause of error."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data['geometry']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 1000
    },
    "colab_type": "code",
    "id": "8N1mjrnqkrZY",
    "outputId": "6947c138-a01f-408e-ce54-c57fcd27ab63"
   },
   "outputs": [],
   "source": [
    "# geopandas adds useful attributes to the geodataframe, such as the ability to get bounds\n",
    "# of all the geometry data\n",
    "data.bounds"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 1000
    },
    "colab_type": "code",
    "id": "aoPpzP7BlBso",
    "outputId": "f0b529ae-6b9a-477b-9187-717914df4171"
   },
   "outputs": [],
   "source": [
    "# similary, we can get attributes such as boundary\n",
    "data.boundary"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Coordinate reference systems\n",
    "\n",
    "There are many different coordinate reference systems (CRS), which refer to different ways of indicating where on the earth you are referring to when you give a coordinate. Different CRS use different models of the earth's surface, map projections, units, and origin points (where 0,0 is). The discussion of the specifics is beyond the scope of this course. \n",
    "\n",
    "For the purposes of this course, we will primarily use the two following:\n",
    "\n",
    "### WGS 84: https://epsg.io/4326\n",
    "```\n",
    "The CRS used by the GPS system\n",
    "units: degrees\n",
    "0,0 is the intersection of greenwich meridian and equator\n",
    "epsg code: 4326\n",
    "```\n",
    "\n",
    "### Web Mercator: https://epsg.io/3857\n",
    "```\n",
    "The CRS used by most web maps, such as Google maps, OSM, Bing, etc.\n",
    "Not accurate at high latitudes >85 degrees, <-85 degrees\n",
    "units: meters\n",
    "0,0 is intersection of greensich meridian and equator\n",
    "epsg code: 3857\n",
    "```\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data.crs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 1000
    },
    "colab_type": "code",
    "id": "UTexsk7ElMV6",
    "outputId": "c0b7e2c3-b710-4c09-919b-0dc2906b8866"
   },
   "outputs": [],
   "source": [
    "# area will warn you if you're trying to do area calculations in geographic CRS\n",
    "data.area"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_in_3857 = data.to_crs('epsg:3857')\n",
    "data_in_3857.area"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "d5qkO9zEry4x"
   },
   "source": [
    "## Exercises\n",
    "Using the polygon objects in the `geometry` column of the data frame:\n",
    "- create a new column called `area` which represent the areas of each row in the shapefile\n",
    "- What are the max, min, median, and quartiles values of the areas?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "weV_oh8rmMfy"
   },
   "source": [
    "## Plotting\n",
    "Geopandas provides a useful `.plot()`  function which creates a matplotlib figure and returns an axes object.\n",
    "\n",
    "There's a ton of additional libraries that provide more plotting functionality, and we'll explore a few of them here. There's no \"correct\" set of libraries to use for GIS in python, and it's up to you to figure out which ones fit the best into your workflow.\n",
    "\n",
    "The `cmap` option to the `.plot()` function allows you to pass in a [matplotlib colormap name](https://matplotlib.org/gallery/color/colormap_reference.html), which are collections of colors used to visualize data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 224
    },
    "colab_type": "code",
    "id": "R_x4EWFgp3fH",
    "outputId": "879774cd-ef72-458d-bc54-f92a73fd9f03"
   },
   "outputs": [],
   "source": [
    "# we can use the built-in geopandas plot function to visualize\n",
    "ax = data.plot(figsize=(10,5), alpha=0.6, cmap='Set2')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "wE97qBNR29Em"
   },
   "source": [
    "currently the colors are assigned arbitrarily. However, we can also use colors to encode information. \n",
    "\n",
    "Let's first use colors to categorize by endangerment status. To do so, we pass the `column` argument to `plot()`. For reference, we also set `legend=True`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 224
    },
    "colab_type": "code",
    "id": "m0ybkz_d23BY",
    "outputId": "d6005970-8615-404c-a820-f2e04030c5db"
   },
   "outputs": [],
   "source": [
    "ax = data.plot(figsize=(10,10), alpha=0.6, cmap='Set2', column='FloodText', legend=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "BVQ8tF1q3d7L"
   },
   "source": [
    "Another common use of colors to encode data is to represent numerical data in an area with colors. This is known as a [choropleth](https://en.wikipedia.org/wiki/Choropleth_map).\n",
    "\n",
    "Let's use this to encode the areas of each region"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 312
    },
    "colab_type": "code",
    "id": "ptYUTftF3cHz",
    "outputId": "ac65bfd7-aae2-4d2d-a3db-8d768dc79372"
   },
   "outputs": [],
   "source": [
    "#then pass the area column as an argument\n",
    "ax = data.plot(figsize=(10,5), alpha=0.6, cmap='Reds', column='FloodClass', legend=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "7nNv0NL2jXLf"
   },
   "source": [
    "The data by itself looks just like a bunch of blobs. Let's put it on a map for context\n",
    "\n",
    "[Contextily](https://github.com/geopandas/contextily) is a library for creating basemaps. It pulls data from a host of different basemap providers - see [documentation](https://contextily.readthedocs.io/en/latest/) for more details.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 34
    },
    "colab_type": "code",
    "id": "mt9Avb0q0IuC",
    "outputId": "de741a50-5aea-4049-d92c-8422298b30f4"
   },
   "outputs": [],
   "source": [
    "# the data is currently in WGS84 (epsg:4326)\n",
    "data.crs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 364
    },
    "colab_type": "code",
    "id": "DCYVcUxTh3F-",
    "outputId": "1b8cb404-7783-4b23-c84e-ec88ed62a354"
   },
   "outputs": [],
   "source": [
    "ax = data.plot(figsize=(10,8), column='FloodText', legend=True)\n",
    "# now we add a basemap. ctx finds a basemap for a background from\n",
    "# an online repository.\n",
    "# It assumes the data is in web mercator (epsg:3857) unless you specify otherwise\n",
    "ctx.add_basemap(ax, crs=data.crs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# we can set bounds using matplotlib\n",
    "ax = data.plot(figsize=(10,5), cmap='Set2', column='FloodText')\n",
    "ax.set_xlim([115,130])\n",
    "ax.set_ylim([-0,25])\n",
    "ctx.add_basemap(ax, crs=data.crs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "QyvB6AAbp3Gu"
   },
   "source": [
    "We can use different background styles:\n",
    "![tile styles](https://contextily.readthedocs.io/en/latest/_images/tiles.png).\n",
    "\n",
    "Note that some styles only contain labels or lines."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 260
    },
    "colab_type": "code",
    "id": "7HGDXxwgp2Fm",
    "outputId": "f0ba1fc7-c2a0-47d2-edf4-bbd69a2cd786",
    "tags": []
   },
   "outputs": [],
   "source": [
    "# to look at all of the different providers, check:\n",
    "ctx.providers"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "previews of the different basemap styles can be viewed at: http://leaflet-extras.github.io/leaflet-providers/preview/ "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 381
    },
    "colab_type": "code",
    "id": "pknCTknLucvf",
    "outputId": "03bd2a95-42ee-4c63-f158-ebdf76e28131"
   },
   "outputs": [],
   "source": [
    "ax = data.plot(figsize=(10,5), alpha=0.6, column='FloodText', legend=True)\n",
    "# to specify the type of basemap, specify the source argument\n",
    "# the syntax is ctx.providers.{provider name}.{provider style}\n",
    "ctx.add_basemap(ax, crs=data.crs, source=ctx.providers.Stamen.Watercolor)\n",
    "# you can add labels independently of the background\n",
    "ctx.add_basemap(ax, crs=data.crs, source=ctx.providers.CartoDB.DarkMatterOnlyLabels)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 54
    },
    "colab_type": "code",
    "id": "Ilaceib1vFFh",
    "outputId": "23c7fd83-5e0e-4a02-96d2-b04a73e7d75b"
   },
   "outputs": [],
   "source": [
    "# we can download background tiles as images for quicker loading (don't need to keep redownloading)\n",
    "# let's use the bounds of the dataframe\n",
    "w,s,e,n = data.total_bounds\n",
    "data.total_bounds"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "x2kZllDjv-iv"
   },
   "source": [
    "the function bounds2img takes coordinates and [zoom level](https://wiki.openstreetmap.org/wiki/Zoom_levels) and downloads the corresponding tiles of the map as images"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 317
    },
    "colab_type": "code",
    "id": "7_uaWWzFvuC1",
    "outputId": "12fcf339-0a80-4994-d66c-d9d218d5c1b2"
   },
   "outputs": [],
   "source": [
    "img, ext = ctx.bounds2img(w, s, e, n, 6, ll=True, source=ctx.providers.Stamen.Watercolor) #ll means coordinates are in lat-lon\n",
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(1,1,1)\n",
    "ax.imshow(img, extent=ext)\n",
    "# bounds2img returns things in epsg:3857, so we need to plot the data in the same crs\n",
    "data.to_crs(epsg=3857).plot(ax=ax, cmap='Set3', alpha=0.8)\n",
    "ax_bounds = data.to_crs(epsg=3857).total_bounds\n",
    "ax.set(xlim=[ax_bounds[0], ax_bounds[2]],ylim=[ax_bounds[1], ax_bounds[3]])\n",
    "plt.axis('off')\n",
    "plt.savefig('watercolor_example.png')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Explore Interface\n",
    "Geopandas also as an interactive interface for plotting maps using the `.explore()` function. The interface uses [folium](https://python-visualization.github.io/folium/) to plot an interactive map.\n",
    "\n",
    "For full documentation on the `.explore()` function, see https://geopandas.org/en/stable/docs/reference/api/geopandas.GeoDataFrame.explore.html"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# choose a column to visualize by passing in the column name\n",
    "m = data.explore('FloodClass')\n",
    "m # folium map object"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can set the basemap using the `tiles` keyword argument and providing an [xyzproviders](https://xyzservices.readthedocs.io/en/stable/api.html#xyzservices.TileProvider) TileProvider. \n",
    "\n",
    "See the gallery preview of all available providers [here](https://xyzservices.readthedocs.io/en/stable/gallery.html). Note that some require an API key with the given provider"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import xyzservices.providers as xyz\n",
    "xyz"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data.explore('FloodClass', tiles=xyz.CartoDB.Voyager)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can plot multiple layers of data onto the same map by passing the `folium` object into the `m=` keyword argument of `.explore()`. We will use another dataframe of health site locations in the philippines as our second layer."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "health_sites_gdf = gpd.read_file('philippines_healthsites.geojson', driver='GeoJSON')\n",
    "health_sites_gdf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# note that .explore() doesn't like timestamp objects\n",
    "health_sites_gdf.explore()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "health_sites_gdf.dtypes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# we can cast that column to string to get around it\n",
    "health_sites_gdf['changeset_timestamp'] = health_sites_gdf['changeset_timestamp'].astype(str)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# plotting both on the same folium map\n",
    "m = data.explore('FloodClass', tiles=xyz.CartoDB.Voyager)\n",
    "# you can select just a subset of columns to include\n",
    "m = health_sites_gdf[['name',\n",
    "                      'amenity',\n",
    "                      'addr_street',\n",
    "                      'addr_city',\n",
    "                      'addr_postcode',\n",
    "                      'healthcare',\n",
    "                      'geometry']].explore(m=m)\n",
    "m"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# you can set the dimensions of a map by creating a \n",
    "# folium object to initially draw on\n",
    "data_centroid = data.dissolve().centroid.values[0] # dissolve combines all the geometries into one collection\n",
    "m = folium.Map(location=(data_centroid.y, data_centroid.x),\n",
    "               zoom_start=5,\n",
    "                height=600, width=500) # in pixels\n",
    "m = data.explore('FloodClass',\n",
    "                 tiles=xyz.CartoDB.Voyager,\n",
    "                 m=m,\n",
    "                 name='Flood Class')\n",
    "# you can select just a subset of columns to include\n",
    "m = health_sites_gdf[['name',\n",
    "                      'amenity',\n",
    "                      'addr_street',\n",
    "                      'addr_city',\n",
    "                      'addr_postcode',\n",
    "                      'healthcare',\n",
    "                      'geometry']].explore(m=m,\n",
    "                                          name='healthcare') # name of layer\n",
    "# add a layer control menu\n",
    "folium.LayerControl().add_to(m)\n",
    "m"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "AxVlRjR8qoHm"
   },
   "source": [
    "## Writing to a different file\n",
    "\n",
    "First we'll make a directory for outputting data to. We use the `mkdir` command which makes an empty folder. The `-p` option will skip it if the directory already exists"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "mTujtJbZqTB2"
   },
   "outputs": [],
   "source": [
    "!mkdir output_data -p"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "ElvvoSQHql4T"
   },
   "outputs": [],
   "source": [
    "# let's write the first 20 rows of the shapefile to a new file\n",
    "outfp = \"output_data/data_selection.json\"\n",
    "\n",
    "# Select first 20 rows\n",
    "selection = data[0:20]\n",
    "\n",
    "# Write those rows into a new file - we will use the GeoJSON file type\n",
    "selection.to_file(outfp, driver='GeoJSON')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "3n7QVWo9tAuV"
   },
   "source": [
    "## Converting shapes to GeoDataFrames\n",
    "You can use Shapely geometric objects to create a GeoDataFrame from scratch. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 69
    },
    "colab_type": "code",
    "id": "r6xTzsiltNiR",
    "outputId": "7e686899-7e11-4812-9a42-fe6a14778341"
   },
   "outputs": [],
   "source": [
    "# Create an empty geopandas GeoDataFrame\n",
    "newdata = gpd.GeoDataFrame()\n",
    "\n",
    "# add a geometry column (necessary for shapefile)\n",
    "newdata['geometry'] = None\n",
    "\n",
    "# Let's see what we have at the moment\n",
    "print(newdata)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 121
    },
    "colab_type": "code",
    "id": "c9fD5asutOlQ",
    "outputId": "5957bc05-4dde-4ef1-eee9-8323a169b0c9"
   },
   "outputs": [],
   "source": [
    "# Coordinates of the MIT main campus in Decimal Degrees\n",
    "coordinates = [(-71.092562, 42.357602), ( -71.080155, 42.361553), ( -71.089817, 42.362584), (-71.094688, 42.360198)]\n",
    "\n",
    "# Create a Shapely polygon from the coordinate-tuple list\n",
    "poly = Polygon(coordinates)\n",
    "\n",
    "# Let's see what we have\n",
    "poly"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 81
    },
    "colab_type": "code",
    "id": "B3-61iLVt_N5",
    "outputId": "52ff4839-275c-4cd4-a588-0a64a1a1e889"
   },
   "outputs": [],
   "source": [
    "# Insert the polygon into 'geometry' -column at index 0\n",
    "newdata.loc[0, 'geometry'] = poly\n",
    "newdata"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 81
    },
    "colab_type": "code",
    "id": "D-5voz3YwU7x",
    "outputId": "58fcce49-588b-4404-ca43-6239e38e37b6"
   },
   "outputs": [],
   "source": [
    "newdata.loc[0, 'location'] = 'MIT main campus'\n",
    "newdata"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "WnCD6scxweE9"
   },
   "source": [
    "Before exporting the data it is necessary to set the coordinate reference system (projection) for the GeoDataFrame. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 34
    },
    "colab_type": "code",
    "id": "tflAvQZ8w-Cr",
    "outputId": "7b9be928-f905-4894-86f2-5c636d0b8f1d"
   },
   "outputs": [],
   "source": [
    "# Set the GeoDataFrame's coordinate system to WGS84 (i.e. epsg code 4326)\n",
    "newdata = newdata.set_crs('epsg:4326')\n",
    "\n",
    "# Let's see how the crs definition looks like\n",
    "newdata.crs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "el0ujZVAxhp6"
   },
   "outputs": [],
   "source": [
    "outfp = \"output_data/MIT_campus.shp\"\n",
    "\n",
    "# Write the data into that Shapefile\n",
    "newdata.to_file(outfp)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 343
    },
    "colab_type": "code",
    "id": "6K9ob7oAyeIo",
    "outputId": "01887b54-cd30-4d4b-842c-30b7096fee39"
   },
   "outputs": [],
   "source": [
    "# Let's plot it\n",
    "ax = newdata.to_crs(epsg=3857).plot(figsize=(10,5),alpha = 0.5, color='#FF55FF')\n",
    "ctx.add_basemap(ax)\n",
    "ax.set_axis_off() # remove the x-y axes\n",
    "plt.savefig('MIT_main_campus_poly.png')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "MCNXf-5K2DDx"
   },
   "source": [
    "# Exercise\n",
    "Find an interesting GIS dataset and:\n",
    "- visualize some raw data\n",
    "- ask an interesting analysis question about it:\n",
    "  - intersections, sizes, quantities\n",
    "  - relationships\n",
    "  - e.g. which latitudes contain the most endangered species? what countries have the most ports per km of coastline?\n",
    "- Visualize some of your analysis\n",
    "You can use the location you've chosen for your location fan-cam as a place of interest! \n",
    "\n",
    "Note that since geopandas is built on pandas, all of your knowledge from pandas should also carry over!\n",
    "\n",
    "As per usual, we'll ask a few volunteers to present their results.\n",
    "\n",
    "Here are some resources to look for GIS datasets:\n",
    "- Cambridge, MA GIS data: http://cambridgegis.github.io/gisdata.html\n",
    "- Awesome GIS data: https://github.com/sshuair/awesome-gis#data\n",
    "- Humanitarian Data Exchange: https://data.humdata.org/\n",
    "- Data.gov: https://www.data.gov/\n",
    "Search for GeoJSON and/or Shapefile file types.\n",
    "\n",
    "One tool to help draw GIS polygons is https://geojson.io; you can export your polygon as a geojson and upload it to your jupyter instance to access from jupyter."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "RhKrcnw27-p8"
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "gfJNj0A9Kd2C"
   },
   "source": [
    "### Covid-19 Resources\n",
    "For those interested in how GIS can be used to analyze the pandemic, here are some ideas and data:\n",
    "\n",
    "*   Visualize raw data collected from sources around the world about the state of the pandemic\n",
    "*   Explore connections between various factors and come up with a hypothesis for your research. Some ideas could be connecting COVID data in different counties to socioeconomy, age, or building architecture data. Remember, mapping data speaks louder than graphs or datasets.\n",
    "*   Present your findings to the rest of the class and come up with a possible solution to the problem or connection that you explored\n",
    "\n",
    "COVID-19 Datasets:\n",
    "* COVID-19 Dataset (Kaggle): www.kaggle.com/imdevskp/corona-virus-report\n",
    "* New York Times Dataset: https://github.com/nytimes/covid-19-data\n",
    "* JHU Dataset: https://github.com/CSSEGISandData/COVID-19/tree/master/csse_covid_19_data\n",
    "*  Feel free to explore more area specific datasets or datasets which outline other conditions. These are just suggestions.\n",
    "\n",
    "To make your research connections, be sure to explore population and demographic datasets of different counties around the country. Be creative with your research!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "colab": {
   "collapsed_sections": [],
   "name": "03_Intro_to_GIS.ipynb",
   "provenance": []
  },
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "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
}