Merge pull request #78 from GeoStat-Framework/vario_fit_update

Variogram fitting update

.travis.yml

@@ -1,3 +1,7 @@
+# do not build pull request from base repo twice (only build PRs of forks)
+if: type != pull_request OR fork = true
+
+# language is python
 language: python
 python: 3.8
 

examples/02_cov_model/00_intro.py

@@ -18,7 +18,7 @@
 # use CovModel as the base-class
 class Gau(gs.CovModel):
     def cor(self, h):
-        return np.exp(-h ** 2)
+        return np.exp(-(h ** 2))
 
 
 ###############################################################################

examples/02_cov_model/05_additional_para.py

@@ -19,7 +19,7 @@ def default_opt_arg(self):
         return {"alpha": 1.5}
 
     def cor(self, h):
-        return np.exp(-h ** self.alpha)
+        return np.exp(-(h ** self.alpha))
 
 
 ###############################################################################

examples/02_cov_model/06_fitting_para_ranges.py

@@ -15,7 +15,7 @@ def default_opt_arg(self):
         return {"alpha": 1.5}
 
     def cor(self, h):
-        return np.exp(-h ** self.alpha)
+        return np.exp(-(h ** self.alpha))
 
 
 # Exemplary variogram data (e.g. estimated from field observations)

examples/03_variogram/01_variogram_estimation.py

@@ -230,7 +230,7 @@ def generate_transmissivity():
 # dashed lines.
 
 plt.figure()  # new figure
-line, = plt.plot(bin_center, gamma, label="estimated variogram (isotropic)")
+(line,) = plt.plot(bin_center, gamma, label="estimated variogram (isotropic)")
 plt.plot(
     bin_center,
     fit_model.variogram(bin_center),
@@ -239,7 +239,7 @@ def generate_transmissivity():
     label="exp. variogram (isotropic)",
 )
 
-line, = plt.plot(x_plot, gamma_x[:21], label="estimated variogram in x-dir")
+(line,) = plt.plot(x_plot, gamma_x[:21], label="estimated variogram in x-dir")
 plt.plot(
     x_plot,
     fit_model_x.variogram(x_plot),
@@ -248,7 +248,7 @@ def generate_transmissivity():
     label="exp. variogram in x-dir",
 )
 
-line, = plt.plot(y_plot, gamma_y[:21], label="estimated variogram in y-dir")
+(line,) = plt.plot(y_plot, gamma_y[:21], label="estimated variogram in y-dir")
 plt.plot(
     y_plot,
     fit_model_y.variogram(y_plot),

examples/03_variogram/02_find_best_model.py

@@ -0,0 +1,64 @@
+"""
+Finding the best fitting variogram model
+----------------------------------------
+"""
+import numpy as np
+import gstools as gs
+from matplotlib import pyplot as plt
+
+###############################################################################
+# Generate a synthetic field with an exponential model.
+
+x = np.random.RandomState(19970221).rand(1000) * 100.0
+y = np.random.RandomState(20011012).rand(1000) * 100.0
+model = gs.Exponential(dim=2, var=2, len_scale=8)
+srf = gs.SRF(model, mean=0, seed=19970221)
+field = srf((x, y))
+
+###############################################################################
+# Estimate the variogram of the field with 40 bins and plot the result.
+
+bins = np.arange(40)
+bin_center, gamma = gs.vario_estimate_unstructured((x, y), field, bins)
+
+###############################################################################
+# Define a set of models to test.
+
+models = {
+    "gaussian": gs.Gaussian,
+    "exponential": gs.Exponential,
+    "matern": gs.Matern,
+    "stable": gs.Stable,
+    "rational": gs.Rational,
+    "linear": gs.Linear,
+    "circular": gs.Circular,
+    "spherical": gs.Spherical,
+}
+scores = {}
+
+###############################################################################
+# Iterate over all models, fit their variogram and calculate the r2 score.
+
+# plot the estimated variogram
+plt.scatter(bin_center, gamma, label="data")
+ax = plt.gca()
+
+# fit all models to the estimated variogram
+for model in models:
+    fit_model = models[model](dim=2)
+    para, pcov, r2 = fit_model.fit_variogram(bin_center, gamma, return_r2=True)
+    fit_model.plot(x_max=40, ax=ax)
+    scores[model] = r2
+
+###############################################################################
+# Create a ranking based on the score and determine the best models
+
+ranking = [
+    (k, v)
+    for k, v in sorted(scores.items(), key=lambda item: item[1], reverse=True)
+]
+print("RANKING")
+for i, (model, score) in enumerate(ranking, 1):
+    print(i, model, score)
+
+plt.show()