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
You can find my data here:https://wetransfer.com/downloads/9c02e4fc1c2da765d5ee9137e6d7df4920220618071144/8f450e
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
def load_data():
df = pd.read_csv(r"File")
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=["functional"],
functional_drift=[func],
)
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
color_bar = basemap.colorbar(colormesh,location='bottom',pad="10%")
df = load_data()
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
conclusion: I don't know how to define the ["external_drift"] to do what I want (this is based on me being so inexperienced in coding in general).
I've been trying to solve these problems for 3 weeks now, but I'm really getting nowhere.
