\u001b[0m\n",
" \u001b[36m │ └ \u001b[0m\u001b[36m\u001b[1mSpectrogram:\u001b[0m\n",
" \u001b[36m │ \u001b[0m\u001b[36m\u001b[1m3 harmonics, 0.0078125Hz spaced \u001b[0m\n",
" \u001b[36m │ \u001b[0m\u001b[36m\u001b[1m from 0.03125 to 0.046875 Hz.\u001b[0m\n",
" \u001b[36m │ \u001b[0m\u001b[36m\u001b[1m\u001b[0m\n",
" \u001b[36m │ \u001b[0m\u001b[36m\u001b[1m416 Time observations \u001b[0m\n",
" \u001b[36m │ \u001b[0m\u001b[36m\u001b[1mStart: 1980-01-01T00:00:0...\u001b[0m\n",
" \u001b[36m └ \u001b[0m\u001b[36m\u001b[1m{\u001b[0m\n",
" \u001b[36m \u001b[0m\u001b[36m\u001b[1m \"band\": {\u001b[0m\n",
" \u001b[36m \u001b[0m\u001b[36m\u001b[1m \"center_averaging_type\": \"geometric\",\u001b[0m\n",
" \u001b[36m \u001b[0m\u001b[36m\u001b[1m \"closed\": \"left\",\u001b[0m\n",
" \u001b[36m \u001b[0m\u001b[36m\u001b[1m \"decimation_level\": 0,\u001b[0m\n",
" \u001b[36m \u001b[0m\u001b[36m\u001b[1m ...\u001b[0m\n",
"\n",
"\u001b[31m\u001b[1mAttributeError\u001b[0m:\u001b[1m 'Band' object has no attribute 'bands'\u001b[0m\n"
]
}
],
"source": [
"show_plot = True\n",
"tf_cls = process_mth5(config,\n",
" kernel_dataset,\n",
" units=\"MT\",\n",
" show_plot=show_plot,\n",
" z_file_path=None,\n",
" )"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Export TF to a file "
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"ename": "AttributeError",
"evalue": "'NoneType' object has no attribute 'write'",
"output_type": "error",
"traceback": [
"\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
"\u001b[0;31mAttributeError\u001b[0m Traceback (most recent call last)",
"Cell \u001b[0;32mIn[10], line 2\u001b[0m\n\u001b[1;32m 1\u001b[0m xml_file_base \u001b[38;5;241m=\u001b[39m \u001b[38;5;124mf\u001b[39m\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124msynthetic_test1.xml\u001b[39m\u001b[38;5;124m\"\u001b[39m\n\u001b[0;32m----> 2\u001b[0m \u001b[43mtf_cls\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mwrite\u001b[49m(fn\u001b[38;5;241m=\u001b[39mxml_file_base, file_type\u001b[38;5;241m=\u001b[39m\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124memtfxml\u001b[39m\u001b[38;5;124m\"\u001b[39m)\n",
"\u001b[0;31mAttributeError\u001b[0m: 'NoneType' object has no attribute 'write'"
]
}
],
"source": [
"xml_file_base = f\"synthetic_test1.xml\"\n",
"tf_cls.write(fn=xml_file_base, file_type=\"emtfxml\")\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"edi_file_base = f\"synthetic_test1.edi\"\n",
"tf_cls.write(fn=edi_file_base, file_type=\"edi\")\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
" NOTE Fourier coefficient section below here is a work in progress.
"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Fourier coefficient storage in MTH5\n",
"\n",
"The capability to store Fourier coeffficients (FCs) in MTH5 is now available. This will enable some different approaches to processing and data quality control (QC). The data QC tools are a work in progress. The following examples show how to add some FC levels to an MTH5, providing a starting point for processing or feature extraction from these data.\n",
"\n",
"There are currently two main ways to add FCs to the MTH5.\n",
"\n",
"1. Store on the fly while processing with Aurora\n",
"2. Use a dedicated method to make FCs only."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Storing while processing with Aurora\n",
"\n",
"We can store FCs while processing by changing the processing config's `save_fcs` attribute to True \n",
"\n",
"If we do that with the current processing config however we encouter a warning:\n",
"\n",
"```\n",
"Saving FCs for remote reference processing is not supported \n",
" - To save FCs, process as single station, then you can use the FCs for RR processing \n",
" - See issue #319 for details \n",
" - forcing processing config save_fcs = False \n",
"```\n",
"\n",
"There are two workarounds for this:\n",
"1. Process each station as a single station, and then the FCs (TODO: Add an example of this)\n",
"2. Explicitly call `add_fcs` function.\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"file_size_before_adding_fcs = mth5_path.stat().st_size\n",
"print(f\"file_size_before_adding_fcs: {file_size_before_adding_fcs}\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Using the \"single-station processing method\" of FC generation"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# set kernel_dataset to work on one station only\n",
"kernel_dataset = KernelDataset()\n",
"kernel_dataset.from_run_summary(run_summary, \"test1\", None) \n",
"\n",
"# Update the processing config\n",
"\n",
"cc = ConfigCreator()\n",
"config = cc.create_from_kernel_dataset(kernel_dataset)\n",
"\n",
"# tell the config to save the FCs\n",
"for dec in config.decimations:\n",
" dec.save_fcs = True\n",
" dec.save_fcs_type = \"h5\"\n",
"\n",
"# you can export the config to a json by uncommenting the following line\n",
"# cfg_json = config.to_json()\n",
"kernel_dataset.mini_summary"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"show_plot = False\n",
"tf_cls = process_mth5(config,\n",
" kernel_dataset,\n",
" units=\"MT\",\n",
" show_plot=show_plot,\n",
" z_file_path=None,\n",
" )"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"file_size_after_adding_fcs_station_1 = mth5_path.stat().st_size\n",
"print(f\"file_size_after_adding_fcs_station_1: {file_size_after_adding_fcs_station_1}\")"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Now the other station\n",
"\n",
"# set kernel_dataset to work on one station only\n",
"kernel_dataset = KernelDataset()\n",
"kernel_dataset.from_run_summary(run_summary, \"test2\", None) \n",
"\n",
"# Update the processing config\n",
"\n",
"cc = ConfigCreator()\n",
"config = cc.create_from_kernel_dataset(kernel_dataset)\n",
"\n",
"# tell the config to save the FCs\n",
"for dec in config.decimations:\n",
" dec.save_fcs = True\n",
" dec.save_fcs_type = \"h5\"\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"show_plot = False\n",
"tf_cls = process_mth5(config,\n",
" kernel_dataset,\n",
" units=\"MT\",\n",
" show_plot=show_plot,\n",
" z_file_path=None,\n",
" )"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"Now, the FCs are there (for the specific processing configuration).\n",
"\n",
"If you were to re-process the data, the FCs are already there, "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"file_size_after_adding_fcs_station_2 = mth5_path.stat().st_size\n",
"print(f\"file_size_after_adding_fcs_station_2: {file_size_after_adding_fcs_station_2}\")"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# set kernel_dataset to work on both stations again\n",
"kernel_dataset = KernelDataset()\n",
"kernel_dataset.from_run_summary(run_summary, \"test1\", \"test2\") \n",
"\n",
"# Update the processing config\n",
"\n",
"cc = ConfigCreator()\n",
"config = cc.create_from_kernel_dataset(kernel_dataset)\n",
"\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"\n",
"from aurora.pipelines.transfer_function_kernel import TransferFunctionKernel\n",
"tfk = TransferFunctionKernel(dataset=kernel_dataset, config=config)\n",
"tfk.update_processing_summary()\n",
"tfk.check_if_fcs_already_exist()\n",
"tfk.dataset_df"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"show_plot = False\n",
"tf_cls = process_mth5(config,\n",
" kernel_dataset,\n",
" units=\"MT\",\n",
" show_plot=show_plot,\n",
" z_file_path=None,\n",
" )"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Using Dedicated FC Generation (`add_fcs` function)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"from mth5.timeseries.spectre.helpers import add_fcs_to_mth5\n",
"from mth5.timeseries.spectre.helpers import fc_decimations_creator\n",
"from mth5.timeseries.spectre.helpers import read_back_fcs"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Build a fresh copy of the synthetic data file"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"mth5_path.unlink()\n",
"create_test12rr_h5(target_folder=target_folder)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In this case we need to prescribe the decimation configuration.\n",
"\n",
"If we dont want to decimate, we can just pass \n",
"`fc_decimations == \"degenerate\"` to add_fcs\n",
"\n",
"Here is how the decimations can be created"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"sample_rate = 1.0\n",
"fc_decimations = fc_decimations_creator(initial_sample_rate=sample_rate, time_period=None)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"And here is what they look like, ... a list, with each element specifying the\n",
"needed info for decimation and FFT."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"fc_decimations\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now add the decimation levels:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"add_fcs_to_mth5(mth5_path, fc_decimations=fc_decimations)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now take a look at the mth5"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"import mth5\n",
"m = mth5.mth5.MTH5(mth5_path)\n",
"m.open_mth5()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"m.close_mth5()\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Alternatively, we could have translated the Aurora processing config to an \n",
"mt_metadata FC decimation and gnerated FCs with the same function."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"mth5_path.unlink()\n",
"create_test12rr_h5(target_folder=target_folder)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Extract FC decimations from processing config and build the layer\n",
"\n",
"fc_decimations = [\n",
" x.to_fc_decimation() for x in config.decimations\n",
"]\n",
"\n",
"add_fcs_to_mth5(mth5_path, fc_decimations=fc_decimations)\n",
" "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"\n",
"read_back_fcs(mth5_path)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"And we can see that when processing, aurora detects that the FC layer is already there, and skips creating it."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"show_plot = False\n",
"tf_cls = process_mth5(config,\n",
" kernel_dataset,\n",
" units=\"MT\",\n",
" show_plot=show_plot,\n",
" z_file_path=None,\n",
" )"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
" NOTE Feature storage section below is a work in progress.
"
]
},
{
"cell_type": "markdown",
"metadata": {
"tags": []
},
"source": [
"\n",
"# Feature Storage **Experimental work in progress**\n",
"\n",
"### Make an example of a feature embedded in mth5 file from the FCs.\n",
"\n",
"Any spectrogram (FC array at some decimation level) can be passed to a feature extraction method.\n",
"These features can, in turn be stored in the mth5.\n",
"\n",
"This will have the following steps:\n",
"\n",
"0. select the source FC decimation level (Will use \"0\" to start, but also test with \"1\",)\n",
"\n",
"1. Select the frequency bands (probably make these wider than normal so we have plenty of FCs per time window)\n",
"\n",
"- Note these are `feature_extraction_bands`, which are possibly different from TF estimation bands.\n",
"\n",
"To do this we will want to look at the frequencies associated with each FCGroup, and each decimation level.\n",
"- Group Runs on Sample Rate\n",
" - tabulate decimation levels\n",
" - group decimation levels by id (and validate that the decimation_level.metadata.sample_rate_decimations agree)\n",
" - for each decimation level, we should get the frequencies associated -- this can be deduced from the shape of the dataset, and\n",
" the min/max frequecies of the metadata, but for now we can just unpack the xarray and take the frequncey axis)\n",
"\n",
"2. loop over the bands (for each band:)\n",
"\n",
"- 2a. extract the band from the FClayer in the mth5 (shape will be nch x ntime x nharmonics)\n",
"\n",
"3. loop over time window (for each time):\n",
"\n",
"- 3a) extract feature\n",
"\n",
"As a placeholder, we will compute the channel crosspowers, and or channel energy.\n",
"\n",
"- 3b) place the feature into an xarray\n",
"\n",
"4) save teh xarray as \"feature_cross_power\"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Imports \n",
"\n",
"from mth5.mth5 import MTH5\n",
"\n",
"import numpy as np\n",
"import pandas as pd"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# # Selection of frequency bands \n",
"\n",
"# active_decimation_level = \"0\"\n",
"\n",
"# m = MTH5(mth5_path)\n",
"# m.open_mth5()\n",
"# run_summary_df = m.run_summary # This fails on current main but works on add_aurora_tools branch.\n",
"# run_summary_df"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"It would be nice if we could get a summary of the fcs available in the mth5. here is a prototype for how that could be done.\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"tags": []
},
"outputs": [],
"source": [
"# # Prototype for fc_summary\n",
"\n",
"# fc_summary = None\n",
"\n",
"# for i, row in run_summary_df.iterrows():\n",
"# print(row.station, row.run)\n",
"# station_obj = m.get_station(row.station)\n",
"# run_fc_group = station_obj.fourier_coefficients_group.get_fc_group(row.run)\n",
"# if fc_summary is None:\n",
"# fc_summary = run_fc_group.decimation_level_summary\n",
"# fc_summary[\"station\"] = row.station\n",
"# fc_summary[\"run\"] = row.run\n",
"# fc_summary[\"sample_rate_time_series\"] = 0.0\n",
"# fc_summary[\"window_num_samples\"] = 0\n",
"# fc_summary[\"harmonic_indices\"] = None\n",
"# else:\n",
"# tmp = run_fc_group.decimation_level_summary\n",
"# tmp[\"station\"] = row.station\n",
"# tmp[\"run\"] = row.run\n",
"# tmp[\"sample_rate_time_series\"] = 0.0\n",
"# tmp[\"window_num_samples\"] = 0\n",
"# tmp[\"harmonic_indices\"] = None\n",
"# fc_summary = pd.concat([fc_summary, tmp], ignore_index=True)\n",
" \n",
"# #print(run_fc_group)\n",
"# for dec_level_id in run_fc_group.groups_list:\n",
"# dec_level = run_fc_group.get_decimation_level(dec_level_id)\n",
"# cond1 = fc_summary.station==row.station\n",
"# cond2 = fc_summary.run==row.run\n",
"# cond3 = fc_summary.component==dec_level_id\n",
"# fc_summary.loc[cond1 & cond2 & cond3] # this should be a unique row\n",
"# ndx = fc_summary.loc[cond1 & cond2 & cond3].index\n",
"# assert(len(ndx))==1\n",
"# #fc_summary[\"sample_rate_time_series\"].at[ndx] = dec_level.metadata.i\n",
"# fc_summary[\"window_num_samples\"].at[ndx[0]] = dec_level.metadata.window.num_samples\n",
"# fc_summary[\"sample_rate_time_series\"].at[ndx[0]] = dec_level.metadata.sample_rate\n",
"# fc_summary[\"harmonic_indices\"].at[ndx[0]] = dec_level.metadata.harmonic_indices\n",
" \n",
" \n",
"# #print(dec_level_id, dec_level)\n",
"# xr_dec_level = dec_level.to_xarray()\n",
"# print(\"station\", row.station, \"run\", row.run, \"decimation level\", dec_level_id, \"SAMPLE rate\", dec_level.metadata.sample_rate) # xr_dec_level.frequency)\n",
" "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# fc_summary[\"delta_f\"] = fc_summary[\"sample_rate_time_series\"] / fc_summary[\"window_num_samples\"]\n",
"# fc_summary"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The above table tells us what we need to know to get frequency axes .. and actually, we can make `_frequencies()` a method of `FC_dec_level`.\n",
"\n",
"Note that when harmonic_indices is [-1,], that means all the (1-sided) FCs are kept.\n",
"- this can be sanity checked by asserting that the dataset size is `window_num_samples`//2\n",
"\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"So now pick a decimation level, ... In this case, the `component` column is telling us the decimation level.\n",
"- It would seem sensible to make sure that the sample rates, (and delta_f) values are same as well.\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# cond0 = fc_summary[\"component\"] == active_decimation_level\n",
"# fc_runs_to_featureize_df = fc_summary[cond0]\n",
"# fc_runs_to_featureize_df"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# frqs = np.fft.fftfreq?"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now pick a scheme for defining the feature frequency bands\n",
"- we prefer not to use the lowest few harmonics ... this is like setting the minimum number of cycles per window, say 10\n",
"- we also prefer not to take harmonics very close to the Nyquist frequency as these maybe attenuated by AAF, say we go up to 80% Nyq\n",
"- For TF estimation, bands may be only a few FCs wide,but for feature extraction we make them much wider.\n",
"- lets try with around half an octave\n",
"- how about 2 bands per octave, mean"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# min_cycles = 10\n",
"# nyquist_fraction = 0.8\n",
"# n_feature_bands_per_octave = 2"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# for i, row in fc_runs_to_featureize_df.iterrows():\n",
"# frequencies = np.fft.rfftfreq(row.window_num_samples, row.sample_rate_time_series)\n",
"# lower_index = int(min_cycles)\n",
"# upper_index = int(nyquist_fraction*(row.window_num_samples//2))\n",
"# print(f\" lower_index={lower_index}, upper_index={upper_index}\")\n",
"# freq_min = frequencies[lower_index]\n",
"# freq_max = frequencies[upper_index]\n",
"# print(f\"Available band: {freq_min:3f} to {freq_max:3f}\")\n",
"# frequency_ratio = freq_max/freq_min\n",
"# available_decades = np.log10(frequency_ratio)\n",
"# available_ocataves = np.log2(frequency_ratio)\n",
"# print(f\"Available decades: {available_decades:3f}\")\n",
"# print(f\"Available octaves: {available_ocataves:3f}\")\n",
"# num_feature_bands = available_ocataves * n_feature_bands_per_octave"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# dec_level.metadata\n",
"# dec_level.to_xarray()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# cond1 = fc_summary.station==\"test1\"\n",
"# cond2 = fc_summary.run==\"001\"\n",
"# cond3 = fc_summary.component==\"0\"\n",
"# fc_summary.loc[cond1 & cond2 & cond3].index"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# fc_summary[\"sample_rate_time_series\"].at[0]"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# m.close_mth5()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "aurora-test",
"language": "python",
"name": "aurora-test"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.9.19"
}
},
"nbformat": 4,
"nbformat_minor": 4
}