.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "content/examples/20-published/plot_load_booky.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_content_examples_20-published_plot_load_booky.py: Heagy et al., 2017 Load and Plot Bookpurnong Data ================================================= In this example, we load and plot the SkyTEM (2006) and RESOLVE (2008) Bookpurnong data, available at `https://storage.googleapis.com/simpeg/bookpurnong/bookpurnong.tar.gz `_ This is published in Lindsey J. Heagy, Rowan Cockett, Seogi Kang, Gudni K. Rosenkjaer, Douglas W. Oldenburg, A framework for simulation and inversion in electromagnetics, Computers & Geosciences, Volume 107, 2017, Pages 1-19, ISSN 0098-3004, http://dx.doi.org/10.1016/j.cageo.2017.06.018. The script and figures are also on figshare: https://doi.org/10.6084/m9.figshare.5107711 This example was updated for SimPEG 0.14.0 on January 31st, 2020 by Joseph Capriotti .. GENERATED FROM PYTHON SOURCE LINES 21-257 .. image-sg:: /content/examples/20-published/images/sphx_glr_plot_load_booky_001.png :alt: Resolve In-phase 400 Hz, SkyTEM High moment 156 $\mu$s :srcset: /content/examples/20-published/images/sphx_glr_plot_load_booky_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none Downloading https://storage.googleapis.com/simpeg/bookpurnong/bookpurnong_inversion.tar.gz saved to: /home/vsts/work/1/s/examples/20-published/bookpurnong_inversion.tar.gz Download completed! | .. code-block:: Python import numpy as np import matplotlib.pyplot as plt import h5py import tarfile import os import shutil from simpeg import utils import discretize def download_and_unzip_data( url="https://storage.googleapis.com/simpeg/bookpurnong/bookpurnong_inversion.tar.gz", ): """ Download the data from the storage bucket, unzip the tar file, return the directory where the data are """ # download the data downloads = utils.download(url) # directory where the downloaded files are directory = downloads.split(".")[0] # unzip the tarfile tar = tarfile.open(downloads, "r") tar.extractall() tar.close() return downloads, directory def save_dict_to_hdf5(fname, dictionary): """ Save a dictionary to hdf5 """ f = h5py.File(fname, "w") for key in dictionary.keys(): f.create_dataset(key, data=dictionary[key]) f.close() def run(plotIt=True, saveIt=False, saveFig=False, cleanup=True): """ Download and plot the Bookpurnong data. Here, we parse the data into a dictionary that can be easily saved and loaded into other worflows (for later on when we are doing the inversion) :param bool plotIt: show the Figures? :param bool saveIt: re-save the parsed data? :param bool saveFig: save the matplotlib figures? :param bool cleanUp: remove the downloaded and saved data? """ downloads, directory = download_and_unzip_data() # data are in a directory inside bookpurnong_inversion data_directory = os.path.sep.join([directory, "bookpurnong_data"]) # Load RESOLVE (2008) header_resolve = "Survey Date Flight fid utctime helicopter_easting helicopter_northing gps_height bird_easting bird_northing bird_gpsheight elevation bird_height bird_roll bird_pitch bird_yaw em[0] em[1] em[2] em[3] em[4] em[5] em[6] em[7] em[8] em[9] em[10] em[11] Line " header_resolve = header_resolve.split() resolve = np.loadtxt( os.path.sep.join([data_directory, "Bookpurnong_Resolve_Exported.XYZ"]), skiprows=8, ) # Omit the cross lines resolve = resolve[(resolve[:, -1] > 30002) & (resolve[:, -1] < 38000), :] dat_header_resolve = "CPI400_F CPQ400_F CPI1800_F CPQ1800_F CXI3300_F CXQ3300_F CPI8200_F CPQ8200_F CPI40k_F CPQ40k_F CPI140k_F CPQ140k_F " dat_header_resolve = dat_header_resolve.split() xyz_resolve = resolve[:, 8:11] data_resolve = resolve[:, 16:-1] # Load SkyTEM (2006) fid = open( os.path.sep.join( [data_directory, "SK655CS_Bookpurnong_ZX_HM_TxInc_newDTM.txt"] ), "rb", ) lines = fid.readlines() fid.close() header_skytem = lines[0].split() info_skytem = [] data_skytem = [] for line in lines[1:]: if len(line.split()) != 65: info_skytem.append(np.array(line.split()[:16], dtype="O")) data_skytem.append(np.array(line.split()[16 : 16 + 24], dtype="float")) else: info_skytem.append(np.array(line.split()[:16], dtype="O")) data_skytem.append(np.array(line.split()[17 : 17 + 24], dtype="float")) info_skytem = np.vstack(info_skytem) data_skytem = np.vstack(data_skytem) lines_skytem = info_skytem[:, 1].astype(float) inds = lines_skytem < 2026 info_skytem = info_skytem[inds, :] data_skytem = data_skytem[inds, :].astype(float) xyz_skytem = info_skytem[:, [13, 12]].astype(float) lines_skytem = info_skytem[:, 1].astype(float) # Load path of Murray River river_path = np.loadtxt(os.path.sep.join([directory, "MurrayRiver.txt"])) # Plot the data nskip = 40 fig = plt.figure(figsize=(12 * 0.8, 6 * 0.8)) title = ["Resolve In-phase 400 Hz", r"SkyTEM High moment 156 $\mu$s"] ax1 = plt.subplot(121) ax2 = plt.subplot(122) axs = [ax1, ax2] out_re = utils.plot2Ddata( xyz_resolve[::nskip, :2], data_resolve[::nskip, 0], ncontour=100, contourOpts={"cmap": "viridis"}, ax=ax1, ) vmin, vmax = out_re[0].get_clim() cb_re = plt.colorbar( out_re[0], ticks=np.linspace(vmin, vmax, 3), ax=ax1, fraction=0.046, pad=0.04 ) temp_skytem = data_skytem[:, 5].copy() temp_skytem[data_skytem[:, 5] > 7e-10] = 7e-10 out_sky = utils.plot2Ddata( xyz_skytem[:, :2], temp_skytem, ncontour=100, contourOpts={"cmap": "viridis", "vmax": 7e-10}, ax=ax2, ) vmin, vmax = out_sky[0].get_clim() cb_sky = plt.colorbar( out_sky[0], ticks=np.linspace(vmin, vmax, 3), ax=ax2, format="%.1e", fraction=0.046, pad=0.04, ) cb_re.set_label("Bz (ppm)") cb_sky.set_label("Voltage (V/A-m$^4$)") for i, ax in enumerate(axs): xticks = [460000, 463000] yticks = [6195000, 6198000, 6201000] ax.set_xticks(xticks) ax.set_yticks(yticks) ax.plot(river_path[:, 0], river_path[:, 1], "k", lw=0.5) ax.set_aspect("equal") if i == 1: ax.plot(xyz_skytem[:, 0], xyz_skytem[:, 1], "k.", alpha=0.02, ms=1) ax.set_yticklabels([str(" ") for f in yticks]) else: ax.plot(xyz_resolve[:, 0], xyz_resolve[:, 1], "k.", alpha=0.02, ms=1) ax.set_yticklabels([str(f) for f in yticks]) ax.set_ylabel("Northing (m)") ax.set_xlabel("Easting (m)") ax.set_title(title[i]) plt.tight_layout() if plotIt: plt.show() if saveFig: fig.savefig("bookpurnong_data.png") cs, ncx, npad = 1.0, 10.0, 20 hx = [(cs, ncx), (cs, npad, 1.3)] npad = 12 temp = np.logspace(np.log10(1.0), np.log10(12.0), 19) temp_pad = temp[-1] * 1.3 ** np.arange(npad) hz = np.r_[temp_pad[::-1], temp[::-1], temp, temp_pad] mesh = discretize.CylindricalMesh([hx, 1, hz], "00C") active = mesh.cell_centers_z < 0.0 dobs_re = np.load(os.path.sep.join([directory, "dobs_re_final.npy"])) dpred_re = np.load(os.path.sep.join([directory, "dpred_re_final.npy"])) mopt_re = np.load(os.path.sep.join([directory, "mopt_re_final.npy"])) # Down sample resolve data nskip = 40 inds_resolve = np.r_[np.array(range(0, data_resolve.shape[0] - 1, nskip)), 16730] booky_resolve = { "data": data_resolve[inds_resolve, :], "data_header": dat_header_resolve, "line": resolve[:, -1][inds_resolve], "xy": xyz_resolve[:, :2][inds_resolve], "src_elevation": resolve[:, 12][inds_resolve], "ground_elevation": resolve[:, 11][inds_resolve], "dobs": dobs_re, "dpred": dpred_re, "mopt": mopt_re, "z": mesh.cell_centers_z[active], "frequency_cp": np.r_[382, 1822, 7970, 35920, 130100], "frequency_cx": np.r_[3258.0], "river_path": river_path, } area = 314.0 waveform = np.loadtxt(os.path.sep.join([directory, "skytem_hm.wf"])) times = np.loadtxt(os.path.sep.join([directory, "skytem_hm.tc"])) booky_skytem = { "data": data_skytem, "data_header": header_skytem[17 : 17 + 24], "line": lines_skytem, "xy": xyz_skytem, "src_elevation": info_skytem[:, 10].astype(float), "ground_elevation": info_skytem[:, 15].astype(float), "area": area, "radius": np.sqrt(area / np.pi), "t0": 0.01004, "waveform": waveform, "times": times, } if saveIt: save_dict_to_hdf5( os.path.sep.join([directory, "booky_resolve.hdf5"]), booky_resolve ) save_dict_to_hdf5( os.path.sep.join([directory, "booky_skytem.hdf5"]), booky_skytem ) if cleanup: os.remove(downloads) shutil.rmtree(directory) if __name__ == "__main__": run(plotIt=True, saveIt=False, cleanup=False) .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 1.632 seconds) **Estimated memory usage:** 9 MB .. _sphx_glr_download_content_examples_20-published_plot_load_booky.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: plot_load_booky.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_load_booky.py ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_