Why OrdinaryKriging predict the same value for different points?

[lzeee]

I have some air quality data, which have latitude and longitude and a quality value.
Records are at the same time but some records lack the quality value. So I want to use kriging to predict.
Here is my code:

                OK = OrdinaryKriging(have_value_lon_list,
                                     have_value_lat_list,
                                     have_value_value_list,
                                     variogram_model='linear', 
                                     verbose=False,
                                     enable_plotting=False,
                                     coordinates_type='geographic')
                predict_value_list, ss = OK.execute('points',no_value_lon_list,no_value_lat_list)
                print(predict_value_list)

And my problem is that, for different points, why the predict values are very close or even the same?
Just like:
[16.97586206896553 16.97586206896553 16.97586206896553 16.97586206896553
16.97586206896553 16.97586206896553 16.97586206896553 16.97586206896553
16.975862068965533 16.975862068965533]

I have tried different variogram_function and it seems the 'spherical' is the best but just like:
[14.86176987024119 12.349359239110434 14.86176987024119 10.695659607721035
13.847852015338699 14.588090815442332 14.349646719219496
30.20493268831776 12.513013744442874 13.676770856568298]

Any help will be appreciated :)

[mjziebarth]

Dear @lzeee,
copying the answer from #113:
This may be related to the error mentioned in the pull request I referenced. In short, unfortunately the Kriging matrix does not properly use the great circle distance while the rest of the code does. I have not done any tests but intuitively, this should cause big errors at higher absolute latitudes.
May I ask where that region of yours is located and how big it is? Until the bug is fixed, I would suggest using 2d planar Kriging in some error-optimized projected coordinates. That should work well if that area you investigate is small. Have you tried that? If the problem persists, it is likely caused by something else.
Regards,
Malte

[lzeee]

Dear @lzeee,
copying the answer from #113:
This may be related to the error mentioned in the pull request I referenced. In short, unfortunately the Kriging matrix does not properly use the great circle distance while the rest of the code does. I have not done any tests but intuitively, this should cause big errors at higher absolute latitudes.
May I ask where that region of yours is located and how big it is? Until the bug is fixed, I would suggest using 2d planar Kriging in some error-optimized projected coordinates. That should work well if that area you investigate is small. Have you tried that? If the problem persists, it is likely caused by something else.
Regards,
Malte

thx for your reply :)
My data is some air quality data from the Los angeles, the locations can be seen directly in this picture.
The y-axis is latitude and x-axis is longitude.
The yellow point is the predict points and they all have the same value like 17.3090909.
Hope my information would help and I will try to project the point first.
Thx for your reply~~~

[lzeee]

Dear @lzeee,
copying the answer from #113:
This may be related to the error mentioned in the pull request I referenced. In short, unfortunately the Kriging matrix does not properly use the great circle distance while the rest of the code does. I have not done any tests but intuitively, this should cause big errors at higher absolute latitudes.
May I ask where that region of yours is located and how big it is? Until the bug is fixed, I would suggest using 2d planar Kriging in some error-optimized projected coordinates. That should work well if that area you investigate is small. Have you tried that? If the problem persists, it is likely caused by something else.
Regards,
Malte

Dear @mjziebarth

