Converting counts in acceleration

Dear Tech team,

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

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')

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

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.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()

``````

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.