%% MATLAB example code
% To accompany the November 8, 2012 ECI Webinar
% updated June 21, 2013 for irisFetch v2.0.1
%
% http://www.iris.edu/hq/webinar/
% Celso Reyes
% IRIS Data Management Center
%% Setup Notes for the first time user
% First, grab irisFetch.m from http://www.iris.edu/manuals/javawslibrary/matlab/
% or, by searching for matlab in http://seiscode.iris.washington.edu
% - [ open browser ]
% - [ click on the IRIS SeisCode link from the irisFetch bookmark ]
% - [ type " irisfetch " into the search field ]
% place it somewhere matlab can find it.
%% Retrieve a seismic trace - and getting started notes.
tr = irisFetch.Traces('IU','ANMO','00','BHZ','2010-02-27 06:30:00','2010-02-27 10:30:00')
% Calling format: NETWORK, STATION, LOCATION, CHANNEL, START_TIME, END_TIME)
% This calls the TRACES routine contained within the irisFetch.m file.
% If A warning message showed up, it was likely because MATLAB didn't yet know
% where the java library (IRIS-WS.jar) lives. However, this routine was smart
% enough to connect to an online version.
% Even so, you should download the library.
% ... Conveniently, the error message tells you where you can find it.
% and then let matlab know where it lives by using the javaddpath routine.
% replace 'PATH_TO_THE_JAR' with the actual path. Also, note the name might be different
% if you are working with a later version. (This example points to version 1.5.2)
javaaddpath('PATH_TO_THE_JAR/IRIS-WS-2.0.2.jar'); %change this to a real path...
% Note: This need only be done when you first start matlab.
%% Retrieve a seismic trace - Describing the returned structure
tr = irisFetch.Traces('IU','ANMO','00','BHZ','2010-02-27 06:30:00','2010-02-27 10:30:00')
% Notice: The retrieved information is put into a matlab struct.
% A structs fields are accessible using the convention VARIABLE.FIELDNAME
% for example:
tr.station % get the station name
tr.startTime % get the start time for this trace
% The above just gave me a number, which is how MATLAB stores dates internally.
% To convert this into a readable text format, use MATLAB's DATESTR function
datestr(tr.startTime)
%% Wildcards work, too
% Now, let's modify the original query slightly. Instead of looking for
% the vertical channel, let's request all BH channels.
% To do this, I use a question mark. (the irisFetch.Traces routine
% understands the asterix as well.
tr = irisFetch.Traces('IU','ANMO','00','BH?','2010-02-27 06:30:00','2010-02-27 10:30:00')
% This retrieves three traces, containing the three BH channels.
% To access each individual trace you specify it by index.
tr(1)
tr(2)
% etc...
% One useful feature of storing data as an array of structs, I can access
% ALL the information of a particular type at once.
% get a list of all the sample rates
allSamples = [tr.sampleRate]
% or azimuths
[tr.azimuth]
% The brackets tell MATLAB to concatenate the asked-for information into an
% array. The brackets are sensitive to the size of each item, and generally
% only work when each item has the same number of elements.
allData = [tr.data];
% we just made a single array from all three components
whos allData
% For items of differing lengths, or for text, (such as station names, comments, or instrument names)
% I use curly brackets. This tells MATLAB to store these in a special type
% of array called a "cell".
% For example, to get a list of all channels.
allChannels = {tr.channel}
%% create a graph with a legend of all 3 channels
% first, plot the data
plot(allData)
% now, add a legend using matlab's "LEGEND" function.
legend(allChannels)
% So, with just a couple lines of MATLAB code, we were able to access data from across the web and plot it.
%% a note about using HELP
% Additional features are available, to learn more about them, consult the help.
%you can get general help for irisFetch by merely typing "help irisFetch"
help irisFetch
%you can also get detailed help for any of the other routines by specifying them. For example:
help irisFetch.Traces
%% EVENT examples
% IrisFetch can also retrieve data from the event webservice.
% The usage is irisFetch.Events ( parameterName1, value1 [, parameterName2, value2[, ... etc])
% For this example, retrieve all events with M >= 6.0 since the beginning of 2010.
ev = irisFetch.Events('StartTime', '1/1/2010', 'MinimumMagnitude', 6.0)
%Again, ev is an array of structs.
% Here is the first one:
ev(1)
%The Second one:
ev(2)
% and so on...
% These have fields that are themselves structs. Never fear, they can be accessed by
% adding another dot and a field.
% So, using the preferredOrigins as an example:
ev(1).PreferredOrigin
ev(1).PreferredOrigin.Latitude
ev(1).PreferredOrigin.Longitude
% In order to use them in bulk, you can create an origin list just like we did for the trace's data
origs = [ev.PreferredOrigin]
% starting with irisFetch version 1.4.0, some of the more accessed "Preferred" fields have been
% copied into the main structure, making them easier to work with.
% Now, we can combine some of the interesting and plottable data...
lats = [ev.PreferredLatitude];
lons = [ev.PreferredLongitude];
% I'm going to see how far away each event is from station ANMO. I only need the BHZ channel so...
% keep only the BHZ trace
tr = tr( ismember({tr.channel},'BHZ') );
% Now using the "distance" function from MATLAB's mapping toolbox, I'm going to get the
% arc length distances. along with a reasonable subsample of these events.
[arclen, az] = distance([tr.latitude],[tr.longitude],[lats],[lons]);
% There are lots of events, but I only want one event for every 2 degrees of distance so:
binned_at_2deg = fix(arclen / 2); %put into 2degree bins
[b,i,j] = unique(binned_at_2deg,'first'); %keep first event in each bin
%% OK, now we're ready for a more complex example...
% some of the plotting examples require MATLAB's mapping toolbox.
%% prepare the map for plotting
%
% Generate some sample plots:
% MAP: show the location of the station [Black square], as well as the
% locations of all earthquakes [gray dots], and the locations of earthquakes used in
% the "traveltime" chart. The most recently plotted event is in RED.
%
% TTCHART: Plot the trace for each event, in order of distance from the station. Each
% trace is scaled and potted in its position at X degrees away.
%
% SPECTROGRAM: The spectral plot for each event.
figure(1); clf;
set(1,'Name','Map and TT plot');
% prepare a world map [requires the mapping toolbox]
subplot(2,1,1);
myMap = worldmap('World');
coast=load('coast');
plotm(coast.lat,coast.long);
% plot the location of the station
plotm(tr.latitude,tr.longitude,'ks','linewidth',2);
% plot all of the earthquake epicenters, as gray dots
alleq = plotm( lats , lons ,'.','markersize',9,'color',[.3 .3 .3]);
% subset of earthquakes that will actually be plotted, as green dots
usedeq = plotm( lats(i) , lons(i) ,'g.','markersize',9)
%currently examined epicenter, in black
h = plotm(tr.latitude,tr.longitude,'k*','linewidth',2);
wiggles = subplot(2,1,2);
% set up spectrogram plot
figure(3);
clf;
set(3,'Name','Spectra');
spectraH = gca;
% Establish some constants for nice readability
HALF_HR = 0.5/24;
WINDOW = 1024;
NOVERLAP = round(WINDOW * 0.8);
NFFT = 1024;
%% plot wiggles at a variety of distances from the chosen station
for idx=i
startNum = datenum(ev(idx).PreferredTime);
endNum = startNum + HALF_HR;
tr = irisFetch.Traces('IU','ANMO','00','BHZ', startNum, endNum);
data = tr.data - mean(tr.data); % demean
scaledData = data / max( data / 2);
timeSteps = [0:(tr.sampleCount-1)] / tr.sampleRate;
plot( wiggles , timeSteps / 60, scaledData + arclen(idx));
xlabel(wiggles , 'Minutes' );
ylabel(wiggles , 'Degrees away' );
hold( wiggles , 'on' );
axes(spectraH);
spectrogram(data,WINDOW,NOVERLAP,NFFT,tr.sampleRate,'yaxis');
ylim([0,6]);
axes(myMap);
delete(h);
h = plotm(lats(idx),lons(idx),'r*','markersize',11,'linewidth',1.5);
end
hold(wiggles,'off');
%% plot earthquakes on map example.
% This example requires the mapping toolbox, and was not demonstrated during the webinar
symbolSize = @(x) ( [x.PreferredMagnitudeValue] - 3) .^ 3 + 5;
colorIndex = @(x) sqrt( [ev.PreferredDepth] );
ev = irisFetch.Events( 'StartTime' , '10/1/2010' , 'MinimumMagnitude' , 5 ,'MaximumMagnitude' , 8 );
[~, I] = sort( [ev.PreferredDepth] , 'descend' );
ev = ev(I);
myMap = worldmap('World');
coast = load('coast');
colormap( (jet + jet + copper) / 3.2 ); % define colors used for events
scatterm([ev.PreferredLatitude] , [ev.PreferredLongitude] , symbolSize(ev) , colorIndex(ev));
plotm( coast.lat , coast.long , 'color' , [0 0.2 0] );