Kriging model seems off

[simonericmoon]

I'm almost done with my Spatial Interpolation, however the result plot isn't what I would like to have. I followed the documentation of scikit-gstools about kriging: https://scikit-gstat.readthedocs.io/en/latest/tutorials/01_getting_started.html

the last plot is what I want to have.

I have a csv consisting of three columnes:
longitude, latiitude and subsidence values.

longitude values are around: 6.68; latitude values are around 50.76; subsidence rates range from +5 to -35

my current code looks like this:
jupyter code

here is the sample csv:

kerpencsv.csv

as you can see I don't get a colored kriging model and a plot which reflects a very high kriging error.
I'm unsure how to deal with this, yet alone why this is.
I'd be greatful for any help!

[mmaelicke]

Hey there,
Thanks for opening this issue. First, I am on vacation right now, so this will be a rather quickly written answer, but I hope it helps anyway.

There are a few things I'd like to comment:

• The Variogramclass is not really made for using WGS84 coords directly. The distance metric becomes quite pointless, then. There was a discussion about implementing the haversine distance in Haversine Distance for Geographic Coordinates #108 , including some examples but also notes (or cases) why (and when) one would not want to use it like this.
• Your variogram models look quite off. None of them is describing the spatial dependence quite well, especially on short distances. This can be a problem, as you have so many observation points at short distances and limit the Kriging algorithm to use 15 neighbors at most. Even with an increased limit, the close neighbors will have substantially larger weights, but are not well supported by the spatial models presented. Especially for the Gaussian and Cubic models I would expect some ridiculous estimations.
• Maybe most important from a technical point of view, I think this line is not doing what you are intending to do:

xx, yy = np.mgrid[6.67:6.69:6.72j, 50.86:50.88:50.92j]

Your xx has an xx.shape of (6,50), which is a quite odd shape for the interpolation grid. You want to end up with something like (500,500) for the data given. I would check out the docs for mgrid (https://numpy.org/doc/stable/reference/generated/numpy.mgrid.html), I guess the complex number has to be different.

Also, I get a warning, that for 275 of the 300 (6 * 50) target locations there are not enough neighbors. Thus, only 25 locations got interpolated.

Finally, concerning the Kriging error, your notebook does not show me any kriging error on the plot when I run it. So I do not really get the point. The associated number for the locations estimated might be high due to the fact that the target grid (xx, yy) is not created as intended.

I hope this helps a bit and feel free to drop another comment when you make some progress or anything of my answer is not clear (might take a time until I respond <- vacation).

Best,
Mirko

[redhog]

As a fairly correct solution to the coordinate problem, use pyproj to transform your coordinates from EPSG 4326 to 32633 (or 25833), then do both the calculation of the variogram, and any kriging of subsidence rates, in that projection. Make sure to keep the original coordinates around so you don't have to transform back and lose precision.

[mmaelicke]

@simonericmoon have you solved your issues and would you mind sharing your solution for reference?
Are still questions left?

[mmaelicke]

I guess this can be closed.