I just tried to project the lat/lon to UTM first using the python lib pyproj and turn to the 'euclidean' OrdinaryKriging. But the result seems to be the same :(

However, the results are not always a same value for different points. I have many data of different time. Some data will get different values for different points. But the 'same' circumstance appears quite frequently.

I wonder if I should do some preprocess to my data or something else before using the kriging or the results are just the correct ones?

[liuchengzimozigreat]

same thing happened to me. I trained the OK with 90 points in an area of 2000 km^2 and the subtraction of extreme lons and lats are all less than 1 degree. All the predicting values are same.
here is my codes:
data = np.array(heights) # heights is a 90×3 dataset, including 3 columns: lons, lats and heights
gridx = np.linspace(min(data[:,0]), max(data[:, 0]), 50)
gridy = np.linspace(min(data[:,1]), max(data[:, 1]), 25)
print(type(gridx))
OK = OrdinaryKriging(data[:, 0], data[:, 1], data[:, 2], variogram_model='linear', verbose=False, enable_plotting=False)
z, ss = OK.execute('grid', gridx, gridy)
print(z)
here is result ( I have tried to adjust number of predicting points, but nothing changed):
[[59.48436956521612 59.48436956521612 59.48436956521612 ...
59.48436956521612 59.48436956521612 59.48436956521612]
[59.48436956521612 59.48436956521612 59.48436956521612 ...
59.48436956521612 59.48436956521612 59.48436956521612]
[59.48436956521612 59.48436956521612 59.48436956521612 ...
59.48436956521612 59.48436956521612 59.48436956521612]
...
[59.48436956521612 59.48436956521612 59.48436956521612 ...
59.48436956521612 59.48436956521612 59.48436956521612]
[59.48436956521612 59.48436956521612 59.48436956521612 ...
59.48436956521612 59.48436956521612 59.48436956521612]
[59.48436956521612 59.48436956521612 59.48436956521612 ...
59.48436956521612 59.48436956521612 59.48436956521482]]

[Vanyary]

same thing happened to me. I trained the OK with 90 points in an area of 2000 km^2 and the subtraction of extreme lons and lats are all less than 1 degree. All the predicting values are same.
here is my codes:
data = np.array(heights) # heights is a 90×3 dataset, including 3 columns: lons, lats and heights
gridx = np.linspace(min(data[:,0]), max(data[:, 0]), 50)
gridy = np.linspace(min(data[:,1]), max(data[:, 1]), 25)
print(type(gridx))
OK = OrdinaryKriging(data[:, 0], data[:, 1], data[:, 2], variogram_model='linear', verbose=False, enable_plotting=False)
z, ss = OK.execute('grid', gridx, gridy)
print(z)
here is result ( I have tried to adjust number of predicting points, but nothing changed):
[[59.48436956521612 59.48436956521612 59.48436956521612 ...
59.48436956521612 59.48436956521612 59.48436956521612]
[59.48436956521612 59.48436956521612 59.48436956521612 ...
59.48436956521612 59.48436956521612 59.48436956521612]
[59.48436956521612 59.48436956521612 59.48436956521612 ...
59.48436956521612 59.48436956521612 59.48436956521612]
...
[59.48436956521612 59.48436956521612 59.48436956521612 ...
59.48436956521612 59.48436956521612 59.48436956521612]
[59.48436956521612 59.48436956521612 59.48436956521612 ...
59.48436956521612 59.48436956521612 59.48436956521612]
[59.48436956521612 59.48436956521612 59.48436956521612 ...
59.48436956521612 59.48436956521612 59.48436956521482]]

I have the same problem，did you find the solution to this problem？

[liuchengzimozigreat]

I finally adopted Inverse Distance Weighted (IDW) interpolation in my work. If you insist in using Kriging method in python, then you may try the related module in ArcGis, which can be transplant into python (I heard it).

[Vanyary]

I finally adopted Inverse Distance Weighted (IDW) interpolation in my work. If you insist in using Kriging method in python, then you may try the related module in ArcGis, which can be transplant into python (I heard it).
Thank you.

[KipCrossing]

I had this issue too. I came up with a crude fix with these assumptions for the variogram:

k_model: OrdinaryKriging = OrdinaryKriging(np.array(x), np.array(y), np.array(z),  
                                       variogram_model="gaussian", 
                                       variogram_parameters={'sill': (max(z)-min(z))*0.9, 
                                                                'range': 0.8*((max(x)-min(x))**2 + (max(y)-min(y))**2)**0.5, 
                                                                'nugget': (max(z)-min(z))*0.05},
                                        exact_values=False)

Note my x and y are converted using:

import utm
utm.from_latlon()

It's not the best, but better than the current.

[KipCrossing]

@MuellerSeb Have you thought about this?

Did my fix give you any clues?

[MuellerSeb]

This is most likely related to the estimated variogram parameters. When a single value is estimated everywhere, this is a sign, that a pure nugget model was estimated and that the standard variogram estimating routine could find any spatial correlation in the data.

You could either use predefined variogram parameters as suggested above, or you could use GSTools to futher analyze the estimated variogram and try find a better setup to estimated the empirical variogram, so the fitting afterwards leads to more useful parameters.

See here for interfacing with GSTools:

• https://geostat-framework.readthedocs.io/projects/gstools/en/stable/examples/05_kriging/02_pykrige_interface.html
• https://geostat-framework.readthedocs.io/projects/gstools/en/stable/examples/03_variogram/index.html