Proof of concept: Using scipy instead evalresp #1718
Conversation
megies
added
.signal
|
+1 ! |
ghost
mentioned this pull request
|
I would like to confirm that this pull request is ready / awaiting review. |
uhm.. it handles all different types of response stages out there already?? |
I believe so. I tested it rather thoroughly with @aringler-usgs in order to make sure it did. (There was a lot of back-and-forth about how to handle digital vs. analog filter stages.) |
|
Removing the response is such an integral part of ObsPy that we absolutely have to make sure that we can explain every difference to the way evalresp does it. When we initially implemented the StationXML response removal and being able to read SEED files we did by comparing the output of raw evalresp to almost every response file in existance: https://github.com/obspy/sandbox/tree/master/the_great_response_test I think something similar would be needed here in addition to all the testing you already did. I can take care of setting up a basic version - should not take too long. That is the most up-to-date branch here, correct? https://github.com/amkearns-usgs/obspy/tree/evalresp_update |
|
Yes. We haven't been rebasing it to any updates on the main branch at all -- I don't think we've touched this since like April or so, but that's where all the updates would happen. |
|
@krischer thanks for setting up an additional location for tests. I am sure there will be a few oddities that come out of this. I will work with @amkearns-usgs to deal with the various issues in the tests. Thanks! |
|
Great job @amkearns-usgs & @aringler-usgs !! This already works really really well and given how much time I spent on this problem in the past I personally really appreciate this. So thanks a lot already. I pushed a first version of the aforementioned test here: https://github.com/obspy/sandbox/tree/master/the_great_response_test_scipy I hope it is easy enough to understand. The main idea is to run this against a significant portion of the data out there and see if the results are identical between evalresp and the new scipy implementation. I'd be fine with some differences as long as they are documented and rational. |
|
I realize that this might be a decent amount of work. But given that response removal is one of the absolute core features of ObsPy we cannot afford to go wrong here. I'm happy to take some of the work. Given that you know the code better than me I think it would be most useful if you try to deal with the initial broader issues and I'll jump it at any time when you need a break. |
|
I also pushed a tiny commit here: https://github.com/amkearns-usgs/obspy/tree/evalresp_update |
|
Let's continue the discussion in #2592 - more useful when the code is directly there. |
Based on the discussion in #1493 it came to me that it would be almost trivial to replace evalresp completely with scipy. I don't know why I did not think of this before. Much simpler than our current wrapping of evalresp.
This PR does this for the PAZ, decimation, and coefficient response stages which are the most common and it also handles the final assembly and sensitivity calculations so it can already do quite a lot of responses. The rest is quite easy to add.
Of course this needs a lot more testing so I'd personally not include this in the upcoming 1.1 but rather wait a release.
The amplitude response matches the evalresp one to a relative error of 1E-10 - the sensitivity scaling makes this a bit different but I think what's done in this PR is correct. The phase differs a bit more for the FIR stage mainly after nyquist but also a tiny bit before. Again I think here scipy does the correct thing as it for example has no Gibb's artifacts which evalresp has.