pfg5-build-kriging.py

# 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)

pfg5-create-grid.py

# 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))

pfg5-display-kriging-and-error.py

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()

pfg5-export-geodf.py

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()

pfg5-export-geotiff.py

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)

pfg5-fit-sigma.py

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()

pfg5-folium1.py

# 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))

pfg5-get-gw-body.py

url = "https://github.com/guiattard/geoscience-notebooks/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()

pfg5-import-dataset.py

url = "https://github.com/guiattard/geoscience-notebooks/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()

pfg5-import-skgstat.py

from skgstat import Variogram, OrdinaryKriging, plotting

# to have nicer plots, we switch to plotly backend
plotting.backend('plotly')

pfg5-import.py

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

pfg5-install-option1

pip install scikit-gstat

pfg5-install-option2

git clone https://github.com/mmaelicke/scikit-gstat.git
cd scikit-gstat
pip install -r requirements.txt
python setup.py install

pfg5-install-other-lib

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

pfg5-interpolate-griddata.py

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')

pfg5-join-gw-body1.py

# 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()

pfg5-join-gw-body2.py

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()

pfg5-mask-outside-doi.py

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)

pfg5-plot-both-techniques.py

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()

pfg5-plot-final-result.py

# 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()

pfg5-plot-hist.py

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()

pfg5-plotly1.py

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))

pfg5-plt1.py

# 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()

pfg5-save-geodf.py

export.to_file('my_kriging_export.shp')

pfg5-variogam1.py

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()