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"]}

# 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)

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)

b8468fbb-3e5b-4b88-b268-2dd8d80798fa /survey
df49f8a4-c5d3-48da-965e-d36f7d219104 /survey/nominal_system
06ce1c1c-c25c-42e3-9179-924bf2e7a7d8 /survey/magnetic_system
9584e8b0-2120-479b-b9ce-1d823ff86920 /survey/data/raw_data
b9987735-d1a4-420b-a695-d79385997cf3 /survey/data/processed_data
df49f8a4-c5d3-48da-965e-d36f7d219104 /survey/data/raw_data/nominal_system
df49f8a4-c5d3-48da-965e-d36f7d219104 /survey/data/processed_data/skytem_system
06ce1c1c-c25c-42e3-9179-924bf2e7a7d8 /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)

b8468fbb-3e5b-4b88-b268-2dd8d80798fa /survey
df49f8a4-c5d3-48da-965e-d36f7d219104 /survey/nominal_system
06ce1c1c-c25c-42e3-9179-924bf2e7a7d8 /survey/magnetic_system
9584e8b0-2120-479b-b9ce-1d823ff86920 /survey/data/raw_data
b9987735-d1a4-420b-a695-d79385997cf3 /survey/data/processed_data
ec0cee42-b500-4882-845a-ba43d78df136 /survey/data/depth_to_bedrock
eb9e16e0-bf2d-4470-8473-b0f4ebfcb4cd /survey/models/inverted_models
dbfb5f8c-bd97-4e75-89ba-6c9e633f210d /survey/derived_products/maps
df49f8a4-c5d3-48da-965e-d36f7d219104 /survey/data/raw_data/nominal_system
df49f8a4-c5d3-48da-965e-d36f7d219104 /survey/data/processed_data/skytem_system
06ce1c1c-c25c-42e3-9179-924bf2e7a7d8 /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.067 seconds)

Gallery generated by Sphinx-Gallery