PySAL - Spatial Autocorrelation - LISA - Start with this script.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "# Spatial Autocorrelation with PySAL"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "## Local Moran's I exercise (Local Indicators of Spatial Association, or LISAs) \n",
    "\n",
    "Using this Jupyter notebook, you will code the same process that we used in GeoDa to run the Univariate Local Moran's I to identify local hot/cold spots in the data.\n",
    "\n",
    "Adapated from https://pysal.readthedocs.org/en/latest/users/tutorials/autocorrelation.html#local-moran-s-i\n",
    "Information also available on pg 39 in PySAL_Documentation.pdf included with the SDS Bootcamp materials.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "### Task #1: read data and create variables"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "#### First, import the needed Python components - all scripting in Python begins with import"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "# PySAL and Numpy the only ones needed to actually run the spatial autocorrelation analysis\n",
    "import pysal\n",
    "import numpy as np\n",
    "\n",
    "#Folium is used to create some map visualizations\n",
    "import folium\n",
    "\n",
    "# These others are to handle, query, and plot data\n",
    "import os\n",
    "import pandas as pd\n",
    "import geopandas as gpd\n",
    "import json\n",
    "import simpledbf\n",
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# This message below will print after the commands above are successfully completed\n",
    "print ('All requested Python libraries imported successfully')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "#### To calculate Moran's I, we need to give it a list of neighbors. We can do this by reading in a spatial weights matrix.\n",
    "Remember that PySAL likes the GAL file format, which we created in the GeoDa exercise. \n",
    "\n",
    "This file can be converted to an ArcGIS Spatial Weights Matrix (SWM) file, and vice versa.\n",
    "PySAL can read the GAL file as follows: counties = pysal.open('path/to/file/called/file.gal').read()\n",
    "\n",
    "#### Alternatively, and much easier, PySAL can read the neighbors directly from a shapefile"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "# Instead of reading in the .gal file we created in GeoDa, we will ask PySAL to create one from the shapefile.\n",
    "counties = pysal.queen_from_shapefile('/home/ubuntu/Documents/Counties/cnty_Lyme_disease.shp')\n",
    "\n",
    "# This message below will print after the command above is successfully completed\n",
    "print ('New weights file successfully created')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "#### The queen_from_shapefile function has defined neighbors using the the queen weights criteria, which defines a location's neighbors as those areas with at least one shared corner\n",
    "\n",
    "#### Other options include a rook weights matrix, in which neighbors need to share an entire border (i.e. a line of two connected vertices)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "# Next, let's read in the dbf that contains data for the counties\n",
    "# http://www.pysal.org/users/tutorials/fileio.html\n",
    "\n",
    "table = pysal.open('/home/ubuntu/Documents/Counties/cnty_Lyme_disease.dbf')\n",
    "\n",
    "# This message below will print after the command above is successfully completed\n",
    "print ('DBF file successfully imported')\n",
    "\n",
    "# Return the column headers from dbf\n",
    "table.header"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "# Next, specify which column contain the variable of interest\n",
    "# Notice that we are using the array function from numpy, which we named np during our import\n",
    "# This array will contain the data from the column called 2005\n",
    "\n",
    "lymecases = np.array(table.by_col['2005'])\n",
    "\n",
    "# This message below will print after the command above is successfully completed\n",
    "print ('A variable called lymecases successfully created')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "### Task #2: complete a single run of Local Moran's I to identify local hot/cold spots in the data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "# Using the functions examples below, update the parameters to run Local Moran's I on the year 2005 data\n",
    "\n",
    "# In the online tutorial, the function reads as follows: lm = pysal.Moran_Local(y,w)\n",
    "# In their example, y = array containing homicide rates and w = spatial weights variable for the neighbors\n",
    "\n",
    "# Another example could be something like: lm = pysal.Moran_Local(crimeindex, blockgroups)\n",
    "# In this example, crimeindex is the array containing a crime index and blockgroups is the spatial weights variable \n",
    "\n",
    "# We also want to run multiple permutations of the random distribution\n",
    "# Hint: lm = pysal.Moran_Local(y,w,permutations = 1)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "# Again, examine the help to learn more about the outputs from your function\n",
    "#help(lm)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "# Which attribute of lm can you use to see the actual observed Local Moran's I values (LISAs)?\n",
    "# Hint: lm.attributename\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "# Which attribute of lm can you use to see the statistical significance of difference between I and simulated I values (LISAs)?\n",
    "# Hint: lm.attributename\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Which attribute of lm can you use to check the cluster type for each LISA?\n",
    "# Hint: lm.attributename\n",
    "# Remember: 1 = High/High 2 = Low/High 3 = Low/Low 4 = High/Low \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "### Task #3: check for significant Local Moran's I values for the LISAs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "# Hint: create a new variable and make it equal to a boolean statement regarding the p value being less than a certain value\n",
    "# Hint: variablename = lm.attributename<value\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "### Task #4: create a csv file of Local Moran's I output results for all counties"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Create a variable that will hold the county identifier, similar to how you made an array of the 2005 values\n",
    "# Hint: variablename = np.array(f.by_col['NameofColumninDBF'])\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "# Next, export results using numpy export to save the LISAs to a csv with the County list\n",
    "# This csv can be joined to the county shapefile to map the LISAs and highlight hot/cold spots \n",
    "# Uncomment following lines and fill in as needed to export to csv\n",
    "# Be sure to include all the output variables you explore previously; there should be 3 plus one you created based on output values\n",
    "\n",
    "# Hint: np.savetxt(\"filename.csv\", np.column_stack((countylist, lm.attributename, lm.attributename, lm.attributename, variableyoucreated)), delimiter=' , ', fmt=\"%s\")\n",
    "# print 'a message to yourself that the file is ready'\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "## Bonus Exercises for Local Moran's I exercise"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "### Bonus #1: add the original lyme disease counts to the csv file of Local Moran's I output results for all counties"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "# Hint: use np.savetxt from Task #4 and add a new parameter to the column_stack function\n",
    "# Hint: examine the help to learn more about the outputs from your function\n",
    "# Hint: help(lm)  --> which output contains the original data? \n",
    "# Hint: you also previously created a variable that containts this info, too\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "### Bonus #2: add a header to the csv file of Local Moran's I output results for all counties"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "# Hint: http://docs.scipy.org/doc/numpy-1.10.0/reference/generated/numpy.savetxt.html\n",
    "# Hint: np.savetxt(\"filename.csv\", np.column_stack((countylist, lm.attributename, lm.attributename, lm.attributename, sig)), delimiter=' , ', fmt=\"%s\", header='FirstColumnName, SecondColumnName, etc', comments = '')\n",
    "# Hint: leave the comments parameter empty (Or, check out what happens if you don't)\n",
    "# print 'a message to yourself that the file is ready'\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "## Advanced Bonus Exercises for Local Moran's I exercise"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "### Advanced Bonus #1: use a loop to run Local Moran's I for other years of data (2006 to 2014)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "# Note: these advanced bonus exercises may be challenging for new Python programmers \n",
    "# Hint: review Advanced Bonus Exercise #1 for Global Moran's I\n",
    "# Hint: create a list of the years to begin\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "# Hint: check out Lists as an iterable on: https://wiki.python.org/moin/ForLoop \n",
    "# Hint: review Advanced Bonus Exercise #1 for Global Moran's I\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "### Advanced Bonus #2: expand your loop to create a different csv for each year"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "# Note: these advanced bonus exercises are challenging for new Python programmers\n",
    "# Hint: you will need to combine Task #4 (or Bonus #1/#2) and Advanced Bonus #1\n",
    "# Hint: your loop will need to include the csv name based on the year of analysis\n",
    "# Hint: create a variable for the csvname and append the year to the name in the loop\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "## End of Local Moran's I exercise\n",
    "\n",
    "### Other options to continue your exploration:\n",
    "#### You can run Global or Local Moran's I on other datasets for crime in the Bonus Data folder\n",
    "#### You can explore the new csv in QGIS and join with the county shapefile to map out the Local Moran's I values. If you need help doing a table join in QGIS, check out this easy tutorial: http://www.qgistutorials.com/en/docs/performing_table_joins.html \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "#### Advanced Python users: maybe continue to explore pandas. See if you can figure out how to combine your multi-year output to a single pandas dataframe:  http://pandas.pydata.org/pandas-docs/stable/10min.html "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "# Hint: check out Lists as an iterable on: https://wiki.python.org/moin/ForLoop \n",
    "# Hint: review Advanced Bonus Exercise #1 for Global Moran's I\n",
    "# Hint: we can create a separate data frame for each year and then concatenate all the years into one data frame\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "# We can play with this data frame a little bit\n",
    "\n",
    "# For example: which counties in 2012 have significant LISA?\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "### Or just for fun: explore some advanced visualizations below"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Create a map of the P Value for each LISA value for 2005\n",
    "\n",
    "# To use Folium for map visualizations, we convert the shapefile to GeoJSON\n",
    "\n",
    "# We can use GeoPandas for this conversion\n",
    "# First, read in the shapefile to geopandas\n",
    "shapefile = gpd.read_file('/home/ubuntu/Documents/Counties/viz/cnty_Lyme_disease_WGS84.shp').set_index('NAME_PCASE')\n",
    "\n",
    "# Next, save the file out to GeoJSON\n",
    "with open('/home/ubuntu/Documents/Counties/viz/cnty_Lyme_disease_WGS84.geojson', 'w') as f:\n",
    "    f.write(shapefile.to_json())\n",
    "\n",
    "# This message below will print after the commands above are successfully completed\n",
    "print ('Successfully converted shapefile to geojson')\n",
    "\n",
    "# Use simpledbf to read the attribute table (dbf file)\n",
    "dbf = simpledbf.Dbf5('/home/ubuntu/Documents/Counties/cnty_Lyme_disease.dbf')\n",
    "\n",
    "# Read in the GeoJSON file created in previous step\n",
    "counties = '/home/ubuntu/Documents/Counties/viz/cnty_Lyme_disease_WGS84.geojson'\n",
    "\n",
    "# Convert dbf file to a pandas data frame\n",
    "df = dbf.to_dataframe()\n",
    "\n",
    "# Store the quadrant values in a variable called q\n",
    "q = lm.q\n",
    "\n",
    "# Store the p values in a variable called q\n",
    "df[\"p\"] = lm.p_sim\n",
    "\n",
    "# Change the quadrant values of the counties whose p values are not significant to 0 for convenience of graphing\n",
    "q = [q[i] if lm.p_sim[i] < 0.0501 else 0 for i in range(0,58)]\n",
    "\n",
    "# Update the dataframe with the new quadrant values \n",
    "df[\"q\"] = q\n",
    "\n",
    "# load GeoJSON data\n",
    "geo_json_data = json.load(open(counties))\n",
    "\n",
    "# Assign the p values and q values in to the GeoJSON object\n",
    "for i in range(0,58):\n",
    "    geo_json_data[\"features\"][i][\"properties\"][\"q\"] = str(df[\"q\"][i])\n",
    "    geo_json_data[\"features\"][i][\"properties\"][\"p\"] = str(df[\"p\"][i])\n",
    "\n",
    "# Create an empty Folium map\n",
    "m = folium.Map([37, -122], zoom_start=6)\n",
    "\n",
    "# Create a map based on significance(p-values)\n",
    "folium.GeoJson(\n",
    "    geo_json_data,\n",
    "    style_function=lambda feature: {\n",
    "        'fillColor': '#006400' if float(feature['properties']['p']) <= 0.0001 else '#008000' if (float(feature['properties']['p']) > 0.0001 and \n",
    "float(feature['properties']['p']) <= 0.001) else '#228B22' if (float(feature['properties']['p']) > 0.001 and \n",
    "float(feature['properties']['p']) <= 0.01) else \"#00FF00\" if (float(feature['properties']['p']) > 0.01 and \n",
    "float(feature['properties']['p']) < 0.05001) else '#808080',\n",
    "        'color': 'black',\n",
    "        'fill_opacity': 1.0,\n",
    "        'weight': 2,\n",
    "        'dashArray': '5, 5',\n",
    "        'fillOpacity': 0.6\n",
    "    }\n",
    ").add_to(m)\n",
    "\n",
    "# Save the map as html\n",
    "m.save(os.path.join(\"/home/ubuntu/Documents/Counties/viz/p_map.html\"))\n",
    "\n",
    "# Print the map to screen\n",
    "m"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "# We can extend the code in the previous cell to create a second map for 2005 \n",
    "# of the quadrants types based on significance from p-values\n",
    "\n",
    "# Counties with non significant p values are in grey\n",
    "# Counties with significant p values are in different colors depending on quadrant type\n",
    "# 1 = High/High in Red\n",
    "# 2 = Low/High in Light Blue\n",
    "# 3 = Low/Low in Dark Blue\n",
    "# 4 = High/Low in Pink\n",
    "\n",
    "# Create an empty Folium map\n",
    "n = folium.Map([37, -122], zoom_start=6)\n",
    "\n",
    "# Create a map of the quadrant labels based on significance (p-values)\n",
    "folium.GeoJson(\n",
    "    geo_json_data,\n",
    "    style_function=lambda feature: {\n",
    "        'fillColor': 'darkblue' if int(feature['properties']['q']) == 3 else 'lightblue' if int(feature['properties']['q']) == 2\n",
    "    else 'darkred' if int(feature['properties']['q']) == 1 else \"pink\" if int(feature['properties']['q'])== 4 else 'gray',\n",
    "        'color': 'black',\n",
    "        'weight': 2,\n",
    "        'dashArray': '5, 5',\n",
    "        'fillOpacity': 0.6\n",
    "    }\n",
    ").add_to(n)\n",
    "\n",
    "# Save the map as html\n",
    "n.save(os.path.join(\"/home/ubuntu/Documents/Counties/viz/q_map.html\"))\n",
    "\n",
    "# Print the map to screen\n",
    "n"
   ]
  }
 ],
 "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.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}