GPS Tool Box

MATLAB Tools for viewing GPS

velocities and time series
Thomas Herring

Abstract Over the past decade, many Global Posi-

Introduction
tioning System (GPS) networks have been installed to The generation of tectonic velocity fields with large numbers
monitor tectonic motions around the world. Some of of GPS sites and complicated time series is becoming com-
these networks contain hundreds of sites spread mon. The complications in the time series are often associ-
across active tectonic margins where the differences ated with geophysical effects such as offsets due to
in velocities across the network can be 50–100 mm/ earthquakes, postseismic transient behavior after earth-
year. For networks that have been running for a quakes, and more noise-like phenomena such as the effects
number of years, the uncertainty in the velocity of groundwater changes. Complications in time series can
estimates can be less than 1 mm/year. In some cases also arise from offsets due to GPS antenna and receiver
the vertical motions can also be significant and of changes and the addition or removal of antenna raydomes
importance. Often, the time series of the motions of from a GPS site. In some cases, GPS receivers fail in modes
the GPS sites show complex non-linear behavior, and that generate reasonable carrier-phase measurements but
in all cases the statistical model of the time series is corrupt position estimates. The GGMatlab tools are designed
more complex than simple white noise. In this article, to allow both the generators and users of large GPS velocity
we describe a set of Matlab tools developed for use fields to examine in detail the velocity fields and the time
with the GAMIT/GLOBK GPS data analysis system series from which the velocities were generated. The tools are
(King 2002; King and Herring 2002) that allow inter- documented and available from http://www-gpsg.mit.edu/
active viewing and manipulation of GPS velocities and tah/GGMatlab. The current version is 1.02 but the tools are
time series with a Matlab-based graphical user still being developed and future versions are expected.
interface (GUI). The formats of the data files used by The GGMatlab toolbox currently contains two main
the tools are specific to GAMIT/GLOBK, but they are components: velview, which allows viewing and analysis of
simple ASCII files that can be generated from other velocity fields; and tsview, which allows viewing and
file formats. The tools are referred to as GGMatlab. manipulation of time series. The time-series viewing tool
can be invoked from within velview by clicking on a
plotted velocity vector. Within velview, profiles of veloci-
The GPS Toolbox is a column dedicated to highlighting algorithms ties across regions can be generated and functions fit to the
and source code utilized by GPS Engineers and scientists. If you velocity changes along the profiles. Velview allows two
have an interesting program or software package you would like different velocity field files to be overlain and the statistics
to share with our readers, please pass it along; e-mail it to us at of the differences and biases between them calculated.
[email protected]/. To comment on any of the source All of the details of the options and installation of
code discussed here, or to download source code, visit our website GGMatlab are given on the Web page. Here, we describe
at http://www.ngs.noaa.gov/gps-toolbox. This column is edited
by Stephen Hilla, National Geodetic Survey, NOAA, Silver Spring, some of main features and the typical scenarios for using
Maryland, and Mike Craymer, Geodetic Survey Division, Natural the packages. The package is designed to work with Matlab
Resources Canada, Ottawa, Ontario, Canada. For the sidebar, see release 12 and greater. Linux executable versions (kernel
the Volume 6, Number 4, 2003 issue of the GPS Toolbox column. version 2.4.7 and greater) are also available. The execut-
able versions do not support all the features available when
the package is run from within Matlab.
Received: 10 August 2003 / Accepted: 18 August 2003
Published online: 1 October 2003
ª Springer-Verlag 2003
Velview
T. Herring The primary aim of velview is to allow plotting, manipu-
Department of Earth, Atmospheric and Planetary Sciences, lation and comparison of velocity fields generated by
Massachusetts Institute of Technology, 77 Massachusetts Avenue, GLOBK or other programs if the appropriate format files
Room 54-618, Cambridge, MA 02139, USA
E-mail: [email protected] are generated. An example of a velview screen view is
Tel.: +1-617-2535941 shown in Fig. 1 with some of the features invoked. The
Fax: +1-617-2531699 example shown is taken from the combined analysis of the

194 GPS Solutions (2003) 7:194–199 DOI 10.1007/s10291-003-0068-0

 GPS Tool Box

Fig. 1 difference columns are not used, but in future versions

Example of a velview screen image. In this example, two velocity fields they are likely to be used to display differences with re-
have been loaded from the SCIGN analysis. The primary field is
displayed in red, and the secondary in blue. The ‘‘Feature File’’ option spect to models rather than total motions. For GAMIT/
has been used to add the California coast (black), the major faults in GLOBK users, the velocity field files can be extracted from
the area (green lines) and the sesismicity (brown dots, earthquake the GLOBK output files.
magnitudes greater than 2 since 1996). The velocity vector of The main features of the velview program allow the
LDES_CHT has been labeled on the figure in a yellow box (box will be
removed if it is clicked on), and the vector options are shown for the
detailed analysis of velocity fields. Two velocity field files,
site ISLK_GHT (grey box; the grey box will disappear once an option referred to as the primary and secondary velocity fields,
is selected or the site name clicked.) can be displayed and compared. A simple map projection
is used to display the site locations. The longitudes are
GPS data from Southern California Integrated GPS scaled to linear units by multiplying by the cosine of the
Network (SCIGN) (Watkins et al. 1997) (http:// latitude in the center of the displayed region. (In polar
www.scign.org). The velocities and time series shown here regions, this projection style will be problematic.) There is
are available from http://www-gpsg.mit.edu/tah/ no wrapping of the longitudes across the Greenwich
SCIGN_MIT. (This web page is regularly updated as new meridian. For work in this area, negative longitudes should
data are added to the analysis). The input velocity field be used west of the Greenwich meridian. When velocities
files for velview are simple ASCII files containing column are loaded, limits can be placed on the sigmas of the dis-
data in free format. Data lines start with at least one blank; played data so that vectors with large uncertainties are not
non-blank characters in column 1 denote comments. The displayed. If vertical rates are also being displayed, the
entries on each line are site longitude (positive east, vertical rate sigma will also be checked. The vertical rates
decimal degrees), latitude (decimal degrees), east velocity, are displayed as north or south lines (positive velocity
north velocity, the differences between the east and north north) without arrowheads and with error bars. The scale
velocity and model values, standard deviations for the east of the vectors can be set separately for horizontal and
and north velocities, correlation between the east and vertical velocities and the velocity scale is displayed in the
north velocities, the vertical velocity, difference between lower left hand corner of the figure window. The nominal
the vertical velocity and a model, the standard deviation of scale of 1 will display the vectors at real scale on most
the vertical velocity, and the 8-character site name. All computer screens (this relationship is not exact because it
units for the velocities are mm/year. There are 13 columns depends on screen resolution and physical size). The
in the format. In the current version of velview, the

GPS Solutions (2003) 7:194–199 195

 GPS Tool Box

