Converting counts in acceleration

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

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

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:

image

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.

Thanks for this information. We however are using the UPD to get the counts and calculate in another computer the velocity and acceleration in real time continously.

Therefore we prefer to start from the counts and obtain the acceleration and the velocity using the relation between counts and acceleration or velocity

In the Technical specification we have found the following

Sensitivity: V6: 3.9965E+08/meter/second

Is it correct that in reality should be
V6=3.9965E8 / (meter/second) ?

So that if we have 34560 counts, it corresponds to:

v=34560/3.9965E08 m/s=8.64e-5 m/s

Similarly, as you pointed out in another post, for the acceleration (ENE, etc), it would be if the counts are 67234, the acceleration would be:

acc=67234/2.845e5= 0.167 m/s2 -> 0.017 g

I would really appreciate your comment
v=

hello alessandro,

instead of providing just the formula for this type of conversion, i would like to direct you to the following online resource which provides a complete explanation of what you are looking for:

cheers,

richard

Thanks Richard, very useful. My assumptions were correct

2 Likes