<jupyter_start><jupyter_text>Spatial Visualization  Table of Contents1. [Simple maps](Simple-maps)2. [Choropleth](Choropleths)2. [Choropleth of rate](Choropleth-of-rate)3. [Choropleth multiple years](Choropleth-multiple-years)[Addendum](Addendum)1. Overview of [spatial data types](Spatial-data-types)2. [Datasets](Datasets) in this Notebook2. [Coordinate Reference Systems (aka projections)](Coordinate-Reference-Systems) Reference material**NOTE: much of the below material combines using SQL via Pandas and Geopandas to get the data in the format we want for spatial analyses or visualization in Python, As a reminder, the "[Spatial Queries in PostGIS](../spatial_notebooks/Spatial-Queries-in-PostGIS.ipynb)" has more detailed explanations of PostGIS and will help you better understand the underlying spatial SQL in this notebook.**  Simple maps- Back to [Table of Contents](Table-of-Contents)In this notebook we will make use of the 2014Q3 subset of IDES establishment point locations created for the "spatial queries in PostGIS" notebook linked above.<jupyter_code>## data manipulation libraries ##
# Pandas for generic manipulation
import pandas as pd
# GeoPandas for spatial data manipulation
import geopandas as gpd

# PySAL for spatial statistics
import pysal as ps

# shapely for specific spatial data tasks (GeoPandas uses Shapely objects)
from shapely.geometry import Point, LineString, Polygon

# SQLAlchemy to get some data from the database
from sqlalchemy import create_engine

# improve control of visualizations
import matplotlib.pyplot as plt
%matplotlib inline
# set up connection paramaters for the 'appliedda' database
host = '10.10.2.10'
db = 'appliedda'

# create DB connection 
# (connection string is "<database-type>://<username>@<host-url>:<port>/<database-name>")
engine = create_engine("postgresql://{host_str}/{db_str}".format(host_str = host, db_str=db))
# read in point locations - the full file has ~455k locations (~208k of which are in Cook County)
# so here we are going to grab a sample of 500 from within Cook County
query = """
SELECT e.*
FROM class2.il_des_subset_2014q3 e
JOIN tl_2016_us_county c
ON ST_Within(e.geom, c.geom)
WHERE c.geoid = '17031'
LIMIT 500;
"""
points = gpd.read_postgis(query, engine, geom_col='geom',
                          crs={'init': u'epsg:4269'})
# what info is contained in this file?
points.info()
# what does the first row look like?
points.head(1)
# what does a simple map of those locations look like? 
# Note that you can pass matplotlib keywords to (geo)pandas, like figsize
points.plot(figsize=(4,6))
plt.annotate('Source: IDES', xy=(0.8,-0.10), xycoords="axes fraction");<jupyter_output><empty_output><jupyter_text>> There's some dots. Not super useful, but we can see by the longitude (x axis) and latitude (y axis) they're where we'd expect for Chicago. Let's add Chicago zip codes too to give some context, to do so we'll pull the polygons from the database.<jupyter_code># create SQL query
sql = """
SELECT z.* 
FROM tl_2016_us_zcta510 z
JOIN tl_2016_us_county c
ON ST_Within(z.geom, c.geom)
WHERE c.geoid = '17031';
"""

# get data from the database
polys = gpd.read_postgis(sql, engine, geom_col='geom', 
                               index_col='gid', crs={'init': u'epsg:4269'})
# note we need to include the "geom_col" and "crs" parameters
# so that geopandas knows which column containes the geographic info
# and the correct Coordinate Reference System of the data


# see info of new geodataframe
polys.info()
# plot counties - using all defaults in pandas's plot() function
# adds colors essentially randomly based on order of shapes in the dataframe
polys.plot()
# now let's try putting the point locations on top of counties, 
# and let's make counties grey

# create a map of IL counties colored grey
ax = polys.plot(color='grey', figsize=(6,8))
# note we assign counties to the object "ax" so we can 
# overlay the points on the same "matplotlib axis"

# use the same "ax" object to plot the prisons on top of the county map, 
# plus resize the markers and remove their outlines
points.plot(ax=ax, markersize=5, markeredgewidth=0); 
plt.annotate('Source: Census TIGERLines & IDES', xy=(0.8,-0.10), xycoords="axes fraction")
plt.title('Business locations');
# pro tip: adding this semi-colon at the end stops Jupyter 
# from printing out the "<matplotlib.axes....>" line<jupyter_output><empty_output><jupyter_text>Choropleths- Back to [Table of Contents](Table-of-Contents)Choropleths are very useful because they can quickly convey how values compare across the study area. Let's start wtih a simple example of the land area of each county. (Note much of the code below comes from Sergio and Dani's Geovisualization notebook -> source: http://darribas.org/gds_scipy16/ipynb_md/02_geovisualization.html)Let's start with a simple dataset of number businesses by Chicago neighborhood. SQL used to create the "chi_nhood_estab_2014q3" table:    SELECT n.neighborhood, count(e.*) tot_biz,         sum(e.total_wages) wages_tot, sum(e.total_jobs) jobs_tot, n.geom     INTO class2.chi_nhood_estab_2014q3    FROM chicago_nhoods AS n    JOIN class2.il_des_subset_2014q3 e    ON ST_Within(e.geom, n.geom)    GROUP BY n.neighborhood, n.geom    ORDER BY n.neighborhood;        ALTER TABLE class2.chi_nhood_estab_2014q3 ADD PRIMARY KEY (neighborhood);    CREATE INDEX chi_nhood_estab_2014q3_geom_gist ON class2.chi_nhood_estab_2014q3 USING gist (geom);<jupyter_code># Get data 
qry = "SELECT * FROM class2.chi_nhood_estab_2014q3"

nhoods = gpd.read_postgis(qry, engine, geom_col='geom', 
                          crs='+init=epsg:4269')
nhoods.info()
# what are the descriptive stats?
nhoods.describe()
# we'll create our matplotlib figure and axis first so we have more control over it later
f, ax = plt.subplots(figsize=(6,8))

# we'll pass the geopandas plot() function the column, scheme (calculation method), 
# number of groups to calculate (k)
# colormap to use, linewidge (to make the edges less noticeable), and 
# the axis object created above
nhoods.plot(column='tot_biz', scheme='QUANTILES', k=10, 
               cmap='plasma', linewidth=0.1, ax=ax)

# and this time we'll turn off the coordinates
ax.set_axis_off();<jupyter_output><empty_output><jupyter_text>> As you can see, Geopandas only allows using the "quantiles" (or any other [scheme supported by PySAL](http://pysal.readthedocs.io/en/latest/library/esda/mapclassify.html)) with between 2 and 9 and if you try soemthing different, it resets to 5. So here is how you can use more categories for your choropleths: create a new column with the appropriate PySAL function and map that, as follows. By "appropriate PySAL function" we mean whichever classification algorithm you select, above we used Quantiles and below we are using the Fisher Jenks algorithm. The PySAL documentation is quite good so we encourage you check there (as well as the notebook from <NAME>'s colleagues referenced above).<jupyter_code># let's try the 'Fisher_Jenks' scheme:
fj10 = ps.Fisher_Jenks(nhoods.tot_biz,k=10)

# the ps.<scheme> function returns two things, the bins used for the cutoffs:
print('bins:')
print(fj10.bins)
# and the assigned bin number to use:
print('\nbin number:')
print(fj10.yb)
# now we can use the new categories to color the choropleth of land area into 10 buckets
# notice the couple new keywords we include

# again we'll create the matplotlib figure and axis objects first
f, ax = plt.subplots(figsize=(6, 8))

# then create our choropleth, 
# the "assign" function adds our Fisher Jenks buckets as a new column to map
## Note with this formulation 'cl' is a temporary column
# the 'categorical=True' tells geopandas we want a different color for each category
nhoods.assign(cl=fj10.yb).plot(column='cl', categorical=True, \
        cmap='viridis', linewidth=0.1, ax=ax, \
        edgecolor='grey', legend=True)

# turn off the latitude/longitude axes
ax.set_axis_off();<jupyter_output><empty_output><jupyter_text>The above map suggests a few neighborhoods have more businesses than others, some of which appear to be in Chicago, however this brings us to an issue central to much of spatial analysis: do we actually care about the _total number of businesses_? Possibly. Often, however, we care more about **where there are _higher concentrations_ of our study variable as compared to some base value**, so next we will consider the rate of IDHS household cases by housing units for 2015.  Choropleth of rate- Back to [Table of Contents](Table-of-Contents)**NOTE: BELOW THIS POINT IS YET TO BE UPDATED**The "chi_nhood_hhcase2015" table in the class2 schema was created with the following SQL:    SELECT n.neighborhood, sum(q.hh_counts) hhcase_total, n.geom    INTO class2.chi_nhood_hhcase2015    FROM chicago_nhoods AS n    JOIN (SELECT h.hh_counts, g.geom_2163         FROM idhs.hh_indcase_spells h        JOIN idhs.case_geocode g        ON h.ch_dpa_caseid = g.ch_dpa_caseid        WHERE h.start_date >= '2015-01-01' AND h.start_date < '2016-01-01') q    ON ST_Within(q.geom_2163, n.geom_2163)    GROUP BY n.neighborhood, n.geom    ORDER BY n.neighborhood;                ALTER TABLE class2.chi_nhood_estab_2014q3 ADD PRIMARY KEY (neighborhood);    CREATE INDEX chi_nhood_estab_2014q3_geom_gist ON class2.chi_nhood_estab_2014q3 USING gist (geom);<jupyter_code># get data
qry ="""
"""

# get new data
nhoods_hh = gpd.read_postgis(qry, engine, geom_col='geom', index_col='gid')

# view info of returned GeoDataFrame
nhoods_hh.info()
# what are the descriptive stats of the rate column?
nhoods_hh.describe()
# replace nulls with 0 as we assume there are no ex-offenders who were released in 2015
nhoods_hh[''].fillna(0, inplace=True)
# again we'll use fisher-jenks with 10 classes
fj10 = ps.Fisher_Jenks(nhoods_hh[''],k=10)

# again we'll create the matplotlib figure and axis objects first
f, ax = plt.subplots(figsize=(6, 8))

# then create our choropleth, the "assign" function adds our 
# Fisher Jenks buckets as a new column to map
## Note as before with this formulation 'cl' is a temporary column
# the 'categorical=True' tells geopandas we want a different color for each category
nhoods_hh.assign(cl=fj10.yb).plot(column='cl', categorical=True, \
        cmap='viridis_r', linewidth=0.1, ax=ax, \
        edgecolor='grey', legend=True)

# let's add a title to describe what we're plotting, the '\n' adds a new line
my_title = ''
ax.set_title(my_title)

# turn off the latitude/longitude axes
ax.set_axis_off();<jupyter_output><empty_output><jupyter_text>Choropleth multiple years- Back to [Table of Contents](Table-of-Contents)For something more interesting let's visualize...<jupyter_code># query for list of zipcodes:
query = """"""

# get the zipcodes as a pandas dataframe
zipcodes = pd.read_sql(query, engine)
# how many zipcodes is this?
len(zipcodes)
# now convert the zipcode dataframe into a string 
# with single quotes and commas separating each
# to select out just those zipcodes in the following SQL query
zipcodes = ','.join(["'"+z+"'" for z in zipcodes.zipcode])
# almost our same query from above Choropleth of rate, with two differences
# 1. get year from the exit_date and limit to our most recent 10 years
# 2. we add year as a GROUP BY in the subquery
qry = """z.geoid10 IN ({});""".format(zipcodes)

# get new data
nhoods_hh_10yr = gpd.read_postgis(qry, engine, geom_col='geom', index_col='gid')

# view info of returned GeoDataFrame
nhoods_hh_10yr.info()
# view summary stats of full dataset
nhoods_hh_10yr.describe()
# let's make a new column of quantiles to make comparing our years easier
nhoods_hh_10yr['quant5'] = ps.Quantiles(ex_off_rate_10yr.exoffrate,k=5).yb
# for each year, how many records are in each of our quantile buckets?
grouped_counts = nhoods_hh_10yr.groupby(['year', 'quant5'])['zip5'].count().reset_index()

# what does this dataframe look like?
grouped_counts.info()
# let's see it as a 10x5 table
grouped_counts.pivot_table(index='year', columns='quant5', values='')
# create 2d list of years
years = [range(2006, 2011, 1), range(2011, 2016, 1)]
years<jupyter_output><empty_output><jupyter_text>**Warning**: following cell takes some time to run so it would be good if people attempt to coordinate in their groups to not all run at the same time<jupyter_code># now we'll make a choropleth for each year all on the same figure
# using mulitple subplots
f, axes = plt.subplots(2, 5, sharex=True, sharey=True, figsize=(10,6))

# do two for loops to plot each year in position
for row in xrange(2):
    for col in xrange(5):
        temp_gdf = nhoods_hh_10yr[nhoods_hh_10yr.year==years[row][col]]
#         print('temp gdf %s %s created' % (row, col))
        temp_gdf.plot('quant5', categorical=True, \
            cmap='viridis_r', linewidth=0.1, ax=axes[row,col], \
            edgecolor='grey')
        # also let's turn off the axis lables as they're not so useful here
        axes[row,col].set_yticklabels([])
        axes[row,col].set_xticklabels([])
        # and set the title column as the year
        axes[row,col].set_title(str(years[row][col]))
# add a title for the whole figure
f.suptitle('');<jupyter_output><empty_output><jupyter_text>Addendum- Back to [Table of Contents](Table-of-Contents) Addendum contents1. [Spatial data types](Spatial-data-types)2. [Datasets](Datasets) in this notebook2. Projections and [Coordinate Reference Systems (CRS)](Coordinate-Reference-Systems)2. [Spatial info of tables in PostGIS](Spatial-info-of-tables-in-PostGIS)  Spatial data types- Back to [Table of Contents](Table-of-Contents)- Back to [Addendum contents](Addendum-contents)There are two generic spatial data types:1. **Vector** - discrete data (usually), represented by points, lines, and polygons2. **Raster** - continuous data (usually), generally represented as square pixels where each pixel (or "grid cell") has some value. Examples of raster data - link to "big data"  * Imagery data (satellite, Google SteetView, traffic cameras, Placemeter)  * Surface data (collected at monitoring stations then interpolated to a 'surface' - eg Array of Things, weather data)  However, raster data is commonly used in few social science contexts, so the below image (courtesy of [Data Science for Social Good](https://github.com/geebioso/postgis-workshop/blob/master/tutorial.org)) is probably sufficiet discussion about rasters for now:![raster](../../data/sample_data/raster_example.png)> Notice the pesky _"usually"_ next to both vector and raster datatypes? Technically any data **_could_** be represented as either vector or raster, but it would be computationally inefficient to create a raster layer of rivers or roads across Illinois because 1. All the non-road and non-river locations would still have some value and 2. You would have to pick a cell size which may not well represent the actual course of a given river (as opposed to a vector - line - that follows a path and could have some value for width of the river or road)  Datasets- Back to [Table of Contents](Table-of-Contents)- Back to [Addendum contents](Addendum-contents)Datasets used in this Notebook1. US Counties - polygons of counties from US Census's TIGER\Lines product (source-> https://www.census.gov/geo/maps-data/data/tiger-line.html)3. Chicago neighborhoods - polygons of Chicago's 88 neighborhoods created using the [Tract to Neighborhood](https://deepdish.adrf.info/detail/adrf-000073) lookup table to create polygons based on underlying Census Blocks. Code to create the neighborhood table is in [this notebook](https://jupyterhub.adrf.info:8000/user/crh278/notebooks/Projects/class2/user/crh278/class2_materials/data/Create-Chicago-neighborhood-spatial-table.ipynb)4. IDES addendum to employer locations, [link to ADRF Explorer page](https://deepdish.adrf.info/detail/adrf-000035)Considerations for spatial data+ Data collection - at what spatial scale were data collected?+ Data aggregation issue: aggregating to different spatial units could give different results due to something known as the "modifiable areal unit problem" (MAUP). MAUP is a statistical bias from summarizing point data to areas as the value of a given area depends on where the boundary is drawn (description paraphrased from Wikipedia, which has more info and links to related issues) One common type of vector data is called a "shapefile". Shapefiles are actually a collection of files and there are:1. Required files  + .dbf - contains attributes (ie variables that describe the object, like number of households in a zip code)  + .shp - the geographic information of the object  + .shx - an index to facilitate searching within data2. Additional / optional files  + .prj - the projetion or coordinate reference sysstem (CRS); technically optional but **highly** recommended as it tells you how the data relates to real-world locations  + .sbn / .sbx - other indexes used by specific software  + .shp.xml - metadata information  > Above description pulled mostly from the Wikipedia article on [shapefiles](https://en.wikipedia.org/wiki/Shapefile) which also has much more information.Other common vector formats:1. GeoJSON - Geospatial standards in the JavaScript Object Notation (JSON) text format2. KML / KMZ - Keyhole Markup Language popularized by Google3. GDB - Esri's Geodatabase, a proprietary data format used in ArcGIS software products  Coordinate Reference Systems- Back to [Table of Contents](Table-of-Contents)- Back to [Addendum contents](Addendum-contents)Coordinate Reference Systems (aka projections) are basically math that (1) describes how information in a given dataset relates to the rest of the world and (2) usually creates a 'flat' surface on which data can be analyzed using more common algorithms (eg Euclidean geometry). >Why do we care?1. Distance / area measurements2. Spatial join - won't work with different CRSAs an example of point 2, consider the distance between two points: Euclidean distance (aka pythagorean theorem) provides an easy way to calculate distance so long as we know the difference in **_x_** (longitude) and **_y_** (latitude) between two points:$$Distance   = \sqrt(({x}_1-{x}_2)^2 + ({y}_1-{y}_2)^2)$$This works fine on **_correctly projected_** data, but **_does not work_** on unprojected data. For one thing the result would be in degrees and degrees are a different distance apart (in terms of meters or miles) at different points on the Earth's surface.All this is to say: if you do a calculation with geographic data and the numbers don't make sense, check the projection. Let's do an example with the IL county areas.<jupyter_code># create DB connection (connection string is "<database-type>://<username>@<host-url>:<port>/<database-name>")
engine = create_engine("postgresql://10.10.2.10:5432/appliedda")

# create SQL query - limit to just IL by using the state FIPS code of 17
sql = "SELECT * FROM tl_2016_us_county WHERE statefp = '17';"

# get data from the database
il_counties = gpd.read_postgis(sql, engine, geom_col='geom', index_col='gid')
# print out the CRS of IL counties:
print(il_counties.crs)<jupyter_output><empty_output><jupyter_text>so first, it needs to be set (I'd guess GeoPandas will appropriately set from a database in the future). If we look it up in the database we'll see that it's NAD83 (North American Datum 1983), which has the [EPSG](www.epsg.org) (European Petroleum Survey Group) code of 4269.<jupyter_code># set the counties crs to 4269
il_counties.crs = {'init': u'epsg:4269'}

# print it out
print(il_counties.crs)
# let's check out the area calculated using Pandas with NAD83
il_counties['area_nad83'] = il_counties.geom.area
# view the first 5 records' aland and calculated area with NAD83:
il_counties.loc[:,('aland', 'area_nad83')].head()<jupyter_output><empty_output><jupyter_text>Clearly not the same. We can look for other projections at a number of websites, but a good one is [epsg.io](www.epsg.io). let's use the US National Atlas Equal Area projection (epsg=2163), which is a meters based equal area projection.<jupyter_code># transform aka re-project the data (use the "inplace=True" flag to perform the operation on this Geodataframe)
il_counties.to_crs(epsg=2163, inplace=True)

# print out the CRS to see it worked
print(il_counties.crs)
# and let's calculate the area with the new CRS
il_counties['area_2163'] = il_counties.geom.area

# and again check the head() of the data, with all 3 area columns:
il_counties.loc[:,('aland', 'area_nad83', 'area_2163')].head()
# let's check if those small differences are just because we're only looking at land area, not full county area:
il_counties['total_area'] = il_counties.aland + il_counties.awater

# and recheck areas against total:
il_counties.loc[:,('total_area', 'area_nad83', 'area_2163')].head()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.09045395465267414, 'fraction_numerical': 0.020218552929283203, 'mean_word_length': 4.876661951909477, 'pattern_counts': {'":': 0, '<': 55, '<?xml version=': 0, '>': 62, 'https://': 6, 'lorem ipsum': 0, 'www.': 3, 'xml': 1}, 'pii_count': 2, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '10', 'dir': 'jupyter-structured-clean-dedup', 'id': '305355', 'license': 'CC0-1.0', 'n_tokens_mistral': 6784, 'n_tokens_neox': 6271, 'n_words': 2512, 'path': 'notebooks/session_04/spatial_visualizations.ipynb', 'repo_name': 'Coleridge-Initiative/ada-2017-welfare', 'size': '20783'}