Variograma Experimental em Python

A seguir será apresentado o código, pyscripts, do cálculo do variograma experimental. A finalidade desde "tutorial" é o desenvolvimento de pyscripts no Qgis e Plugins de Geoestátistica. Os dados de teste podem ser obtidos no seguinte link.

Para iniciantes em python podem ir no tutorial Instalação de Python e Pycharm, e para o entendimento de forma gráfica e complementar a implementação do variograma experimental, ir no tutorial Variograma Experimental en R.

import pandas as pd import numpy as np import math import matplotlib.pylab as plt # Importacion de datos data = pd.read_csv("data.csv", delimiter=";") data = np.array(data) # Parametros de tolerancias lagtol = 1.25 lag = 2.5 numlags = 12 # Parametros del vector direcional ndir = 2 azimuth = [45, 135] * np.repeat(math.pi / 180, ndir) tolerance = np.repeat(45 * math.pi / 180, ndir) bandwidth = np.repeat(500, ndir) # variable de interes y sus respectivas coordenadas coord = data[:, [0, 1]] variable = data[:, 2] # Definiendo el tamanho de las matrices npairs = np.zeros((ndir, numlags), dtype=int) mcorr = np.zeros((ndir, numlags), dtype=float) lags = lag + np.arange(0, 12) * lag for dir in range(ndir): if tolerance[dir] >= (math.pi / 2): coneh = -1 bandwidthsq = -1 else: coneh = bandwidth[dir] / math.tan(tolerance[dir]) bandwidthsq = coneh * math.tan(tolerance[dir]) bandwidthsq = bandwidthsq * bandwidthsq # Sistema de coordenadas polares y el sistema acimutal dirx = math.cos((math.pi / 2) - azimuth[dir]) diry = math.sin((math.pi / 2) - azimuth[dir]) nd = len(data) for i in range(0, nd): for j in range(i + 1, nd): # print(i, j) v = coord[i, :] - coord[j, :] if tolerance[dir] >= (math.pi / 2): bool = True else: h = dirx * v[0] + diry * v[1] d2 = sum(v ** 2) - h * h if h > coneh: bool = d2 < bandwidthsq else: dmax = h * math.tan(tolerance[dir]) bool = d2 < (dmax * dmax) if ~bool: continue d = math.sqrt(sum(v ** 2)) lagid = -1 for k in range(numlags): if (lags[k] - lagtol) <= d and (d < (lags[k] + lagtol)): lagid = k if lagid < 0: continue hpropi = variable[i] hpropj = variable[j] tpropi = variable[i] tpropj = variable[j] mcorr[dir, lagid] = mcorr[dir, lagid] + 0.5 * (hpropi - hpropj) * (tpropi - tpropj) npairs[dir, lagid] = npairs[dir, lagid] + 1 mcorr[dir, :] = mcorr[dir, :] / npairs[dir, :] plt.plot(lags, mcorr[0,:], '.r', lags, mcorr[1,:], '.b') plt.axis([0, 30, 0, max(mcorr[0,:])]) plt.show()

#python #variograma #geoestatistica