Murdock

Published:

murdock
In [2]:
#%% regular scientific python imports
import os
import sys
import glob
import numpy as np
import pandas as pd
import statsmodels.api as sm
import statsmodels.formula.api as smf
# %matplotlib inline
# run for jupyter notebook
from IPython.core.interactiveshell import InteractiveShell
InteractiveShell.ast_node_interactivity = 'all'
#%%
/home/alal/anaconda3/envs/gds/lib/python3.5/site-packages/statsmodels/compat/pandas.py:56: FutureWarning: The pandas.core.datetools module is deprecated and will be removed in a future version. Please use the pandas.tseries module instead.
  from pandas.core import datetools
In [3]:
# spatial libraries
import matplotlib.pyplot as plt
import seaborn as sns
from seaborn import palplot
%matplotlib inline
import geopandas as gpd
import rasterio
import pysal as ps
from pysal.contrib.viz import mapping as maps
In [4]:
# national boundaries shapefile can be downloaded 
# from http://www.naturalearthdata.com/downloads/50m-cultural-vectors/50m-admin-0-countries-2/
world_boundaries = '/home/alal/Desktop/data/spatial/_country-map/ne_50m_admin_0_countries.shp'

# murdock map from nunn's digitised version - https://scholar.harvard.edu/nunn/pages/data-0
nunn_murdock = '/home/alal/Desktop/data/spatial/murdock_shapefile1_0/borders_tribes.shp'
In [5]:
murdock    = gpd.read_file(nunn_murdock)
boundaries = gpd.read_file(world_boundaries)
In [6]:
boundaries['REGION_UN'].value_counts()
Out[6]:
Africa                     57
Asia                       53
Americas                   51
Europe                     50
Oceania                    24
Seven seas (open ocean)     5
Antarctica                  1
Name: REGION_UN, dtype: int64
In [7]:
afr = boundaries.loc[boundaries['REGION_UN']=='Africa',:]
boundaries.shape
afr.shape
Out[7]:
(241, 72)
Out[7]:
(57, 72)

Plot Africa country boundaries + Murdock Atlas

In [8]:
f, ax = plt.subplots(1, figsize=(15, 15))
afr.plot(edgecolor='black',facecolor='green',linewidth=0.5,ax=ax)
murdock.plot(edgecolor='red',facecolor='none',linewidth=0.5,ax=ax)
ax.set_axis_off()
Out[8]:
<matplotlib.axes._subplots.AxesSubplot at 0x7f9dbbee6748>
Out[8]:
<matplotlib.axes._subplots.AxesSubplot at 0x7f9dbbee6748>

zoom to any country

In [9]:
afr.SOVEREIGNT.unique()
Out[9]:
array(['Angola', 'Burundi', 'Benin', 'Burkina Faso', 'Botswana',
       'Central African Republic', 'Ivory Coast', 'Cameroon',
       'Democratic Republic of the Congo', 'Republic of the Congo',
       'Comoros', 'Cabo Verde', 'Djibouti', 'Algeria', 'Egypt', 'Eritrea',
       'Ethiopia', 'Gabon', 'Ghana', 'Guinea', 'Gambia', 'Guinea-Bissau',
       'Equatorial Guinea', 'Kenya', 'Liberia', 'Libya', 'Lesotho',
       'Morocco', 'Madagascar', 'Mali', 'Mozambique', 'Mauritania',
       'Mauritius', 'Malawi', 'Namibia', 'Niger', 'Nigeria', 'Rwanda',
       'Western Sahara', 'Sudan', 'South Sudan', 'Senegal',
       'United Kingdom', 'Sierra Leone', 'Somaliland', 'Somalia',
       'Sao Tome and Principe', 'Swaziland', 'Seychelles', 'Chad', 'Togo',
       'Tunisia', 'United Republic of Tanzania', 'Uganda', 'South Africa',
       'Zambia', 'Zimbabwe'], dtype=object)
In [10]:
def murdock_map(country):
    # subset africa shapefile
    outline = afr.loc[afr['SOVEREIGNT']==country,:]
    # overlay (not real spatial merge because no common keys exist)
    tribes = gpd.overlay(outline, murdock, how='intersection')
    # plot
    f, ax = plt.subplots(1, figsize=(10, 10))
    tribes.plot(edgecolor='red',ax=ax)
    ax.set_axis_off()
    f.suptitle('Ethnic homelands in {0}'.format(country))
In [15]:
murdock_map('Democratic Republic of the Congo')

Modern day light-levels

In [12]:
from rasterstats import zonal_stats
import rasterio.plot as rioplot
In [14]:
dmsp_1992 = '/home/alal/tmp/gis-sandbox/dmsp/F101992.v4/F101992.v4b_web.stable_lights.avg_vis.tif'
In [15]:
lights_by_country = gpd.GeoDataFrame.from_features(zonal_stats(afr,dmsp_1992,prefix='lights_',geojson_out=True))
/home/alal/anaconda3/envs/gds/lib/python3.5/site-packages/rasterstats/main.py:145: FutureWarning: The value of this property will change in version 1.0. Please see https://github.com/mapbox/rasterio/issues/86 for details.
  with Raster(raster, affine, nodata, band) as rast:
/home/alal/anaconda3/envs/gds/lib/python3.5/site-packages/rasterstats/io.py:242: FutureWarning: GDAL-style transforms are deprecated and will not be supported in Rasterio 1.0.
  self.affine = guard_transform(self.src.transform)
/home/alal/anaconda3/envs/gds/lib/python3.5/site-packages/rasterstats/io.py:294: UserWarning: Setting nodata to -999; specify nodata explicitly
  warnings.warn("Setting nodata to -999; specify nodata explicitly")
/home/alal/anaconda3/envs/gds/lib/python3.5/site-packages/rasterstats/main.py:165: FutureWarning: Conversion of the second argument of issubdtype from `float` to `np.floating` is deprecated. In future, it will be treated as `np.float64 == np.dtype(float).type`.
  np.issubdtype(fsrc.array.dtype, float)
In [16]:
lights_by_homeland = gpd.GeoDataFrame.from_features(zonal_stats(murdock,dmsp_1992,prefix='lights_',geojson_out=True))
/home/alal/anaconda3/envs/gds/lib/python3.5/site-packages/rasterstats/main.py:145: FutureWarning: The value of this property will change in version 1.0. Please see https://github.com/mapbox/rasterio/issues/86 for details.
  with Raster(raster, affine, nodata, band) as rast:
/home/alal/anaconda3/envs/gds/lib/python3.5/site-packages/rasterstats/io.py:242: FutureWarning: GDAL-style transforms are deprecated and will not be supported in Rasterio 1.0.
  self.affine = guard_transform(self.src.transform)
/home/alal/anaconda3/envs/gds/lib/python3.5/site-packages/rasterstats/io.py:294: UserWarning: Setting nodata to -999; specify nodata explicitly
  warnings.warn("Setting nodata to -999; specify nodata explicitly")
/home/alal/anaconda3/envs/gds/lib/python3.5/site-packages/rasterstats/main.py:165: FutureWarning: Conversion of the second argument of issubdtype from `float` to `np.floating` is deprecated. In future, it will be treated as `np.float64 == np.dtype(float).type`.
  np.issubdtype(fsrc.array.dtype, float)
In [17]:
f, ax = plt.subplots(1, figsize=(11, 13))
ax.set_title('Average lights observed in country')
lights_by_country.plot(column='lights_mean', scheme='quantiles', k=10, 
                   cmap='viridis', legend=True, linewidth=1, ax=ax)
plt.axis('off');
In [18]:
f, ax = plt.subplots(1, figsize=(13, 13))
ax.set_title('Average lights observed in country')
lights_by_homeland.plot(column='lights_mean', scheme='quantiles', k=10, 
                   cmap='viridis', legend=True, linewidth=1, ax=ax)
plt.axis('off');

Homeland-level lights by country

In [21]:
# subset to country and plot lights
def murdock_map_lights(country):
    # subset africa shapefile
    outline = afr.loc[afr['SOVEREIGNT']==country,:]
    # rasterstats merge
    tribes = gpd.overlay(outline, murdock, how='intersection')
    lights_by_homeland = gpd.GeoDataFrame.from_features(zonal_stats(tribes,dmsp_1992,prefix='lights_',geojson_out=True))
    # plot
    f, ax = plt.subplots(1, figsize=(10, 10))
    lights_by_homeland.plot(column='lights_mean', scheme='quantiles', k=10, 
                           cmap='viridis', legend=True, linewidth=1, ax=ax)
    ax.set_axis_off()
    f.suptitle('Observed light in ethnic homelands in {0}'.format(country))
In [20]:
murdock_map_lights('Nigeria')
/home/alal/anaconda3/envs/gds/lib/python3.5/site-packages/rasterstats/main.py:145: FutureWarning: The value of this property will change in version 1.0. Please see https://github.com/mapbox/rasterio/issues/86 for details.
  with Raster(raster, affine, nodata, band) as rast:
/home/alal/anaconda3/envs/gds/lib/python3.5/site-packages/rasterstats/io.py:242: FutureWarning: GDAL-style transforms are deprecated and will not be supported in Rasterio 1.0.
  self.affine = guard_transform(self.src.transform)
/home/alal/anaconda3/envs/gds/lib/python3.5/site-packages/rasterstats/io.py:294: UserWarning: Setting nodata to -999; specify nodata explicitly
  warnings.warn("Setting nodata to -999; specify nodata explicitly")
/home/alal/anaconda3/envs/gds/lib/python3.5/site-packages/rasterstats/main.py:165: FutureWarning: Conversion of the second argument of issubdtype from `float` to `np.floating` is deprecated. In future, it will be treated as `np.float64 == np.dtype(float).type`.
  np.issubdtype(fsrc.array.dtype, float)
/home/alal/anaconda3/envs/gds/lib/python3.5/site-packages/pysal/esda/mapclassify.py:267: RuntimeWarning: invalid value encountered in greater
  binIds += (x > l) * (x <= r) * k
/home/alal/anaconda3/envs/gds/lib/python3.5/site-packages/pysal/esda/mapclassify.py:267: RuntimeWarning: invalid value encountered in less_equal
  binIds += (x > l) * (x <= r) * k
/home/alal/anaconda3/envs/gds/lib/python3.5/site-packages/numpy/lib/function_base.py:4033: RuntimeWarning: Invalid value encountered in median
  r = func(a, **kwargs)