Commit 281fd7c3 authored by doetschj's avatar doetschj
Browse files
parents cbe28a9f 82ce7057
......@@ -4,20 +4,28 @@ Create a new file when needed.
import time
import logging
import os
import io
from ctypes import c_int32
from math import floor
from obspy.core import Stream, Trace, UTCDateTime
from obspy.core.inventory import Longitude, Latitude
import pyasdf
from dug_seis.acquisition.flat_response_stationxml import get_flat_response_inventory
class DataToASDF:
def __init__(self, param):
self.folder = param['Acquisition']['asdf_settings']['data_folder']
self.folder_tmp = param['Acquisition']['asdf_settings']['data_folder_tmp']
if self.folder[len(self.folder)-1] is not "/":
self.folder += "/"
if self.folder_tmp[len(self.folder_tmp)-1] is not "/":
self.folder_tmp += "/"
self.filename = param['General']['project_name']
self.compression = param['Acquisition']['asdf_settings']['compression']
self.file_length_sec = param['Acquisition']['asdf_settings']['file_length_sec']
......@@ -35,6 +43,8 @@ class DataToASDF:
self._sampling_rate = param['Acquisition']['hardware_settings']['sampling_frequency']
self._gain_selection = param['Acquisition']['hardware_settings']['gain_selection']
self._sensor_coords = param['General']['sensor_coords']
self._sensor_count = param['General']['sensor_count']
self._nr_of_datapoints = floor(self.l_notify_size.value / 16 / 2) # nr of channels & 16 bit = 2 bytes
self._file_handle = None
self._time_age_of_file = 0 # keeps track internally how old the file is
......@@ -57,6 +67,7 @@ class DataToASDF:
file_name = "{0}_{1}.h5".format(UTCDateTime(), self.filename)
file_name = file_name.replace(":", "_")
folder_file_name = "{0}{1}".format(self.folder_tmp, file_name)
# print("_create_new_file with folder_file_name = {0}".format(folder_file_name))
logging.info("_create_new_file with folder_file_name = {0}".format(folder_file_name))
self._time_age_of_file = time.time()
......@@ -71,6 +82,8 @@ class DataToASDF:
os.rename(self.folder_tmp + self._last_used_file_name, self.folder + self._last_used_file_name)
self._last_used_file_name = file_name
self._add_all_station_xml_s(self._file_handle)
def _create_new_file_if_needed(self):
# check if file_handle exists
if self._file_handle is None:
......@@ -83,6 +96,40 @@ class DataToASDF:
# logging.info("file too old, creating new file.\n")
self._create_new_file()
def _add_all_station_xml_s(self, ds):
for i in range(self._sensor_count):
ds.add_stationxml(self._create_station_xml(i))
def _create_station_xml(self, channel_nr):
# print("channel_nr = {}, len(self._sensor_coords) = {}".format(channel_nr, len(self._sensor_coords)))
inv = get_flat_response_inventory(
sensitivity_value=self._gain_selection[channel_nr] * 2 / 65536, # conversion 16bit int to mV
sensitivity_frequency=1.0,
input_units="M/S", # ?
output_units="Counts", # ?
creation_date=UTCDateTime(self._time_age_of_file),
network_code=self.stats['network'],
station_code=self._get_station_name(channel_nr),
location_code=self.stats['location'],
channel_code=self.stats['channel'],
sampling_rate=self.stats['sampling_rate'],
latitude=Latitude(self._sensor_coords[channel_nr][1]), # from here on ?
longitude=Longitude(self._sensor_coords[channel_nr][0]),
depth=self._sensor_coords[channel_nr][2],
elevation=0.0,
azimuth=0.0,
dip=0.0)
# Test if the response makes up a valid StationXML file.
with io.BytesIO() as buf:
inv.write(buf, format="stationxml", validate=True)
return inv
def _get_station_name(self, channel_nr):
return str(self.station_naming[channel_nr]).zfill(3)
def data_to_asdf(self, np_data_list):
time_start_buffer_data_available = time.time()
......@@ -94,8 +141,7 @@ class DataToASDF:
time_tmp = time.time()
for i in range(16):
self.stats['station'] = str(self.station_naming[i + 16*card_nr]).zfill(3)
self.stats['gain'] = self._gain_selection[i + 16*card_nr]
self.stats['station'] = self._get_station_name(i + 16 * card_nr)
stream += Trace(np_data[i], header=self.stats) # without transpose np_data[:, i]
# logging.info("{}, {}\n".format(self.stats['station'], self.stats['starttime']))
......
# example from lion HowTo generate a StationXML file
#
# ObsPy StationXML: https://docs.obspy.org/tutorial/code_snippets/stationxml_file_from_scratch.html
# pyasdf adding StationXML: http://seismicdata.github.io/pyasdf/tutorial.html#adding-station-information
# pyasdf reading StationXML: http://seismicdata.github.io/pyasdf/tutorial.html#reading-waveforms-and-stationxml
import obspy
from obspy.core.inventory import (
Inventory, Network, Station, Channel, Site, Response,
InstrumentSensitivity, PolesZerosResponseStage, Longitude, Latitude)
def get_flat_response_inventory(
sensitivity_value: float, sensitivity_frequency: float,
input_units: str, output_units: str,
sampling_rate: float,
creation_date: obspy.UTCDateTime,
network_code: str, station_code: str,
location_code: str, channel_code: str,
latitude: Latitude, longitude: Longitude, elevation: float, depth: float,
azimuth: float, dip: float,
site_name: str = "site",
source_str: str = "self"):
_inv = Inventory(
networks=[],
source=source_str)
net = Network(
code=network_code,
stations=[])
sta = Station(
code=station_code,
latitude=latitude,
longitude=longitude,
elevation=elevation,
site=Site(name=site_name),
creation_date=creation_date)
cha = Channel(
code=channel_code,
location_code=location_code,
latitude=latitude,
longitude=longitude,
elevation=elevation,
depth=depth,
azimuth=azimuth,
dip=dip,
sample_rate=sampling_rate)
response = Response(
instrument_sensitivity=InstrumentSensitivity(
value=sensitivity_value, frequency=sensitivity_frequency,
input_units=input_units, output_units=output_units),
response_stages=[PolesZerosResponseStage(
stage_sequence_number=1, stage_gain=sensitivity_value,
stage_gain_frequency=sensitivity_frequency,
input_units=input_units, output_units=output_units,
pz_transfer_function_type="LAPLACE (RADIANS/SECOND)",
normalization_frequency=sensitivity_frequency,
zeros=[], poles=[]
)])
# Now tie it all together.
cha.response = response
sta.channels.append(cha)
net.stations.append(sta)
_inv.networks.append(net)
return _inv
#!/usr/bin/env python
import yaml
with open("config/dug-seis.yaml", 'r') as stream:
try:
param = yaml.safe_load(stream)
print(param)
except yaml.YAMLError as exc:
print(exc)
# example from lion HowTo generate a StationXML file
#
# ObsPy StationXML: https://docs.obspy.org/tutorial/code_snippets/stationxml_file_from_scratch.html
# pyasdf adding StationXML: http://seismicdata.github.io/pyasdf/tutorial.html#adding-station-information
# pyasdf reading StationXML: http://seismicdata.github.io/pyasdf/tutorial.html#reading-waveforms-and-stationxml
import io
import obspy
from obspy.core.inventory import (Longitude, Latitude)
from dug_seis.acquisition.flat_response_stationxml import get_flat_response_inventory
tr = obspy.read()[0]
inv = get_flat_response_inventory(
sensitivity_value=200.0,
sensitivity_frequency=1.0,
input_units="M/S",
output_units="Counts",
creation_date=obspy.UTCDateTime(2012, 1, 1),
network_code=tr.stats.network,
station_code=tr.stats.station,
location_code=tr.stats.location,
channel_code=tr.stats.channel,
sampling_rate=tr.stats.sampling_rate,
latitude=Latitude(0.0),
longitude=Longitude(0.0),
depth=0.0,
elevation=0.0,
azimuth=0.0,
dip=0.0)
# Test if the response makes up a valid StationXML file.
with io.BytesIO() as buf:
inv.write(buf, format="stationxml", validate=True)
inv.plot_response(0.001)
# Now plot the trace once before and once after
# removing the response.
tr.plot()
tr.remove_response(inventory=inv)
tr.plot()
# end lion
# test file generation
import pyasdf
ds = pyasdf.ASDFDataSet("test_file.h5", compression="gzip-3")
print(ds)
ds.add_waveforms(tr, tag="raw_recording")
print(ds)
ds.add_stationxml(inv)
print(ds.waveforms.BW_RJOB)
print(ds.waveforms.BW_RJOB.raw_recording)
print(ds.waveforms.BW_RJOB.StationXML)
# ds.waveforms.XB_001.StationXML.plot_response(1)
\ No newline at end of file
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment