""" Daugaard: probabilistic raw-material assessment with INTEGRATE ============================================================== GEUS previously assessed raw-material (sand/gravel) potential at three target sub-areas in Daugaard ("Delområde 1/2/3") using a **sequential, deterministic** workflow: invert the tTEM data to a single best-fit resistivity model, interpret lithology, and hand-estimate overburden and raw-material volumes per sub-area, constrained to what raw-material boreholes could document. Those results live in ``ReferenceProjects/Sdr Felding og Daugard til Integrate.pptx`` and the target-area outlines in ``ReferenceProjects/Fokusområder_Daugaard/``. This notebook is organised in three parts: **Part A -- Resistivity only: probabilistic vs. deterministic.** Invert the tTEM data with a *generic* resistivity-only prior (a plain layered-Earth model, no lithology classes) and compare the resulting harmonic-mean / std resistivity section to an externally computed **deterministic** WorkBench least-squares (LSQ) inversion. This part is purely a demonstration of how probabilistic resistivity-only inversion and deterministic inversion relate to each other -- it is not used for the raw-material estimate. **Part B -- The full INTEGRATE workflow.** Build an *informed* prior from two geological scenarios (``daugaard_standard.xlsx`` + ``daugaard_valley.xlsx``, built with ``geoprior1d``, merged 50/50), forward-model the prior tTEM data, load and add the borehole data, run the joint (tTEM + borehole) rejection inversion, demonstrate the ``ig.query`` tool, grow a data-driven raw-material region from the posterior probability map (B8), and turn it into a low / median / high volume estimate (B9), built from per-sounding thickness percentiles times the Voronoi-cell area around each sounding. Part B is purely about the probabilistic result and ends with a bar chart of the grown region's raw-material volume. **Part C -- Comparison to Mette's original estimate.** Everything involving Mette's hand-drawn target polygons lives here: their own low/median/high volumes, and the region-vs-polygon comparison next to the single-number estimates from the old sequential/deterministic assessment. One important caveat, from GEUS's own notes on the old assessment (see ``ReferenceProjects/EmailFromMette.md``): the deterministic raw-material estimate was **capped at the depth documented by raw-material boreholes**, even where the tTEM inversion suggested raw material could extend deeper. The probabilistic estimate here is *not* capped this way, so a meaningfully larger estimate is expected -- a genuine methodological difference, not an error in either approach. Run this notebook from within the ``examples/`` directory (as with the other INTEGRATE examples) so that the case data downloads and the ``integrate_rawmaterial_utils`` helper module are found correctly. """ # %% try: get_*ipython() get_ipython().run_line_magic('load_ext', 'autoreload') get_ipython().run_line_magic('autoreload', '2') except Exception: pass # %% import os # Trim JAX/XLA compile time on the first run: # - xla_gpu_autotune_level=1 skips exhaustive GEMM/conv autotuning # - xla_backend_optimization_level=1 lowers XLA/LLVM optimisation effort # Neither affects the ptxas (PTX -> SASS) step. If compilation still takes # many minutes with a `ptxas` process at 100 % CPU, XLA is using the wrong # (system) ptxas -- see _ensure_bundled_ptxas in integrate_rejection_jax.py. # NOTE: XLA_FLAGS is read when jax is first imported (via `import integrate` # below), so in a notebook this cell must run before that -- restart the # kernel after editing it. os.environ["XLA_FLAGS"] = "--xla_gpu_autotune_level=1 --xla_backend_optimization_level=1" import h5py import numpy as np import matplotlib.pyplot as plt import time import integrate as ig import integrate_rawmaterial_utils as rmu hardcopy = True # Choose the forward model (the default is ga-aem) os.environ["EM_FORWARD_METHOD"] = "anemone" os.environ["EM_FORWARD_DEVICE"] = "cuda" #os.environ["EM_FORWARD_DEVICE"] = "cpu" #os.environ["EM_FORWARD_METHOD"] = "ga-aem" t_start = time.time() # %% # 0. Settings and fixed file names # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ # # Every expensive step (prior sampling, forward modelling, rejection) writes # to a **fixed** path defined below and is skipped with ``os.path.exists`` when # that file already exists. This makes it possible to run any one section on # its own: as long as the file-name variables in this cell are defined, a # later section can just load ``f_post_h5`` (etc.) without re-running the # earlier ones -- e.g. run Part C directly once ``f_post_h5`` exists. # # Every generated file name is prefixed with ``N_`` (and the working data # file also encodes ``inflateNoise``), so a different ``N`` writes to a separate # set of files / figures. If you change ``inflateNoise``, delete the affected # ``N*_DAUGAARD_*.h5`` files (or bump ``SUFFIX``) so they are regenerated. # %% # --- run-size settings ------------------------------------------------- N = 1_000_001 # production-scale #N = 10_000 # demo-scale; increase for a production-quality run N = 10_005 N = 100_006+1 # Prior size used everywhere: the generic prior (Part A) and each of the two # geological-scenario priors merged into the informed prior (Part B, N // 2 # realizations each). inflateNoise = 2 # multiply the tTEM data std by this (0 = use as-is) SUFFIX = '_N%d_iN%d' % (N, inflateNoise) # optional tag appended to every generated file name # Prefix prepended to every generated file name (h5 + figures) so that a # given choice of N writes to its own set of files, e.g. 'N12000_iN2_'. PREFIX = '' # --- fixed output file names ----------------------------------------- # working data file (noise-inflated copy; == downloaded file if inflateNoise == 0) f_data_work_h5 = '%sDAUGAARD_AVG_gf%g%s.h5' % (PREFIX, inflateNoise, SUFFIX) # Part A -- resistivity-only, generic prior f_prior_generic_h5 = '%sDAUGAARD_PRIOR_GENERIC%s.h5' % (PREFIX, SUFFIX) f_prior_generic_data_h5 = '%sDAUGAARD_PRIOR_GENERIC_DATA%s.h5' % (PREFIX, SUFFIX) f_post_generic_h5 = '%sDAUGAARD_POSTERIOR_GENERIC%s.h5' % (PREFIX, SUFFIX) # Part B -- informed prior + boreholes f_prior_scenario_h5_map = { 'daugaard_standard': '%sDAUGAARD_PRIOR_STANDARD%s.h5' % (PREFIX, SUFFIX), 'daugaard_valley': '%sDAUGAARD_PRIOR_VALLEY%s.h5' % (PREFIX, SUFFIX), } f_prior_merged_h5 = '%sDAUGAARD_PRIOR_MERGED%s.h5' % (PREFIX, SUFFIX) f_prior_data_h5 = '%sDAUGAARD_PRIOR_MERGED_DATA%s.h5' % (PREFIX, SUFFIX) f_prior_data_bh_h5 = '%sDAUGAARD_PRIOR_MERGED_DATA_BH%s.h5' % (PREFIX, SUFFIX) f_post_h5 = '%sDAUGAARD_POSTERIOR%s.h5' % (PREFIX, SUFFIX) f_post_jax_h5 = '%sDAUGAARD_POSTERIOR%s_jax.h5' % (PREFIX, SUFFIX) # %% # 1. Load the data # ~~~~~~~~~~~~~~~~ # # 1a. tTEM data and GEX system file # ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ # %% case = 'DAUGAARD' files = ig.get_case_data(case=case, showInfo=1) 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) X, Y, LINE, ELEVATION = ig.get_geometry(f_data_h5) ig.plot_geometry(f_data_h5, pl='LINE') ig.plot_data(f_data_h5, hardcopy=hardcopy) ig.plot_data_xy(f_data_h5, data_channel=15, cmap='jet') # %% if inflateNoise != 0: print("=" * 60) print("Using tTEM data with noise (std) inflated by a factor of %d" % inflateNoise) print("=" * 60) f_data_src_h5 = f_data_h5 f_data_h5 = f_data_work_h5 if not os.path.exists(f_data_h5): D = ig.load_data(f_data_src_h5) ig.copy_hdf5_file(f_data_src_h5, f_data_h5) ig.save_data_gaussian(D['d_obs'][0], D_std=D['d_std'][0] * inflateNoise, f_data_h5=f_data_h5, file_gex=file_gex) ig.plot_data(f_data_h5, useLog=0, hardcopy=hardcopy) plt.show() # %% # 1b. Target-area polygons and the west-to-east profile line # ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ # # The Daugaard focus-area shapefile (the areas Mette selected from the # deterministic interpretation) has an almost-empty attribute table, so the 3 # Delområder are matched to their names by polygon area, which is reported # (and matches exactly) in the old-approach PowerPoint. The polygons are # loaded here because the profile line used for the resistivity sections in # Parts A and B runs through their centroids; the polygons themselves are # used for the volume analysis in Part C. # %% DAUGAARD_AREA_ID_MAP = { 308812.2: 'Delområde 1', # DAU_02, DAU_06 106727.7: 'Delområde 2', # DAU_07 196598.0: 'Delområde 3', # DAU_05, DAU_03 } f_shp_daugaard = os.path.join('ReferenceProjects', 'Fokusområder_Daugaard', 'Fokusområder_polygon.shp') polygons = rmu.load_target_polygons(f_shp_daugaard, area_id_map=DAUGAARD_AREA_ID_MAP) # %% # Build a profile line through the centroid of each target polygon, ordered # from west (smallest X) to east (largest X), and find the sounding indices # along that line with ``ig.find_points_along_line_segments`` -- the same # function used for profile selection in ``integrate_workflow.py`` and # ``integrate_profiles.py``. These indices (``id_line``) are reused in Parts A # and B to plot resistivity/lithology sections along this single line. # %% centroids = {name: (geom.centroid.x, geom.centroid.y) for name, geom in polygons.items()} centroids_sorted = sorted(centroids.items(), key=lambda kv: kv[1][0]) # west (small X) -> east (large X) for name, (cx, cy) in centroids_sorted: print(" %-15s centroid = (%9.1f, %9.1f)" % (name, cx, cy)) Xl = np.array([cx for _, (cx, cy) in centroids_sorted]) Yl = np.array([cy for _, (cx, cy) in centroids_sorted]) Yl[0]=Yl[0]-200 buffer = 10.0 # Widen the profile at both ends so it also picks up every sounding # strictly west of the westernmost Delområde center and strictly east of # the easternmost one, not just soundings near the 3 centroids themselves: # extend the line with one extra waypoint beyond each end of the data's # X-range (at the Y of the nearest centroid), then re-search along the # now-longer line. pad = buffer Xl_wide = np.concatenate(([X.min() - pad], Xl, [X.max() + pad])) Yl_wide = np.concatenate(([Yl[0]], Yl, [Yl[-1]])) indices, distances, segment_ids = ig.find_points_along_line_segments(X, Y, Xl_wide, Yl_wide, tolerance=buffer) id_line = indices print("Found %d soundings within %.0f m of the widened west-to-east center profile" % (len(id_line), buffer)) pl_data = [id_line[10], id_line[100]] # Sanity check: are the tTEM soundings and the target polygons in the same # coordinate frame? (Both are UTM32N, but this is worth confirming visually # rather than assuming.) The selected profile line and its soundings are # overlaid on the same plot. rmu.plot_polygons_over_points(X, Y, polygons, title='Daugaard', hardcopy=hardcopy, profile_xy=(Xl_wide, Yl_wide), profile_idx=id_line, suffix=SUFFIX) # %% PLot data locations as well as the selected profile fig, ax, sc = ig.plot_xy(ELEVATION, X=X, Y=Y, cmap='jet', colorbar_label='Elevation [m]', title='Data locations and selected profile', s=3) ax.plot(X[id_line], Y[id_line], 'k.', label='profile soundings (id_line)') j=0 for i in pl_data: j=j+1 ax.plot(X[i], Y[i], 'wo', ms=16) ax.plot(X[i], Y[i], 'k*', ms=6+j*5, label='example data location %d (id=%d)' % (j,i)) ax.text(X[i] + 25, Y[i] + 25, 'd%d (%d)' % (j, i), color='k', fontsize=18) ax.legend(loc='best') if hardcopy: fig.savefig(PREFIX + 'daugaard_data_locations_profile' + SUFFIX + '.png', dpi=200, bbox_inches='tight') plt.show() # %% # Part A -- Resistivity only: probabilistic (generic prior) vs. deterministic # --------------------------------------------------------------------------- # # Invert the tTEM data with a *generic* resistivity-only prior (a plain # layered-Earth model, ``lay_dist='chi2'``, no lithology classes, as in # ``integrate_getting_started.py``). Plot the posterior harmonic-mean and std # resistivity along the west-to-east profile, then load the externally # computed deterministic WorkBench LSQ inversion and plot it the same way. # The point of this part is only to show how the probabilistic # resistivity-only result and the deterministic result relate to each other. t_generic_start = time.time() # %% if not os.path.exists(f_prior_generic_h5): ig.prior_model_layered( N=N, lay_dist='chi2', NLAY_deg=3, RHO_min=1, RHO_max=3000, f_prior_h5=f_prior_generic_h5, showInfo=1) ig.plot_prior_stats(f_prior_generic_h5, hardcopy=hardcopy) # %% if not os.path.exists(f_prior_generic_data_h5): ig.prior_data_em(f_prior_generic_h5, file_gex, doMakePriorCopy=True, f_prior_data_h5=f_prior_generic_data_h5) # %% if not os.path.exists(f_post_generic_h5): ig.integrate_rejection( f_prior_generic_data_h5, f_data_h5, f_post_h5=f_post_generic_h5, N_use=N, id_use=[1], autoT=1, T_base=1, showInfo=0, updatePostStat=True, backend='jax') t_generic_end = time.time() ig.plot_T_EV(f_post_generic_h5, pl='CHI2', hardcopy=hardcopy) # %% # Probabilistic resistivity-only profile: harmonic mean and std (log space) # ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ # %% ig.plot_profile(f_post_generic_h5, ii=id_line, im=1, panels=['harmonicmean', 'std'], xaxis='x', gap_threshold=100, hardcopy=hardcopy, txt='probabilistic_generic') # %% same but with alpha scaled by std # alpha>0 fades the harmonic-mean panel by uncertainty, normalised in the same # log10(std) space shown in the 'std' panel (default logstd_min=0.5, logstd_max=1.0): # opaque at LogStd<=0.5, fully transparent at LogStd>=1.0. ig.plot_profile(f_post_generic_h5, ii=id_line, im=1, panels=['harmonicmean', 'std'], xaxis='x', gap_threshold=100, hardcopy=hardcopy, alpha=.95, logstd_min=0.1, logstd_max=1.0, txt='probabilistic_generic_alpha') # %% # Deterministic (WorkBench LSQ) inversion # ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ # # Download the WorkBench LSQ result (smooth and sharp variants) and reproject # it onto the same depth grid as the probabilistic posterior, following the # pattern in ``integrate_paper_daugaard_supp.py``, then plot it as a normal # profile (mean + std, log space) along the same ``id_line``. # %% import libaarhusxyz ig.get_case_data(case=case, loadType='WB_smooth') ig.get_case_data(case=case, loadType='WB_sharp') f_xyz_list = { 'smooth': 'SCI7_40_ml_Daugaard_I01_MOD_inv.xyz', 'sharp': 'SCI7_40_ml_sharp2_I02_MOD_inv.xyz', } f_lsq_h5_list = {} for label, f_xyz in f_xyz_list.items(): f_lsq_h5 = os.path.splitext(f_xyz)[0] + '.h5' ig.copy_hdf5_file(f_post_generic_h5, f_lsq_h5) Xp, Yp, _, _ = ig.get_geometry(f_lsq_h5) with h5py.File(f_lsq_h5, 'r') as f: M1_median = f['/M1/Median'][:] with h5py.File(f_prior_generic_data_h5, 'r') as f: z = f['/M1'].attrs['x'][:].flatten() model = libaarhusxyz.XYZ(f_xyz) rho = model.layer_data['rho'].values rho_std = model.layer_data['rho_std'].values dep_top = model.layer_data['dep_top'].values utmx = model.flightlines['utmx'].values utmy = model.flightlines['utmy'].values M1_mean = np.full_like(M1_median, np.nan) M1_stdf = np.full_like(M1_median, np.nan) # raw WorkBench rho_std (uncertainty factor, >= 1) for i in range(len(utmx)): idx = np.argmin((Xp - utmx[i]) ** 2 + (Yp - utmy[i]) ** 2) j_idx = np.searchsorted(dep_top[i], z, side='right') - 1 valid = j_idx >= 0 M1_mean[idx, valid] = rho[i, j_idx[valid]] M1_stdf[idx, valid] = rho_std[i, j_idx[valid]] # Aarhus Workbench RHO_STD is the *standard-deviation factor* (STDF): the # parameter is estimated in log space, so it is multiplicative -- # 68% CI = [rho / RHO_STD, rho * RHO_STD] # std(log10 rho) = log10(RHO_STD) <- this is what the std panel plots # std(ln rho) = ln(RHO_STD) # RHO_STD == 1 -> perfectly resolved; large -> unresolved (Workbench caps # it at 99). It is NOT additive and NOT a percentage, so (RHO_STD - 1) is # only valid when RHO_STD is close to 1 (see the smooth vs. sharp header: # /MODEL TYPE, /MODEL UNIT, /INVERSION DATA SPACE = Logarithmic). M1_logstd = np.log10(M1_stdf) # std of log10(rho); values > 1 = unresolved M1_std = M1_mean * np.log(M1_stdf) # approx linear-space std [ohm-m] (rho * sigma_ln) # f_lsq_h5 was copied from the probabilistic posterior, so it still holds # that run's /M1 statistics. Overwrite EVERY statistic plot_profile might # read, or a stale probabilistic value shows through. In particular the # std panel plots /M1/LogStd. The LSQ result is a single model, so # mean = median = harmonic mean = the model itself. with h5py.File(f_lsq_h5, 'a') as f: for key in ['/M1/Mean', '/M1/LogMean', '/M1/Median', '/M1/HarmonicMean', '/M1/Std', '/M1/LogStd']: if key in f: del f[key] f['/M1/Mean'] = M1_mean f['/M1/Median'] = M1_mean f['/M1/HarmonicMean'] = M1_mean f['/M1/LogMean'] = np.log10(M1_mean) f['/M1/Std'] = M1_std # linear-space std [ohm-m] f['/M1/LogStd'] = M1_logstd # std of log10(rho) -- what the std panel plots f_lsq_h5_list[label] = f_lsq_h5 # %% # Same west-to-east profile (`id_line`) as the probabilistic plot above. for label, f_lsq_h5 in f_lsq_h5_list.items(): ig.plot_profile(f_lsq_h5, ii=id_line, im=1, panels=['mean', 'std'], xaxis='x', gap_threshold=100, hardcopy=hardcopy, txt='WB_%s' % label) ig.plot_profile(f_post_generic_h5, ii=id_line, im=1, panels=['harmonicmean', 'std'], xaxis='x', gap_threshold=100, hardcopy=hardcopy, txt='probabilistic_generic') # %% # **Takeaway**: the probabilistic resistivity-only inversion recovers a # section broadly consistent with the deterministic WorkBench LSQ result, but # additionally quantifies the uncertainty (std) at every point rather than # delivering only a single smooth/sharp model. This resistivity-only # comparison says nothing about raw-material volumes -- that needs the # informed lithology prior in Part B. # %% # Part B -- The full INTEGRATE workflow (informed prior + boreholes) # ------------------------------------------------------------------ t_full_start = time.time() # %% # B1. Informed prior from two geological scenarios # ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ # # The prior for the real workflow is built with ``geoprior1d`` # (https://github.com/GEUSjesper/geoprior1d) from Excel specifications of the # expected layer sequences, thicknesses, and resistivities for the Daugaard # geology -- this is the step GEUS owns for a new target area. # # Two scenarios are available: ``daugaard_standard.xlsx`` (background/regional # geology) and ``daugaard_valley.xlsx`` (geology inside a buried meltwater # valley). Since it is not known in advance whether a given sounding sits # inside or outside a buried valley, ``N // 2`` realizations are generated # independently from *each* Excel specification and combined into a single # prior with ``ig.merge_prior`` (50% from each scenario), letting the data # indicate, location by location, which scenario it is more consistent with # -- the same construction used in ``integrate_workflow.py`` and # ``integrate_daugaard_multi_prior.py``. # %% t_full_start_prior = time.time() from geoprior1d import geoprior1d ig.get_case_data(case=case, filelist=['daugaard_standard.xlsx', 'daugaard_valley.xlsx']) # N_prior_each realizations per geological scenario; the merged prior has 2x this. N_prior_each = N // 2 f_xlsx_files = ['daugaard_standard.xlsx', 'daugaard_valley.xlsx'] f_prior_h5_list = [] for file_xlsx in f_xlsx_files: fname = file_xlsx.split('.')[0] f_scen_h5 = f_prior_scenario_h5_map[fname] if not os.path.exists(f_scen_h5): geoprior1d(file_xlsx, Nreals=N_prior_each, dz=1, dmax=90, output_file=f_scen_h5) else: print("Using existing prior realizations: %s" % f_scen_h5) f_prior_h5_list.append(f_scen_h5) if not os.path.exists(f_prior_merged_h5): ig.merge_prior(f_prior_h5_list, f_prior_merged_h5=f_prior_merged_h5) else: print("Using existing merged prior: %s" % f_prior_merged_h5) f_prior_h5 = f_prior_merged_h5 t_full_end_prior = time.time() ig.plot_prior_stats(f_prior_h5, hardcopy=hardcopy) ig.prior_describe(f_prior_h5) # %% # B2. Boreholes # ^^^^^^^^^^^^^ # %% BHOLES = ig.read_borehole('daugaard_12boreholes.json', showInfo=1) ig.plot_boreholes(BHOLES, f_prior_h5, fontsize=17, hardcopy=hardcopy); # %% # B3. Prior tTEM data # ^^^^^^^^^^^^^^^^^^^ # # Forward-model the merged prior once (this is the expensive step; the result # is reused for the borehole prior data below). # %% t_full_start_prior_data_em = time.time() if not os.path.exists(f_prior_data_h5): ig.prior_data_em(f_prior_h5, file_gex, doMakePriorCopy=True, f_prior_data_h5=f_prior_data_h5) t_full_end_prior_data_em = time.time() # %% # B4. Borehole prior data # ^^^^^^^^^^^^^^^^^^^^^^^ # # Add the boreholes as extra, jointly inverted data types onto a fresh copy # of the tTEM prior-data file. ``ig.save_borehole_data`` takes the whole # ``BHOLES`` list in one call (no Python loop) and returns one prior/data ``/D`` # index per borehole, exactly as in ``integrate_workflow.py``. Boreholes are # appended right after the tTEM data (``/D1``), so on the cached path the data # IDs are the deterministic ``2 .. 1+len(BHOLES)``. Skipped entirely if the # posterior already exists. # %% t_full_start_prior_data_bd = time.time() im_prior = 2 id_borehole_list = list(range(2, 2 + len(BHOLES))) # deterministic /D indices (cached path) if not os.path.exists(f_post_h5) and not os.path.exists(f_prior_data_bh_h5): ig.copy_hdf5_file(f_prior_data_h5, f_prior_data_bh_h5) id_prior_list, id_borehole_list = ig.save_borehole_data( f_prior_data_bh_h5, f_data_h5, BHOLES, im_prior=im_prior, range_xyz=100, range_data=1, range_data_i_use=[17,18,19,20], doPlot=True, showInfo=0) else: print("Skipping borehole prior-data build (posterior or %s already exists)." % f_prior_data_bh_h5) t_full_end_prior_data_bd = time.time() # %% # ig.plot_discrete_data_entropy(f_data_h5, id_list=[2,3],plotPoints=True, ) # Entropy map over ALL multinomial (borehole) datasets — no id_list needed fig, ax, sc = ig.plot_discrete_data_entropy(f_data_h5, cmap = 'hot', plotPoints=True, plotPoints_size=1, plotPoints_color='k', plotPoints_alpha=0.1, hardcopy=False) # Overlay the borehole collar locations from the BHOLES list bx = [bh['X'] for bh in BHOLES] by = [bh['Y'] for bh in BHOLES] ax.plot(bx, by, 'rx', ms=9, mew=1.5, label='boreholes') for bh in BHOLES: ax.annotate(bh['name'], (bh['X'], bh['Y']), xytext=(4, 4), textcoords='offset points', fontsize=7) ax.legend(loc='best') fig.savefig('DAUGAARD_entropy_with_boreholes_N%d.png' % (N), dpi=150, bbox_inches='tight') # %% # B5. Joint (tTEM + borehole) rejection inversion # ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ # %% t_full_start_rejection = time.time() N_use = N # subset of the merged prior used in the rejection sampler if not os.path.exists(f_post_h5): id_use = [1] + id_borehole_list # tTEM (id 1) jointly with all borehole logs ig.integrate_rejection( f_prior_data_bh_h5, f_data_h5, f_post_h5=f_post_h5, N_use=N_use, id_use=id_use, nr=1000, T_N_above=50, T_P_acc_level=0.2, autoT=1, showInfo=1, updatePostStat=True, backend ='numpy') id_use = [1] ig.integrate_rejection( f_prior_data_bh_h5, f_data_h5, f_post_h5=f_post_jax_h5, N_use=N_use, id_use=id_use, nr=1000, T_N_above=50, T_P_acc_level=0.2, autoT=1, showInfo=1, updatePostStat=True, backend ='jax') else: print("Using existing posterior: %s" % f_post_h5) t_full_end_rejection = time.time() ig.plot_T_EV(f_post_h5, pl='CHI2', hardcopy=hardcopy) # %% # B6. Inspect the posterior # ^^^^^^^^^^^^^^^^^^^^^^^^^ # # Resistivity (harmonic mean, ``im=1``) and mode lithology (``im=2``) along the # west-to-east profile, then depth-slice maps. # %% ig.plot_profile(f_post_h5, im=1, ii=id_line, key='HarmonicMean', gap_threshold=100, xaxis='x', hardcopy=hardcopy) ig.plot_profile(f_post_h5, ii=id_line, im=1, xaxis='x', gap_threshold=100, hardcopy=hardcopy, alpha=.95, logstd_min=0.1, logstd_max=1.0, txt='_alpha' ) # %% ig.plot_profile(f_post_h5, im=2, ii=id_line, gap_threshold=100, xaxis='x', hardcopy=hardcopy) ig.plot_profile(f_post_h5, im=2, ii=id_line, gap_threshold=100, xaxis='x', alpha=0.95, entropy_min=0.1, entropy_max=1.0, txt='_alpha', hardcopy=hardcopy) # %% plot posterior models at selected data locations (similar the prior realisationz) for i_plot in pl_data: ig.plot_post_stats(f_post_h5, i_plot=i_plot, im=1, hardcopy=hardcopy) ig.plot_post_stats(f_post_h5, i_plot=i_plot, im=2, hardcopy=hardcopy) # %% for elevation in [40, 30, 20, 10]: ig.plot_feature_2d(f_post_h5, im=1, key='HarmonicMean', elevation=elevation, uselog=1, s=2, hardcopy=hardcopy) plt.show() for elevation in [40, 30, 20, 10]: ig.plot_feature_2d(f_post_h5, im=2, key='Mode', elevation=elevation, s=2, hardcopy=hardcopy) plt.show() # %% Plot prior data and posterior data fo a few select lociations # Prior (wheat/gray), posterior (black) and observed (red) tTEM data at a # couple of specific (X, Y) locations, using `ig.plot_data_prior_post`, which # takes a sounding *index* (`i_plot`) -- so each location is first snapped to # its nearest sounding via `argmin` on the same (X, Y) loaded in section 1a # (the same pattern used for the WorkBench reprojection in Part A and the # region-growing seed in B8). for i_plot in pl_data: #print("Location (%.1f, %.1f) -> nearest sounding #%d at (%.1f, %.1f)" # % (X_loc, Y_loc, i_plot, X[i_plot], Y[i_plot])) ig.plot_data_prior_post(f_post_h5, i_plot=i_plot, hardcopy=hardcopy) ig.plot_data_prior(f_post_h5, f_data_h5, hardcopy=hardcopy) # %% # B7. The ``ig.query`` tool # ^^^^^^^^^^^^^^^^^^^^^^^^^ # # We need to find WHERE we potentially have raw material worth producing! # %% # %% useLLM = False if 'ANTHROPIC_API_KEY' in os.environ: useLLM = True if useLLM: #os.environ['ANTHROPIC_API_KEY']='sk-ant-XXXXX' query_raw, interp_raw, prompt_raw = ig.query_from_text( 'What is the probbability that the cumulative thickness af raw materials (any sand and gravel) is greater than 10 m within the top 30 m, and where the oberburden (the top layer of non-raw material) is no more than 3 meters thick?', f_prior_h5=f_prior_h5, api_key=os.environ.get('ANTHROPIC_API_KEY')) else: query_raw = {'constraints': [{'im': 2, 'classes': [2, 5, 6], 'thickness_mode': 'cumulative', 'thickness_comparison': '>', 'thickness_threshold': 10.0, 'depth_min': 0.0, 'depth_max': 30.0, 'negate': False}, {'im': 2, 'classes': [1, 3, 4, 7, 8], 'thickness_mode': 'first_occurrence', 'thickness_comparison': '<=', 'thickness_threshold': 3.0, 'depth_min': 0.0, 'depth_max': 30.0, 'negate': False}]} # %% # perform the que raw_classes = query_raw['constraints'][0]['classes'] # sand + gravel: coarse raw material fine_classes = query_raw['constraints'][1]['classes'] # everything else: fine, non-raw material (overburden) t_query1_start = time.time() P_raw, meta_raw = ig.query(f_post_h5, query_raw) ig.query_plot(P_raw, meta_raw, query_text="P(raw material)", text_panel=True, hardcopy=PREFIX + 'daugaard_P_raw' + SUFFIX if hardcopy else False) t_query1_end = time.time() # %% # B8. Automatic region search: grow an area of potential raw material # ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ # # Delineate an area of potential raw material by letting a region grow itself # outward from the most promising spot: start at the sounding with the highest # raw-material probability ``P_raw``, then repeatedly add the neighbouring # sounding with the next-highest ``P_raw`` -- as long as it is still >= ``P_MIN``. # The region simply keeps spreading into adjacent high-probability ground and # stops on its own when the surrounding soundings are no longer likely enough. # # ``P_MIN`` is the main knob (the polygon "width"/area): "include every # connected sounding where the probability of raw material is at least # this". P_MIN = 0.5 -> the region is the contiguous patch (around the seed) # where raw material is more likely than not. Lower it to grow the region, # raise it to shrink it. ``MAX_AREA_M2`` is an optional hard cap (None = no # cap). # # Neighbours come from a Voronoi tessellation of the sounding locations (two # soundings are neighbours iff their Voronoi cells share an edge); each cell # is clipped to a *concave* hull of the soundings so its area is meaningful. # The region area is then the sum of its soundings' cell areas, and its # raw-material volume is sum(cell_area * thickness_percentile) (see B9). # # Before growing, **edge-affected soundings are dropped** (``flag_edge_cells``): # a sounding on the sparse rim of the survey has a huge, badly constrained # Voronoi cell that would otherwise let the region balloon into empty ground. # # **Organization.** This block first defines the reusable region-analysis # functions (below), then B9 runs the clean end-to-end pipeline (query -> # pick point(s) -> grow -> percentile -> volume). # # Two entry points:: # # * ``find_coherent_area(X, Y, P, p_min, ...)`` -- the full pipeline in one # call: builds the Voronoi adjacency + per-sounding cells, applies the # edge filter, grows a connected region from a seed, and returns the # region indices, its boundary polygon, AND the Voronoi scaffold (graph / # cells / areas / outline / keep-mask) for reuse. These now live in the # ``integrate`` package (integrate.integrate_query), called as ``ig.``. # # * ``grow_connected_region(P, neighbors, cell_area, p_min, ...)`` -- the # pure graph-grow step on a pre-built scaffold (no Voronoi/shapely). # ``find_coherent_area`` uses it internally. # %% # The region-search functions (voronoi_graph, voronoi_cells_ordered, # cells_to_polygon, flag_edge_cells, grow_connected_region, find_coherent_area) # are part of the `integrate` package -- call them as ig.(...). # See the B8 markdown above for what each does. POS_CENTER = [] POS_CENTER.append((543039.3,6175596.0)) # west end of the profile POS_CENTER.append((544500.0,6175800.0)) # west end of the profile AREA_LIST = [] for icenter in range(len(POS_CENTER)): # --- central nodel X_center = POS_CENTER[icenter][0] Y_center = POS_CENTER[icenter][1] # --- region-growing knobs (the polygon "width"/area is set by these two) --- P_MIN = 0.2 # keep growing while the next sounding has P_raw >= this MAX_AREA_M2 = None # optional hard cap on the region area [m^2]; None = no cap # --- edge-affected-cell filter (soundings on the sparse rim of the survey) --- HULL_RATIO = 0.10 # concave-hull tightness for the survey outline (0 = tight, 1 = convex) EDGE_BUFFER = None # m; a sounding this close to the outline is "near the edge". # None -> auto = 2 * sqrt(interior-median cell area) CELL_AREA_K = 6.0 # a cell is "oversized" if larger than K x the interior-median cell area ELONG_MAX = 4.0 # drop near-edge cells more elongated than this # (perimeter^2 / (4 pi area); 1 = disc); None to skip this test CELL_SIZE = 5.0 # m; raster resolution for the polygon-clipped cell areas used in B9 # Grow the region at the single highest-P_raw sounding. `AREA` (idx / mask / # polygon / vor / neighbors / cells / cell_area / boundary / good / ...) is # used directly everywhere below, including in B9. AREA = ig.find_coherent_area(X, Y, P_raw, P_MIN, max_area_m2=MAX_AREA_M2, X_center = X_center, Y_center = Y_center, hull_ratio=HULL_RATIO, edge_buffer=EDGE_BUFFER, cell_area_k=CELL_AREA_K, elong_max=ELONG_MAX) edge_affected = ~AREA['good'] P_raw_eff = np.where(AREA['good'], P_raw, -np.inf) # edge-affected soundings are unreachable print("Edge filter dropped %d / %d soundings; default region seed %d at (%.0f, %.0f), P_raw=%.2f" % (int(edge_affected.sum()), len(X), AREA['seed'], X[AREA['seed']], Y[AREA['seed']], P_raw[AREA['seed']])) plt.figure(figsize=(8, 8)) plt.plot(X, Y,'.', markersize=.1, color='lightgray', label='all soundings') plt.scatter(X, Y, c=P_raw, s=2, cmap='hot_r', vmin=0, vmax=1, label='kept soundings') #plt.plot(X[AREA['idx']], Y[AREA['idx']], 'k.', markersize=1, label='region') plt.plot(*AREA['polygon'].exterior.xy, 'k-', lw=2, label='polygon') plt.plot(AREA['X_center'], AREA['Y_center'], 'ko', ms=22, label='CENTER') plt.legend() plt.xlabel('X') plt.ylabel('Y') plt.title('P_raw with edge-affected soundings dropped') plt.axis('equal') plt.show() AREA_LIST.append(AREA) # %% # B9. Volume: query -> percentile -> volume # ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ # # One ``ig.query`` percentile query per quantity -- raw material (``raw_classes``: # sand+gravel) and overburden (``fine_classes``: every other, non-raw class) -- # run once over ALL soundings, exactly like the percentile-query examples in # ``integrate_query.py``. Volume at each percentile = sum(percentile_thickness * # cell_area) over the grown region's soundings (``AREA['mask']``, from B8), via # ``ig.region_volumes(pct, area)``. # %% # ---- A. one percentile query per quantity (run ONCE, shared by every AREA) - PCT = [5, 50, 95] # -> low / median / high query_overburden = {"im": 2, "classes": fine_classes, "thickness_mode": "cumulative", "depth_min": 0.0} query_raw_material = {"im": 2, "classes": raw_classes, "thickness_mode": "cumulative", "depth_min": 0.0} t_query2_start = time.time() pct_overburden, _ = ig.query(f_post_h5, {"metric": query_overburden, "percentiles": PCT}) pct_raw_material, _ = ig.query(f_post_h5, {"metric": query_raw_material, "percentiles": PCT}) # each: (N_sounding, len(PCT)) posterior thickness percentiles [m], one row per sounding t_query2_end = time.time() # ---- B. volume per grown area (one entry per AREA in AREA_LIST) ----------- AREA_NAMES = ['Area %d' % i for i in range(len(AREA_LIST))] V_overburden_list = [] V_raw_material_list = [] for iarea, AREA in enumerate(AREA_LIST): V_overburden = ig.region_volumes(pct_overburden, AREA) V_raw_material = ig.region_volumes(pct_raw_material, AREA) V_overburden_list.append(V_overburden) V_raw_material_list.append(V_raw_material) print("%-8s overburden PCT%s = %s m^3" % (AREA_NAMES[iarea], PCT, np.round(V_overburden).astype(int))) print("%-8s raw material PCT%s = %s m^3" % (AREA_NAMES[iarea], PCT, np.round(V_raw_material).astype(int))) # Voronoi cells coloured by P_raw; edge-affected cells hatched grey; this # area's grown region outlined. fig, ax, mappable = ig.plot_voronoi_cells(AREA, P=P_raw, vmin=0, vmax=1) ax.plot(*AREA['polygon'].exterior.xy, color='k', lw=2.5, label='grown region') ax.plot(AREA['X_center'], AREA['Y_center'], 'ko', ms=25, label='center location') ax.set_xlabel('UTM X (m)') ax.set_ylabel('UTM Y (m)') ax.set_title('%s -- Voronoi cells coloured by P_raw (grey hatch = edge-affected, dropped)' % AREA_NAMES[iarea]) fig.colorbar(mappable, ax=ax, label='P_raw') ax.legend(fontsize=8) if hardcopy: fig.savefig('%sdaugaard_voronoi_cells_area%d%s.png' % (PREFIX, iarea, SUFFIX), dpi=200, bbox_inches='tight') plt.show() # ---- C. raw-material volume, one bar per grown area ---------------------------- # Purely probabilistic -- comparison to Mette's polygons is in Part C. Vg = np.vstack(V_raw_material_list) xg = np.arange(len(AREA_LIST)) fig, ax = plt.subplots(figsize=(3 + 1.5 * len(AREA_LIST), 5)) ax.bar(xg, Vg[:, 1], yerr=[Vg[:, 1] - Vg[:, 0], Vg[:, 2] - Vg[:, 1]], capsize=5, color='C0') ax.set_xticks(xg) ax.set_xticklabels(AREA_NAMES) ax.set_ylabel('Raw-material volume (m$^3$)') ax.set_title('Raw-material volume per grown area (bar = P50, whiskers = P5-P95)') ax.grid(True, axis='y', ls='--', alpha=0.4) if hardcopy: fig.savefig(PREFIX + 'daugaard_rawmat_volume_B' + SUFFIX + '.png', dpi=200, bbox_inches='tight') plt.show() # %% # Part C -- Comparison to Mette's original (deterministic) estimate # ----------------------------------------------------------------- # # Everything that involves Mette's hand-drawn target polygons lives here -- # Part B is about the probabilistic result only. Per-polygon low / median / # high volumes reuse the same ``ig.region_volumes(pct, area)`` helper and the # ``pct_*`` arrays queried in B9, on an AREA-like dict built from the # per-sounding area inside each of Mette's polygons (rasterized at # ``CELL_SIZE`` resolution via ``rmu.compute_point_footprint_area``). Reference # numbers are transcribed from # ``ReferenceProjects/Sdr Felding og Daugard til Integrate.pptx`` # (cross-checked against the shapefile polygon areas). # # The grown region is no longer a single area: every area in ``AREA_LIST`` # (grown from its own centre in B8) is carried through the comparison, so # C2's bar chart shows one bar per grown area alongside one per Mette # polygon, and C3 overlays both sets on the Voronoi tessellation. # %% # C1. Low / median / high volumes for each of Mette's target polygons. prob_results = {} for name, polygon in polygons.items(): footprint_area = rmu.compute_point_footprint_area(X, Y, polygon, cell_size=CELL_SIZE) poly_area = {'mask': footprint_area > 0, 'cell_area': footprint_area} prob_results[name] = {'overburden': ig.region_volumes(pct_overburden, poly_area), 'raw_material': ig.region_volumes(pct_raw_material, poly_area)} print("%-15s overburden PCT%s = %s m^3" % (name, PCT, np.round(prob_results[name]['overburden']).astype(int))) print("%-15s raw material PCT%s = %s m^3" % (name, PCT, np.round(prob_results[name]['raw_material']).astype(int))) # %% # C2. Raw-material volume, every grown area vs. Mette's polygons (bar = P50, # whiskers = P5-P95), then the side-by-side tables from `rmu.compare_to_reference`. names = AREA_NAMES + list(prob_results) V = np.vstack(V_raw_material_list + [prob_results[n]['raw_material'] for n in prob_results]) x = np.arange(len(names)) fig, ax = plt.subplots(figsize=(5 + 1.2 * len(names), 5)) ax.bar(x, V[:, 1], yerr=[V[:, 1] - V[:, 0], V[:, 2] - V[:, 1]], capsize=5, color=['C0'] * len(AREA_NAMES) + ['0.7'] * len(prob_results)) ax.set_xticks(x) ax.set_xticklabels(names, rotation=30, ha='right') ax.set_ylabel('Raw-material volume (m$^3$)') ax.set_title('Raw-material volume (bar = P50, whiskers = P5-P95)') ax.grid(True, axis='y', ls='--', alpha=0.4) if hardcopy: fig.savefig(PREFIX + 'daugaard_rawmat_volume_C' + SUFFIX + '.png', dpi=200, bbox_inches='tight') plt.show() # %% # C3. Voronoi overview: cells coloured by P_raw, with Mette's target polygons # (shades of blue, thin solid lines) and the grown areas (black->light-grey, # thicker solid lines) drawn on top. The Voronoi scaffold (cells / edge filter # / survey outline) is identical for every AREA, so AREA_LIST[0] is the backdrop. fig, ax, mappable = ig.plot_voronoi_cells(AREA_LIST[0], P=P_raw, vmin=0, vmax=1) def _grey_shades(n, lo=0.0, hi=0.75): """n greyscale levels from black (`lo`) to light grey (`hi`).""" return [str(lo + (hi - lo) * (i / max(n - 1, 1))) for i in range(n)] area_greys = _grey_shades(len(AREA_LIST)) for i, AREA in enumerate(AREA_LIST): shade = area_greys[i] ax.plot(*AREA['polygon'].exterior.xy, color=shade, lw=3.0, zorder=3, label=AREA_NAMES[i]) ax.plot(AREA['X_center'], AREA['Y_center'], marker='o', color=shade, ms=12, mec='k', zorder=5) mette_blues = [plt.cm.Blues(0.45 + 0.5 * (i / max(len(polygons) - 1, 1))) for i in range(len(polygons))] for i, (name, polygon) in enumerate(polygons.items()): ax.plot(*polygon.exterior.xy, color=mette_blues[i], lw=1.0, zorder=4, label="Mette: %s" % name) ax.set_xlabel('UTM X (m)') ax.set_ylabel('UTM Y (m)') ax.set_title("Voronoi cells (P_raw) -- Mette's polygons vs. the grown areas") fig.colorbar(mappable, ax=ax, label='P_raw') ax.legend(fontsize=8, loc='best') if hardcopy: fig.savefig(PREFIX + 'daugaard_voronoi_areas_vs_mette' + SUFFIX + '.png', dpi=200, bbox_inches='tight') plt.show() DAUGAARD_REFERENCE = { # 'coarser' (gravel-only) is kept for documentation but not compared below -- # this workflow only distinguishes raw material (sand+gravel) vs. overburden # (everything else), not a further gravel-only split. 'Delområde 1': {'overburden': 865_000, 'raw_material': 4_500_000, 'coarser': 1_700_000}, 'Delområde 2': {'overburden': 566_000, 'raw_material': 1_700_000, 'coarser': 1_000_000}, 'Delområde 3': {'overburden': 452_000, 'raw_material': 3_300_000, 'coarser': 2_200_000}, } rmu.compare_to_reference(prob_results, DAUGAARD_REFERENCE, quantities=('overburden', 'raw_material'), hardcopy=hardcopy, f_name=PREFIX + 'daugaard_rawmaterial_comparison' + SUFFIX) # %% # Discussion # ^^^^^^^^^^ # # The probabilistic estimates come with an explicit min/median/max range # instead of a single number, which is the key practical addition of this # workflow: a decision-maker can see not just a central raw-material volume # estimate, but how confident that estimate is, sub-area by sub-area. # # If the probabilistic median raw-material volume is noticeably larger than # the old deterministic estimate, recall the caveat from GEUS's notes on the # old assessment (``ReferenceProjects/EmailFromMette.md``): the deterministic # estimate was capped at the depth documented by raw-material boreholes, even # where the tTEM inversion suggested raw material could extend deeper. The # probabilistic model is not capped this way, so part of any gap reflects # that methodological difference rather than a disagreement about the shallow # geology. Whether deep/older units in ``raw_classes`` (e.g. Miocene sand) # should count as producible raw material at all is a # geological judgement call that should be reviewed before treating this as a # strict apples-to-apples comparison (see Part B's printed class list). t_end = time.time() t_full_end = time.time() # %% TIMING, print timings for each major step in the workflow # total running time print("Total running time: %.1f s" % (t_end - t_start)) print(32*"*") # generic print("Generic posterior: %.1f s" % (t_generic_end - t_generic_start)) print(32*"*") print("Total running time (full workflow): %.1f s" % (t_full_end - t_full_start)) # prior data (tTEM) print("Prior data (tTEM): %.1f s" % (t_full_end_prior - t_full_start_prior)) # prior data (boreholes) print("Prior data (boreholes): %.1f s" % (t_full_end_prior_data_bd - t_full_start_prior_data_bd)) # rejection inversion print("Rejection inversion: %.1f s" % (t_full_end_rejection - t_full_start_rejection)) # query print("Query: %.1f s" % (t_query1_end - t_query1_start)) # percentile query print("Percentile query: %.1f s" % (t_query2_end - t_query2_start))