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 with exact_values when ordinary kriging with backend = 'C' #255

Closed JuHolland closed 2 years ago

JuHolland commented 2 years ago

Hello,

It seems that the exact_values boolean parameter is inverted when doing Ordinary Kriging with the backend = 'C'.

When setting exact_values = True:

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],
                                  exact_values = True)

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

init = np.full((20,20), np.nan)
init[y,x] = z

fig, (ax0, ax1, ax2) = plt.subplots(1, 3, figsize = (15,5))
kriged_result_C = krig.execute(style = 'grid', xpoints = xgrid, ypoints = ygrid, backend = 'C')[0].data
kriged_result_C = np.asarray(kriged_result_C)
kriged_result_loop = krig.execute(style = 'grid', xpoints = xgrid, ypoints = ygrid, backend = 'loop')[0].data
kriged_result_loop = np.asarray(kriged_result_loop)

im0 = ax0.imshow(init, vmin=20, vmax=40)
ax0.set_title('No kriging')
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 = [ax0,ax1, ax2])

Output:

image

And when setting exact_values = False:

image

MuellerSeb commented 2 years ago

Good catch!

Problem is here: https://github.com/GeoStat-Framework/PyKrige/blob/4b2ce01d1fd4b769f640b912a69ed2b513d6e1e4/src/pykrige/lib/cok.pyx#L69 The not is wrong.

And here, the condition check is even missing: https://github.com/GeoStat-Framework/PyKrige/blob/4b2ce01d1fd4b769f640b912a69ed2b513d6e1e4/src/pykrige/lib/cok.pyx#L163

MuellerSeb commented 2 years ago

Bug was introduced here: https://github.com/GeoStat-Framework/PyKrige/commit/ffbd472996a9c1da10e5d78423f848555338efaf

PR #256 should fix this.