2-GetCounties-Documented.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## GeoPandas Demo: Get Counties\n",
    "This example demonstrates how to grab data from an ArcGIS MapService and pull it into a GeoPandas data frame."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import requests\n",
    "import pandas as pd\n",
    "import geopandas as gpd\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Fetching some data\n",
    "We'll tap into a NOAA map server to pull some state boundary features...\n",
    "* Build the request\n",
    "* Send the request, receive the response"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Build the request and parameters to fetch county features\n",
    "#  from the NOAA ArcGIS map server end point\n",
    "stateFIPS = '37' #This is NC\n",
    "\n",
    "url = 'https://nowcoast.noaa.gov/arcgis/rest/services/nowcoast/mapoverlays_political/MapServer/find'\n",
    "params = {'searchText':stateFIPS,\n",
    "          'contains':'true',\n",
    "          'searchFields':'STATEFP',\n",
    "          'sr':'',\n",
    "          'layers':'2',\n",
    "          'layerDefs':'',\n",
    "          'returnGeometry':'true',\n",
    "          'maxAllowableOffset':'',\n",
    "          'geometryPrecision':'',\n",
    "          'dynamicLayers':'',\n",
    "          'returnZ':'false',\n",
    "          'returnM':'false',\n",
    "          'gdbVersion':'',\n",
    "          'returnUnformattedValues':'false',\n",
    "          'returnFieldName':'false',\n",
    "          'datumTransformations':'',\n",
    "          'layerParameterValues':'',\n",
    "          'mapRangeValues':'',\n",
    "          'layerRangeValues':'',\n",
    "          'f':'json'}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Fetch the data\n",
    "response = requests.get(url,params)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Examining the response\n",
    "* Convert the response to a JSON object\n",
    "* Examine its structure\n",
    "* Extract the `attributes` and `geometry` elements."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Convert to a JSON object (i.e. a dictionary)\n",
    "respons_js = response.json()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#The 'results' object contains a record for each county returned, i.e., a feature\n",
    "results = respons_js['results']\n",
    "len(results)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Within each item in the results object are the following items\n",
    "results[0].keys()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#The 'attributes' item contains the feature attributes\n",
    "results[0]['attributes']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#And the geometry object contains the shape\n",
    "results[0]['geometry']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Convert the elements to dataFrames\n",
    "* Creating a dataFrame from the Results object\n",
    "* \"Exploding\" the dictionary values in the `attributes` and `geometry` columns\n",
    "* Concatenating dataFrames lengthwise (adding columns)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Create a dataFrame from the results, \n",
    "#  keeping just the attributes and geometry objects\n",
    "df = pd.DataFrame(results,columns=('attributes','geometry'))\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Explode the dictionary values into fields\n",
    "dfCounties = df['attributes'].apply(pd.Series)\n",
    "dfGeom = df['geometry'].apply(pd.Series)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Combine the two\n",
    "dfAll = pd.concat((dfCounties,dfGeom),axis='columns')\n",
    "dfAll.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Converting the [ESRI] geometry coordinates to a [shapely] geometric feature\n",
    "The `dfAll` dataframe now has all feature attributes and the geometry object stored in the `rings` column. \n",
    "* Exploring the 'rings' object\n",
    "* Exploring the `shapely` package: rings, polygons, and multipolygons\n",
    "* Using shapely to create features\n",
    "* Converting the dataFrame to geodataFrame\n",
    "* Plotting the output"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Explore the values in the \"ring\" column, looking at the first row of data\n",
    "rings = dfAll['rings'][0]\n",
    "print (\"There is/are {} ring(s) in the record\".format(len(rings)))\n",
    "print (\"There are {} vertices in the first ring\".format(len(rings[0])))\n",
    "print (\"The first vertex is at {}\".format(rings[0][0]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So, the \"ring\" value in each row of our dataframe contains a *list* of rings, with each ring being a list of coordinates defining the vertices of our polyon. Usually the list of rings only includes one ring, the outer boundary of a single polygon. However, it's possible it contains more than one, e.g. the boundary of Hawaii. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we'll extract the first ring object from the ring list of the first record in our dataframe and convert it to a Shapely polygon object. To do this we need to import a few Shapely geometry class objects."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Import the shapely objects we'll need\n",
    "from shapely.geometry import LinearRing\n",
    "from shapely.geometry import Polygon"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Create a shapely polygon from the first ring \n",
    "ring = rings[0]        # Get the outer ring, in coordinates\n",
    "r = LinearRing(ring)   # Convert coordinates to shapely ring object\n",
    "s = Polygon(r)         # Convert shapely ring object to shapely polygon object\n",
    "s.area                 # Show the area of the polygon"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we've seen the proof of concept, we'll form a Python function that \n",
    "* takes a list of rings (i.e., the value of one row's `rings` field), \n",
    "* converts each ring item in this ring list into a Shapely LinearRing object, \n",
    "* converts *that* into a Shapely polygon object, adding each these polygons to a list, \n",
    "* and then constructs a Shapely MultiPolygon object from the list of polygons"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#A function to convert all rings into a Shapely multipolygon object\n",
    "def polyFromRings(rings):\n",
    "    #Import necessary Shapely classes\n",
    "    from shapely.geometry import LinearRing, Polygon, MultiPolygon\n",
    "    #Construct an empty list of polygons\n",
    "    polyList = []\n",
    "    #Compile a list of shapely ring objects and convert to polygons\n",
    "    for ring in rings:\n",
    "        #Construct a ring from the ring coordinates\n",
    "        r = LinearRing(ring)\n",
    "        #Convert the ring to a shapely polygon\n",
    "        s = Polygon(r)\n",
    "        #Add the polygon to the polyList\n",
    "        polyList.append(s)\n",
    "    #Convert the list of polyongs to a multipolygon object\n",
    "    multiPoly = MultiPolygon(polyList)\n",
    "    return multiPoly"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we use Panda's `apply` method to apply the \"polyFromRings\" function above to each row's \"ring\" values. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Apply the function to each item in the geometry column\n",
    "dfAll['geometry']=dfAll.apply(lambda x: polyFromRings(x.rings),axis='columns')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Convert dataframe to a *geo*dataframe\n",
    "With the rings successfuly converted to Shapely geometry objects, we can now \"upgrade\" our Pandas dataframe to a GeoPandas dataframe, capable of spatial analysis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Create a geodataframe from our pandas dataframe (the geometry column must exist)\n",
    "gdf=gpd.GeoDataFrame(dfAll)\n",
    "gdf.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Set the projection (obtained from spatialReference column)\n",
    "gdf.crs = {'init': 'epsg:3857'}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Check the data types; note some should be fixed!\n",
    "gdf.dtypes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Convert the `ALAND` and `AWATER` to floating point values\n",
    "gdf['ALAND']=gdf['ALAND'].astype('double')\n",
    "gdf['AWATER']=gdf['AWATER'].astype('double')\n",
    "gdf.dtypes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Use familiar Pandas operation to select a feature and \n",
    "gdf[gdf['NAME'] == 'Durham'].plot();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Saving the data\n",
    "We can save the attributes to CSV file or save the feature class to a shapefile"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Save our attribute data to a shapefile\n",
    "gdf.to_csv(\"counties_{}.csv\".format(stateFIPS))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Saving data to a shapefile is a bit more finicky. In particular, we need to remove the old \"rings\" field in the geodataframe."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Delete the 'rings' column from our geodataframe\n",
    "gdf.drop('rings',axis='columns',inplace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Write the data to a file\n",
    "gdf.to_file(driver='ESRI Shapefile',filename=\"./data/NC_Counties.shp\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Recap\n",
    "So here we've imported a layer from an ESRI Map Service and done the necessary conversions to get this into a GeoPandas dataframe -- and also export it. This reveals a bit about the requirements of a geodataframe, namely the structure of the geometry column and how the Shapely package helps with that. \n",
    "\n",
    "From here, we can explore more about the cool things we can do with a GeoPandas dataframe."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}