GeoStat-Framework / PyKrige

Kriging Toolkit for Python
https://pykrige.readthedocs.io
BSD 3-Clause "New" or "Revised" License
741 stars 186 forks source link

No specified drift in loop-based calculation of Universal Kriging ? #173

Closed codename5281 closed 1 year ago

codename5281 commented 3 years ago

Hello and thanks for your tool, I am trying to perform a UK with specified drift on a rather large matrix, thus i activated the backend="loop" argument. However the returned output corresponds to a UK without drift. Is that a behavior to be expected ? Thanks in advance, Best, J.

MuellerSeb commented 3 years ago

I need to take a closer look at this. Sounds strange to me.

codename5281 commented 3 years ago

Ok thanks, let me know if you need more info about how this issue is happening for us.

MuellerSeb commented 3 years ago

@codename5281 : a minimal example to reproduce your problem would be good either way!

codename5281 commented 3 years ago

Sure. Here is the code i am using :

data_x=list(gt["x_absolute_conv"])
data_y=list(gt["y_absolute_conv"])

UK = UniversalKriging(
    data_x,
    data_y,
    list(gt["Argile (%)"]),
    variogram_model="spherical",
    verbose=True,
    enable_plotting=False,
    variogram_parameters={'sill': V.describe()["sill"],
                          'range': V.describe()["effective_range"],
                          'nugget': V.describe()["nugget"]},
    drift_terms = 'specified',
    specified_drift=[gt["soilgrids_argile"].to_numpy()]                   
)

coeff=1
gridx = xnew
gridy = ynew

pivot=downscaled.to_numpy()
z, ss = UK.execute("masked",gridx,gridy,specified_drift_arrays =[pivot],mask=(pivot==0))

kt.write_asc_grid(gridx, gridy, z, filename="output.asc")
plt.figure(num=None, figsize=(20, 15), dpi=150, facecolor='w', edgecolor='k')
plt.imshow(z, origin = "lower")
plt.colorbar()

plt.show()
plt.close()

This returns the following, which is correct with what is expected:

Adjusting data for anisotropy...
Initializing variogram model...
Using 'spherical' Variogram Model
Partial Sill: 2542728.6074524377
Full Sill: 2542728.60745248
Range: 1903784396.4456012
Nugget: 4.2317871705984184e-08 

Calculating statistics on variogram model fit...
Q1 = 2.1155914729822847
Q2 = 422.3268706579135
cR = 362.1010268831601 

Initializing drift terms...
Executing Universal Kriging...

image

However, expliciting the argument backend="vectorize" throws the following error : ValueError: Specified backend vectorize is not supported for 2D universal kriging.. Plus, setting vectorize="loop" returns the following, which corresponds to the output i obtain when performing UK without indicating drift terms. image

MuellerSeb commented 1 year ago

I am a bit puzzled.

  1. it should be backend="vectorized" (with a "d" at the end)
  2. where did you set vectorize="loop"? There is no option vectorize
  3. From the code I at least think that backend="loop" should respect the specified drift.