1up noob: 1月 2018

地球統計学での内挿法の一つであるクリギング法では、データが距離と方向にどのような関係を持つかを表すバリオグラムを元に重み付けを行います。バリオグラムやクリギング法については下記の書籍を参考に、バリオグラムを計算するプログラムをPythonで書きました。

x,y,valueという並び順のcsvで与えたデータに対して、バリオグラム（等方性）を計算します。
バリオグラムを計算したのち、指定した理論バリオグラムの形にフィットする値を計算します。

実行環境

Python 3.6.1

	# coding:utf-8
	import numpy as np
	import matplotlib.pyplot as plt
	from scipy import stats
	from scipy import optimize as opt
	from scipy.spatial.distance import pdist, squareform

	# 元データを読み込み、配列に格納する
	source_arr = np.genfromtxt('source.csv', delimiter=',', skip_header=1)

	# セミバリオグラムを計算する
	def variogram_2d(xyv_array):
	xy_dist = squareform(pdist(xyv_array[:, 0:2], 'euclidean'))
	s_vario = squareform(pdist(xyv_array[:, 2:3], 'euclidean')**2 / 2)
	xy_dist[np.tri(xyv_array.shape[0], dtype=bool)] = np.nan
	s_vario[np.tri(xyv_array.shape[0], dtype=bool)] = np.nan
	return [xy_dist, s_vario]
	z_vario = variogram_2d(source_arr)

	# 経験バリオグラム(階級幅の中でのセミバリオグラムの平均値)を計算する
	def emp_variogram(z_vario, lag_h):
	bin_h = (z_vario[0][~np.isnan(z_vario[0])].flatten()/lag_h).astype(int)
	num_rank = int(np.max(z_vario[0][~np.isnan(z_vario[0])]) / lag_h) + 1
	mean_h, _, _ = stats.binned_statistic(bin_h, z_vario[0][~np.isnan(z_vario[0])].flatten(), statistic='mean', bins=num_rank)
	mean_vario, _, _ = stats.binned_statistic(bin_h, z_vario[1][~np.isnan(z_vario[0])].flatten(), statistic='mean', bins=num_rank)
	e_vario = np.stack([mean_h[~np.isnan(mean_h)], mean_vario[~np.isnan(mean_vario)]], axis=0)
	e_vario = np.delete(e_vario, np.where(e_vario[1] <= 0)[0], axis=1)
	#e_vario = np.round(e_vario, 5)
	return e_vario
	# lag_hには階級幅を設定する
	e_vario = emp_variogram(z_vario, 6)

	# 理論バリオグラムモデル
	def liner_model(x, a, b):
	return a + b * x
	def gaussian_model(x, a, b, c):
	return a + b * (1 - np.exp(-(x / c)**2))
	def exponential_model(x, a, b, c):
	return a + b * (1 - np.exp(-(x / c)))
	def spherical_model(x, a, b, c):
	cond = [x < c, x > c]
	func = [lambda x : a + (b / 2) * (3 * (x / c) - (x / c)**3), lambda x : a + b]
	return np.piecewise(x, cond, func)

	# モデルへのフィッティング関数
	def auto_fit(e_vario, fitting_range, selected_model):
	# フィッティングレンジまで標本バリオグラムを削る
	data = np.delete(e_vario, np.where(e_vario[0]>fitting_range)[0], axis=1)
	# 選択したモデルに対してナゲットが0以下にならない様にフィッティングさせる
	if selected_model == 0:
	param, _ = opt.curve_fit(liner_model, data[0], data[1], bounds=(0, np.inf))
	elif selected_model == 1:
	param, _ = opt.curve_fit(gaussian_model, data[0], data[1], [0, 0, fitting_range], bounds=(0, np.inf))
	elif selected_model == 2:
	param, _ = opt.curve_fit(exponential_model, data[0], data[1], [0, 0, fitting_range], bounds=(0, np.inf))
	elif selected_model == 3:
	param, _ = opt.curve_fit(spherical_model, data[0], data[1], [0, 0, fitting_range], bounds=(0, np.inf))
	param = np.insert(param, 0, [selected_model, fitting_range])
	return param
	fitting_range = 60
	selected_model = 2
	param = auto_fit(e_vario, fitting_range, selected_model)
	#np.savetxt("param.csv", param, fmt="%.10f", delimiter=',')

	# グラフ作成
	fig = plt.figure()
	ax = plt.subplot(111)
	ax.plot(e_vario[0], e_vario[1], 'o')
	xlim_arr = np.arange(0, np.max(e_vario[0]), lag_h)
	if (param[0] == 0):
	ax.plot(xlim_arr, liner_model(xlim_arr, param[2], param[3]), 'r-')
	print(param[2], param[3])
	elif (param[0] == 1):
	ax.plot(xlim_arr, gaussian_model(xlim_arr, param[2], param[3], param[4]), 'r-')
	print(xlim_arr, param[3], param[4])
	elif (param[0] == 2):
	ax.plot(xlim_arr, exponential_model(xlim_arr, param[2], param[3], param[4]), 'r-')
	print(param[2], param[3], param[4])
	elif (param[0] == 3):
	ax.plot(xlim_arr, spherical_model(xlim_arr, param[2], param[3], param[4]), 'r-')
	print(param[2], param[3], param[4])
	# グラフのタイトルの設定
	ax.set_title('Semivariogram')
	# 軸ラベルの設定
	ax.set_xlim([0, np.max(e_vario[0])])
	ax.set_ylim([0, np.max(e_vario[1])])
	ax.set_xlabel('Distance [m]')
	ax.set_ylabel('Semivariance')
	# グラフの縦横比を調整
	aspect = 0.8 * (ax.get_xlim()[1] - ax.get_xlim()[0]) / (ax.get_ylim()[1] - ax.get_ylim()[0])
	ax.set_aspect(aspect)
	# グラフの描画
	plt.show()

view raw vario.py hosted with ❤ by GitHub

デモファイルに対する実行結果

【関連記事】

1up noob

2018/01/08

【Python】バリオグラムを計算する