from pykrige.uk import UniversalKriging # selecting the dates we want to interpolate date = '1990-01-03' data_date = data.loc[data["date"]==date].dropna() # extent in which we want the interpolation to be performed lat_sw, lat_ne, lon_sw, lon_ne = (14.4, 14.9, -61.3, -60.8) # resolution of the interpolation, °/px res= 0.005 # loading the DEM with xarray dem = xr.open_rasterio("./output_SRTMGL3.tif") # clipping the DEM according to the extent dem = dem.sel(x=slice(lon_sw, lon_ne), y=slice(lat_ne, lat_sw)) # applying a buffer to the interpolation extent to mitigate the errors related to the # extent of the DEM being bigger than the interpolation extent buffer = 0.005 gridy = np.arange(lat_ne-buffer, lat_sw+buffer, -res) gridx = np.arange(lon_sw+buffer, lon_ne-buffer, res) # return coordinates of DEM as one array similar to the pos,field formalism dem_y, dem_x = np.meshgrid(dem.y.values, dem.x.values) dem_pos = np.array([dem_y.flatten(), dem_x.flatten()]) dem_field = dem.values.flatten() # defining a loop to try out different variogram models for model in ["power"]: # ["linear", "power", "gaussian", "spherical", "exponential", "hole-effect"]: UK = UniversalKriging( data_date["lon"], data_date["lat"], data_date["TN"], variogram_model=model, drift_terms=["external_Z"], external_drift=dem[0].T.values, external_drift_x=dem.x.values, external_drift_y=dem.y.values, verbose=True, pseudo_inv=True, enable_plotting=True, exact_values=True, anisotropy_scaling=1, ) z, ss = UK.execute("grid", gridx, gridy) plt.imshow(z, extent=(lon_sw, lon_ne, lat_sw, lat_ne)) plt.colorbar() plt.show()