Analysing GPS timeseries
by Andrew Newman
Caveat Emptor! : I will be updating these through the class so that they become more coherent, and a
better overall resource guide for future GIPSY 6 processing.
We have quite a few tools available, some that both Lujia Feng and I have developed. Hers are more
extensively tested, but Matlab based.
Before we begin
- Getting your data back on the GT unix system
% # on a local machine
% cd ~/GPS/MGM/my_project # or some similar directory that you want to store your data in
% sftp 'username@fembot1.rsmas.miami.edu' # insert your password on request
sftp> cd GIPS/RINEX/gdfiles
sftp> mget *.GD # answer 'a' if prompted, to get all data without further prompt
sftp> mget *.GD.ps # answer 'a' if prompted, to get all data without further prompt
- Go ahead and look at your .ps files. They may or may not be of use to you right now. In my experience
these automated plots work quite well, when the timeseries is for a signal that is highly linear, and of sufficient
duration (usually >~1 year).`
% gs *.ps # space to tab through them, or 'quit' to quit
# or
% for file in *.ps ; do gv $file ; done # 'q' to quit individual plots.
Lujia's codes
- The code library is available in ~lfeng/GPS and can either be copied to a local directory, or called directly.
If you haven't done this before, use the addpath command in matlab.
You can additionally but an addpath line in your matlab startup file.
% echo "path('~lfeng/GPS',path);" >> ~/matlab/startup.m
- HELP: Lujia wrote some well documented codes, and hence if you'd like you can look directly at them. Alternatively, you
can usually type help 'commandname' within matlab for basic usage information.
- Convert GD files to N, E, and Z: The program GPS_GD2XYZ2rneu.m will convert GD files to more human readable North, East, and Vertical format in meters.
As well, the program converts the data into a 'YYYYMMDD' and 'YYYY.decimal' format. All position information is relative now to the first data point.
% matlabcj
>> GPS_GD2XYZ2rneu('STAT.GD') % for a single station
>> % or
>> GPS_GD2XYZ2rneu('ALL.GD') % for a files in the current directory.
>> quit
- Viewing the new format
% head STAT.rneu
# rneu converted from KERA.GD
# origin is lon=25.348112716 lat=36.417845361 height=244.5785
# 1 2 3 4 5 6 7 8
# YEARMODY YEAR.DCML NORTH EAST VERT N_err E_err V_err [m]
20060507 2006.345205 0.00000 0.00000 0.00000 0.00320 0.00260 0.00980
20060508 2006.347945 -0.00133 -0.00036 0.00010 0.00180 0.00150 0.00580
20060509 2006.350685 -0.00178 -0.00045 -0.00260 0.00180 0.00160 0.00590
20060510 2006.353425 -0.00166 0.00269 0.00580 0.00180 0.00150 0.00580
20060511 2006.356164 -0.00266 -0.00135 -0.00250 0.00230 0.00200 0.00760
20060512 2006.358904 0.00377 -0.00735 0.00190 0.00320 0.00260 0.01020
20060513 2006.361644 -0.00044 -0.00260 -0.00340 0.00210 0.00180 0.00680
20060514 2006.364384 -0.00355 -0.00036 0.00210 0.00180 0.00160 0.00580
20060515 2006.367123 -0.00277 -0.00179 0.00580 0.00180 0.00150 0.00580
20060516 2006.369863 -0.00666 0.00242 0.00130 0.00180 0.00160 0.00600
20060517 2006.372603 -0.00433 -0.00063 0.00490 0.00180 0.00160 0.00600
20060518 2006.375342 -0.00433 -0.00009 0.00010 0.00180 0.00160 0.00580
% gnuplot
gnuplot> plot 'STAT.rneu' u 2:4:7 w e # East data with errors (w e)
gnuplot> plot 'STAT.rneu' u 2:3:6 w e # North data with errors
gnuplot> plot 'STAT.rneu' u 2:5:8 w e # Vertical data with errors
gnuplot> splot 'STAT.rneu' u 4:3:5 w p pt 7 # 3D position change (without time)
- Calculating Velocities: The program GPS_rneu2rate.m performs linear fitting of the time series, and has a lot more control
than the ones created at Miami. I'll repeat her help information hear for ease:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
GPS_rneu2rate.m %
linear fitting long-term trend through SITE.rneu data %
outliers are removed via two steps %
If fin_name = 'all.rneu', run this script for all *.rneu files in the current dir %
Note: 'all.rneu' is not case sensitive %
%
INPUT: %
(1) fin_name - *.rneu all in [m] %
1 2 3 4 5 6 7 8 %
YEARMODY YEAR.DCML NORTH EAST VERT N_err E_err V_err %
Note: fin_name examples %
'all.rneu'/'ALL.rneu' or 'ACOS.rneu' or 'B*.rneu' %
Do not use ? in fin_name %
%
(2) period - time period used %
period = [ sta end ] in [decimal year] %
if period = [], all available data will be used %
%
(3) feulerin_name - a file storing an euler pole %
Two formats are supported %
a. euler pole without covariance %
euler = [ pole_lat pole_lon omega ] %
omega - angular velocity [deg/Myr] %
eg. CA-IGSb00 [ 37.6 -99.3 0.247 ] from Lopez_etal_GRL_2006 %
%
b. euler pole with covariance (5x3 matrix) %
euler = [ pole_lat pole_lon omega ] [deg/Myr] Note: Units different %
[ Cxx Cxy Cxz ] %
Cw = [ Cyx Cyy Cyz ] [rad^2/Myr^2] %
[ Czx Czy Czz ] %
Trans = [ Tx Ty Tz ] [mm/yr] %
Cw is the covariance matrix for angular velocity in xyz cartesian coordinate %
x - direction of (0N,0E) Greenwich %
y - direction of (0N,90E) %
z - direction of 90N %
Cxx, Cyy, Czz - varicance [rad^2/Myr^2] %
Cxy, Cxz, Cyz - covariance [rad^2/Myr^2] %
The ITRF geocentre CM (mass center including atmosphere, oceans, and ice sheets) %
is translating relative to the solid Earth centre CE %
if feulerin_name = '', euler.vel will not be output %
%
(4) subplot_var - parameters for subplot %
= [ subplot_row subplot_col] %
if subplot_var = [], do not do subplotting %
%
OUTPUT: %
(1) *.vel for models %
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 %
Sta Lon Lat Height VE VN VU ErE ErN ErU Weight T0 T0D T1 T1D Cne %
Note: Cne is correlation coefficient between North and East uncertainties! %
%
(2) *.fitsum for fit summary %
1 2 3 4 5 6 %
Sta Lon Lat Height TT Days %
1 2 3 4 5 6 7 8 %
Sta North VN ErN w_amp f_amp wrms rchi2 %
1 2 3 4 5 6 7 8 %
Sta East VE ErE w_amp f_amp wrms rchi2 %
1 2 3 4 5 6 7 8 %
Sta Up VU ErU w_amp f_amp wrms rchi2 %
%
(3) new rneu file after %
a. removing outlies and %
b. removing an euler motion if an euler pole is provided %
It does not change erros for each data point %
%
(4) plots for individual sations %
a. raw data with outliers and b. new rneu data with outliers %
%
(5) plots for all stations %
a. raw data without outliers and b. new rneu data without outliers %
%
%
first created by lfeng Thu Oct 21 22:23:46 EDT 2010 %
considered error from euler vectors lfeng Thu Nov 4 02:24:21 EDT 2010 %
corrected a bug with site name like "BALL" lfeng Mon Nov 8 15:36:44 EST 2010 %
* could be used in fin_name %
incorporated covariance propagation lfeng Mon Feb 28 16:56:05 EST 2011 %
added the ITRF translation relative to centre-of-mass lfeng Thu Mar3 21:04:38 EST 2011%
added fitsum and subplot lfeng Sat Mar 5 12:27:36 EST 2011 %
corrected the wrong alphabetic order lfeng Sat Mar 5 22:46:46 EST 2011 %
added saving new rneu lfeng Fri Apr 1 21:21:34 EDT 2011 %
removed Ceu Cnu lfeng Sun Apr 3 20:34:00 EDT 2011 %
last modified by lfeng Sun Apr 3 20:57:17 EDT 2011 %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- A simple example to create a local velocity for a defined window without correcting for local plate motion.
>> GPS_rneu2rate('NOMI.rneu',[2010.9 2012.1],[],[])
>> quit
- An example correcting for a defined euler plate (note that your current output is probably in ITRF2008 which will be slightly
different than ITRF2005, and so the correct terrestrial reference frame is important.
If you are just looking at differential within a region, this usually isn't a problem.
I've found a very recent poster that gives new estimates of plate motions relative to ITRF2008. I haven't tested these, so be cautious/
>> GPS_rneu2rate('NOMI.rneu',[2010.9 2012.1],'~lfeng/GPS/euler_pole/EU-ITRF2005.euler',[])
>> quit
The output files are well explained by the help file. Not only do they report velocities, and formal errors
but also the white noise (w_amp), flicker (f_amp), wrms and reduced-chi2 errors, and finally the component coveriances.
More information about the errors sources can be found in Dixon et al., [2000].
- Calculating Displacements:
If you are not intersted in velocity fields but are instead trying to determine displacement vectors, the above program is not
going to be the most useful for you. You would normally want to do this if you are looking either at a discrete offset (like an earthquake slip), or
if the data is highly non-linear (such as in the case of post-seismic relaxation or volcano deformation). In these instances, it is
instead, preferential to subset the data into descrete before and after windows for analysis.
I have a program that is specific for one such instance, however you may want to modify this, or make your own for your own situation.
Because my program is rather specific, I'll post the code rather than have you mistakenly run my code for data that is not appropriate.
#!/bin/bash
# rneu_disps
# calculate appropriately weighted displacement fields using GPS data
# currently best suited for continuous GPS data
# check to make certain variable was cleared
if [ ${#} -lt "4" ] ; then
printf "Usage: `basename $0` 'rneu_file' 'Tstart' 'Tend' 'DT'
rneu_file is the GPS timeseries file of interest (usually created from GPS_GD2XYZ2rneu.m)
Tstart is the start time for calculating displacements (in decimal year)
Tend is the end time for calculating displacements (in decimal year)
DT is the window after/before Tstart/Tend that will be use for determining displacments (in decimal year)
\nFor example: `basename $0` NOMI.rneu 2010.45 2011.68 0.06 >> outfile \n" ; exit 1
fi
INFILE=$1
Tstart=$2
Tend=$3
DT=$4
printf "`basename $INFILE .rneu.rmean` $Tstart $Tend $DT "
awk -F"=| " '$2~"origin"{printf "%8.4f %9.4f %8.2f ", $5,$7,$9}' $INFILE
awk 'BEGIN{BNWGT=0; BNNUM=0; BEWGT=0; BENUM=0; BUWGT=0; BUNUM=0; BCOUNT=0; Tstart='$Tstart';DT='$DT';
ENWGT=0; ENNUM=0; EEWGT=0; EENUM=0; EUWGT=0; EUNUM=0; ECOUNT=0; Tend='$Tend'; DT='$DT'} $1!~"#"{
if($2>=Tstart && $2<=(Tstart+DT)) {
BNWGT=BNWGT+1/$6/$6; BNNUM=BNNUM+$3/$6/$6
BEWGT=BEWGT+1/$7/$7; BENUM=BENUM+$4/$7/$7
BUWGT=BUWGT+1/$8/$8; BUNUM=BUNUM+$5/$8/$8
BCOUNT=BCOUNT+1;
}
if($2<=Tend && $2>=(Tend-DT)) {
ENWGT=ENWGT+1/$6/$6; ENNUM=ENNUM+$3/$6/$6
EEWGT=EEWGT+1/$7/$7; EENUM=EENUM+$4/$7/$7
EUWGT=EUWGT+1/$8/$8; EUNUM=EUNUM+$5/$8/$8
ECOUNT=ECOUNT+1;
}
} END{
BN=BNNUM/BNWGT
BE=BENUM/BEWGT
BU=BUNUM/BUWGT
BNe=(1/BNWGT)**.5
BEe=(1/BEWGT)**.5
BUe=(1/BUWGT)**.5
EN=ENNUM/ENWGT
EE=EENUM/EEWGT
EU=EUNUM/EUWGT
ENe=(1/ENWGT)**.5
EEe=(1/EEWGT)**.5
EUe=(1/EUWGT)**.5
Ne=BNe+ENe # summation or difference of independent indeterminent errors results in the straight addition of the errors
Ee=BEe+EEe
Ue=BUe+EUe
printf "%8.5f %8.5f %8.5f %7.5f %7.5f %7.5f %4d %4d\n", EN-BN,EE-BE,EU-BU,Ne,Ee,Ue,BCOUNT,ECOUNT }' $INFILE
exit 0
The program will report its help if you use no variables. As an example:
% ~/bin/GPS/rneu_disps NOMI.rneu 2010.5 2011.5 .05
NOMI.rneu 2010.5 2011.5 .05 25.4286 36.4217 309.64 -0.01129 0.03419 0.02528 0.00082 0.00071 0.00266 18 18
#STAT Sdate Edate DT Longitude Latitude height Dnorth Deast Dvert Enorth Eeast Evert #strt #end
You may want to set up a script to call a number of stations, and
pipe the output from a number of these calls into a single file.
|
gatech.edu | Updated:
Mon Feb 13 14:22:43 EST 2012