length of the arrowheads is specified in millimeters, which program. If the time series files are not in the same directory
is useful for judging the size of very small velocities. The from which velview was run, the name of the directory can
number of sites being displayed for each of the velocity be specified in the box next to the ‘‘Time Series’’ label at the
fields is shown at the top of the figure. top of the figure. Tsview is discussed below.
Once the velocities fields are loaded, there are options for Another feature of velview is the ability to select a transect
viewing and manipulating the fields. There are two types of across the velocity field for profiling where the velocities are
zoom available. The ‘‘Vec Zoomin’’ feature allows a click- shown as a function of distance along the profile and are
and-drag selection of an area to be viewed in detail while projected into directions based on the azimuth given in the
the vectors maintain their original scaling. The ‘‘Zoom In’’ box next to the ‘‘AZ’’ button. The azimuth of the projection
feature rescales the vectors during the zoom so that their can be either typed into box, graphically selected by clicking
lengths increase with zooming (this is the standard zoom the ‘‘AZ’’ button when the box is empty, or default to along
feature in Matlab). Both zoom features maintain a stack of the profile direction and normal to the profile direction. The
previous zooms in the zoom out buttons so that previous profile is selected by first drawing the profile centerline and
views can be recovered by successive zoom outs. (If both then by graphical selecting its width. The selected profile
types of zooms are used, the zoom out feature should be box is displayed on the figure and the sites inside the box are
used in the reverse order of the zoom ins.) At anytime, marked with yellow circles. The profile box information is
numerical values can be entered into the arrow length placed on a stack as profiles are selected. The boxes and site
boxes, and the redraw button used to return the vectors to markings on the figure can be removed from the top of the
well-defined lengths. The ‘‘Pan’’ box arrows and minus stack using the ‘‘POP’’ button (most recent removed) or
sign can be used to move the view by half the current from the bottom of the stack (earliest removed) using the
width and zoom out by 50% of the current size. ‘‘PIP’’ button. Only results from the primary velocity field
Manipulation of the velocity fields includes the ability to are shown on the profiles. (The ‘‘SWAP’’ button at the bot-
apply offsets to the motions either by typing the values into tom of the velview window can be used to exchange the
the box next to ‘‘Offset’’ or, when this box is empty, to click- velocity fields thus allowing profiles for each field to be
and-drag a region in which the weighted mean velocities will plotted in separate windows). An example of the profile plot
be zero after the offsets are applied. The offsets are com- is shown in Fig. 2. Part (a) of the figure shows the velview
puted from the primary velocity field and applied to both window used to select the profile and part (b) shows the
primary and secondary (if being displayed) fields. A dif- profile figure itself. For this figure, the azimuth for the
ferential offset between the two velocity fields can be cal- projection of the velocities along the profile was graphically
culated from a click-and-drag region or specified directly in selected to be parallel to the CICE velocity vector.
the box next to the ‘‘Align’’ button. Both the offset and align In Fig. 2b, we have used the ‘‘Fitting Functions’’ window to fit
options display the results in boxes on the figure and as text a function to the profile velocities based on the co-seismic
in the Matlab or terminal window. The boxes on the figure strain accumulation for a locked strike-slip fault (the SS
can be removed by clicking on them. The offsets are applied function) (Okada 1985). When velview is run in Matlab, any
as simple additions to the velocity components. For large valid Matlab function, defined to be a function of distance
areas, these offsets cannot adequately represent the effects of along the profile, can be used in the fitting functions. In the
rotations and translations between two reference frames. In executable version, only specific functions can be used with
this case, the GLOBK utility programs cvframe and velrot or the pre-defined SS function being available in both the
the userÕs own tools should be used to define reference Matlab and executable versions. The green shaded boxes on
frames and to align velocity fields in different frames. the left of the figure show the results of the fit in the direction
Other manipulations available in velview include the removal of the azimuth (labeled ‘‘along track’’), normal to this
of duplicate sites with exactly the same velocities and sigmas direction (labeled ‘‘cross track’’) and for the vertical velocity.
and within the specified distance of each other (given in the Because space is limited in these text boxes, the units are not
box next to the ‘‘RmDups’’ and ‘‘Stats’’ boxes). The removal given. For the weighted root mean square (WRMS) scatter
of duplicates is useful for GLOBK analyses where sites over and the amplitudes of the estimated coefficients (FCN0 and
some intervals of time may have been renamed because of FCN1 in the figure), they are in mm/year. The FCN0 value is
possible changes in position, but the velocities of the multiple the offset of the profile velocities. The normalized root-
site names have been forced to be equal. The ‘‘Stats’’ button mean-square (NRMS) scatter is dimensionless and is the
computes the weighted root mean square (WRMS) differ- ratio of the scatter to the expected scatter based on the
ences in velocities between sites within each velocity field velocity error bars (square root of v2-per-degree-of-free-
that are separated by less than the distance given in the dom). A single fault with a deep locking depth has been used
box next to the button. The ‘‘ID Points’’ list all sites in a to fit the profile in the example but a better fit can be obtained
click-and-drag region and is useful when multiple sites are by fitting to multiple faults with shallower locking depths.
overlaid on each other and clicking on a specific vector is Individual sites may be deleted before the fitting functions
difficult. The list of sites appears both in a box at the top of are estimated and they appear with a cross over the point.
figure (removed by clicking on it) and in the terminal window. The results of the fit and the profile information can be
The velocity vectors plotted by velview can be selected and saved to a file by putting the file name in the yellow ‘‘File
the grey popup menu (example shown in Fig. 1) can be used box’’ and clicking the save button. Information may be
to obtain information about the vector. If position time extracted from this file so that it can be plotted in another
series are available, they can be plotted using the tsview graphics program if desired.

196 GPS Solutions (2003) 7:194–199

 GPS Tool Box

Fig. 2 coordinate does not change as a site moves by small

a The velview window used to select the profile sites shown in Fig. 2b. amounts in latitude (the files are generated with the GAMIT/
The solid straight-line shows the center of the profile and the dotted
lines mark the region of the sites to be included in the profile. The GLOBK program multibase). An example of the tsview
profile starts in the northeast and ends in the southwest. b The profile window is shown in Fig. 3 after detrending the time series
window generated by velview when the profile button is selected. The and removing some outliers.
SS function in the ‘‘Fitting Functions’’ models the behavior of an When tsview is run, it checks the directory in which is was
infinite strike-slip fault. In the example, the fault is locked to a depth
of 20 km and crosses the profile 200 km from the start of the profile.
run and looks for files with names of the form
The blue line in each frame shows the model values and the red line in mb_XXXXXXXX.dat1 where XXXXXXX is an 8-character
each frame is the mean value of the velocities. The green boxes on the GPS site name. This file naming convention is the standard
left show the results of the fit and can be removed from the figure by name given to the north coordinate time series files
clicking on them. All error bars are one standard deviation generated by the suite of GLOBK programs. The east and
height components have the same style names and end
Tsview with the numerals 2 and 3. If tsview is invoked from within
velview, the same list is formed and tsview tries to match
Tsview is the time series viewing and manipulation module the site name passed to it from velview. The time series
in the GGMatlab toolbox. Tsview can be invoked by select- files are simple ASCII files with three header records (not
ing ‘‘TimeSeries’’ in the popup menu that appears when a used by tsview) and data records that consist of a time in
velocity vector is clicked in velview or it can be run as a decimal years, and a NEU coordinate and sigma in meters.
stand-alone module. Tsview plots times series of GLOBK When invoked in stand-alone mode, a site name is selected
north, east and height (NEU) coordinate estimates from the list on the left of the screen and ‘‘load’’ button is
contained in GLOBK time series files. The north coordinates used to load the time series. After a series is loaded, data
are the geodetic latitudes of the sites multiplied by the from another site name may be appended. The append
WGS84 semi major axis, the east coordinates are the option can also be used when tsview is invoked from
distances from the Greenwich meridian along the small velview. In standard GLOBK GPS processing, breaks in a
circles at (quantized) latitudes of the sites, and the heights time series, due, for example to earthquakes or equipment
are WGS84 ellipsoidal heights. The latitudes of the sites are changes are accounted for by changing the last 2 or 3
quantized for the east coordinate calculation so that the east characters of eight-character site name. Segmented time

GPS Solutions (2003) 7:194–199 197

 GPS Tool Box

Fig. 2 b the time series and the statistical properties of the time series
(Contd.) residuals with either white or correlated noise assumptions.
If the ‘‘RealSigma’’ (meaning realistic sigma) button is
series of this type can be shown as single time series in unchecked, a white noise model is used. The white noise
tsview. Once a time series is loaded, there are numerous assumption almost invariably generates very optimistic
options for manipulating the time series. estimates of the uncertainties of the parameter estimates.
The primary aims of tsview are to assess the quality of time When the RealSigma box is checked, a time-correlated noise
series and to generate files that can be used to control the model is used to estimate the parameter uncertainties. With
treatment of data from specific sites at specific times when this model, a correlation time of the residuals for each
data are combined in GLOBK. The center set of buttons in the coordinate component is estimated by computing the in-
tsview window control editing and fitting of functions to the crease in the chi-squared-per-degree of freedom of succes-
time series. The function boxes in the lower left hand part of sively longer time averages of the residuals. For a white noise
the window control mainly the output files from tsview. error model, the chi-squared-per-degree of freedom would
Breaks in time series can be treated in multiple ways in not depend on averaging time. With temporal correlations
tsview. The break can be a simple offset in the times series in the time series, chi-squared-per-degree of freedom
or it can be an offset followed by either an exponential or increases as the residuals are averaged over successively
logarithmic functions to represent the time evolution after longer time intervals. The character of the averaged residual
the break. When exponential and logarithmic functions can be viewed in tsview using the ‘‘Average’’ button. The
are used, a change in the secular rate can also be estimated value in the box below the button gives the length of time in
after the break. The main use of these different function days that will be averaged. The statistics of the averaged
forms is to study the behavior of time series after earth- residuals are shown in bold on the display. For the example
quakes. Annual and semiannual components in the time shown in Fig. 3, the realistic sigma estimates of the secular
series can also be estimated. By selecting the ‘‘Linear Only velocity in each coordinate of this site are 3 to 5 times larger
Residual’’ option, the residuals to the fit can be displayed than the white noise estimates. Thirty-day averages of the
with only the linear motions and offsets removed. In this coordinate residuals have scatters of 0.4 mm in the
mode, a thick black line shows the time variations of the horizontal components and 0.8 mm for the vertical
estimated model parameters and the red model uncer- component. These scatters are only a factor 2 smaller for
tainty lines are placed plus and minus one standard horizontal coordinates and a factor 4 smaller for the vertical
deviation about this black line. coordinate than the scatters of daily values. If the residuals
The uncertainties of the parameter estimates in tsview are were white noise, the scatters should have reduced by a
calculated using the sigmas of the coordinate estimates in factor of 5.5 (square root of 30).

