Gis Notes I

Published:

murdock

Misc spatial operations (spatial overlay, zonal statistics) on the Murdock map. For reference, esp for tedious tasks like clipping a raster.

In [1]:
#%% 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 [2]:
# spatial libraries
import geopandas as gpd
import rasterio
import pysal as ps
from rasterstats import zonal_stats
from pysal.contrib.viz import mapping as maps
In [3]:
# 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 [4]:
murdock    = gpd.read_file(nunn_murdock)
boundaries = gpd.read_file(world_boundaries)
In [5]:
boundaries['REGION_UN'].value_counts()
Out[5]:
Africa                     57
Asia                       53
Americas                   51
Europe                     50
Oceania                    24
Seven seas (open ocean)     5
Antarctica                  1
Name: REGION_UN, dtype: int64
In [6]:
afr = boundaries.loc[boundaries['REGION_UN']=='Africa',:]
boundaries.shape
afr.shape
Out[6]:
(241, 72)
Out[6]:
(57, 72)

Plot Africa country boundaries + Murdock Atlas

In [7]:
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[7]:
<matplotlib.axes._subplots.AxesSubplot at 0x7f4b00bd0eb8>
Out[7]:
<matplotlib.axes._subplots.AxesSubplot at 0x7f4b00bd0eb8>

zoom to any country

In [8]:
afr.SOVEREIGNT.unique()
Out[8]:
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 [9]:
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 [10]:
murdock_map('Democratic Republic of the Congo')

Modern day light-levels

In [11]:
from rasterstats import zonal_stats
import rasterio.plot as rioplot
In [12]:
dmsp_1992 = '/home/alal/Desktop/data/spatial/Lights/DMSP/F101992.v4/F101992.v4b_web.stable_lights.avg_vis.tif'
In [13]:
source = rasterio.open(dmsp_1992)
bounds = (source.bounds.left, source.bounds.right, \
          source.bounds.bottom, source.bounds.top)
#%%
In [14]:
w,s,e,n = afr.total_bounds
#%%
f, ax = plt.subplots(1, figsize=(13,13))
ax.imshow(source.read(1), extent=bounds)
ax.set_xlim((w, e))
ax.set_ylim((s, n))
afr.plot(color='green', alpha=0.1, linewidth=1, edgecolor='w',ax=ax)
ax.set_axis_off()
Out[14]:
<matplotlib.image.AxesImage at 0x7f4b00b57780>
Out[14]:
(-25.341552734374943, 57.7919921875)
Out[14]:
(-46.96289062500003, 37.34038085937502)
Out[14]:
<matplotlib.axes._subplots.AxesSubplot at 0x7f4afc709e10>

Clip Raster

In [15]:
import fiona
from rasterio.tools.mask import mask
from shapely.geometry import box
In [16]:
w,s,e,n = afr.total_bounds
bbox = box(w,s,e,n)
geo = gpd.GeoDataFrame({'geometry': bbox},index=[0], crs=fiona.crs.from_epsg(4326))
In [17]:
f, ax = plt.subplots(1)
geo.plot(ax=ax)
afr.plot(edgecolor='k',ax=ax)
Out[17]:
<matplotlib.axes._subplots.AxesSubplot at 0x7f4afc6a3908>
Out[17]:
<matplotlib.axes._subplots.AxesSubplot at 0x7f4afc6a3908>
In [18]:
def getFeatures(gdf):
    """Function to parse features from GeoDataFrame in such a manner that rasterio wants them"""
    import json
    return [json.loads(gdf.to_json())['features'][0]['geometry']]
In [19]:
coords = getFeatures(geo)
coords
Out[19]:
[{'coordinates': [[[57.7919921875, -46.96289062500003],
    [57.7919921875, 37.34038085937502],
    [-25.341552734374943, 37.34038085937502],
    [-25.341552734374943, -46.96289062500003],
    [57.7919921875, -46.96289062500003]]],
  'type': 'Polygon'}]
In [20]:
source = rasterio.open(dmsp_1992)
In [21]:
out_image, out_transform = mask(raster=source, shapes=coords, crop=True)
out_meta = source.meta.copy()
    
out_meta.update({
    "driver": "GTiff",
     "height": out_image.shape[1],
     "width": out_image.shape[2],
     "transform": out_transform
})
In [22]:
clipped_file = 'dmsp_africa_clipped.tif'
with rasterio.open(clipped_file, 'w', **out_meta) as dest:
    dest.write(out_image)
In [23]:
new = rasterio.open(clipped_file)
bounds = (new.bounds.left, new.bounds.right, \
          new.bounds.bottom, new.bounds.top)
In [24]:
plt.imshow(new.read(1), cmap='viridis')
Out[24]:
<matplotlib.image.AxesImage at 0x7f4b00a9b668>
In [25]:
def dmsp_country_zoom(country):
    outline = afr.loc[afr['SOVEREIGNT']==country,:]
    w,s,e,n = outline.total_bounds
    f, ax = plt.subplots(1, figsize=(13,13))
    ax.imshow(new.read(1), extent=bounds)
    ax.set_xlim((w, e))
    ax.set_ylim((s, n))
    outline.plot(color='green', alpha=0.1, linewidth=1, edgecolor='w',ax=ax)
    ax.set_axis_off()
In [26]:
dmsp_country_zoom('South Africa')

Aggregate to polygons

In [27]:
lights_by_country = gpd.GeoDataFrame.from_features(zonal_stats(afr,clipped_file,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 [28]:
lights_by_homeland = gpd.GeoDataFrame.from_features(zonal_stats(murdock,clipped_file,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 [29]:
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 [30]:
f, ax = plt.subplots(1, figsize=(13, 13))
ax.set_title('Average lights observed in country')
lights_by_homeland.plot(column='lights_max', scheme='quantiles', k=10, 
                   cmap='viridis', legend=True, linewidth=1, ax=ax)
plt.axis('off');