Spatial interpolation methods are included here.

Using IDW

```python from math import sqrt import sys sys.path.append('/Users/xiao/programs/lib/gisalgs')

from geom.point import * from interpolation.idw import * from interpolation.prepare_interpolation_data import *

fname = '/Users/xiao/lib/gisalgs/data/necoldem.dat'

f = open(fname, 'r') Z = f.readlines() Z = [x.strip().split() for x in Z] Z = [ [ float(x[0]),float(x[1]),float(x[2])] for x in Z]

x = Point(337000, 4440911) N = 10

Z1 = prepare_interpolation_data(x, Z, N)[0]

print 'power=0.0:', IDW(Z1, 0) print 'power=0.5:', IDW(Z1, 0.5) print 'power=1.0:', IDW(Z1, 1.0) print 'power=1.5:', IDW(Z1, 1.5) print 'power=2.0:', IDW(Z1, 2.0) ```

Cross validation on different power values:

```python N = len(Z) numNeighbors = 10 mask = [True for i in range(N)] powers = [0, 0.5, 1, 1.5, 2, 2.5, 3, 3.5, 4, 4.5, 5] test_results = [] for i in range(N): mask[i] = False x = Point(Z[i][0], Z[i][1]) P = [ Z[j] for j in range(N) if mask[j] == True] P1 = prepare_interpolation_data(x, P, numNeighbors)[0] diff = [] for n in powers: zz = IDW(P1, n) diff.append(zz-Z[i][2]) test_results.append(diff) mask[i] = True

for i in range(len(powers)): rmse = sqrt(sum([r[i]**2 for r in test_results])/len(test_results)) print rmse, '[', powers[i], ']' ```

Kriging

```python import sys sys.path.append('/Users/xiao/lib/gisalgs')

from math import sqrt import numpy as np

from geom.point import * from indexing.kdtree1 import * from indexing.kdtree3 import *

from interpolation.fitsemivariance import * from interpolation.semivariance import * from interpolation.covariance import * from interpolation.read_data import * from interpolation.prepare_interpolation_data import * from interpolation.okriging import * from interpolation.skriging import *

fname = '/Users/Xiao/lib/gisalgs/data/necoldem.dat'

Z = read_data(fname)

hh = 50 lags = range(0, 3000, hh) gamma = semivar(Z, lags, hh) covariance = covar(Z, lags, hh)

Z1 = np.array(Z) semivariogram = fitsemivariogram(Z1, gamma, spherical)

x = Point(337000, 4440911) P1 = prepare_interpolation_data(x, Z1)[0] print okriging(np.array(P1), semivariogram)[0:2]

x = Point(Z1[0,0], Z1[0,1]) P1 = prepare_interpolation_data(x, Z1)[0] print okriging(np.array(P1), semivariogram)[0:2] ```

Features in the kriging functions are adopted from the geostatsmodels by cjohnson318.

Issues

fitsemivariance.py: nugget if statement

opened on 2018-04-29 14:05:23 by SimonMolinsky

When we calculate theoretical semivariance then we must set a nugget for each model. As I understand from the book (and the code itself) a nugget value should be set to zero (?) even if our first lag distance from empirical semivariance is not 0.0. But I do not understand the else statement in a condition when we set a nugget. Why do we choose the second value of lag if the first is equal to 0? Then we will start with nugget different than 0. Is this a correct way of setting our nugget value?

Code related to the question:

python if s[0][0] is not 0.0: # c0, nugget c0 = 0.0 else: c0 = s[0][1]

fitsemivariance.py: typo in spherical function's docstring

opened on 2018-04-29 13:49:12 by SimonMolinsky

Function spherical has a small typo in docstring. Variable c0 is described as a sill and variable c is described as a nugget but they should be described inversely:

python def spherical(h, c0, c, a): """ Input h: distance c0: nugget c: sill a: range """

Update idw.py

opened on 2017-10-09 14:38:01 by JialinLiOSU None