Interface with GSTools

[MuellerSeb]

Hey there,

We are working on a project for geostatistical simulaton (https://github.com/GeoStat-Framework), where we, among other things, provide tools to generate spatial random fields based on variogram methods.

Now we wanted to provide conditioned random fields by incorporating data via Kriging and stumbled over your neat project! Your input information is totally similar to our approach, so I wrote a small interface for pykrige, which looks and works like this:

import numpy as np
from gstools import Gaussian
from pykrige.ok import OrdinaryKriging

data = np.array([[0.3, 1.2, 0.47],
                 [1.9, 0.6, 0.56],
                 [1.1, 3.2, 0.74],
                 [3.3, 4.4, 1.47],
                 [4.7, 3.8, 1.74]])

gridx = np.arange(0.0, 5.5, 0.1)
gridy = np.arange(0.0, 5.5, 0.1)

# a GSTools based covariance model
model = Gaussian(dim=2, len_scale=1, anis=0.2, angles=0.5, var=0.5, nugget=0.1)

OK1 = OrdinaryKriging(
    data[:, 0], data[:, 1], data[:, 2],
    **model.pykrige_kwargs,  # magic
)
z1, ss1 = OK1.execute('grid', gridx, gridy)

We provide a flexible Covariance-Class that contains all informations about the variogram and has a huge range of methods (like spectrum, spectral density, integral scale, variogram fitting, anisotropy, rotation... and so on) You can also create user defined models by just defining a correlation function (See here: https://github.com/GeoStat-Framework/GSTools#user-defined-covariance-models)

This class now also contains a property called pykrige_kwargs which can be used with your project as shown above. Only difference:
The anisotropy ratios are calculated inversely in both packages. So anis=0.2 in GSTools will result in anisotropy_scaling=1.0/0.2 in pykrige.
The angles are also calculated differently. Your orientation is the opposite of ours and you use deg instead of rad. So angles=0.5 in GSTools will result in anisotropy_angle=-np.rad2deg(0.5) in pykrige.

Same in 3D. pykrige_kwargs takes care of that depending on the dimension.

But still, the provided functionality overlaps perfectly! And I think, If you are still working on refactoring your code, It could be a cool thing working together!

What do you think about that? Maybe we can join forces! I'd love to have a chat about that.

Cheers,
Sebastian

[MuellerSeb]

We just updated our rotation behavior in GSTools (GeoStat-Framework/GSTools@0dfaf5f), so that the orientation is the same as in PyKrige.
Also we now provide simple and ordinary kriging in GSTools and recommend your package for the more fancier stuff.
The above mentioned Interface to PyKrige will be included in our next release GSTools 1.1.0. I hope you are OK with that.

[mjziebarth]

Hey Sebastian,

please excuse the late response. Personally, I am quite busy right now and seeing that there has been no response the last 2 weeks, I suppose it's similar for the others. Looking at the code you provided, it seems like an excellent integration and I can't see any reason against including it in your release. While I can't speak for the others, I suspect they feel similarly.

Regarding working together, do you have anything specific in mind? I have a feeling that the project has become kind of a side project for everyone involved so that fresh air is probably helpful :-)
What do you say, @bsmurphy, @rth, @basaks ?

Cheers,
Malte

[MuellerSeb]

Hey Malte,

thanks for the nice words. One thing I have in mind is, that your routines could be able to directly take a covariance model from gstools directly, since the presented interface with a provided kwargs dictionary is a bit of a hack. I think it would be nice, if one could provide an instance of our covariance model as a variogram model to your routines, like:

model = Gaussian(dim=2, len_scale=1, anis=0.2, angles=0.5, var=0.5, nugget=0.1)

OK1 = OrdinaryKriging(
    data[:, 0], data[:, 1], data[:, 2],
    variogram_model=model,
)

Then there could be a check, for example, at https://github.com/bsmurphy/PyKrige/blob/a4db3003b0b5688658c12faeb95a5a8b2b14b433/pykrige/ok.py#L216 if the given model is an instance of our model and everything from the input can be overwritten with information from the model and the variogram_model can be set internally as "custom". I think this could be quite convenient.

This should work for ordinary, and universal kriging in your case. The regression kriging in your package seems to require an internal model.

Cheers,
Sebastian

[MuellerSeb]

I created a pull request providing the above mentioned Interface: #125

[rth]

Thanks for your PR, it is quite nice!

And I think, If you are still working on refactoring your code, It could be a cool thing working together!

We have somewhat a lack of maintainer resources at the moment.

If you have an interest and the availability to improve PyKrige that would be very much appreciated! Personnally, I won't be able to contribute much beyound reviewing some PRs. I got into pykrige in 2015 for a side project, and haven't realy used in my work/outside of work projects since.

[MuellerSeb]

@rth and @bsmurphy : If there is interest in joining forces, one idea could be to integrate PyKrige into the GeoStat-Framework. Then we could care about PyKrige together and since the Framework is a GitHub organization, it would be easier for newbies to join the project and create a community maintained package base.

One step in that direction is of course #125. (Sorry for the long delay!)

Maybe we could have a conversation about that. What do you think?

Cheers, Sebastian

[MuellerSeb]

First step done with #125