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
Retrieve the study area of this interactive spatial interpolation jupyter notebook
# 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)
0 POLYGON ((-87.94010 42.00093, -87.94003 41.998...
place_name bbox_north bbox_south \
0 Chicago, Cook County, Illinois, United States ... 42.02304 41.644531
bbox_east bbox_west
0 -87.524081 -87.940101
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
'sensors.csv',
delim_whitespace=False, header=None,
names=["ID","Lat", "Lon", "Z"])
#data = pd.DataFrame(dd)
data.head(len(data))
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)
[ 28.24 19.83 26.17 45.7 35.07 36.47 22.45 125.01 19.82 26.21
22.04 20.2 20.5 42.15 21.67 25.14 125.01 21.16 19.5 21.61
32.03 28.3 20.11 40.6 42.8 31.46 21.35 21.99 21.62 20.88
20.55 20.41]
[28.24 19.83 26.17 45.7 35.07 36.47 22.45 45. 19.82 26.21 22.04 20.2
20.5 42.15 21.67 25.14 45. 21.16 19.5 21.61 32.03 28.3 20.11 40.6
42.8 31.46 21.35 21.99 21.62 20.88 20.55 20.41]
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
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)
Adjusting data for anisotropy...
Initializing variogram model...
Coordinates type: 'euclidean'
Using 'gaussian' Variogram Model
Partial Sill: 23.13242836937993
Full Sill: 93.76725570798001
Range: 0.3088207463280863
Nugget: 70.63482733860009
Calculating statistics on variogram model fit...
Executing Ordinary Kriging...
[[27.92963843313738 27.8858039970517 27.840219117484956 ...
31.386443563996654 31.40455626777043 31.415938741488716]
[27.905896385139272 27.860683189061522 27.81366895016103 ...
31.44165040499307 31.460519582840632 31.47249374835975]
[27.88211218938647 27.835521908251973 27.787080234206595 ...
31.49439048187477 31.514039150198606 31.52663138120343]
[29.130760923267935 29.155208360026453 29.180115169970442 ...
29.049409504791104 29.041974883402496 29.034401750480836]
[29.135870577895176 29.160612631119065 29.185827333210476 ...
29.039409247482 29.03163868914854 29.023773683305]
[29.13975307065343 29.164724932520965 29.190180905287033 ...
29.029451254242556 29.021385177488398 29.0132656101344]]
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')
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)
Adjusting data for anisotropy...
Initializing variogram model...
Coordinates type: 'euclidean'
Using 'linear' Variogram Model
Slope: 78.08683547092697
Nugget: 69.86456553372653
Calculating statistics on variogram model fit...
Executing Ordinary Kriging...
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')
