GeoStat-Framework / PyKrige

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

Bug: Problem when ordinary kriging with backend = 'C' #219

Closed JuHolland closed 2 years ago

JuHolland commented 2 years ago

Hi all,

I get very different results when doing Ordinary Kriging with the backend = 'loop' and backend = 'C':

x = [2, 7, 15]
y = [10,18,4]
z = [10, 40, 30]
krig = pykrige.ok.OrdinaryKriging(x,y,z, variogram_model = 'exponential', variogram_parameters = [1,5,0.6])

xgrid = np.arange(0.0,20.0)
ygrid = np.arange(0.0,20.0)

fig, (ax1, ax2) = plt.subplots(1, 2)
kriged_result_C = krig.execute(style = 'grid', xpoints = xgrid, ypoints = ygrid, backend = 'C')[0].data
kriged_result_loop = krig.execute(style = 'grid', xpoints = xgrid, ypoints = ygrid, backend = 'loop')[0].data

im1 = ax1.imshow(kriged_result_C, vmin=20, vmax=40)
ax1.set_title('backend = C')
im2 = ax2.imshow(kriged_result_loop, vmin=20, vmax=40)
ax2.set_title('backend = loop')
fig.colorbar(im2, ax = (ax1, ax2))

Output:

image

There seems to be a problem with the variogram sill when doing ordinary kriging with an exponantial variogram and backend = 'C'. At instantiation of the OrdinaryKriging class, the sill gets converted to a partial sill and is then kept as an instance variable with the other model parameters.

When the .execute function is called with backend = 'C', the parameters stored with the instance are used. However the c_exec_loop_moving_window funtion converts the sill to partial sill a second time:

https://github.com/GeoStat-Framework/PyKrige/blob/f5b39204b18a552f5f162a3c07caeec78c2b8520/pykrige/lib/variogram_models.pyx#L59-L66

MuellerSeb commented 2 years ago

Thanks for pointing this out.

Would you mind creating a PR to fix this?