Multi-dataset Survey

This example demonstrates the typical workflow for creating a GS file for an AEM survey in its entirety, i.e., the NetCDF file contains all related datasets together, e.g., raw data, processed data, inverted models, and derivative products. Specifically, this survey contains:

Minimally processed (raw) AEM data and raw/processed magnetic data provided by SkyTEM
Fully processed AEM data used as input to inversion
Laterally constrained inverted resistivity models
Point-data estimates of bedrock depth derived from the AEM models
Interpolated magnetic and bedrock depth grids

Note: To make the size of this example more managable, some of the input datasets have been downsampled relative to the source files in the data release referenced below.

Dataset Reference: Minsley, B.J, Bloss, B.R., Hart, D.J., Fitzpatrick, W., Muldoon, M.A., Stewart, E.K., Hunt, R.J., James, S.R., Foks, N.L., and Komiskey, M.J., 2022, Airborne electromagnetic and magnetic survey data, northeast Wisconsin (ver. 1.1, June 2022): U.S. Geological Survey data release, https://doi.org/10.5066/P93SY9LI.

import matplotlib.pyplot as plt
from os.path import join
import numpy as np
import gspy
from gspy import Survey
import xarray as xr
from pprint import pprint

Convert the Skytem csv data to NetCDF

Initialize the Survey

# Path to example files
data_path = '..//data_files//skytem_csv'

# Survey metadata file
metadata = join(data_path, "data//skytem_survey.yml")

# Establish the Survey
survey = Survey.from_dict(metadata)

data_container = survey.gs.add_container('data', **dict(content = "raw and processed data",
                                                        comment = "This is a test"))

1 - Raw Data - Import raw AEM data from CSV-format. Define input data file and associated metadata file

d_data1 = join(data_path, 'data//skytem_contractor_data.csv')
d_supp1 = join(data_path, 'data//skytem_contractor_data.yml')

# Add the raw AEM data as a tabular dataset
data_container.gs.add(key='raw_data', data_filename=d_data1, metadata_file=d_supp1, system=survey.nominal_system)

2 - Processed Data - Import processed AEM data from CSV-format. Define input data file and associated metadata file

d_data2 = join(data_path, 'data//skytem_processed_data.csv')
d_supp2 = join(data_path, 'data//skytem_processed_data.yml')

system = {"skytem_system" : survey["nominal_system"].isel(lm_gate_times=np.s_[1:], hm_gate_times=np.s_[10:]),
          "magnetic_system" : survey["magnetic_system"]}

sys =  system["skytem_system"].to_dataset()
print(sys.isel({"n_receiver": 0}))

# Add the processed AEM data as a tabular dataset
pd = data_container.gs.add(key='processed_data', data_filename=d_data2, metadata_file=d_supp2, system=system)

<xarray.Dataset> Size: 4kB
Dimensions:                                      (gate_times: 22, nv: 2,
                                                  lm_gate_times: 27,
                                                  hm_gate_times: 22,
                                                  n_loop_vertices: 8, xyz: 3,
                                                  n_transmitter: 2,
                                                  transmitter_lm_waveform_time: 21,
                                                  transmitter_hm_waveform_time: 36,
                                                  n_component: 4)
Coordinates:
  * gate_times                                   (gate_times) float64 176B 5....
  * nv                                           (nv) int64 16B 0 1
  * lm_gate_times                                (lm_gate_times) float64 216B ...
  * hm_gate_times                                (hm_gate_times) float64 176B ...
  * n_loop_vertices                              (n_loop_vertices) int64 64B ...
  * xyz                                          (xyz) int64 24B 0 1 2
  * n_transmitter                                (n_transmitter) int64 16B 0 1
  * transmitter_lm_waveform_time                 (transmitter_lm_waveform_time) float64 168B ...
  * transmitter_hm_waveform_time                 (transmitter_hm_waveform_time) float64 288B ...
    n_receiver                                   int64 8B 0
  * n_component                                  (n_component) int64 32B 0 1 2 3
