0
$\begingroup$

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=["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
    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:

enter image description here

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:

-0.2 (O18)/100 m(altitude)

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?

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

0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy