mseed读取[2]
折腾了好久,总算可以读出来了。
今天解决的
(1)批量转换mseed格式为sac格式
(2)找了比较适合的读取sac格式的文件rdsac.m
(3)matlab下批量读取sac文件并绘图
file_read=dir('E:\***\');
dirs=dir('E:\***\SACDATA\20171109\lhz'); % 用你需要的目录以及文件扩展名替换。读取某个目录的指定类型文件列表,返回结构数组。
dircell=struct2cell(dirs)' ; % 结构体(struct)转换成元胞类型(cell),转置一下是让文件名按列排列。
filenames=dircell(:,1) % 第一列是文件名
wjl=length(filenames);
qq=0;
d=[];t=[];
% for wj=3:wjl
for wj=3:wjl
yyy=filenames(wj);
zz=char(yyy);
f = fullfile('E:\***\SACDATA\20171109\lhz',zz);
fdata=rdsac(f);
data=fdata.d;tm=fdata.t;
d=[d;data];t=[t;tm];
% fid = fopen(f,'rb','ieee-be');
end
plot(t,d)
(4)老刘说重力与强震仪的波形应该是反项的,不懂原因。下周做看看。
============
https://cn.mathworks.com/matlabcentral/fileexchange/46356-rdsac-and-mksac--read-and-write-sac-seismic-data-file
function varargout=rdsac(varargin) %RDSAC Read SAC data file. % X=RDSAC(FILE) reads the Seismic Analysis Code (SAC) FILE and returns a % structure X containing the following fields: % t: time vector (DATENUM format) % d: data vector (double) % HEADER: header sub-structure (as defined in the IRIS/SAC format). % % [D,T0,H]=RDSAC(FILE) returns data vector D (single), origin time T0 as % a scalar (DATENUM format) and optional header as structure H. % % RDSAC without input argument will open a file browser window. % % [...]=RDSAC(...,'plot') or RDSAC(...) without output argument will plot % the data in a new figure. % % RDSAC(...,'enumerated') returns original integer values for enumerated % header fields (name start with an I), instead of descriptive string. % % Notes: % - RDSAC tries to detect automatically byte ordering of the file; % - time is corrected from B value; % % Acknowledgments: Arnesha Threatt, Rall Walsh % Reference: http://www.iris.edu/files/sac-manual/ % % Author: F. Beauducel <beauducel@ipgp.fr> % Created: 2014-04-01 % Updated: 2016-03-05 % Release history: % [2016-03-05] v1.2 % - adds NZDTTM header vector. % - enumerated header fields are now replaced by descriptive strings. % Use the new 'enumerated' option for backward compatibility. % % [2015-11-11] v1.1 % - output alternative returning data, origin time and header. % - undefined header fields are not returned. % % Copyright (c) 2016, Fran鏾is Beauducel, covered by BSD License. % All rights reserved. % % Redistribution and use in source and binary forms, with or without % modification, are permitted provided that the following conditions are % met: % % * Redistributions of source code must retain the above copyright % notice, this list of conditions and the following disclaimer. % * Redistributions in binary form must reproduce the above copyright % notice, this list of conditions and the following disclaimer in % the documentation and/or other materials provided with the % distribution % % THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS % IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED % TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A % PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT % OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, % SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT % LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, % DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY % THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT % (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE % OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. if nargin == 0 [filename,pathname] = uigetfile('*.SAC;*.sac','Select a SAC file'); f = [pathname,filename]; if filename == 0 error('Please select a SAC file or use function arguments.'); end end if nargin > 0 f = varargin{1}; if ~ischar(f) || ~exist(f,'file') error('FILENAME must be a valid file name.') end end fid = fopen(f, 'rb', 'ieee-le'); if fid == -1 error('Cannot open input data file %s',f); end [H,t0] = readheader(fid,varargin); % inconsistent header content might be due to big-endian byte ordering... if isnan(t0) fclose(fid); fid = fopen(f, 'rb', 'ieee-be'); % closes and re-open [H,t0] = readheader(fid); end d = fread(fid,H.NPTS,'*float32'); % imports data as single class fclose(fid); if length(d) ~= H.NPTS || isnan(t0) warning('Inconsistent data header: may be not a SAC file.'); end % makes time vector (using sampling interval DELTA and time correction B) t = t0 + (H.B + (0:H.DELTA:(H.NPTS - 1)*H.DELTA)')/86400; if nargout == 1 varargout{1} = struct('t',t,'d',double(d),'HEADER',H); elseif nargout > 1 varargout{1} = d; varargout{2} = t(1); varargout{3} = H; end % plots the data if (nargout == 0 || any(strcmpi(varargin,'plot'))) && length(d) == length(t) figure plot(t,d) xlim = [min(t),max(t)]; set(gca,'XLim',xlim) datetick('x','keeplimits'); xlabel(sprintf('%s to %s',datestr(xlim(1)),datestr(xlim(2)))) ylabel('Count') if isfield(H,'KNETWK') net = H.KNETWK; else net = '?'; end if isfield(H,'KSTNM') stn = H.KSTNM; else stn = '?'; end if isfield(H,'KCMPNM') cmp = H.KCMPNM; else cmp = '?'; end title(sprintf('[%s:%s:%s] %s',net,stn,cmp,f)) end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function [H,t0] = readheader(fid,vararg) novalue = -12345; hn = [fread(fid,[5,14],'float32'),fread(fid,[5,8],'int32')]; hs = fread(fid,[8,24],'*char')'; % --- classifies header fields % numerical variables v = { ... 'DELTA', 'DEPMIN', 'DEPMAX', 'SCALE', 'ODELTA'; 'B', 'E', 'O', 'A', 'INTERNAL'; 'T0', 'T1', 'T2', 'T3', 'T4'; 'T5', 'T6', 'T7', 'T8', 'T9'; 'F', 'RESP0', 'RESP1', 'RESP2', 'RESP3'; 'RESP4', 'RESP5', 'RESP6', 'RESP7', 'RESP8'; 'RESP9', 'STLA', 'STLO', 'STEL', 'STDP'; 'EVLA', 'EVLO', 'EVEL', 'EVDP', 'MAG'; 'USER0', 'USER1', 'USER2', 'USER3', 'USER4'; 'USER5', 'USER6', 'USER7', 'USER8', 'USER9'; 'DIST', 'AZ', 'BAZ', 'GCARC', 'INTERNAL'; 'INTERNAL','DEPMEN', 'CMPAZ', 'CMPINC', 'XMINIMUM'; 'XMAXIMUM','YMINIMUM','YMAXIMUM','UNUSED', 'UNUSED'; 'UNUSED', 'UNUSED', 'UNUSED', 'UNUSED', 'UNUSED'; 'NZYEAR', 'NZJDAY', 'NZHOUR', 'NZMIN', 'NZSEC'; 'NZMSEC', 'NVHDR', 'NORID', 'NEVID', 'NPTS'; 'INTERNAL','NWFID', 'NXSIZE', 'NYSIZE', 'UNUSED'; 'IFTYPE', 'IDEP', 'IZTYPE', 'UNUSED', 'IINST'; 'ISTREG', 'IEVREG', 'IEVTYP', 'IQUAL', 'ISYNTH'; 'IMAGTYP', 'IMAGSRC', 'UNUSED', 'UNUSED', 'UNUSED'; 'UNUSED', 'UNUSED', 'UNUSED', 'UNUSED', 'UNUSED'; 'LEVEN', 'LPSPOL', 'LOVROK', 'LCALDA', 'UNUSED'; }'; for n = 1:numel(v) if ~strcmp(v(n),'UNUSED') && hn(n) ~= novalue H.(v{n}) = hn(n); end end % string variables v = { ... 'KSTNM', 'KEVNM0', 'KEVNM1'; 'KHOLE', 'KO', 'KA'; 'KT0', 'KT1', 'KT2'; 'KT3', 'KT4', 'KT5'; 'KT6', 'KT7', 'KT8'; 'KT9', 'KF', 'KUSER0'; 'KUSER1', 'KUSER2', 'KCMPNM'; 'KNETWK', 'KDATRD', 'KINST'; }'; for n = 1:numel(v) s = deblank(hs(n,:)); if ~strcmp(s,num2str(novalue)) && ~isempty(s) H.(v{n}) = s; end end % concatenates KEVNM (to respect exactly the IRIS format) if isfield(H,'KEVNM0') && isfield(H,'KEVNM1') H.KEVNM = [H.KEVNM0,H.KEVNM1]; H = rmfield(H,v(2:3)); end % checks the origin time validity t0 = NaN; if H.NZYEAR >= novalue ... && (H.NZJDAY >= 1 && H.NZJDAY <= 366 || H.NZJDAY == novalue) ... && (H.NZHOUR >= 0 && H.NZHOUR < 24 || H.NZHOUR == novalue) ... && (H.NZMIN >= 0 && H.NZMIN < 60 || H.NZMIN == novalue) ... && (H.NZSEC >= 0 && H.NZSEC < 60 || H.NZSEC == novalue) t0 = datenum(H.NZYEAR,1,H.NZJDAY,H.NZHOUR,H.NZMIN,H.NZSEC + H.NZMSEC/1e3); % readable origin time H.NZDTTM = [H.NZYEAR,H.NZJDAY,H.NZHOUR,H.NZMIN,H.NZSEC,H.NZMSEC]; H.KZDATE = datestr(t0,sprintf('mmm dd (%03d) yyyy',H.NZJDAY)); H.KZTIME = datestr(t0,'HH:MM:SS.FFF'); end % replaces enumerated values by their explicite description (string) if ~any(strcmpi(vararg,'enumerated')) fields = fieldnames(H); enum = fields(strncmpi(fields,'I',1)); for n = 1:length(enum) E = enumheader(H.(enum{n})); if ~isempty(E.code) H.(enum{n}) = sprintf('%s {%g}',upper(E.description),H.(enum{n})); end end end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function E = enumheader(n) v = { ... 01 , 'ITIME' , 'Time series file'; 02 , 'IRLIM' , 'Spectral file - real and imaginary'; 03 , 'IAMPH' , 'Spectral file - amplitude and phase'; 04 , 'IXY' , 'General x versus y data'; 05 , 'IUNKN' , 'Unknown'; 06 , 'IDISP' , 'Displacement in nm'; 07 , 'IVEL' , 'Velocity in nm/s'; 08 , 'IACC' , 'Acceleration in nm/s/s'; 09 , 'IB' , 'Begin time'; 10 , 'IDAY' , 'GMT day'; 11 , 'IO' , 'Event origin time'; 12 , 'IA' , 'First arrival time'; 13 , 'IT0' , 'User defined time pick 0'; 14 , 'IT1' , 'User defined time pick 1'; 15 , 'IT2' , 'User defined time pick 2'; 16 , 'IT3' , 'User defined time pick 3'; 17 , 'IT4' , 'User defined time pick 4'; 18 , 'IT5' , 'User defined time pick 5'; 19 , 'IT6' , 'User defined time pick 6'; 20 , 'IT7' , 'User defined time pick 7'; 21 , 'IT8' , 'User defined time pick 8'; 22 , 'IT9' , 'User defined time pick 9'; 23 , 'IRADNV' , ''; 24 , 'ITANNV' , ''; 25 , 'IRADEV' , ''; 26 , 'ITANEV' , ''; 27 , 'INORTH' , ''; 28 , 'IEAST' , ''; 29 , 'IHORZA' , ''; 30 , 'IDOWN' , ''; 31 , 'IUP' , ''; 32 , 'ILLLBB' , ''; 33 , 'IWWSN1' , ''; 34 , 'IWWSN2' , ''; 35 , 'IHGLP' , ''; 36 , 'ISRO' , ''; 37 , 'INUCL' , 'Nuclear event'; 38 , 'IPREN' , 'Nuclear pre-shot event'; 39 , 'IPOSTN' , 'Nuclear post-shot event'; 40 , 'IQUAKE' , 'Earthquake'; 41 , 'IPREQ' , 'Foreshock'; 42 , 'IPOSTQ' , 'Aftershock'; 43 , 'ICHEM' , 'Chemical explosion'; 44 , 'IOTHER' , 'Other'; 45 , 'IGOOD' , 'Good data'; 46 , 'IGLCH' , 'Glitches'; 47 , 'IDROP' , 'Dropouts'; 48 , 'ILOWSN' , 'Low signal to noise ratio'; 49 , 'IRLDTA' , 'Real data'; 50 , 'IVOLTS' , 'Velocity in V'; 52 , 'IMB' , 'Bodywave Magnitude'; 53 , 'IMS' , 'Surfacewave Magnitude'; 54 , 'IML' , 'Local Magnitude'; 55 , 'IMW' , 'Moment Magnitude'; 56 , 'IMD' , 'Duration Magnitude'; 57 , 'IMX' , 'User Defined Magnitude'; 58 , 'INEIC' , 'National Earthquake Information Center'; 59 , 'IPDEQ' , ''; 60 , 'IPDEW' , ''; 61 , 'IPDE' , 'Preliminary Determination of Epicenter'; 62 , 'IISC' , 'Internation Seismological Centre'; 63 , 'IREB' , 'Reviewed Event Bulletin'; 64 , 'IUSGS' , 'US Geological Survey'; 65 , 'IBRK' , 'UC Berkeley'; 66 , 'ICALTECH' , 'California Institute of Technology'; 67 , 'ILLNL' , 'Lawrence Livermore National Laboratory'; 68 , 'IEVLOC' , 'Event Location (computer program)'; 69 , 'IJSOP' , 'Joint Seismic Observation Program'; 70 , 'IUSER' , 'The individual using SAC2000'; 71 , 'IUNKNOWN' , 'Unknown'; 72 , 'IQB' , 'Quarry or mine blast confirmed by quarry'; 73 , 'IQB1' , 'Quarry/mine blast with designed shot info-ripple fired'; 74 , 'IQB2' , 'Quarry/mine blast with observed shot info-ripple fired'; 75 , 'IQBX' , 'Quarry or mine blast - single shot'; 76 , 'IQMT' , 'Quarry/mining-induced events: tremors and rockbursts'; 77 , 'IEQ' , 'Earthquake'; 78 , 'IEQ1' , 'Earthquakes in a swarm or aftershock sequence'; 79 , 'IEQ2' , 'Felt earthquake'; 80 , 'IME' , 'Marine explosion'; 81 , 'IEX' , 'Other explosion'; 82 , 'INU' , 'Nuclear explosion'; 83 , 'INC' , 'Nuclear cavity collapse'; 84 , 'IO_' , 'Other source of known origin'; 85 , 'IL' , 'Local event of unknown origin'; 86 , 'IR' , 'Regional event of unknown origin'; 87 , 'IT' , 'Teleseismic event of unknown origin'; 88 , 'IU' , 'Undetermined or conflicting information'; 89 , 'IEQ3' , ''; 90 , 'IEQ0' , ''; 91 , 'IEX0' , ''; 92 , 'IQC' , ''; 93 , 'IQB0' , ''; 94 , 'IGEY' , ''; 95 , 'ILIT' , ''; 96 , 'IMET' , ''; 97 , 'IODOR' , ''; 103 , 'IOS' , ''; }; v(strcmp(v(:,3),''),3) = v(strcmp(v(:,3),''),2); k = find(cat(1,v{:,1})==n); if isempty(k) E = struct('code',[],'description',[]); else E.code = v{k,2}; E.description = v{k,3}; end