{"id":3035,"date":"2025-01-28T10:02:57","date_gmt":"2025-01-28T15:02:57","guid":{"rendered":"https:\/\/blogs.mathworks.com\/matlab\/?p=3035"},"modified":"2025-01-28T10:23:08","modified_gmt":"2025-01-28T15:23:08","slug":"visualizing-yesterdays-earthquake-in-boston-using-matlab-and-a-raspberry-shake","status":"publish","type":"post","link":"https:\/\/blogs.mathworks.com\/matlab\/2025\/01\/28\/visualizing-yesterdays-earthquake-in-boston-using-matlab-and-a-raspberry-shake\/","title":{"rendered":"Visualizing yesterday’s earthquake in Boston using MATLAB and a Raspberry Shake"},"content":{"rendered":"
During an online meeting with my colleagues in Boston yesterday, some of them were rather startled when their offices started violently shaking for a brief time! We quickly hopped onto BlueSky (Feel free to <\/span>follow me there<\/span><\/a>) where we saw thousands of messages referring to what everyone suspected was an earthquake. We now know that it was indeed an earthquake that had a magnitude of 3.8 on the Richter scale. <\/span><\/div>

Learning that MathWorks has a seismograph<\/span><\/h2>
\"You know we have a seismograph at MathWorks right?\" one of my colleagues asked me. Another colleague, Christopher Stapels, <\/span>wrote about it in 2022<\/span><\/a> when he used it to study traffic passing by MathWorks' headquarters. The seismograph is a <\/span>Raspberry Shake<\/span><\/a> and is constantly sending data to the cloud. <\/span><\/div>
Our shake has an ID of RF23B and you can visualize its location, along with all of the other shakes online in the world, at <\/span>stationview.raspberryshake.org<\/span><\/a>. <\/span><\/div>
Now you know the ID of our device, you can download data from it as a <\/span>miniSeed<\/span><\/a> file. The API call below requests the data between 3pm and 4pm GMT on 27th January for the device RF23B and saves it as the file <\/span>waveform.mseed<\/span>.<\/span><\/div>
urlQuery =<\/span>...<\/span><\/span><\/div><\/div>
\"https:\/\/data.raspberryshake.org\/fdsnws\/dataselect\/1\/query?starttime=2025-01-27T15:00:00&endtime=2025-01-27T16:00:00&network=AM&station=RF23B\"<\/span>;<\/span><\/span><\/div><\/div>
websave(<\/span>\"waveform.mseed\"<\/span>,urlQuery);<\/span><\/span><\/div><\/div><\/div>

Loading the seismograph data into MATLAB<\/span><\/h2>
In order to load the miniSeed file I made use of the <\/span>RDMSEED and MKMSEED<\/span><\/a> package, written by <\/span>Fran\u00e7ois Beauducel<\/span><\/a>, on File Exchange. You'll need to install this before proceeding. <\/span><\/div>
Once this is installed we can load and visualize the <\/span>waveform.seed<\/span> file we downloaded earlier<\/span><\/div>
waveform = rdmseed(<\/span>'waveform.mseed'<\/span>,<\/span>'simple'<\/span>)<\/span><\/span><\/div>
waveform = struct with fields:<\/span><\/span><\/div>
t: [360454\u00d71 double]\n d: [360454\u00d71 double]\n name: 'AM:RF23B:00:EHZ'\n<\/div><\/div><\/div><\/div><\/div><\/div>
The <\/span>t<\/span> field of this structure represents time and is an array of numbers that records time relative to some point in the past. Here is the first one<\/span><\/div>
waveform.t(1)<\/span><\/span><\/div>
ans = <\/span>7.3964e+05<\/div><\/div><\/div><\/div><\/div><\/div>
I'd prefer to have the times represented as <\/span>datetime<\/span><\/a> values so need to figure out how to convert them. There are a <\/span>range of conversion possibilities<\/span><\/a> available via the <\/span>'ConvertFrom'<\/span> option of <\/span>datetime<\/span>. For example '<\/span>datenum'<\/span> assumes everything is relative to 0-Jan-0000 while <\/span>'excel'<\/span> assumes everything is relative to 0-Jan-1900. The correct offset is probably discussed in the documentation somewhere but I'm lazy and just guessed.<\/span><\/div>
datetime(waveform.t(1),<\/span>\"ConvertFrom\"<\/span>,<\/span>\"datenum\"<\/span>)<\/span><\/span><\/div>
ans = datetime<\/span><\/span><\/div>
27-Jan-2025 14:59:58\n<\/div><\/div><\/div><\/div><\/div><\/div>
Looks right! I Let's plot<\/span><\/div>
waveform.u = datetime(waveform.t,<\/span>'ConvertFrom'<\/span>,<\/span>'datenum'<\/span>);<\/span><\/span><\/div><\/div>
plot(waveform.u,waveform.d)<\/span><\/span><\/div><\/div>
axis <\/span>tight<\/span><\/span><\/div>
<\/div>
<\/div>
<\/div>
<\/div>
<\/div><\/div><\/div><\/div><\/div>
I'm not sure what the y-axis represents yet but I imagine this to be like a pen on an old-school seismograph until I learn better. Something big happened between 15:15 and 15:30, let's see when this was exactly by finding at what time the maximum amplitude occurred. <\/span><\/div>
[maxD,Idx] = max(abs(waveform.d));<\/span><\/span><\/div><\/div>
waveform.u(Idx)<\/span><\/span><\/div>
ans = datetime<\/span><\/span><\/div>
27-Jan-2025 15:23:29\n<\/div><\/div><\/div><\/div><\/div><\/div>
Looking online, the earthquake is reported to have occurred eight miles underground near York Harbor at 15:22 GMT - a minute and a half from where I found the maximum amplitude. I know next to nothing about seismograph analysis and so even this discrepancy is confusing to me. Several things occur to me that might explain it but I'm not sure what the answer is. <\/span><\/div>
I'm also wondering what else we could learn about the earthquake from this dataset. Do any of you have any insight here? Let me know in the comments section, or on via any of my social accounts below<\/span><\/div>