Data variables: (12/31)
    gate_times_bnds                              (gate_times, nv) float64 352B ...
    lm_gate_times_bnds                           (lm_gate_times, nv) float64 432B ...
    hm_gate_times_bnds                           (hm_gate_times, nv) float64 352B ...
    n_loop_vertices_bnds                         (n_loop_vertices, nv) float64 128B ...
    xyz_bnds                                     (xyz, nv) float64 48B -0.5 ....
    transmitter_label                            (n_transmitter) <U2 16B 'LM'...
    ...                                           ...
    component_receivers                          (n_component) <U1 16B 'z' .....
    component_txrx_dx                            (n_component) float64 32B -1...
    component_txrx_dy                            (n_component) float64 32B 0....
    component_txrx_dz                            (n_component) float64 32B -2...
    component_data_type                          (n_component) <U4 64B 'dBdt'...
    component_gate_times                         (n_component) <U13 208B 'lm_...
Attributes:
    type:                      system
    mode:                      airborne
    method:                    electromagnetic
    submethod:                 time domain
    instrument:                skytem 304M
    name:                      nominal_system
    data_normalized:           True
    skytem_skb_gex_available:  True
    reference_frame:           right-handed positive down
    coil_orientations:         X, Z

3 - Inverted Models -

# Create a new container for models
model_container = survey.gs.add_container('models', **dict(content = "Inverted models",
                                                          comment = "This is a test"))

# pprint(survey.gs.get_all_attr('standard_name'))
print(survey.gs.tree)

# Import inverted AEM models from CSV-format.
# Define input data file and associated metadata file
m_data3 = join(data_path, 'model//skytem_inverted_models.csv')
m_supp3 = join(data_path, 'model//skytem_inverted_models.yml')

# Add the inverted AEM models as a tabular dataset
model_container.gs.add(key='inverted_models', data_filename=m_data3, metadata_file=m_supp3)

/survey
/survey/nominal_system
/survey/magnetic_system
/survey/data
/survey/models
/survey/data/raw_data
/survey/data/processed_data
/survey/data/raw_data/nominal_system
/survey/data/processed_data/skytem_system
/survey/data/processed_data/magnetic_system

4 - Bedrock Picks - Import AEM-based estimated of depth to bedrock from CSV-format. Define input data file and associated metadata file

d_data4 = join(data_path, 'data//top_dolomite_blocky_lidar.csv')
d_supp4 = join(data_path, 'data//bedrock_picks.yml')

# Add the AEM-based estimated of depth to bedrock as a tabular dataset
data_container.gs.add(key='depth_to_bedrock', data_filename=d_data4, metadata_file=d_supp4)

5 - Derivative Maps -

# We can add arbitrarily named containers to the survey
derived_products = survey.gs.add_container('derived_products', **dict(content = "products derived from other data and models"))

# Import interpolated bedrock and magnetic maps from TIF-format.
# Define input metadata file (which contains the TIF filenames linked to variable names)
m_supp5 = join(data_path, 'data//magnetics_bedrock_picks.yml')

# Add the interpolated maps as a raster dataset
derived_products.gs.add(key='maps', metadata_file=m_supp5)

print(survey.gs.tree)

/survey
/survey/nominal_system
/survey/magnetic_system
/survey/data
/survey/models
/survey/derived_products
/survey/data/raw_data
/survey/data/processed_data
/survey/data/depth_to_bedrock
/survey/models/inverted_models
/survey/derived_products/maps
/survey/data/raw_data/nominal_system
/survey/data/processed_data/skytem_system
/survey/data/processed_data/magnetic_system

Save to NetCDF file

d_out = join(data_path, 'skytem.nc')
survey.gs.to_netcdf(d_out)

The gspy goal is to have the complete survey in a single file. However, we can also save containers or datasets separately.

data_container.gs.to_netcdf(join(data_path, 'test_datacontainer.nc'))

Reading back in

new_survey = gspy.open_datatree(d_out)['survey']

# print(new_survey)

# #%%
# # Plotting
# plt.figure()
# new_survey['data']['raw_data']['height'].plot()
# plt.tight_layout()

# pd = new_survey['data']['processed_data']
# plt.figure()
# pd['elevation'].plot()
# plt.tight_layout()

# m = new_survey['derived_products']['maps']
# plt.figure()
# m['magnetic_tmi'].plot(cmap='jet')
# plt.tight_layout()

# plt.show()

Total running time of the script: (0 minutes 1.105 seconds)

Gallery generated by Sphinx-Gallery