Ordinary kriging on a moving window

[rth]

This should be merged after the pull request #4 (as right now this partially redundant with it).

Ordinary kriging on a moving window (see issue #2 ) is implemented , and will be enabled if not Null n_closest_points is passed to the routine,

  z, ss = OK.execute(..., n_closest_points=100)

vectorized backend is not supported in this case, mostly intentionaly, since potentially the A martix can have variable shape, and it is therefore not suitable for vectorization.

In terms of performance, say we take a grid of 100x100 and 100 closest points, bellow is a timing (with the loop backend) as a function of the kriging data size,

time to solution (sec)	n=100	n=200	n=400	n=800	n=1600	n=3200	n=6200
OK	0.63	0.798	1.156	3.76	15.8	67.49	289.6
OK moving window	3.63	3.56	3.68	3.72	4.23	5.32	12.03

so for small data sizes, this is actually slower. However as the data size, increases, OK on moving window execution time remains almost constant, and there is simply no other way to do kriging of say 10k points.
As to the results, they look reasonable enough, although this migh need more testing. Another thing is scipy.spatial.KDTree is used to get the nearest neighbours, which is faster then calculating all the distances between the mesh and the data points.

This also adds a cython backend for both OK, and OK with kriging window, all styles (masked, grid, points) are supported, that can be called with,

   z, ss = OK.execute(..., backend='C')

To give an idea of the speedup, assuming, that we take 50 closest points, the speed up between the Cython version and the C version is given below (the format is OK speed up / OK moving window speed up),

relative speed up	n=100	n=200	n=400
100x100	3.67 / 3.68	2.87 / 3.80	1.88 / 3.54
400x400	1.35/ 3.25	1.63 / 3.57	1.65 / 3.64

it is not that large, but there still might be some room for improvement.

[bsmurphy]

Awesome. Thanks for putting all this together, and sorry I haven't gotten back to it sooner. I'll get working on a PyPI release.

[dshean]

This looks like the solution I need, as I am working with (lots of) 3D points. Any plans to add this functionality to the 3D kriging options?

[bsmurphy]

Sure, why not. It needs to be done eventually. It shouldn't be too difficult to extend @rth's pure Python implementation to 3D. Might be a few weeks, but I'll work on implementing this capability. @rth - I've also been thinking recently about adding an 'adaptive' capability, so that if the spatial covariance isn't stationary (changes across the dataset) you can automatically re-estimate the variogram model on subsets of points. I haven't thought too much about this, but I think it should fit relatively easily into your moving window implementation.

[dshean]

Sounds great - thanks for considering. The adaptive capability you mention would also be useful.

[bsmurphy]

Sorry it took me a while to get to this. I've made a quick modification to the 3D ordinary kriging code to include the moving window option. I've put this update on the development branch, as I've realized there are a lot of changes I'd like to make and I'm going to keep tinkering around on this development version for a little while. The OK3D moving window change does seem to work, so give it a try. Extending to UK and UK3D is the next step.