Spatial Interpolation¶
Spatial interpolation is used to predicts values for cells in a raster from a limited number of sample data points around it. We are studying streaming high frequency temperature data in Chicago retrieved from Array of Thing (AoT).
Kriging is a family of estimators used to interpolate spatial data. This family includes ordinary kriging, universal kriging, indicator kriging, co-kriging and others (Taken from Lefohn et al., 2005). The choice of which kriging to use depends on the characteristics of the data and the type of spatial model desired. The most commonly used method is ordinary kriging, which was selected for this study. Reference:
Lefohn, Allen S. ; Knudsen, H. Peter; and Shadwick, Douglas S. 2005. Using Ordinary Kriging to Estimate the Seasonal W126, and N100 24-h Concentrations for the Year 2000 and 2003. A.S.L. & Associates, 111 North Last Chance Gulch Suite 4A Helena , Montana 59601. contractor_2000_2003.pdf
Notebook Outline¶
Setup¶
Retrieve the study area of this interactive spatial interpolation jupyter notebook
1) setting the environment, import the library
#!pip install pykrige
import osmnx as ox
%matplotlib inline
ox.config(log_console=True, use_cache=True)
#ox.__version__
2) we choose the Chicago as study area, download the distrct using osmnx, and save the dataset as shapefile
# from some place name, create a GeoDataFrame containing the geometry of the place
city = ox.gdf_from_place('Chicago, IL')
print (city)
# save the retrieved data as a shapefile
ox.save_gdf_shapefile(city)
3) plot the Chicago city
fig, ax = ox.plot_shape(city)
Spatial Interpolation¶
The pykrige is a Kriging Toolkit for Python. The code supports 2D and 3D ordinary and universal kriging. Standard variogram models (linear, power, spherical, gaussian, exponential) are built in.
import numpy as np
import pandas as pd
import glob
from pykrige.ok import OrdinaryKriging
from pykrige.kriging_tools import write_asc_grid
import pykrige.kriging_tools as kt
import matplotlib.pyplot as plt
#from mpl_toolkits.basemap import Basemap
from matplotlib.colors import LinearSegmentedColormap
from matplotlib.patches import Path, PathPatch
Show the points available in Chicago with accurate data¶
Read sensor data in the CSV file, including the sensor ID, latitude, longitude and tempreture.
# uncomment to get from CSV
data = pd.read_csv(
'sensors.csv',
delim_whitespace=False, header=None,
names=["ID","Lat", "Lon", "Z"])
#data = pd.DataFrame(dd)
data.head(len(data))
2） Data processing part, if the tempreture is greater than 45, then set the data be 45
lons=np.array(data['Lon'])
lats=np.array(data['Lat'])
zdata=np.array(data['Z'])
print (zdata)
#If some data are greate than 50, then
for r in range(len(zdata)):
if zdata[r]>50:
zdata[r] = 45
print (zdata)
Use ordinary kriging to do the spatial interpolation¶
In order to run spatial interpolation, we should define the boundary for the Chicago. Get the bounday value from the shapefile.
import geopandas as gpd
Chicago_Boundary_Shapefile = './data/il-chicago/il-chicago.shp'
boundary = gpd.read_file(Chicago_Boundary_Shapefile)
# get the boundary of Chicago
xmin, ymin, xmax, ymax = boundary.total_bounds
Generate the grids for longitude and lantitude, the number of bins is 100¶
xmin = xmin-0.01
xmax = xmax+0.01
ymin = ymin-0.01
ymax = ymax+0.01
grid_lon = np.linspace(xmin, xmax, 100)
grid_lat = np.linspace(ymin, ymax, 100)
Run the OrdinaryKriging method, the variogram_model is gaussian¶
OK = OrdinaryKriging(lons, lats, zdata, variogram_model='gaussian', verbose=True, enable_plotting=False,nlags=20)
z1, ss1 = OK.execute('grid', grid_lon, grid_lat)
print (z1)
Plot the spatial interpolation result with ordinary kriging using 'gaussian' variogram model¶
Generate the result and the legend
xintrp, yintrp = np.meshgrid(grid_lon, grid_lat)
fig, ax = plt.subplots(figsize=(30,30))
#ax.scatter(lons, lats, s=len(lons), label='Input data')
boundarygeom = boundary.geometry
contour = plt.contourf(xintrp, yintrp, z1,len(z1),cmap=plt.cm.jet,alpha = 0.8)
plt.colorbar(contour)
boundary.plot(ax=ax, color='white', alpha = 0.2, linewidth=5.5, edgecolor='black', zorder = 5)
npts = len(lons)
plt.scatter(lons, lats,marker='o',c='b',s=npts)
#plt.xlim(xmin,xmax)
#plt.ylim(ymin,ymax)
plt.xticks(fontsize = 30, rotation=60)
plt.yticks(fontsize = 30)
#Tempreture
plt.title('Spatial interpolation from temperature with gaussian (%d points)' % npts,fontsize = 40)
plt.show()
#ax.plot(grid_lon, grid_lat, label='Predicted values')
Run the OrdinaryKriging method, the variogram_model is gaussian¶
OK = OrdinaryKriging(lons, lats, zdata, variogram_model='linear', verbose=True, enable_plotting=False,nlags=20)
z2, ss1 = OK.execute('grid', grid_lon, grid_lat)
#print (z2)
Plot the spatial interpolation result with ordinary kriging using 'linear' variogram model¶
Generate the result and the legend
xintrp, yintrp = np.meshgrid(grid_lon, grid_lat)
fig, ax = plt.subplots(figsize=(30,30))
#ax.scatter(lons, lats, s=len(lons), label='Input data')
boundarygeom = boundary.geometry
contour = plt.contourf(xintrp, yintrp, z2,len(z2),cmap=plt.cm.jet,alpha = 0.8)
plt.colorbar(contour)
boundary.plot(ax=ax, color='white', alpha = 0.2, linewidth=5.5, edgecolor='black', zorder = 5)
npts = len(lons)
plt.scatter(lons, lats,marker='o',c='b',s=npts)
#plt.xlim(xmin,xmax)
#plt.ylim(ymin,ymax)
plt.xticks(fontsize = 30, rotation=60)
plt.yticks(fontsize = 30)
#Tempreture
plt.title('Spatial interpolation from temperature with linear function (%d points)' % npts,fontsize = 40)
plt.show()
#ax.plot(grid_lon, grid_lat, label='Predicted values')
