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()
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:
@MuellerSeb