Dear Tech team,
How do I convert from counts in acceleration in a Raspberry shake 4D?
Thank you in advance,
edited 2020-07-08: refined code and updated quake
Hi Timo,
Welcome to the community forum! Using Python you can easily convert from counts to acceleration (this is also called “removing the instrument response”) by installing obspy and running the following two pieces of code:
from obspy.clients.fdsn import Client
from obspy.core import UTCDateTime, Stream
rs = Client('RASPISHAKE')
# set the station name and download the response information
stn = 'R24FA' ### <<< replace this with your station name ###
inv = rs.get_stations(network='AM', station=stn, level='RESP')
# set data start/end times
start = UTCDateTime(2020, 3, 7, 2, 40, 45) # (YYYY, m, d, H, M, S)
end = UTCDateTime(2020, 3, 7, 2, 42, 0) # (YYYY, m, d, H, M, S)
# set the FDSN server location and channel names
channels = ['EHZ', 'ENE', 'ENN', 'ENZ'] # ENx = accelerometer channels; EHx or SHZ = geophone channels
# get waveforms and put them all into one Stream
stream = Stream()
for ch in channels:
trace = rs.get_waveforms('AM', stn, '00', ch, start, end)
stream += trace
# here you will plot raw counts for all four channels
stream.plot()
Plot of raw data in counts:
# attach the response info
stream.attach_response(inv)
# remove the response and plot
resp_removed = stream.remove_response(output='ACC') # convert to acceleration in M/S
# now plot all channels in acceleration units (m/s/s)
resp_removed.plot()
Plot of acceleration in m/s2:
Note that the ground motion for this quake near my house in Panama is just barely enough to register on the accelerometers—since the accelerometer channels are not very sensitive to weak motion—so they are more noisy. Keep in mind that most of the time, when the ground is not moving enough for humans to feel, the accelerometer channels will be too noisy to get anything good out of.
There are other ways to do it, but as far as I am aware, Python/obspy this is the most reliable (free) way. Others may be able to comment on how to do it in MATLAB or R. I think it is also possible in SWARM, but only if the server provides metadata for the station (not all servers do).
Ian
Dear Ian,
Thanks for the code.
When I wrote my station’s name, I get the following error:
HTTPError: 204 HTTP Error: No Content for url: https://fdsnws.raspberryshakedata.com/fdsnws/station/1/query?network=AM&station=REE93&level=resp&format=xml
My station has no internet connection, so I downloaded the data in “mseed” files for each component.
Kind regards,
Timo
Timo:
The convenience of the web services is only available to Shakers who publicly share their data in real-time.
In your case, as this application is offline/ datalogger mode, you would have to generate the inventory-xml instrument response file manually using one of the templates available at: http://manual.raspberryshake.org/metadata.html#metadata. Then you should be good to go.
Yours, branden
Timo
With Obspy and R4D connected to the PC
from obspy.core import UTCDateTime
from obspy.clients.earthworm import Client
start=UTCDateTime("2019-01-20 1:35:00") # another way to place the date
end=UTCDateTime("2019-01-20 1:40:00") # another way to place the date
client=Client("IP of your R4D", 16032)
st=client.get_waveforms('AM',"Name of your station",'00','EN?',start,end)
st.plot()
st=st.detrend('constant')
r1=st[2].data # ENE
r2=st[1].data # ENN
r3=st[0].data # ENZ
for i in range(len(r1)):
r1[i]=r1[i]*(1/3.85e5)/9.81 # (from Tech. Spec R4D V5 1/3.85e5 count/m/seg2 )
r2[i]=r2[i]*(1/3.85e5)/9.81
r3[i]=r3[i]*(1/3.85e5)/9.81
s1='file output in ascii'
f=open(s1,'w')
for i in range(len(r1)):
s=('{0:20.17f} {1:20.17f} {2:20.17f} \n'.format(r1[i],r2[i],r3[i]))
f.write(s)
f.close()
Alejandro
1 Like
Thank you Branden and Alejandro. Finally, I could include the instrument response.
1 Like
Thanks Alejanro for your effective small code.
However how can I continously evaluate the PGA instead of accessing a stored waveform time period ?
Thanks
Alessandro
Alessandro
With the code you can not get in real time PGA with Obspy for the delay that they have spoken several times in this forum.
Using the UDP real-time packet output indicated in the manual, I have not been able to obtain the 4 separate channels yet. Requirement to obtain PGA in real time
Alejandro
I have created and attached a dummy RS4D XML file with no location. Right now the station code in this file is Timo’s REE93, but can be changed in a text editor. The XML station code tag to look for is <Station code="REE93" startDate="2019-03-10T19:43:32.072" restrictedStatus="open">
This should allow you to attach to a RS4D stream from a REE93 miniSEED file using the read_inventory() function in obspy. You can then calculate things like ground acceleration on all channels without having to contact the RShake servers.
Let me know if you have questions.
Ian
REE93-no_location-dummyfile.xml (18.3 KB)
dear Iannesbitt,
I don’t know if your answer is for me or for other issues in the same thread.
If yes, could you be more explicit on how I could translate this into code:
You can then calculate things like ground acceleration on all channels without having to contact the RShake servers.
Thanks
Alessandro
Sure. Right now there’s no way to easily do this in real time, although later this year we will release a python library that will be able to receive UDP data directly to an obspy stream. This example will calculate PGA for a MiniSEED from the same station as plotted above, SF4A5. Example files are attached below.
from obspy import read_inventory, read
import os
stream = read(os.path.expanduser('~/Downloads/SF4A5-20190528-124800.mseed'))
inv = read_inventory(os.path.expanduser('~/Downloads/SF4A5-no_location-dummyfile.xml'))
stream.attach_response(inv)
stream.remove_response(output='ACC') # convert to acceleration in M/S
for trace in stream:
print('Channel %s PGA: %.5f' % (trace.stats.channel, max(abs(trace.data))))
Here’s what the output of that print statement looks like:
Channel EHZ PGA: 0.01045
Channel ENE PGA: 0.01438
Channel ENN PGA: 0.02364
Channel ENZ PGA: 0.01375
SF4A5-20190528-124800.mseed (109 KB)
SF4A5-no_location-dummyfile.xml (18.3 KB)
You could also do this with miniSEED files on your Shake, but we don’t recommend installing obspy on the Shake because it uses significant system resources, so you’d have to periodically transfer from the Shake to the computer running obspy somehow.
Let me know if you have questions.
Ian
can we save it to a new miniseed file after removing the response or will there be a dataloss?
stream.write(os.path.expanduser('~/Downloads/SF4A5-respremoved.mseed'))
This line will save a response-removed miniseed file in Downloads. This file will have only ground acceleration, not counts, so in that sense there is data loss. It’s mathematically possible to convert back to raw data but I’m not sure if obspy has a function for it.
If you’re looking to preserve raw data and response removed data, you should write the stream both before and after you remove the response.
Ian
If I use the function below and then use the stream.write() function then the saved stream should not contain acceleration data.
stream.remove_response()
stream.write(os.path.expanduser('~/Downloads/SF4A5-respremoved.mseed'))
but will there be a data loss?
And do you recommend using pre_filter and water_level with the instrument response to remove the noise as explained in obspy because it totally relies on the instrument response.
pre_filt = [0.001, 0.005, 45, 50]
tr.remove_response(inventory=inv, pre_filt=pre_filt, output="VEL", water_level=60, plot=True)
Hi shakeel
I realize now that I wasn’t clear in my previous post. You are correct that this will not contain acceleration, it will contain velocity since obspy’s documentation shows that the default for the remove_response function is velocity.
So for example remove_response() is equivalent to remove_response(output='VEL').
As I say above, applying the inventory function alone with will not result in data loss, but where my answer was incomplete was that it did not mention that the default water_level will cause data loss. obspy currently has no way to reattach the response, partially because of this.
If you use water_level=60 (the default), then there will be unrecoverable data loss since you will be filtering out all response information below a certain power threshold in the frequency domain. A good visual example is shown here:
There’s nothing wrong with this per se, because it will filter out noise in less powerful spectral bands which is often desirable in this situation, but you will definitely lose data from the less powerful frequency bands, which may not be acceptable for your case.
Using pre_filt in the function call is up to you. I have not encountered a scenario where I need it but it seems like it could be useful for removing unwanted noise from certain spectral bands that you know to be noisy. This is where looking at the data in spectrogram or frequency domain plots is often very useful. The obspy documentation states the following about pre_filt:
Apply a bandpass filter in frequency domain to the data before deconvolution. The list or tuple defines the four corner frequencies (f1, f2, f3, f4) of a cosine taper which is one between f2 and f3 and tapers to zero for f1 < f < f2 and f3 < f < f4.
If you are looking to retain as much data as possible without data loss, then I would recommend saving the files as miniSEED without removing the response, then building in a response removal function to your software that gets applied at run time. You can also write your inventory information locally using inv.write(os.path.expanduser('~/Downloads/SF4A5-inventory.txt') so you don’t have to contact the RS servers each time your software runs.
I hope this answer was more clear.
Ian
Thank you for a detailed reply. I think we can avoid water_level by setting it’s value to water_level=None. In this way we can preserve the data. Secondly I agree that pre_filt function is only used for some specific applications.
1 Like
dear Ian,
I see from your code above that you convert also EHZ in PGA. I understood that EHZ should be the velocity component while ENE, ENN and ENZ the acceleration ones. So EHZ should be in m/s and not m/s2. How do I convert from counts to m/s ?
For the other you already indicated ENNcounts/385000/9.81 but what about EHZ ?
Thanks
Alessandro
We are in the end of a project so that this information is very important for us. Can you please indicate how to convert EHZ in cm/s or m/s ?
thanks
Alessandro
Hi Alessandro, see the obspy documentation about remove_response(): https://docs.obspy.org/packages/autogen/obspy.core.trace.Trace.remove_response.html
remove_response(output='VEL') will give you velocity.