#!/usr/bin/env python import h5py import numpy as np import matplotlib.pyplot as plt import matplotlib.colors as mcolors import datetime import pytz import obspy import io import os import sys import glob import pandas as pd from obspy import read, Trace, Stream, UTCDateTime from tqdm import tqdm,trange def file_name(data_path, start): jday = start.strftime('%j') year = start.strftime('%Y') month = start.strftime('%m') day = start.strftime('%d') hour = str(start)[11:13] minute = str(start)[14:16] second = str(start)[17:19] time = year+"-"+month+"-"+day+"T"+hour+minute path = data_path+year+"/"+year+"."+jday+"/MBARI_GL20m_OCP5m_FS200Hz_"+time+"*Z.h5" real_path = str(glob.glob(path))[2:-2] return real_path """ Original code from Peter Hubbard The conversion of optical phase shift is in details described in this SEAFOM report : https://seafom.com/mdocs-posts/seafom_msp_02-august-2018-pdf , as shown in equations 9 and 10 on page 26. """ def optic_phase_to_starin(obspy_st,gauge_length, refractive_index=1.468, optical_wavelength=1550, pe_scaling_factor = 0.78): #-------------------- # This function is in-place operation in order to save memory, make a copy before using # # obspy_st : Obspy Stream object # gauge_length : gauge length in m # refractive_index : refractive index which depends on type of fiber # optical_wavelength : the operational optical (vacuum) wavelength of the DAS system # pe_scaling_factor : the photo-elastic scaling factor for longitudinal strain in isotropic material #-------------------- channels = len(obspy_st) record_length = len(obspy_st[0].data) strain = np.zeros((channels,record_length)) overscale = [] Max_amp = (2147483647/10430.378850470453) for i in trange(channels): data = obspy_st[i].data sig_phi = data #radians if max(sig_phi)-min(sig_phi) > 1.9*Max_amp: ## check overscale.append(i) sign0 = np.sign(sig_phi[0]) asign = (np.sign(sig_phi)-1*sign0)*Max_amp sig_phi = sig_phi + asign # nm epsilon sig_ep2 = optical_wavelength*(10**-9)*sig_phi/(4*np.pi*refractive_index*gauge_length*pe_scaling_factor) obspy_st[i].data = sig_ep2 print(overscale) return obspy_st #========================================================================================= if len(sys.argv) < 3: print(" ") print("USAGE: python h5tom.py [index] [start time \"2022-10-25T18:41:06\"] [duration \"420\"]") print(" ") sys.exit() #========================================================================================= ind = sys.argv[1] start_time = sys.argv[2] duration = sys.argv[3] data_path = "/work/das01/BR/liwei/DATA/DAS_FULL_DATA/" start = pd.to_datetime(start_time) end_time = start + pd.DateOffset(seconds=int(duration)) record_path = ind+"_"+str(start_time)+"_"+str(duration)+"s" if not os.path.isdir(record_path): os.makedirs(record_path) file_list = [] for i in range(0,int(duration)+60,60): file_time = start + pd.DateOffset(seconds=int(i)) h5_name = file_name(data_path, file_time) file_list.append(h5_name) # ## Channel range start_channel = 0 channel_num = 10245 #start_channel = 8000 #channel_num = 10 end_channel = start_channel + channel_num print("# start_channel = ", start_channel) print("# end_channel = ", end_channel) chan_range = [start_channel, end_channel] print("# chan_range = ", chan_range) # ## Raeding h5 file # original code from Peter. added start_unix_time from the first h5 file. this is used to define the startime time of miniseed file. for h5name in sorted(file_list): print("# h5name = ", h5name) f = h5py.File(h5name, 'r') list(f.keys()) RawData = f['Acquisition']['Raw[0]']['RawData'] RawDataTime = f['Acquisition']['Raw[0]']['RawDataTime'] RawData[100][0:10] # Original code from Peter Hubbard row1 = 0 count = 0 for h5name in tqdm(sorted(file_list)): print("# h5name = ", h5name) count += 1 f = h5py.File(h5name, 'r') RawData = f['Acquisition']['Raw[0]']['RawData'] RawDataTime = f['Acquisition']['Raw[0]']['RawDataTime'] if row1 == 0: all_RawData = np.transpose(RawData[:,chan_range[0]:chan_range[1]]) all_RawDataTime = np.array(RawDataTime) row1 = all_RawData.shape[0] start_unix_time = RawDataTime[0] else: all_RawData = np.append(all_RawData,np.transpose(RawData[:,chan_range[0]:chan_range[1]]),axis=1) all_RawDataTime = np.append(all_RawDataTime,RawDataTime,axis=0) row1 = all_RawData.shape[0] f.close() # This depends on the shape of raw data array, which may differ from different set-up. #all_RawData = np.transpose(all_RawData) # if the data is in HDF5 PRODML format, # to obtain phase shift values in “rad” divide each sample value by 10430.378850470453 #all_RawData = all_RawData*np.pi/2**15 all_RawData = all_RawData / 10430.378850470453 channels = all_RawData.shape[0] total_time = (all_RawDataTime[-1] - all_RawDataTime[0])/10**6 print(str(count)+' files imported. There are '+ str(total_time) + ' seconds recorded on '+ str(channels) + ' channels...') #print ("# start_unix_time = ", start_unix_time) time_utc = datetime.datetime.fromtimestamp(start_unix_time/10**6, pytz.timezone('UTC')) fs = 1/(all_RawDataTime[1]-all_RawDataTime[0])*10**6 print("# h5 file raw data sampling rate fs =", fs, "Hz") t = (all_RawDataTime - all_RawDataTime[0])/10**6 # ## Reading data # 1. read each trace, 2. decimate into 100sps data (or 1sps) assuming that the original is 1000sps data. station code=channel, location code = BK, channel/component = EPX, location code = 00 # # for i in trange(all_RawData.shape[0]): data = all_RawData[i,:] tr = Trace(data=data) tr.stats.sampling_rate = fs tr.stats.starttime = time_utc tr.stats.network = "MARS" tr.stats.station = str(start_channel+i).zfill(5) tr.stats.channel = "EPX" tr.stats.location = "00" # remove linear trend #tr.detrend(type='linear') # remove mean #tr.detrend(type='demean') # taper #tr.taper(0.05) # assuming original 1000 Hz data -> decimated into 100Hz #tr_deci = tr_fil.copy() #.decimate(factor=10,,no_filter=True) # 1sps case #tr_deci = tr_fil.copy().decimate(factor=5).decimate(factor=5).decimate(factor=2).decimate(factor=5).decimate(factor=2).decimate(factor=5) #print(tr_deci.stats) if i == 0: st = Stream(traces=tr) else: st += Stream(traces=tr) st = optic_phase_to_starin(st,gauge_length=20) #st.decimate(factor=2,no_filter=True).decimate(factor=5,no_filter=True).decimate(factor=5,no_filter=True) # ## output path #msOUT = record_path+'/'+'strain_'+str(st[0].stats.starttime)+'_'+str(st[0].stats.endtime)+'.mseed' #print("# msOUT = ", msOUT) #st.write(msOUT, format="MSEED") st_filtered = st.copy() st_filtered.detrend(type='linear') st_filtered.detrend(type='demean') st_filtered.taper(0.01) #st_filtered.decimate(factor=2,no_filter=True).decimate(factor=5,no_filter=True) st_filtered.filter('bandpass',freqmin=0.1,freqmax=1,corners=6, zerophase=True) #msOUT = record_path+'/'+'strain_'+str(st[0].stats.starttime)+'_'+str(st[0].stats.endtime)+'.mseed' #print("# msOUT = ", msOUT) #st.write(msOUT, format="MSEED") def plotSpaceTime(someData,minSec,maxSec,minCh,maxCh,title,sampleRate,plot_scale): # Basic error checking if (minSec >= maxSec) or (minSec < 0) or (maxSec*sampleRate > someData.shape[1]): print("ERROR in plotSpaceTime inputs minSec: "+str(minSec)+" or maxSec: "+str(maxSec)) return if (minCh >= maxCh) or (minCh < 0) or (maxCh > someData.shape[0]): print("Error in plotSpaceTime inputs minCh: "+str(minCh)+" or maxCh: "+str(maxCh)+" referring to array with "+str(someData.shape[0])+" channels.") return #norm = abs(someData[:-2000]).mean() #v_min= norm*-2.5 #v_max= norm*2.5 v_min= plot_scale*-0.000000001 v_max= plot_scale*0.000000001 norm = mcolors.TwoSlopeNorm(vmin=v_min, vmax=v_max, vcenter=0) # turn time range (in seconds) to indices minSecID = int(minSec*sampleRate) maxSecID = int(maxSec*sampleRate) # make the plot plt.figure(figsize= (20, 12)) plt.imshow(someData[minCh:maxCh,minSecID:maxSecID],aspect='auto',interpolation='none',cmap='seismic',extent=(minSec,maxSec,maxCh,minCh),norm=norm) plt.xlabel('time (s)') plt.ylabel('channel') plt.title(title) plt.colorbar() plt.savefig(record_path+'/Plot_'+str(ind)+'_'+str(title)+'_sc'+str(plot_scale)+'.png'); someData = np.array(st_filtered) max_sec = int(duration) + 30 plotSpaceTime(someData,minSec=0,maxSec=max_sec,minCh=0,maxCh=10244,title=str(start_time)+'_2-99Hz',sampleRate=200, plot_scale=1.5) plotSpaceTime(someData,minSec=0,maxSec=max_sec,minCh=0,maxCh=10244,title=str(start_time)+'_2-99Hz',sampleRate=200, plot_scale=15)