Note; No special knowledge of Pykrige is needed to answer the question, as I already mention examples in the question!

Hi I would like to use Universal Kriging in my code. For this I have data that is structured as follows:

      Latitude   Longitude  Altitude          H2        O18        Date    Year  month       dates  a_diffO       O18a
0    45.320000  -75.670000     114.0  -77.500000 -11.110000  2004-09-15  2004.0    9.0  2004-09-15   -0.228 -10.882000
1    48.828100    9.200000     314.0  -31.350000  -4.880000  2004-09-15  2004.0    9.0  2004-09-15   -0.628  -4.252000
2    51.930000  -10.250000       9.0  -18.800000  -3.160000  2004-09-15  2004.0    9.0  2004-09-15   -0.018  -3.142000
3    48.248611   16.356389     198.0  -45.000000  -6.920000  2004-09-15  2004.0    9.0  2004-09-15   -0.396  -6.524000
4    50.338100    7.600000      85.0  -19.200000  -3.190000  2004-09-15  2004.0    9.0  2004-09-15   -0.170  -3.020000



I want to interpolate the data (Latitude, Longitude, Altitude and O18) with Universal Kriging and use the height as a drift function.

So far I have programmed this here but I am not getting anywhere, e.g. I don't know how to effectively use the height as a drift function and the information from the Pykrige documentation is of limited help:

from traceback import print_tb
from typing_extensions import Self
import numpy as np
from pykrige.uk import UniversalKriging
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.patches import Path, PathPatch
import pandas as pd
from osgeo import gdal

return(df)

def get_data(df):
return {
"lons": df['Longitude'],
"lats": df['Latitude'],
"values": df['O18'],
}

def extend_data(data):
return {
"lons": np.concatenate([np.array([lon-360 for lon in data["lons"]]), data["lons"], np.array([lon+360 for lon in data["lons"]])]),
"lats": np.concatenate([data["lats"], data["lats"], data["lats"]]),
"values":  np.concatenate([data["values"], data["values"], data["values"]]),
}

def generate_grid(data, basemap, delta=1):
grid = {
'lon': np.arange(-180, 180, delta),
'lat': np.arange(np.amin(data["lats"]), np.amax(data["lats"]), delta)
}
grid["x"], grid["y"] = np.meshgrid(grid["lon"], grid["lat"])
grid["x"], grid["y"] = basemap(grid["x"], grid["y"])
return grid

def interpolate(data, grid):
uk = UniversalKriging(
data["lons"],
data["lats"],
data["values"],
variogram_model='exponential',
verbose=True,
drift_terms=["point_log"],
)
return uk.execute("grid", grid["lon"], grid["lat"])

def prepare_map_plot():
figure, axes = plt.subplots(figsize=(10,10))
basemap = Basemap(projection='robin', lon_0=0, lat_0=0, resolution='h',area_thresh=1000,ax=axes)
return figure, axes, basemap

def plot_mesh_data(interpolation, grid, basemap):
colormesh = basemap.contourf(grid["x"], grid["y"], interpolation,32, cmap='RdBu_r') #plot the data on the map. plt.cm.RdYlBu_r

base_data = get_data(df)
figure, axes, basemap = prepare_map_plot()
grid = generate_grid(base_data, basemap, 90)
extended_data = extend_data(base_data)
interpolation, interpolation_error = interpolate(extended_data, grid)
plot_mesh_data(interpolation, grid,basemap)
plt.show()



I now only use universal kriging and create these images:

I get the expected error: ValueError: Must specify location(s) and strength(s) of point drift terms.

I just know that I have to create a grid with the height, but I don't know how and I don't know how to make the drift dependent on the altitude. The altitude formula is:

where 100 m represents 100 m hight difference.

The interesting thing is that there is this website with examples: however, I am too inexperienced in coding to understand the examples and to transfer them to my example: https://python.hotexamples.com/examples/core/-/calc_cR/python-calc_cr-function-examples.html

I've been trying to solve these problems for 3 weeks now, but I'm really getting nowhere.

The question, how can I create an external drift for the UK interpolation that works as a function of altitude?

• I think that the problem is that the other region (somehow) is within the convex hull of the data and gets interpolated. May 14 at 17:18
• After your edit I don't see a question in your post. Jun 18 at 13:54
• @nicoguaro Ok I fixed it ;) Jun 18 at 14:12
• @nicoguaro If you could help me, it would bring important progress in my project. Jun 18 at 14:14