Note that you DO have to work on redoubt for some things today!
Last week, you wrote some code to estimate time series properties. This week, you'll use some other software to do that. Now that you know how to process data into position estimates and can convert them into, e.g., a local coordinate system, we make the assumption that you can repeat that for every day that you have data. If you concatenate the daily position solutions you get a time series and from that, we can finally start learning something about the long and short term geologic processes that affect the site.
The goal for this week is familiarize yourself a bit with a the program CATS (Create and Analyze Time Series) that analyzes time series by fitting a multi-parameter model to the data.
Change to your working directory and copy the test data:
$> cd $GEOP555
$> mkdir lab11
$> cd lab11
$> cp -r /usr/local/CATS/Bezy_2 ./
A directory listing should look similar to this:
redoubt:lab11> ls Bezy_2/
BEZD.cats BEZR.cats BZ01.cats BZ03.cats BZ05.cats BZ07.cats BZ09.cats KAMD.cats KLU.cats MAYS.cats
BEZH.cats BZ00.cats BZ02.cats BZ04.cats BZ06.cats BZ08.cats ES1.cats KBG.cats KLUC.cats PETS.cats
For now, we will limit ourselves to working only on 1 file - BZ09.cats, which contains the timeseries in local NEU format for a station near Bezymianny Volcano, Kamchatka (beautiful spot!). The contents of the files look like this:
2006.7644 0.00000 0.00000 0.00000 0.00110 0.00140 0.00320
2006.7671 -0.00579 0.00599 -0.00170 0.00100 0.00130 0.00290
2006.7699 -0.00412 0.00231 -0.00530 0.00100 0.00130 0.00290
2006.7726 -0.00468 0.00318 0.00380 0.00100 0.00130 0.00290
2006.7753 -0.00936 0.00212 -0.00100 0.00100 0.00130 0.00290
2006.7781 -0.00590 0.00050 -0.00310 0.00100 0.00130 0.00290
2006.7808 -0.00590 0.00412 0.00060 0.00100 0.00130 0.00290
2006.7836 -0.00479 0.00456 -0.00260 0.00100 0.00130 0.00290
2006.7863 -0.00590 0.00300 0.00160 0.00100 0.00130 0.00290
2006.7890 -0.00991 0.00225 -0.00180 0.00100 0.00130 0.00290
2006.7918 -0.00991 -0.00162 0.00620 0.00100 0.00130 0.00280
...
Where columns are:
decimal year | North | East | Up | North error | East Error | Up Error (all meters)
Plot the north, east and up components of station BZ09 over time in 3 panels in a single figure. Note that this should be a discrete plot, i.e. disconnected dots. With an error bar function, you add the respective errors.
The benefit of CATS is that it allows you to choose between a range of stochastic models for noise approximation. Traditionally, we've been using white noise approximations to estimate our uncertainties in the site velocities. By now we know that the noise characteristics are more complicated than this.
You can find a lot of the detail on the various stochastic models in the manual. In addition to estimating noise parameters, the program also allows you to estimate slopes, and offsets at given times, as well as annual and semiannual parameters. Let's do that!
In your lab directory run:
$> cats Bezy_2/BZ09.cats --model wh: --verbose
This estimates, for this station, the intercept, slope and white noise parameters according to the models in the manual, for all three components:
... +NORT MLE : -2859.806080 +NORT INTER : 2.4353 +- 0.2396 +NORT SLOPE : -13.6079 +- 0.0818 +NORT WHITE NOISE +NORT WH : 2.9409 +- 0.0615 ...
(units are mm) Let's see what happens to these parameters with a more realistic noise model, when we add flicker noise.
$> cats Bezy_2/BZ09.cats --model pl:k-1 --model wh: --verbose --output BZ09.out
The output is being redirected into BZ09.out. You can get the results to the screen with:
$> grep "^+" BZ09.out
Having the slopes and intercepts for both models, add them to your plot using the line equation:
y = a+b*t
Where t is time (set up a vector over the times given in the file or draw a line for 2 points in time), a is intercept and b is slope. Note that the intercept is defined at t0, which is the start of the first year of data (here: 2006.0), so you will have to subtract that value from your times in the above equation.
Some other parameters that we are often interested in are those for annual and semi-annual sinusoids, e.g., seasonal loading:
$> cats Bezy_2/BZ09.cats --model pl:k-1 --model wh: --verbose --sinusoid 1y1 --output BZ09_seasonal.out
The sinusoid parameter allows to add the period, time flag and harmonics. Here, the period is 1 year (1y) and we're interested in the first harmonic (1y1), i.e. the semiannual solution. Note that CATS appends to output files! So either remove the old one or create the new output file as I named it here: BZ09_seasonal.out.
The result file will now contain information on COS/SIN components. You could add these to your plots.
What happens to your velocity estimate if you have an actual offset in the data? Let's add an artificial offset:
$> awk '{if ($1 > 2008.0) {print $1, $2+1, $3-0.75, $4+0.5, $5, $6, $7} else {print $0;}}' Bezy_2/BZ09.cats > Bezy_2/BZ09_steps.cats
Plot the new file BZ09_steps.cats
Now re-estimate the velocity, e.g,:
$> cats Bezy_2/BZ09_steps.cats --model pl:k-1 --model wh: --verbose --output BZ09_steps.out
How has the velocity changed? Is this a realistic estimate for the long-term velocity of the site?
For known steps in time, CATS has a feature to help estimate the offsets and still get a long-term velocity. Create a new file Bezy_2/BZ09_step_fixed.cats and add information on the offset:
$> cp Bezy_2/BZ09_steps.cats Bezy_2/BZ09_steps_fixed.cats
To the very top of Bezy_2/BZ09_steps_fixed.cats add the line:
# offset 2008.0 7
In the manual, find what the 7 means
Now rerun the solution and compare the output. You should find information on the offset in all three components. How do the velocities compare to the "no-earthquake" data and the estimates without offset correction? What is the estimated offset for the components? Why isn't it exactly what we put in?
$> cats Bezy_2/BZ09_steps_fixed.cats --model pl:k-1 --model wh: --verbose --sinusoid 1y1 --output BZ09_fixed.out
rg <at> nmt <dot> edu | Last modified: December 18 2017 17:46.