Support for Cross Variogram

[digital-idiot]

Support for cross variogram would be great.

[mmaelicke]

Hey there, Yes, that would be great. Unfortunately, it would also be a considerable amount of work, that I cannot afford right now. The main reason is, that the Variogram class would need a second input for the co-variable. One of the very early design decisions I made, was to require only the coordinates and one array of values from the user. Breaking this up would mean to change the class all over the place. Thus, from my point of view, it would be easier to introduce a new class for Cross variograms. The other reason is, that the next issue will most likely ask for CoKriging, which would double the amount of work ;)

Nevertheless, implementing cross variograms for just one specific application should not be too complicated. If I can help, I will happily try to give some hints how to do that. You are also very welcome to implement and contribute a CrossVariogram class yourself, and I will happily give you some guidance, how Variogram works and at which places changes would be necessary. Sorry but best,

Mirko

[digital-idiot]

I would very much like to contribute. As you mentioned, a few pointers from you would be help me jump start, otherwise I have to read a huge chunk of code very carefully. Thanks for your help.

[mmaelicke]

Hey, Thanks, sounds great! I am not so sure anymore if it is really that much work. But first, I have to say, that I am not too familiar with cross-variograms. I read the 1998 paper by Cressie about it, but have never applied it so far anywhere.

So my suggestion would be to introduce a new class called CrossVariogram, inheriting from Variogram as a first step. As soon as we see the differences side-by-side, we might back-implement it into the existing Variogram class.

As far as I get it, we can assume, that instead of having a N-dimensional array of coordinates and a one dimensional array of observations (of same length), we have two one dimensional arrays (let's call them v and w) for the same coordinates, right? Basically, what we need to do is to change the way how the Z(v) - Z(v +h) part is calculated. All we need is to change this part to Z(v) - Z(w) for each combination of v,w.

These are the pairwise differences, stored in an array defined here: https://github.com/mmaelicke/scikit-gstat/blob/0257d4ec21f5d54fb5f08092b48c964d4192be3a/skgstat/Variogram.py#L158

The structure of this array is the same as for the distance, which is the row-wise upper triangle of the distance matrix. For _diff, it is just not the * distance* in space, but the difference in value.

You can obtain the content using this property: https://github.com/mmaelicke/scikit-gstat/blob/0257d4ec21f5d54fb5f08092b48c964d4192be3a/skgstat/Variogram.py#L259

These should both be fine. The array is filled by this function:

https://github.com/mmaelicke/scikit-gstat/blob/0257d4ec21f5d54fb5f08092b48c964d4192be3a/skgstat/Variogram.py#L1019

which further uses this function to yield the correct indexes: https://github.com/mmaelicke/scikit-gstat/blob/0257d4ec21f5d54fb5f08092b48c964d4192be3a/skgstat/Variogram.py#L1040 (btw. something that I always wanted to implement more properly... =))

I I am not mistaken, the __vdiff_indexer is also fine, as the the indices of the value are the same as long as the two observation arrays are obtained at the same points, in the same order. If not, we would have to change that. Not sure if that makes sense at all, just thinking loud.

So back into _calc_diff, there are two important lines: https://github.com/mmaelicke/scikit-gstat/blob/0257d4ec21f5d54fb5f08092b48c964d4192be3a/skgstat/Variogram.py#L1031

here, we basically need a v and a w, where w needs to be implemented (which is a bit of work, as self.values is a property, not an attribute). With that change in place, https://github.com/mmaelicke/scikit-gstat/blob/0257d4ec21f5d54fb5f08092b48c964d4192be3a/skgstat/Variogram.py#L1037

would be changed to something like:

for t, k in zip(self.__vdiff_indexer(), range(len(self._diff))):
            self._diff[k] = np.abs(v[t[0]] - w[t[1]])

Finally, we need to implement two value arrays, that hold the two covariables. The way they are handled can be found in the property values: https://github.com/mmaelicke/scikit-gstat/blob/0257d4ec21f5d54fb5f08092b48c964d4192be3a/skgstat/Variogram.py#L236

and its setter: https://github.com/mmaelicke/scikit-gstat/blob/0257d4ec21f5d54fb5f08092b48c964d4192be3a/skgstat/Variogram.py#L255

That's what we need twice. Just a short note. I would not implement the variable and the covariable as a 2d-array, as that is the input format I use in SpaceTimeVariograms to pass time series data. We should use another notion for cross-variograms.

Finally, if I didn't miss anything, it should be possible to implement cross variograms by introducing a cross_variable attribute and at the few places above check if that attribute is set.

If variable and co-variable can be observed at different locations (not sure if that makes sense), we would have to keep the new class and re-implement the distance matrix, to calculate all distances across both sets of coordinates.

[mmaelicke]

Hey there, @digital-idiot , just out of curiosity, is there any progress on this issue? Best, Mirko

[digital-idiot]

Thanks @mmaelicke for finally implementing this. I have been very busy with other things and it skipped my attention.