GeoStat-Framework / PyKrige

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

Find optimal parameters for UK with external drifts #196

Closed jlemond closed 2 years ago

jlemond commented 3 years ago

Hello,

I currently use your very useful kriging package PyKrige for Universal Kriging with 2 external drifts (gridded temperature and gridded elevation) as follow:

UK = UniversalKriging(
    lon,
    lat,
    temp,
    variogram_model=variogram,
    nlags=nlags,
    exact_values=True,
    verbose=True,
    enable_plotting=True,
    drift_terms=['specified'],
    specified_drift=[z,raw],
)

krigresult,var =UK.execute('grid',gridx,gridy,specified_drift_arrays=[dem_indo,raw_temp_fit.values])

where: lon= longitudes of stations, lat= latitudes of stations, temp= the temperature values of stations that I want to krige, z= the elevation of stations used as external drift, raw= the temperature values used as external drift, gridx= vector of longitudes, gridy= vector of latitudes , dem_indo= gridded elevation data , raw_temp_fit.values= gridding temperature data

I would like to use Kriging CV to search optimal parameters amongst different variograms and number of bins (nlags).

I’ve tried following code:

param_dict= {
    "method": ['universal'],
    "variogram_model": ["linear", "power", "gaussian", "spherical", "exponential"],
    "nlags": [6, 10, 15, 20,30],
    "drift_terms":['specified']
}

estimator = GridSearchCV(Krige(), param_dict, verbose=True, return_train_score=True)

X=np.array([lon,lat]); X=X.T
y=val    
estimator.fit(X=X, y=val)

As I don’t find the specified drift argument in the Krige() function, I wonder how to manage the drift values. Could you help me?

Thanks in advance, Julien

MuellerSeb commented 3 years ago

The Krige class doesn't support external drift kriging at the moment.

I extended the original class last year to be able to use universal kriging, but since external drift kriging needs a specified drift at the output locations, I didn't know how to properly pass them.

This was decided here: https://github.com/GeoStat-Framework/PyKrige/pull/158#issuecomment-659968767

We need further discussion about that internally.

MuellerSeb commented 3 years ago

For now you could use ext_drift_grid of the Krige class And use a grid defined with the external drift, where all values are extracted during calculation. But that only works for 2D.

See: https://geostat-framework.readthedocs.io/projects/pykrige/en/stable/generated/pykrige.rk.Krige.html#pykrige.rk.Krige