# Install spatial packages
#!apt install gdal-bin python-gdal python3-gdal --quiet
#!apt install python3-rtree --quiet
#!pip install git+git://github.com/geopandas/geopandas.git --quiet
#!pip install descartes --quiet
# Install Pysal packages
#!pip install splot --quiet
#!pip install libpysal --quiet
#!pip install esda --quiet
import numpy as np
import pandas as pd
import geopandas as gpd
import matplotlib.pyplot as plt
from libpysal import weights
import esda
from esda.moran import Moran, Moran_Local
import splot
from splot.esda import moran_scatterplot, plot_moran, lisa_cluster
!wget https://raw.githubusercontent.com/quarcs-lab/data-quarcs/master/nepal/nepal.geojson
--2021-05-20 17:04:22-- https://raw.githubusercontent.com/quarcs-lab/data-quarcs/master/nepal/nepal.geojson
Resolving raw.githubusercontent.com (raw.githubusercontent.com)... 185.199.111.133, 185.199.108.133, 185.199.110.133, ...
Connecting to raw.githubusercontent.com (raw.githubusercontent.com)|185.199.111.133|:443... connected.
HTTP request sent, awaiting response... 200 OK
Length: 1474734 (1.4M) [text/plain]
Saving to: ‘nepal.geojson’
nepal.geojson 100%[===================>] 1.41M --.-KB/s in 0.01s
2021-05-20 17:04:23 (128 MB/s) - ‘nepal.geojson’ saved [1474734/1474734]
gdf = gpd.read_file("nepal.geojson")
#gdf.head()
gdf.plot("pcinc")
fig, ax = plt.subplots(figsize=(12,10))
gdf.plot(column="pcinc", scheme='Quantiles', k=3, cmap='viridis', legend=True, ax=ax)
plt.tight_layout()
plt.show()
fig, ax = plt.subplots(figsize=(12,10))
gdf.plot(column="pcinc", scheme='FisherJenks', k=3, cmap='viridis', legend=True, ax=ax)
plt.tight_layout()
ax.axis("off")
plt.show()
fig, ax = plt.subplots(figsize=(12,10))
gdf.plot(column="pcinc", scheme='BoxPlot', cmap='viridis', legend=True, ax=ax)
plt.tight_layout()
ax.axis("off")
plt.show()
fig, ax = plt.subplots(figsize=(12,10))
gdf.plot(column="pcinc", scheme='MaxP', cmap='viridis', legend=True, ax=ax)
plt.tight_layout()
ax.axis("off")
plt.show()
# Create W
w = weights.Queen.from_dataframe(gdf, idVariable="id" )
w.transform = "R"
#Spatial lag
gdf["w_pcinc"] = weights.lag_spatial(w, gdf["pcinc"])
gdf[["pcinc","w_pcinc"]]
# Global Moran
y = gdf["pcinc"]
moran = Moran(y, w)
moran.I
# P-value of the Global Moran
moran.p_sim
# Moran scatterplot
fig, ax = moran_scatterplot(moran, aspect_equal=True)
plt.savefig("scatter-moran.png")
plt.show()
# Local Moran
m_local = Moran_Local(y, w)
# Plot local Moran
fig, ax = moran_scatterplot(m_local, p=0.05)
ax.set_xlabel('Per-capita Income')
ax.set_ylabel('Spatial Lag of Per-capita Income')
plt.show()
# Plot LISA map
fig, ax = plt.subplots(figsize=(14,12))
lisa_cluster(m_local, gdf, p=0.01, figsize = (16,12),ax=ax)
plt.show()