.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "examples/03_variogram/04_directional_3d.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. .. rst-class:: sphx-glr-example-title .. _sphx_glr_examples_03_variogram_04_directional_3d.py: Directional variogram estimation and fitting in 3D -------------------------------------------------- In this example, we demonstrate how to estimate a directional variogram by setting the estimation directions in 3D. Afterwards we will fit a model to this estimated variogram and show the result. .. GENERATED FROM PYTHON SOURCE LINES 10-17 .. code-block:: Python import matplotlib.pyplot as plt import numpy as np from mpl_toolkits.mplot3d import Axes3D import gstools as gs .. GENERATED FROM PYTHON SOURCE LINES 18-19 Generating synthetic field with anisotropy and rotation by Tait-Bryan angles. .. GENERATED FROM PYTHON SOURCE LINES 19-29 .. code-block:: Python dim = 3 # rotation around z, y, x angles = [np.deg2rad(90), np.deg2rad(45), np.deg2rad(22.5)] model = gs.Gaussian(dim=3, len_scale=[16, 8, 4], angles=angles) x = y = z = range(50) pos = (x, y, z) srf = gs.SRF(model, seed=1001) field = srf.structured(pos) .. GENERATED FROM PYTHON SOURCE LINES 30-32 Here we generate the axes of the rotated coordinate system to get an impression what the rotation angles do. .. GENERATED FROM PYTHON SOURCE LINES 32-37 .. code-block:: Python # All 3 axes of the rotated coordinate-system main_axes = gs.rotated_main_axes(dim, angles) axis1, axis2, axis3 = main_axes .. GENERATED FROM PYTHON SOURCE LINES 38-42 Now we estimate the variogram along the main axes. When the main axes are unknown, one would need to sample multiple directions and look for the one with the longest correlation length (flattest gradient). Then check the transversal directions and so on. .. GENERATED FROM PYTHON SOURCE LINES 42-54 .. code-block:: Python bin_center, dir_vario, counts = gs.vario_estimate( pos, field, direction=main_axes, bandwidth=10, sampling_size=2000, sampling_seed=1001, mesh_type="structured", return_counts=True, ) .. GENERATED FROM PYTHON SOURCE LINES 55-57 Afterwards we can use the estimated variogram to fit a model to it. Note, that the rotation angles need to be set beforehand. .. GENERATED FROM PYTHON SOURCE LINES 57-64 .. code-block:: Python print("Original:") print(model) model.fit_variogram(bin_center, dir_vario) print("Fitted:") print(model) .. rst-class:: sphx-glr-script-out .. code-block:: none Original: Gaussian(dim=3, var=1.0, len_scale=16.0, anis=[0.5, 0.25], angles=[1.57, 0.785, 0.393]) Fitted: Gaussian(dim=3, var=0.972, len_scale=13.0, nugget=0.0138, anis=[0.542, 0.281], angles=[1.57, 0.785, 0.393]) .. GENERATED FROM PYTHON SOURCE LINES 65-66 Plotting main axes and the fitted directional variogram. .. GENERATED FROM PYTHON SOURCE LINES 66-95 .. code-block:: Python fig = plt.figure(figsize=[10, 5]) ax1 = fig.add_subplot(121, projection=Axes3D.name) ax2 = fig.add_subplot(122) ax1.plot([0, axis1[0]], [0, axis1[1]], [0, axis1[2]], label="0.") ax1.plot([0, axis2[0]], [0, axis2[1]], [0, axis2[2]], label="1.") ax1.plot([0, axis3[0]], [0, axis3[1]], [0, axis3[2]], label="2.") ax1.set_xlim(-1, 1) ax1.set_ylim(-1, 1) ax1.set_zlim(-1, 1) ax1.set_xlabel("X") ax1.set_ylabel("Y") ax1.set_zlabel("Z") ax1.set_title("Tait-Bryan main axis") ax1.legend(loc="lower left") x_max = max(bin_center) ax2.scatter(bin_center, dir_vario[0], label="0. axis") ax2.scatter(bin_center, dir_vario[1], label="1. axis") ax2.scatter(bin_center, dir_vario[2], label="2. axis") model.plot("vario_axis", axis=0, ax=ax2, x_max=x_max, label="fit on axis 0") model.plot("vario_axis", axis=1, ax=ax2, x_max=x_max, label="fit on axis 1") model.plot("vario_axis", axis=2, ax=ax2, x_max=x_max, label="fit on axis 2") ax2.set_title("Fitting an anisotropic model") ax2.legend() plt.show() .. image-sg:: /examples/03_variogram/images/sphx_glr_04_directional_3d_001.png :alt: Tait-Bryan main axis, Fitting an anisotropic model :srcset: /examples/03_variogram/images/sphx_glr_04_directional_3d_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 96-97 Also, let's have a look at the field. .. GENERATED FROM PYTHON SOURCE LINES 97-99 .. code-block:: Python srf.plot() .. image-sg:: /examples/03_variogram/images/sphx_glr_04_directional_3d_002.png :alt: Field 3D structured (50, 50, 50), Plane :srcset: /examples/03_variogram/images/sphx_glr_04_directional_3d_002.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 4.751 seconds) .. _sphx_glr_download_examples_03_variogram_04_directional_3d.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: 04_directional_3d.ipynb <04_directional_3d.ipynb>` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: 04_directional_3d.py <04_directional_3d.py>` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: 04_directional_3d.zip <04_directional_3d.zip>` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_