Krigging for isopleths

krigging_iso.py

import geopandas as gpd
import numpy as np
import matplotlib.pyplot as plt
from pykrige.ok import OrdinaryKriging
from sklearn.linear_model import LinearRegression
import rasterio
from rasterio.transform import from_origin

gdf = gpd.read_file("your_input_points.gpkg")
gdf = gdf.to_crs(2193)
gdf = gdf.dropna(subset=['thickness'])

x = gdf.geometry.x.values
y = gdf.geometry.y.values
log_thick = np.log10(gdf['thickness'].values + 1)

X = np.vstack([x, y]).T
trend_model = LinearRegression().fit(X, log_thick)
trend = trend_model.predict(X)
residuals = log_thick - trend

xi = np.linspace(x.min(), x.max(), 200)
yi = np.linspace(y.max(), y.min(), 200)
grid_x, grid_y = np.meshgrid(xi, yi)

OK = OrdinaryKriging(
    x, y, residuals,
    variogram_model='spherical',
    verbose=False,
    enable_plotting=False
)
res_kriged, ss = OK.execute('grid', xi, yi)

grid_coords = np.vstack([grid_x.ravel(), grid_y.ravel()]).T
grid_trend = trend_model.predict(grid_coords).reshape(grid_x.shape) #Predicted log(thickness) from X, Y trend
grid_log_thickness = grid_trend + res_kriged #Combined prediction in log space
grid_thickness = 10**grid_log_thickness - 1 #Final thickness prediction, removing log
grid_thickness = np.clip(grid_thickness, 0, None) #Final thickness prediction, removing negative values
grid_trend_thickness = 10**grid_trend - 1 #Trend-only thickness (no kriging), removing log
grid_trend_thickness = np.clip(grid_trend_thickness, 0, None)#Trend-only thickness (no kriging),  removing negative values

plt.figure(figsize=(10, 6))
trend_contour = plt.contourf(grid_x, grid_y, grid_trend_thickness, levels=20, cmap='viridis')
plt.colorbar(trend_contour, label='Trend Thickness')
plt.scatter(x, y, c='black', s=10, label='Data points')
plt.title('Trend Surface (Log-linear on X/Y)')
plt.xlabel('Easting (m)')
plt.ylabel('Northing (m)')
plt.axis('equal')
plt.legend()
plt.tight_layout()
plt.show()

plt.figure(figsize=(10, 6))
interp_contour = plt.contourf(grid_x, grid_y, grid_thickness, levels=20, cmap='plasma')
plt.colorbar(interp_contour, label='Interpolated Thickness')
plt.scatter(x, y, c='black', s=10, label='Data points')
plt.title('Full Interpolated Thickness (Trend + Residuals)')
plt.xlabel('Easting (m)')
plt.ylabel('Northing (m)')
plt.axis('equal')
plt.legend()
plt.tight_layout()
plt.show()

res_x = (xi[-1] - xi[0]) / (len(xi) - 1)
res_y = (yi[0] - yi[-1]) / (len(yi) - 1)

transform = from_origin(
    xi[0],
    yi[0],
    res_x,
    res_y
)

with rasterio.open(
    "trend_thickness.tif",
    "w",
    driver="GTiff",
    height=grid_trend_thickness.shape[0],
    width=grid_trend_thickness.shape[1],
    count=1,
    dtype=grid_trend_thickness.dtype,
    crs=gdf.crs,
    transform=transform,
) as dst:
    dst.write(grid_trend_thickness, 1)

with rasterio.open(
    "interpolated_thickness.tif",
    "w",
    driver="GTiff",
    height=grid_thickness.shape[0],
    width=grid_thickness.shape[1],
    count=1,
    dtype=grid_thickness.dtype,
    crs=gdf.crs,
    transform=transform,
) as dst:
    dst.write(grid_thickness, 1)

print("Exported GeoTIFFs: trend_thickness.tif and interpolated_thickness.tif")