Converting counts in acceleration

#1

Dear Tech team,

How do I convert from counts in acceleration in a Raspberry shake 4D?

Thank you in advance,

#2

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 import read_inventory, read
from obspy.clients.fdsn import Client
from obspy.core import UTCDateTime, Stream

# set the station name and download the response information
stn = 'R6A3B'                            # your station name
inv = read_inventory('https://fdsnws.raspberryshakedata.com/fdsnws/station/1/query?network=AM&station=%s&level=resp&format=xml' % (stn))

# set data start/end times
start = UTCDateTime(2019, 5, 28, 12, 48, 0) # (YYYY, m, d, H, M, S)
end = UTCDateTime(2019, 5, 28, 12, 50, 0) # (YYYY, m, d, H, M, S)

# set the FDSN server location and channel names
rs = Client(base_url='https://fdsnws.raspberryshakedata.com/')
channels = ['EHZ', 'ENE', 'ENN', 'ENZ'] # ENx = accelerometer channels; EHx or SHZ = geophone channels

# get waveforms, put them all into one Stream, and attach the response
stream = Stream()
for ch in channels:
    trace = rs.get_waveforms('AM', stn, '00', ch, start, end)
    stream += trace
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
resp_removed.plot()

Plot of acceleration in M/S:

Note that the ground motion for this small quake near @ivor’s house in Panama is just barely enough to register on the accelerometers—since the accelerometer channels are not very sensitive to weak motion—so this plot is very 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

#3

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

#4

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

#5

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
#6

Thank you Branden and Alejandro. Finally, I could include the instrument response.

1 Like
#9

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

#10

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

#11

@Timo @Alejandro @annunal:

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)

#12

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

#13

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