Consider the height when kriging #Pykrige

[Weiss1hexe]

Hi I am creating an isoscape model in which I am plotting 2H and 18O of precipitation. Until now I have always used ordinary kriging for the interpolation, but I am somewhat dissatisfied with the results. Since the data I use for interpolation comes from the precipitation, I assume that the expected value cannot be constant. However, this is assumed with ordinary kriging. Since universal kriging adds a trend to the data, this seems more sensible to me in this case. Is it possible to use the height as a trend or does that make less sense?

my python code:

from traceback import print_tb
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

def load_data():
    df = pd.read_csv(r"C:/Users/Name/Desktop/Bachelor/Programm/Daten/Classified/O18/2012_9_O18.csv")
    return(df)

def get_data(df):
    #todo import data
    return {
        "lons": df['Longitude'],
        "lats": df['Latitude'],
        "values": df['O18'],
        "alts": df['Altitude'],
    }

def extend_data(data):
    # Copy data to the "left" and "right" to allow for interpolation to the edges of the map
    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"]]),
        "alts": np.concatenate([data["alts"], data["alts"], data["alts"]]),
    }

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) # dont extrapolate towards the poles
    }
    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',
        specified_drift = data["alts"], 
        #nlags=100, # TODO select plausible value
    )
    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) 
    basemap.drawcoastlines() 
    basemap.drawparallels(np.arange(-90.,120.,30.))
    basemap.drawmeridians(np.arange(0.,420.,60.))
    return figure, axes, basemap

def plot_mesh_data(interpolation, grid, basemap):
    colormesh = basemap.contourf(grid["x"], grid["y"],  interpolation,50, cmap='jet', ) #plot the data on the map. plt.cm.RdYlBu_r
    color_bar = basemap.colorbar(colormesh,location='bottom',pad="10%") #plot the colorbar on the map


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

@MuellerSeb

[MuellerSeb]

Hi there,

it is a common approach to use elevation as external drift for precipitation. From looking at your code, this is what you already did with specified_drift = data["alts"]. So this should be fine.

Sebastian