#!/usr/bin/env python
# coding: utf-8

# # Reading LSOAs through the ONS API
# 
# - [Dani Arribas-Bel](http://darribas.org/) [[`@darribas`](https://twitter.com/darribas)]

# In[1]:


get_ipython().run_line_magic('matplotlib', 'inline')

import geopandas
import contextily
import requests
from io import BytesIO


# In this short tutorial, we will explore how to pull down data from the ONS' Open Geography Portal. In particular, we will be working with 2011 LSOAs:

# In[2]:


from IPython.display import IFrame

site = ('http://geoportal.statistics.gov.uk/datasets/'\
        'lower-layer-super-output-areas-december-2011-'\
        'generalised-clipped-boundaries-in-england-and-wales/geoservice')
IFrame(site, 600, 400)


# To read off the ONS API, we can use the base URL provided [here](http://geoportal.statistics.gov.uk/datasets/lower-layer-super-output-areas-december-2011-generalised-clipped-boundaries-in-england-and-wales/geoservice):

# In[3]:


base_url = ("https://ons-inspire.esriuk.com/arcgis/rest/services/"\
            "Census_Boundaries/"\
            "Lower_Super_Output_Areas_December_2011_Boundaries/"\
            "MapServer/2/query")


# It is a standard `REST` API, which means we can pass the options as a set of parameters in a dictionary. 
# 
# For example, let us ask for all the fields of the LSOAs with an area larger than 282,152,900 $m^2$, returned in LonLat (`EPSG:4326`) and in the `GeoJSON` format:

# In[4]:


params = {
    'outFields': '*',
    'where': 'st_area(shape) > 282152900',
    'outSR': 4326,
    'f': 'geojson'
}


# **NOTE**: you can get all the details on which arguments are allowed [here](https://developers.arcgis.com/rest/services-reference/query-feature-service-layer-.htm).

# Now, we can pass the base URL with the query parameters to `requests` and let it do the hard work for us:

# In[5]:


r = requests.get(base_url, params=params)


# At this point we have already downloaded the data (you can print it raw if you want calling `r.content`).
# 
# However, this is only half the job. We also need to turn that into a data structure that allows us to work with it in Python. Since the payload is returned as a `GeoJSON`, you can pass that to the file reader in `geopandas` tricking it to believe it's a file instead of a string with `BytesIO`:

# In[6]:


db = geopandas.read_file(BytesIO(r.content))


# This operation builds a `GeoDataFrame`, with which we can work with in Python. for example, we can make a map that combines the polygons returned with a basemap pulled in from `contextily`:

# In[7]:


ax = db.to_crs(epsg=3857)\
       .plot(alpha=0.5, color='red',
             figsize=(12, 12))
contextily.add_basemap(ax);


# Now, just for fun, let us also demonstrate how you can make a slightly more sophisticated query. Exploring the table, we find that the `lsoa11nm` contains the name of the Local Authority District (LAD) each LSOA belongs to:

# In[8]:


db.head()


# We can use this then to query polygons from a single LAD, such as Liverpool (note this uses standard SQL):

# In[9]:


params = {
    'where': "lsoa11nm LIKE '%Liverpool%'",
    'outSR': 4326,
    'f': 'geojson'
}


# And we can now pull the data and turn it into a `GeoDataFrame` in the same way as above:

# In[10]:


r = requests.get(base_url, params=params)

db = geopandas.read_file(BytesIO(r.content))


# Let's check the results:

# In[11]:


ax = db.to_crs(epsg=3857)\
       .plot(alpha=0.5, color='red',
             figsize=(12, 12))
contextily.add_basemap(ax, 
                       url=contextily.sources.ST_TONER_BACKGROUND);


# 🎉

# ## Connecting ONS geometries with IMD

# In[12]:


import pandas

imd_url = ('https://assets.publishing.service.gov.uk'
           '/government/uploads/system/uploads/attachment_data'
           '/file/467764/File_1_ID_2015_Index_of_Multiple_Deprivation.xlsx')
imd = pandas.read_excel(imd_url, sheet_name='IMD 2015')
imd.head()


# In[13]:


db_imd = db.join(imd.set_index('LSOA code (2011)')\
                    .iloc[:, -1:]\
                    .rename(columns={('Index of Multiple Deprivation (IMD) Decile '\
                                      '(where 1 is most deprived 10% of LSOAs)'): 'score'}),
                 on='lsoa11cd')
db_imd.info()


# In[14]:


ax = db_imd.to_crs(epsg=3857)\
       .plot(column='score', scheme='quantiles',
             alpha=0.75, figsize=(12, 12))
contextily.add_basemap(ax, 
                       url=contextily.sources.ST_TONER_BACKGROUND);