198 GPS Solutions (2003) 7:194–199

 GPS Tool Box

Fig. 3 they are used.) The edit file from tsview can be used in
Example of a tsview screen image. The time series shown is for the site GLOBK to remove edited sites from an analysis and to
DHLG which is one of the many sites in Southern California affected
by the October 1999 Hector Mine earthquake (Hudnut et al. 2002). In implement discontinuities in time series through renaming
the GLOBK processing, this site is divided into a pre-Hector mine sites. The lowest box in this part of the figure allows an output
name (DHLG_GLA) and a post-Hector mine name (DHLG_GHT). file from the GLOBK ensum program to be used to sort the
The plot shows positions from both site names appended together list of sites based on the statistics of their fits to linear mo-
with a break in late 1999 (1999/10/16). A logarithmic function is used
to model the behavior of the site position after the earthquake. The
tions with the poorest fit sites appearing at the top of the list.
green horizontal lines show the bounds of 3-times the WRMS scatter
of the residuals. The red, nearly horizontal lines near zero on the plots
show the 1-r limits of the model time series based on the standard
deviations of the parameter estimates with the contribution from the References
overall mean value of the time series removed. (When only a rate is
estimated, this display forms a thin ‘‘x’’ cross shape that crosses zero
Hudnut KW, King NE, Galetzka JE, Stark KF, Behr JA, Aspiotes A,
near the center of the plot.) The full set of parameter estimates and
van Wyk S, Moffitt R, Dockter S, Wyatt F (2002) Continuous
their standard deviations are written to the Matlab workspace and to a
GPS observations of postseismic deformation following the 16
pop-up window in the upper right hand corner of the computer October 1999 Hector Mine, California, earthquake (Mw 7.1),
screen. In this figure we have used the ‘‘Realistic Sigma’’ model Bull Seis Soc Am vol 92, no 4, pp 1403–1422
discussed in the text
King RW (2002) Documentation for the GAMIT GPS analysis
software, MIT Internal Report, 206 pp (http://www-gpsg.mit.-
For users of GAMIT/GLOBK, the results from tsview can be edu/simon/gtgk/GAMIT.pdf)
used in later GLOBK analyses of the data set. The ‘‘Save’’ King RW, Herring TA (2002) Global Kalman filter VLBI and GPS
button in center column will save the results from the time analysis program, MIT Internal Report, 98 pp (http://www-
gpsg.mit.edu/simon/gtgk/GLOBK.pdf)
series being displayed in the files whose names are given in
Okada Y (1985) Surface deformation due to shear and tensile
the bottom part of the left-hand column. (The small W faults in a half space. Bull Seismol Soc Am 75:1135–1154
buttons show the files to be written, and the E buttons Watkins MM, Bock Y, Hudnut KW, Prescott WH (1997) The
indicate which files will be erased before they are written. The Southern California Integrated GPS Network: Status Report, Eos
erase buttons reset automatically to non-erase mode once Trans AGU, AGU Spring Meeting, 29 April 1997, p 105

GPS Solutions (2003) 7:194–199 199