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
# create folder to save to
!mkdir -p data/il-chicago/
# from some place name, create a GeoDataFrame containing the geometry of the place
city = ox.geocode_to_gdf('Chicago, IL')
# project to epsg:
print(city)
# save the retrieved data as a shapefile
city.to_file("./data/il-chicago/il-chicago.shp")
3) plot the Chicago city
ax = ox.project_gdf(city).plot(fc='gray', ec='none')
_ = ax.axis('off')
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¶
In this next several lines of codes, we divide the area of Chicago into multiple rasters by longitude and latitude. And the chicago area is divided into 100*100 subarea based on the longitude and latitude.
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¶
In spatial statistics the theoretical variogram is a function describing the degree of spatial dependence of a spatial random field or stochastic process.
And Ordinary Kriging is a very popular method for spatial interpolation.
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. The red area are places where temperature is high while the blue area are places where temperature is low.
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¶
In this case, we are using ordinary kriging using another variogram model which is linear instead of the gaussian variogram model that was used previously. And you may notice some differences between using these two different models by looking at the following plot.
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')
