.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "examples/generated/scripts_synth_data/three_survey_lines_parameter_estimation.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_generated_scripts_synth_data_three_survey_lines_parameter_estimation.py: Invert three surveys line for a thin plate ========================================== .. raw:: html .. GENERATED FROM PYTHON SOURCE LINES 13-18 |Open In Colab| .. |Open In Colab| image:: https://img.shields.io/badge/open%20in-Colab-b5e2fa?logo=googlecolab&style=flat-square&color=ffd670 :target: https://colab.research.google.com/github/inlab-geo/cofi-examples/blob/main/examples/airborne_em/airborne_em_three_lines_transmitters.ipynb .. GENERATED FROM PYTHON SOURCE LINES 21-40 .. raw:: html .. raw:: html .. If you are running this notebook locally, make sure you’ve followed `steps here `__ to set up the environment. (This `environment.yml `__ file specifies a list of packages required to run the notebooks) .. GENERATED FROM PYTHON SOURCE LINES 40-52 .. code-block:: Python # -------------------------------------------------------- # # # # Uncomment below to set up environment on "colab" # # # # -------------------------------------------------------- # # !pip install -U cofi # !pip install git+https://github.com/JuergHauser/PyP223.git # !git clone https://github.com/inlab-geo/cofi-examples.git # %cd cofi-examples/examples/vtem_max .. GENERATED FROM PYTHON SOURCE LINES 54-59 .. code-block:: Python # If this notebook is run locally PyP223 needs to be installed separately by uncommenting the following line, # that is by removing the # and the white space between it and the exclamation mark. # !pip install git+https://github.com/JuergHauser/PyP223.git .. GENERATED FROM PYTHON SOURCE LINES 61-80 .. code-block:: Python import numpy import matplotlib.pyplot as plt from matplotlib.lines import Line2D import cofi from vtem_max_forward_lib import ( problem_setup, system_spec, survey_setup, ForwardWrapper, plot_predicted_profile, plot_transient, plot_plate_faces, plot_plate_faces_single ) numpy.random.seed(42) .. GENERATED FROM PYTHON SOURCE LINES 85-94 Background ---------- This example inverts a three survey line of VTEM max data using the vertical component for a thin plate target. It thus becomes possible to invert for the easting, depth of the plate reference point, the plate dip and plate azimuth and the plate length. The forward problem and model setup is described `here <./thin_plate_inversion.ipynb>`__. .. GENERATED FROM PYTHON SOURCE LINES 97-100 Problem definition ------------------ .. GENERATED FROM PYTHON SOURCE LINES 100-114 .. code-block:: Python tx_min = 115 tx_max = 281 tx_interval = 15 ty_min = 25 ty_max = 176 ty_interval = 75 tx_points = numpy.arange(tx_min, tx_max, tx_interval) ty_points = numpy.arange(ty_min, ty_max, ty_interval) n_transmitters = len(tx_points) * len(ty_points) tx, ty = numpy.meshgrid(tx_points, ty_points) tx = tx.flatten() ty = ty.flatten() .. GENERATED FROM PYTHON SOURCE LINES 116-123 .. code-block:: Python fiducial_id = numpy.arange(len(tx)) line_id = numpy.zeros(len(tx), dtype=int) line_id[ty==ty_points[0]] = 0 line_id[ty==ty_points[1]] = 1 line_id[ty==ty_points[2]] = 2 .. GENERATED FROM PYTHON SOURCE LINES 125-142 .. code-block:: Python survey_setup = { "tx": tx, # transmitter easting/x-position "ty": ty, # transmitter northing/y-position "tz": numpy.array([50]*n_transmitters), # transmitter height/z-position "tazi": numpy.deg2rad(numpy.array([90]*n_transmitters)), # transmitter azimuth "tincl": numpy.deg2rad(numpy.array([6]*n_transmitters)), # transmitter inclination "rx": tx, # receiver easting/x-position "ry": numpy.array([100]*n_transmitters), # receiver northing/y-position "rz": numpy.array([50]*n_transmitters), # receiver height/z-position "trdx": numpy.array([0]*n_transmitters), # transmitter receiver separation inline "trdy": numpy.array([0]*n_transmitters), # transmitter receiver separation crossline "trdz": numpy.array([0]*n_transmitters), # transmitter receiver separation vertical "fiducial_id": fiducial_id, # unique id for each transmitter "line_id": line_id # id for each line } .. GENERATED FROM PYTHON SOURCE LINES 144-160 .. code-block:: Python true_model = { "res": numpy.array([300, 1000]), "thk": numpy.array([20]), "peast": numpy.array([175]), "pnorth": numpy.array([100]), "ptop": numpy.array([30]), "pres": numpy.array([0.1]), "plngth1": numpy.array([100]), "plngth2": numpy.array([100]), "pwdth1": numpy.array([0.1]), "pwdth2": numpy.array([90]), "pdzm": numpy.array([75]), "pdip": numpy.array([60]) } .. GENERATED FROM PYTHON SOURCE LINES 162-166 .. code-block:: Python forward = ForwardWrapper(true_model, problem_setup, system_spec, survey_setup, ["pdip","pdzm", "peast", "ptop", "pwdth2"], data_returned=["vertical"]) .. rst-class:: sphx-glr-script-out .. code-block:: none ['pdip', 'pdzm', 'peast', 'ptop', 'pwdth2'] .. GENERATED FROM PYTHON SOURCE LINES 168-172 .. code-block:: Python # check the order of parameters in a model vector forward.params_to_invert .. rst-class:: sphx-glr-script-out .. code-block:: none ['pdip', 'pdzm', 'peast', 'ptop', 'pwdth2'] .. GENERATED FROM PYTHON SOURCE LINES 174-177 .. code-block:: Python true_param_value = numpy.array([60,65, 175, 30, 90]) .. GENERATED FROM PYTHON SOURCE LINES 182-185 True model ~~~~~~~~~~ .. GENERATED FROM PYTHON SOURCE LINES 185-203 .. code-block:: Python _, axes = plt.subplots(2, 2) axes[1,1].axis("off") plot_plate_faces( "plate_true", forward, true_param_value, axes[0,0], axes[0,1], axes[1,0], color="purple", label="True model" ) plt.tight_layout() point = Line2D([0], [0], label='Fiducial', marker='o', markersize=5, markeredgecolor='orange', markerfacecolor='orange', linestyle='') handles, labels = axes[1,0].get_legend_handles_labels() handles.extend([point]) axes[1,0].legend(handles=handles,bbox_to_anchor=(1.04, 0), loc="lower left") .. image-sg:: /examples/generated/scripts_synth_data/images/sphx_glr_three_survey_lines_parameter_estimation_001.png :alt: three survey lines parameter estimation :srcset: /examples/generated/scripts_synth_data/images/sphx_glr_three_survey_lines_parameter_estimation_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 208-211 Generate synthetic data ~~~~~~~~~~~~~~~~~~~~~~~ .. GENERATED FROM PYTHON SOURCE LINES 211-223 .. code-block:: Python # The data absolute_noise= 0.05 # create data and ad a realisation of the noise data_pred_true = forward(true_param_value) data_obs = data_pred_true + numpy.random.randn(len(data_pred_true))*absolute_noise # define data covariance matrix sigma=absolute_noise Cdinv=numpy.identity(len(data_obs))*(1.0/(sigma*sigma)) .. GENERATED FROM PYTHON SOURCE LINES 228-231 Parameter estimation -------------------- .. GENERATED FROM PYTHON SOURCE LINES 234-236 **Initialise a model for inversion** .. GENERATED FROM PYTHON SOURCE LINES 236-239 .. code-block:: Python init_param_value = numpy.array([45, 90, 150, 20, 80]) .. GENERATED FROM PYTHON SOURCE LINES 244-246 **Define helper functions for CoFI** .. GENERATED FROM PYTHON SOURCE LINES 246-273 .. code-block:: Python def my_objective(model): dpred = forward(model) residual = dpred - data_obs return residual.T @ Cdinv @ residual def my_gradient(model): dpred = forward(model) jacobian = forward.jacobian(model, relative_step=0.1) residual = dpred - data_obs return jacobian.T @ Cdinv @ residual def my_hessian(model): jacobian = forward.jacobian(model) return jacobian.T @ Cdinv @ jacobian class PerIterationCallbackFunction: def __init__(self): self.x = None self.i = 0 def __call__(self, xk): print(f"Iteration #{self.i+1}") print(f" objective value: {my_problem.objective(xk)}") self.x = xk self.i += 1 .. GENERATED FROM PYTHON SOURCE LINES 278-280 **Define CoFI problem** .. GENERATED FROM PYTHON SOURCE LINES 280-287 .. code-block:: Python my_problem = cofi.BaseProblem() my_problem.set_objective(my_objective) my_problem.set_gradient(my_gradient) my_problem.set_hessian(my_hessian) my_problem.set_initial_model(init_param_value) .. GENERATED FROM PYTHON SOURCE LINES 292-294 **Define CoFI options** .. GENERATED FROM PYTHON SOURCE LINES 294-299 .. code-block:: Python my_options = cofi.InversionOptions() my_options.set_tool("scipy.optimize.minimize") my_options.set_params(method="Newton-CG",callback=PerIterationCallbackFunction()) .. GENERATED FROM PYTHON SOURCE LINES 304-306 **Run CoFI inversion** .. GENERATED FROM PYTHON SOURCE LINES 306-311 .. code-block:: Python my_inversion = cofi.Inversion(my_problem, my_options) my_result = my_inversion.run() print(my_result.model) .. rst-class:: sphx-glr-script-out .. code-block:: none Iteration #1 objective value: 25216.587572321012 Iteration #2 objective value: 7290.678246750972 Iteration #3 objective value: 3101.2129903115647 Iteration #4 objective value: 1982.6413889061635 Iteration #5 objective value: 1684.5275371363728 Iteration #6 objective value: 1585.5627896068984 Iteration #7 objective value: 1569.5279365139831 Iteration #8 objective value: 1553.5849324815508 Iteration #9 objective value: 1552.8316200317665 Iteration #10 objective value: 1551.64208272022 Iteration #11 objective value: 1550.9006994110896 [ 59.89497772 65.24627406 174.87958828 29.9595983 90.03333672] .. GENERATED FROM PYTHON SOURCE LINES 316-319 Plotting -------- .. GENERATED FROM PYTHON SOURCE LINES 322-325 Data - Fiducials ~~~~~~~~~~~~~~~~ .. GENERATED FROM PYTHON SOURCE LINES 325-334 .. code-block:: Python _, ax = plt.subplots() plot_transient(true_param_value, forward, "data from true model", ax, color="purple") plot_transient(init_param_value, forward, "data from init model", ax, color="green", linestyle=":") plot_transient(my_result.model, forward, "data from inverted model", ax, color="red", linestyle="-.") ax.legend(loc="upper center") ax.set_title("vertical") plt.tight_layout() .. image-sg:: /examples/generated/scripts_synth_data/images/sphx_glr_three_survey_lines_parameter_estimation_002.png :alt: vertical :srcset: /examples/generated/scripts_synth_data/images/sphx_glr_three_survey_lines_parameter_estimation_002.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 339-342 Data - Profiles ~~~~~~~~~~~~~~~ .. GENERATED FROM PYTHON SOURCE LINES 342-359 .. code-block:: Python # Select gates to plot idx_to_plot = numpy.arange(8, 30) _, axes = plt.subplots(3, 1, figsize=(12,12)) for i in range(3): plot_predicted_profile(true_param_value, forward, "Data from true model", gate_idx=idx_to_plot, line_id=[i], ax=axes[i], color="purple") plot_predicted_profile(init_param_value, forward, "Data from starting model", gate_idx=idx_to_plot, line_id=[i], ax=axes[i], color="green", linestyle=":") plot_predicted_profile(my_result.model, forward, "Data from MAP model", gate_idx=idx_to_plot, line_id=[i], ax=axes[i], color="red", linestyle="-.") axes[i].set_title("Crossline {} (m)".format(ty_points[i])) axes[i].legend(bbox_to_anchor=(1.04, 0), loc="lower left") plt.tight_layout() .. image-sg:: /examples/generated/scripts_synth_data/images/sphx_glr_three_survey_lines_parameter_estimation_003.png :alt: Crossline 25 (m), Crossline 100 (m), Crossline 175 (m) :srcset: /examples/generated/scripts_synth_data/images/sphx_glr_three_survey_lines_parameter_estimation_003.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 364-374 Model ~~~~~ As a general rule we can only constrain parateters if there are fiducials in the survey that are ssnsitvit tothem and also fiducials that are not sensitive to them. In order words the anomly needs to be closed with respect to the model paratmeter in question. So to be able to constraint plate lenghy we would need survey lines to the north and south that are not overflying the thin plate. .. GENERATED FROM PYTHON SOURCE LINES 374-400 .. code-block:: Python _, axes = plt.subplots(2, 2) axes[1,1].axis("off") plot_plate_faces( "plate_true", forward, true_param_value, axes[0,0], axes[0,1], axes[1,0], color="purple", label="True model" ) plot_plate_faces( "plate_init", forward, init_param_value, axes[0,0], axes[0,1], axes[1,0], color="green", label="Starting model" ) plot_plate_faces( "plate_inverted", forward, my_result.model, axes[0,0], axes[0,1], axes[1,0], color="red", label="MAP model", linestyle="dotted" ) plt.tight_layout() point = Line2D([0], [0], label='Fiducial', marker='o', markersize=5, markeredgecolor='orange', markerfacecolor='orange', linestyle='') handles, labels = axes[1,0].get_legend_handles_labels() handles.extend([point]) axes[1,0].legend(handles=handles,bbox_to_anchor=(1.04, 0), loc="lower left") .. image-sg:: /examples/generated/scripts_synth_data/images/sphx_glr_three_survey_lines_parameter_estimation_004.png :alt: three survey lines parameter estimation :srcset: /examples/generated/scripts_synth_data/images/sphx_glr_three_survey_lines_parameter_estimation_004.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 405-418 -------------- Watermark ========= .. raw:: html .. raw:: html .. GENERATED FROM PYTHON SOURCE LINES 418-424 .. code-block:: Python watermark_list = ["cofi", "numpy", "scipy", "matplotlib","smt"] for pkg in watermark_list: pkg_var = __import__(pkg) print(pkg, getattr(pkg_var, "__version__")) .. rst-class:: sphx-glr-script-out .. code-block:: none cofi 0.2.11 numpy 2.3.5 scipy 1.17.0 matplotlib 3.10.8 smt 2.10.1 .. GENERATED FROM PYTHON SOURCE LINES 425-425 sphinx_gallery_thumbnail_number = -1 .. rst-class:: sphx-glr-timing **Total running time of the script:** (5 minutes 42.213 seconds) .. _sphx_glr_download_examples_generated_scripts_synth_data_three_survey_lines_parameter_estimation.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: three_survey_lines_parameter_estimation.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: three_survey_lines_parameter_estimation.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: three_survey_lines_parameter_estimation.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_