Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
# We build our kriging model:
ok = OrdinaryKriging(V, min_points=2, max_points=10, mode='exact')
# We calculate the hydraulic head on our regular grid,
# and we make the result in a good shape
hh_hat = ok.transform(xx.flatten(), yy.flatten()).reshape(xx.shape)
# We calculate the kriging error on our grid:
s2 = ok.sigma.reshape(xx.shape)
# We determine xmin, xmax, ymin and ymax from our dataset:
xv = gdf['x'].values
yv = gdf['y'].values
xmin, xmax = min(xv), max(xv)
ymin, ymax = min(yv), max(yv)
# We determine the resolution in x direction:
res_x = 100
# We determine the resolution in y direction
# based on res_x to make a regular grid:
res_y = int((ymax - ymin)*res_x/(xmax - xmin))
# We build the grid:
xx,yy = np.mgrid[xmin:xmax:complex(res_x),
ymin:ymax:complex(res_y)]
print('X coordinates:\n', xx[:5,:5].round(1))
print('Y coordinates:\n', yy[:5,:5].round(1))
fig, axs = plt.subplots(1, 2 ,figsize=(15, 10), sharey=True)
ax = axs[0]
# Contour fringes of the kriging process:
ctr_hh = ax.contourf(xx, yy, hh_hat,
range(150,200,5),
cmap = "viridis_r",
alpha = 0.5)
ax.set_title("Ordinary kriging estimation")
divider = make_axes_locatable(ax)
cax = divider.append_axes("right", size="5%", pad=0.05)
plt.colorbar(ctr_hh, cax=cax, label = 'Hydraulic head [m asl]')
# We add the location of points of our dataset
gdf.plot(ax = ax, c = "black", marker= '.', markersize = 2)
ax = axs[1]
# Contour fringes of the kriging error:
ctr_err = ax.contourf(xx, yy, s2,
range(0,12,2),
cmap = "plasma",
alpha = 0.5)
ax.set_title("Kriging error estimation")
# We add the location of points of our dataset
gdf.plot(ax = ax, c = "black", marker= '.', markersize = 2)
divider = make_axes_locatable(ax)
cax = divider.append_axes("right", size="5%", pad=0.05)
plt.colorbar(ctr_err, cax=cax, label = 'Error [m]')
plt.show()
df_ok = pd.DataFrame({'x': xx.flatten(),
'y': yy.flatten(),
'hh_k': hh_ma.flatten(),
's2': s2_ma.flatten()})
# We remove points without values (in the False part of our mask)
df_ok = df_ok.dropna()
# We transform our df into a geodataframe
export = gpd.GeoDataFrame(df_ok,
geometry=gpd.points_from_xy(df_ok.x, df_ok.y))
export.head()
from osgeo import gdal, osr
# we determine the pixel size in x and y directions.
dx = abs(xmax - xmin)/res_x
dy = abs(ymax - ymin)/res_y
# we determine paramaters associated to our image
# top left x
# w-e pixel resolution
# rotation, 0 if image is "north up"
# top left y
# n-s pixel resolution
params =(xmin - dx/2, dx, 0, ymax + dy/2, 0, -dy)
# the name of our geoTIFF
tif_name = "my_kriging_geotiff.tif"
# Create/Open the raster
output_raster = gdal.GetDriverByName('GTiff').Create(tif_name, res_x+1, res_y+1, 1 ,gdal.GDT_Float32)
# Specify its coordinates
output_raster.SetGeoTransform(params)
# Establish its coordinate encoding:
srs = osr.SpatialReference()
# Our projection system is specified:
srs.ImportFromEPSG(2154)
# Exports the coordinate system to the file
output_raster.SetProjection(srs.ExportToWkt())
# Writes my array to the raster after some transformation due to the resulting shape of the kriging:
output_raster.GetRasterBand(1).WriteArray(np.transpose(np.flip(hh_hat[::-1])))
output_raster.GetRasterBand(1).SetNoDataValue(0)
output_raster.FlushCache()
# Uncomment below if you want to preview the result:
#import rasterio
#src = rasterio.open(tif_name)
#plt.imshow(src.read(1), vmin = 150, vmax = 180)
V.fit_sigma = 'linear'
fig1 = V.plot(show=False)
fig1.update_layout(title='Linear Weights', width = 1000)
V.fit_sigma = 'sqrt'
fig2 = V.plot(show=False)
fig2.update_layout(title='Sqrt-decrease Weights', width = 1000)
fig1.show()
fig2.show()
# We build our map focusing on our site and specifying the zoom start
mymap = folium.Map(location=[gdf.lat.mean(), gdf.lon.mean()],
zoom_start=12,
show=True, control_scale = True)
# Add the geological map of France
url = 'http://geoservices.brgm.fr/geologie'
layer1 = 'GEOLOGIE'
folium.WmsTileLayer(url, layer1, attr = 'BRGM', name = 'Geological map', show=True).add_to(mymap)
# Add a feature group to add borehole with hydraulic head measurements
fg_gwt = folium.FeatureGroup(name = 'Hydraulic head', show = True)
mymap.add_child(fg_gwt)
colormap = cm.LinearColormap(colors=['orange', 'yellow', 'green', 'lightblue', 'blue'],
index=[gdf["hh"].quantile(0.1),
gdf["hh"].quantile(0.25),
gdf["hh"].quantile(0.5),
gdf["hh"].quantile(0.75),
gdf["hh"].quantile(0.9)],
vmin=gdf["hh"].quantile(0.1),
vmax=gdf["hh"].quantile(0.9))
# We add the caption of our colormap
colormap.caption = 'Hydraulic head [m asl]'
# We add all points of our dataset
for i, h in zip(gdf.index.values, gdf.hh.values):
dfi = folium.CircleMarker(
location = [gdf.at[i,'lat'], gdf.at[i,'lon']],
popup =str(round(h,2)) + " m asl",
radius= 7,
fill=True,
fill_opacity=0.7,
color = colormap(h),
fill_color = colormap(h))
dfi.add_to(fg_gwt)
mymap.add_child(colormap)
mymap.add_child(folium.map.LayerControl(collapsed=True))
url = "https://github.com/guiattard/PythonForGeosciences/raw/master/geostatistics-applied-to-hydrogeology-with-scikit-gstat/gw-body.zip"
file = "gw-body.zip"
# Command to donwload the file at the given url
r = requests.get(url)
# Then we open the file
open(file, 'wb').write(r.content)
# We extract the content of the .zip file
with zipfile.ZipFile(file, 'r') as unzip: unzip.extractall("./dat")
# we finally read the shapefile
area = gpd.read_file("./dat/gw-body.shp")
area.head()
url = "https://github.com/guiattard/PythonForGeosciences/raw/master/geostatistics-applied-to-hydrogeology-with-scikit-gstat/hydraulic-head-lyon-sample.zip"
file = "hydraulic-head-lyon-sample.zip"
# Command to donwload the file at the given url
r = requests.get(url)
# Then we open the file
open(file, 'wb').write(r.content)
# We extract the content of the .zip file
with zipfile.ZipFile(file, 'r') as unzip: unzip.extractall("./dat")
# we finally read the shapefile and make some cleaning
gdf = gpd.read_file("./dat/hydraulic-head-sample-lyon.shp")
# We rename the hydraulic name column by hh
gdf = gdf.rename(columns = {'hydraulic_' : "hh"})
gdf.head()
from skgstat import Variogram, OrdinaryKriging, plotting
# to have nicer plots, we switch to plotly backend
plotting.backend('plotly')
import pandas as pd
import geopandas as gpd
import matplotlib.pyplot as plt
from mpl_toolkits.axes_grid1 import make_axes_locatable
import numpy as np
# We import the folium library to make interactive maps
import folium
# We need to import branca.colormap to give pretty colors to our points according
# to groundwater table elevation
import branca.colormap as cm
# We need requests to get our dataset and zipfile to unzip our dataset
import requests
import zipfile
%matplotlib inline
pip install scikit-gstat
git clone https://github.com/mmaelicke/scikit-gstat.git
cd scikit-gstat
pip install -r requirements.txt
python setup.py install
apt install gdal-bin python-gdal python3-gdal
pip install fiona shapely pyproj
apt install python3-rtree
pip install git+git://github.com/geopandas/geopandas.git
pip install git+https://github.com/python-visualization/folium
pip install plotly
pip install rasterio
from scipy.interpolate import griddata
# We put the x/y coordinates of our dataset into an array
points = np.asanyarray(gdf[['x', 'y']])
# We put our hydraulic head values of our dataset into list of the same dimension
values = gdf['hh'].values
# We make the linear interpolation
HH_LINEAR = griddata(points, values, (xx, yy), method='linear')
# We make the cubic interpolation
HH_CUBIC = griddata(points, values, (xx, yy), method='cubic')
# Spatial join
gdf = gpd.sjoin(gdf, area, op = "within")
# PLot the results
fig, ax = plt.subplots(figsize=(10,8))
area.plot(ax = ax, alpha = 0.25)
# We add the location of borehole with water level in blue
gdf.plot(ax = ax, cmap = "viridis_r",
column = "hh",
markersize=10,
legend = True,
legend_kwds={'label': "hydraulic head [m a.s.l]"})
plt.show()
gdf = gdf[gdf.hh > 145]
# PLot the results
fig, ax = plt.subplots(figsize=(10,8))
area.plot(ax = ax, alpha = 0.25)
# We add the location of borehole with water level in blue
gdf.plot(ax = ax, cmap = "viridis_r",
column = "hh",
markersize=10,
legend = True,
legend_kwds={'label': "hydraulic head [m a.s.l]"})
plt.show()
from shapely.geometry import Point
# We initiale a mask with 0 everywhere with the similar shape as earlier
maskin = xx - xx
# In this mask, we assign True if the point is inside the groundwater body
# and we assign False if the point is outside the groundwater body
for i in range(maskin.shape[0]):
for j in range(maskin.shape[1]):
xi = xx[i][j]
yj = yy[i][j]
if Point(xi, yj).within(area["geometry"][0]):
maskin[i, j] = 0
else:
maskin[i, j] = 1
# We apply our mask to previous calculated arrays:
hh_ma = np.ma.masked_array(hh_hat, maskin)
s2_ma = np.ma.masked_array(s2, maskin)
fig, axs = plt.subplots(1,2, figsize=(15, 10), sharey=True)
ax =axs[0]
# We plot contour fringes based on our linear griddata interpolation
ctr_linear = ax.contourf(xx,
yy,
HH_LINEAR,
range(150,200,5),
cmap = "viridis_r",
alpha = 0.5)
# Some parameters to make a pretty fig:
divider = make_axes_locatable(ax)
cax = divider.append_axes("right", size="5%", pad=0.05)
plt.colorbar(ctr_linear, cax=cax, label = 'Hydraulic head [m asl]')
ax.set_title("Interpolation with the linear method")
# We add the location of points of our dataset
gdf.plot(ax = ax, c = "black", marker= '.', markersize = 2)
ax =axs[1]
# We plot contour fringes based on our cubic griddata interpolation
ctr_cubic = ax.contourf(xx,
yy,
HH_CUBIC,
range(150,200,5),
cmap = "viridis_r",
alpha = 0.5)
# Some parameters to make a pretty fig:
divider = make_axes_locatable(ax)
cax = divider.append_axes("right", size="5%", pad=0.05)
plt.colorbar(ctr_cubic, cax=cax, label = 'Hydraulic head [m asl]')
ax.set_title("Interpolation with the cubic method")
# We add the location of points of our dataset
gdf.plot(ax = ax, c = "black", marker= '.', markersize = 2)
plt.show()
# We plot the results:
fig, axs = plt.subplots(1, 2 ,figsize=(15, 10), sharey=True)
ax = axs[0]
# Contour fringes of the kriging process:
ctr_hh = ax.contourf(xx, yy, hh_ma, range(150, 200, 5),
cmap = "viridis_r",
alpha = 0.5)
ax.set_title("Ordinary kriging estimation")
divider = make_axes_locatable(ax)
cax = divider.append_axes("right", size="5%", pad=0.05)
plt.colorbar(ctr_hh, cax=cax, label = 'Hydraulic head [m asl]')
# We add the location of points of our dataset
gdf.plot(ax = ax, c = "black", marker= '.', markersize = 2)
ax = axs[1]
# Contour fringes of the kriging error:
ctr_err = ax.contourf(xx, yy, s2_ma,
range(0,12,2),
cmap = "plasma",
alpha = 0.5)
ax.set_title("Kriging error")
# We add the location of points of our dataset
gdf.plot(ax = ax, c = "black", marker= '.', markersize = 2)
divider = make_axes_locatable(ax)
cax = divider.append_axes("right", size="5%", pad=0.05)
plt.colorbar(ctr_err, cax=cax, label = 'Error [m]')
plt.show()
fig, axs = plt.subplots(1,3, figsize=(10,4))
ax = axs[0:2]
gdf[['hh', 'z']].hist(ax = ax, bins=25)
ax = axs[2]
ax.set_ylabel('hydrualic head [m]')
gdf[['hh', 'z']].boxplot(ax = ax)
plt.show()
import plotly.express as px
# create the scatter_mapbox
fig = px.scatter_mapbox(gdf,
lat='lat', lon='lon', zoom=12, height=450,
hover_data=['hh', 'depth'], color='hh'
)
# use a base-style that does not need a MAPBOX token
fig.update_layout(
mapbox_style='stamen-terrain',
mapbox=dict(pitch=45),
)
# get rid of margins
fig.update_layout(margin=dict(l=0, r=0, t=0, b=0))
# We plot the result
fig, ax = plt.subplots(figsize=(10,8))
# We add the location of borehole with hydraulic head
gdf.plot(ax = ax, cmap = "viridis_r",
column = "hh",
markersize=10,
legend = True,
legend_kwds={'label': "hydraulic head [m a.s.l]"})
plt.show()
export.to_file('my_kriging_export.shp')
V = Variogram(list(zip(gdf.x, gdf.y)),
gdf.hh,
normalize = False,
model = "spherical",
use_nugget = True,
n_lags=80,
maxlag=1700)
fig = V.plot()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment