Kriging with external drift for temperature interpolation using elevation... again

[codename5281]

Hello,

I am trying to perform the interpolation of an air temperature field using weather station measurements and a digital elevation model (DEM).

As presented in https://github.com/GeoStat-Framework/PyKrige/issues/155 i used the "external_Z" argument for the drift_terms option, but whether or not i use drift terms, the result remains the same and i don't understand why.

I acknowledged https://github.com/GeoStat-Framework/PyKrige/issues/235 and i am using v1.7.1.

Am i missing something ?

It may be interesting to note that the resolution of my DEM (90m) is finer than the resolution of the interpolation i want to perform (0.005° ~ 200m under these latitudes).

Below is the used code :

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()

Thank you for your package and thanks in advance for your help, Best, J.