""" Getting started with INTEGRATE ============================== This notebook contains a simple example for getting started with INTEGRATE. The workflow follows these main steps: 0. Load electromagnetic data (DAUGAARD.h5) 1a. Set up a prior model, creating the file PRIOR.h5 1b. Update the prior model with forward-modeled data in PRIOR.h5 2. Perform probabilistic inversion using integrate_rejection, creating POST.h5 3. Plot and analyze the results """ # %% import os os.environ["EM_FORWARD_METHOD"] = "anemone" os.environ["EM_FORWARD_DEVICE"] = "cuda" os.environ["REJECTION_BACKEND"] = "jax" import integrate as ig hardcopy = True import matplotlib.pyplot as plt import numpy as np # %% # 0. Get TTEM data # ---------------- # # Several test cases are available in the INTEGRATE package. # To see which cases are available, check the ``get_case_data`` function. # # The code below downloads the file DAUGAARD_AVG.h5 that contains # TTEM (time-domain electromagnetic) soundings from Daugaard, Denmark. # It also downloads the corresponding GEX file, TX07_20231016_2x4_RC20-33.gex, # which contains information about the TTEM system configuration and parameters. # %% case = 'DAUGAARD' files = ig.get_case_data(case=case, showInfo=2) f_data_h5 = files[0] file_gex= ig.get_gex_file_from_data(f_data_h5) print("Using data file: %s" % f_data_h5) print("Using GEX file: %s" % file_gex) # %% # Plot the geometry and data # ~~~~~~~~~~~~~~~~~~~~~~~~~~ # # ``ig.plot_geometry`` plots the spatial geometry of the data (i.e., the locations of the soundings). # ``ig.plot_data`` plots the measured electromagnetic data for each sounding. # %% # The commented line below would plot LINE, ELEVATION and data id as three scatter plots: # ig.plot_geometry(f_data_h5) # Each plot can be generated separately by specifying the 'pl' argument: ig.plot_geometry(f_data_h5, pl='LINE') ig.plot_geometry(f_data_h5, pl='ELEVATION') ig.plot_geometry(f_data_h5, pl='id') # %% useSubset = False X, Y, LINE, ELEVATION = ig.get_geometry(f_data_h5) # Find points within buffer distance Xl = np.array([544000, 543550]) Yl = np.array([6174500, 6176500]) #Xl = np.array([544000, 543550, 543000]) #Yl = np.array([6174500, 6176500, 6176400]) buffer = 10.0 indices, distances, segment_ids = ig.find_points_along_line_segments( X, Y, Xl, Yl, tolerance=buffer ) i_line = indices if useSubset: i_use = i_line plt.figure() ig.plot_geometry(f_data_h5, pl='ELEVATION') plt.plot(Xl, Yl,'ko', markersize=15) plt.plot(X[i_use], Y[i_use], 'ko', markersize=5) plt.title('Profile line') plt.show() else: # Use all data 1:len(X) i_use = np.arange(len(X)) # %% # The electromagnetic data (d_obs and d_std) can be plotted using ig.plot_data: ig.plot_data(f_data_h5, hardcopy=hardcopy) # Plot data channel 15 in an XY grid ig.plot_data_xy(f_data_h5, data_channel=15, cmap='jet'); # %% # 1. Set up the prior model ($\rho(\mathbf{m},\mathbf{d})$) # --------------------------------------------------------- # # In this example, a simple layered prior model will be considered. # The prior represents our initial beliefs about subsurface resistivity structure. # %% # 1a. Generate prior model parameters # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ # # First, we generate a sample of the prior model parameters, $\rho(\mathbf{m})$. # %% # As an example, we choose a simple layered Earth model. # The number of layers follows a chi-squared distribution with 4 degrees of freedom, # and the resistivity in each layer is log-uniformly distributed between [1,3000] Ωm. # # This creates N realizations of 3 types of model parameters:: # # PRIOR:/M1: 1D resistivity values in 1m thick layers down to 90m depth # PRIOR:/M2: 1D resistivity values in discrete parameter sets where the first # Nlayer parameters are resistivities, and the last Nlayer-1 # parameters are depths to the base of each layer # PRIOR:/M3: The number of layers in each model realization # %% # Select how many prior model realizations (N) should be generated N=2_000_000 N=1_000_000 N=100_000 f_prior_h5 = ig.prior_model_layered(N=N,lay_dist='chi2', NLAY_deg=3, RHO_min=1, RHO_max=3000, f_prior_h5='PRIOR_N%d.h5' % N, showInfo=1) #print('%s is used to hold prior realizations' % (f_prior_h5)) # %% # Plot summary statistics of the prior model for quality control of the prior choice ig.plot_prior_stats(f_prior_h5, hardcopy=hardcopy) # %% # 1b. Generate corresponding prior data # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ # # Next, we generate a corresponding sample of $\rho(\mathbf{d})$ (prior data distribution). # # The prior data, corresponding to the prior model parameters, are computed using # the GA-AEM electromagnetic forward modeling code and the GEX system configuration file. # %% # Option 1: Update the existing PRIOR.h5 file with forward-modeled data # set the default hdf file to be the f_prior_data_h5 (without extension) + file_gex (without extension) + .h5 f_prior_data_h5 = '%s_%s_Nh280_Nf12.h5' % (f_prior_h5[:-3], file_gex[:-4]) useExistingData = True if useExistingData: # if exist f_post_data_h5 then used it, otherwise create it import os if not os.path.exists(f_prior_data_h5): f_prior_data_h5 = ig.prior_data_em(f_prior_h5, file_gex, doMakePriorCopy=True, f_prior_data_h5=f_prior_data_h5) else: print('Using existing prior data file: %s' % f_prior_data_h5) #f_prior_data_h5 = ig.prior_data_em(f_prior_h5, file_gex, doMakePriorCopy=False) print('Updated %s to hold prior data (forward-modeled responses)' % (f_prior_data_h5)) # %% D = ig.load_prior_data(f_prior_data_h5)[0][0] # %% # It is useful to compare the prior data to the observed data before inversion. # If there is little to no overlap between observed and prior data, the inversion # is unlikely to succeed, indicating inconsistency between the prior model and observations. # In the figure below, you can see that the observed data (red) falls clearly within # the range of the prior data distribution. # %% ig.plot_data_prior(f_prior_data_h5,f_data_h5,nr=1000,hardcopy=hardcopy) # %% # 2. Sample the posterior distribution $\sigma(\mathbf{m})$ # --------------------------------------------------------- # # The posterior distribution is sampled using the extended rejection sampler. # %% # Rejection sampling of the posterior can be done with default settings using: #f_post_h5 = ig.integrate_rejection(f_prior_h5, f_data_h5) # However, you can control several important options. # You can choose to use only a subset of the prior data. Decreasing the sample # size makes the inversion faster but increasingly approximate. N_use = N # Number of prior samples to use (use all available) T_base = 1 # Base annealing temperature for rejection sampling autoT = 1 # Automatically estimate optimal annealing temperature f_post_h5 = ig.integrate_rejection(f_prior_data_h5, f_data_h5, f_post_h5 = 'POST.h5', N_use = N_use, autoT = autoT, T_base = T_base, showInfo=0, ip_range=i_use) # %% # Posterior statistics computation (typically done after inversion) # This computes summary statistics like mean, median, standard deviation # ig.integrate_posterior_stats(f_post_h5) # ig.integrate_posterior_stats(f_post_h5, ip_range = i_use) # %% # 3. Plot statistics from the posterior $\sigma(\mathbf{m})$ # ---------------------------------------------------------- # # Compare prior and posterior data # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ # # First, compare prior (beige) to posterior (black) data, along with observed data (red), # for specific measurement locations (data IDs). # %% ig.plot_data_prior_post(f_post_h5, i_plot=100,hardcopy=hardcopy) ig.plot_data_prior_post(f_post_h5, i_plot=i_use[0],hardcopy=hardcopy) ig.plot_data_prior_post(f_post_h5, i_plot=i_use[-1],hardcopy=hardcopy) # %% # Evidence and annealing temperature # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ # # The evidence quantifies how well the data fits the model, # while temperature controls the acceptance rate in rejection sampling. # %% # Plot the annealing temperature used for inversion ig.plot_T_EV(f_post_h5, pl='T',hardcopy=hardcopy) # Plot the evidence (log-likelihood) estimated during inversion ig.plot_T_EV(f_post_h5, pl='EV',hardcopy=hardcopy) # Plot the chi-squared data fit (normalized mean-loglikelihood) # Values less than one suggest overfitting # Values above one suggest underfitting ig.plot_T_EV(f_post_h5, pl='CHI2',hardcopy=hardcopy) # %% # Resistivity profiles # ~~~~~~~~~~~~~~~~~~~~ # # Plot a profile showing posterior statistics of model parameter M1 (resistivity) # along a section of the survey line. # %% # Plot resistivity profile for model M1 ig.plot_profile(f_post_h5, im=1, key='HarmonicMean', hardcopy=hardcopy) # Plot resistivity 'Median' profile for model M1 from data point i1 to i2 ig.plot_profile(f_post_h5, i1=1401, i2=2000, im=1, key='Median', hardcopy=hardcopy) # Plot resistivity 'HarmonicMean' profile for model M1 from data point i1 to i2 ig.plot_profile(f_post_h5, i1=1401, i2=2000, im=1, key='HarmonicMean', hardcopy=hardcopy) # Plot resistivity 'Mean' profile for model M1 from data point i1 to i2 ig.plot_profile(f_post_h5, i1=1401, i2=2000, im=1, key='Mean', hardcopy=hardcopy) # %% # Plot resistivity profile for model M1 for specific data points, along 'x', 'y' and 'index' axes ig.plot_profile(f_post_h5, ii=i_line, im=1, hardcopy=hardcopy, xaxis='x', gap_threshold=10) ig.plot_profile(f_post_h5, ii=i_line, im=1, hardcopy=hardcopy, xaxis='y', gap_threshold=10) ig.plot_profile(f_post_h5, ii=i_line, im=1, hardcopy=hardcopy, xaxis='index', gap_threshold=10) # ig.plot_profile(f_post_h5, ii=i_line, im=1, hardcopy=hardcopy, xaxis='index', panels=['Median'], gap_threshold=10) # %% plot random realizations for model M1 # PLot a single posterior realization for model M1 for i in range(1): ig.plot_profile(f_post_h5, ii=i_line, im=1, hardcopy=hardcopy, xaxis='x', gap_threshold=5, panels=['realization'], seed=i, title='Posterior realization', f_png='post_%d.png' % (i)) # Plot a single prior realization for model M1 ig.plot_profile(f_post_h5, ii=i_line, im=1, hardcopy=hardcopy, xaxis='x', gap_threshold=5, panels=['realization'], seed=i, plot_prior=True, title='Prior realization', f_png='prior_%d.png' % (i)) ''' mogrify -trim prior*png mogrify -trim post*png magick -delay 10 -loop 0 prior*.png prior_reals.gif magick -delay 10 -loop 0 post*.png post_reals.gif ''' ig.plot_profile(f_post_h5, ii=i_line, im=1, hardcopy=hardcopy, xaxis='x', gap_threshold=5) # %% # Plot 2D spatial features # ~~~~~~~~~~~~~~~~~~~~~~~~ # # Plot the median resistivity at specific depths (layers 5, 30, and 50) # to show lateral variations in subsurface structure. # %% # Plot 2D features: Resistivity at different depths try: ig.plot_feature_2d(f_post_h5,im=1,iz=5, key='Median', uselog=1, cmap='jet', s=2, hardcopy=hardcopy) plt.show() ig.plot_feature_2d(f_post_h5,im=1,iz=5, key='Mean', uselog=1, cmap='jet', s=2, hardcopy=hardcopy) plt.show() ig.plot_feature_2d(f_post_h5,im=1,iz=5, key='HarmonicMean', uselog=1, cmap='jet', s=2, hardcopy=hardcopy) plt.show() except: pass try: ig.plot_feature_2d(f_post_h5,im=1,iz=30, key='Median', uselog=1, cmap='jet', s=2, hardcopy=hardcopy) plt.show() except: pass try: ig.plot_feature_2d(f_post_h5,im=1,iz=50, key='Median', uselog=1, cmap='jet', s=2,hardcopy=hardcopy) plt.show() except: pass try: # Plot number of layers ig.plot_feature_2d(f_post_h5,im=3,key='Median', uselog=0, cmap='jet', s=2,hardcopy=hardcopy) plt.show() except: pass # %% # Plot 2D features: Resistivity at different elevations try: for ele in np.arange(60,-51,-20): plt.figure() ig.plot_feature_2d(f_post_h5,im=1,elevation=ele, key='Median', uselog=1, clim=[1,2000],s=2, hardcopy=hardcopy) plt.show() except: pass # %% try: # Plot a 2D feature: The estimated number of layers ig.plot_feature_2d(f_post_h5,im=3,iz=0,key='Median', uselog=0, clim=[1,6], cmap='jet', s=2,hardcopy=hardcopy) plt.show() except: pass # %% # Export results to CSV format # ---------------------------- # # Export the posterior results to CSV files for use in GIS software or further analysis. # %% f_csv, f_point_csv = ig.post_to_csv(f_post_h5) # %% # Read the exported CSV file for inspection # Example filename (actual filename will be generated automatically): #f_point_csv = 'POST_DAUGAARD_AVG_PRIOR_CHI2_NF_3_log-uniform_N100000_TX07_20231016_2x4_RC20-33_Nh280_Nf12_Nu100000_aT1_M1_point.csv' import pandas as pd df = pd.read_csv(f_point_csv) df.head() # %% # Optional: Use PyVista for 3D visualization of X,Y,Z coordinates with median resistivity plPyVista = False if plPyVista: import pyvista as pv import numpy as np from pyvista import examples #pv.set_jupyter_backend('client') pv.set_plot_theme("document") p = pv.Plotter(notebook=True) p = pv.Plotter() filtered_df = df[(df['Median'] < 50) | (df['Median'] > 200)] #filtered_df = df[(df['LINE'] > 1000) & (df['LINE'] < 1400) ] points = filtered_df[['X', 'Y', 'Z']].values[:] median = np.log10(filtered_df['Mean'].values[:]) opacity = np.where(filtered_df['Median'].values[:] < 100, 0.5, 1.0) #p.add_points(points, render_points_as_spheres=True, point_size=3, scalars=median, cmap='jet', opacity=opacity) p.add_points(points, render_points_as_spheres=True, point_size=6, scalars=median, cmap='hot') p.show_grid() p.show()