Note that you DO have to work on the Linux cluster today today! THIS WILL NOT WORK ON badger!
In an earlier lab, you went through the exercise of determining a position
from pseudoranges; by hand. Good job! Today we will start using
GIPSYx
for precise positioning. GIPSYx is a collection of many programs
and scripts to analyze GPS data (and other data). You can find them
under /usr/local/GIPSY/GipsyX_current/bin
. All of these
come with a help, that displays when you call the program on the command
line with the -h
flag. The main command we will use is
gd2e.py
.
Today focuses on setting you up for processing of GPS data. The processing will remain quite arbitrary and next week we'll move toward improving on this!
$> cd ~/geodesy
$> mkdir lab04
$> cd lab04
$> source /usr/local/gipsyx/rc_GipsyX.sh
You're now in a directory for this lab. Stay there. The last command sets up the environment for GipsyX. You could add this command as it is to your ~/.bashrc
file if you desired.
Otherwise, just run this every time you open a new terminal in which you'd want to
run GipsyX.
You should now have two new environment variables set. Make sure you get these results:
$> echo $GIPSYX /usr/local/gipsyx/GipsyX_current $> echo $GOA_VAR /gps/goa-var $>$ which gd2e.py /usr/local/gipsyx/GipsyX_current/bin/gd2e.pyIf you're getting any errors, let me know (send me the error messages if you're doing the lab by yourself.)!
The command
$> gd2e.py -h | lessgives you full help and examples for this interface to GIPSYx. (you can quit
less
by
typing q
). Many of the command line parameters are for advanced processing, which
we won't be getting into. It's worthwhile knowing about all these options, though. Note that
the actual processing engine of GIPSYx is rtgx
and gd2e.py
provides
a simple front-end to this more complicated (and more flexible) program. For now,
focus on using gd2e.py
- I am bringing this up, because some of the
output will refer to or label things with rtgx
!
In this lab, we can only cover some basics, I highly recommend that you take some time and explore this tool! It's what we use to make GPS time series!
Initially, we will just run a simple example to get used to this. We first need some
data. GIPSYx comes with a handy tool for download, which we'll use for now: rinexFetch.py
.
If you execute it on the command line with the -h
flag, you'll get some detailed information. After you explore that output, run this (may take a while):
$>rinexFetch.py -outDir . -day -t1 2019-07-03 -t2 2019-07-04 -stns p595 -srv unavcoThe output should be:
1 daily files downloaded of 2 found on unavco
. If you
run ls -R
you should get:
>$ ls -R .: y2019 ./y2019: d184 ./y2019/d184: p5951840.19d.Z
Using the rinexFetch.py
help, explain what this command does. What are the individual
parameters?
We've got some data, you're now ready to run the first solution:
$> gd2e.py -rnxFile y2019/d184/p5951840.19d.Z
The program output should look something like this:
2019-10-02 01:47:44.064 : gd2e.py started by roon on hare.giseis.alaska.edu E::1000011::Reader Event::Read ConstellationInfo file /gps/goa-var/etc/ConstellationInfo 2019-10-02 01:47:48.704 : Starting "gde /home/roon/geodesy/lab04/gde.tree", stdout=editData.log, stderr=editData.err 2019-10-02 01:47:57.509 : Finished "gde /home/roon/geodesy/lab04/gde.tree", stdout=editData.log, stderr=editData.err 2019-10-02 01:48:01.501 : Fetched GNSS products. Found GCORE style. Source: https://sideshow.jpl.nasa.gov/pub/JPL_GNSS_Products/Final. Recovered using files with date labels 2019-07-03 to 2019-07-03. 2019-10-02 01:48:02.510 : Starting "rtgx Trees/ppp_0.tree", stdout=rtgx_ppp_0.tree.log0_0, stderr=rtgx_ppp_0.tree.err0_0 2019-10-02 01:48:23.734 : Finished "rtgx Trees/ppp_0.tree", stdout=rtgx_ppp_0.tree.log0_0, stderr=rtgx_ppp_0.tree.err0_0 2019-10-02 01:48:23.862 : gd2e.py finished Elapsed run time: 0:00:39.797848
You should actually read through this output it tells you a few important things. (1) There's some
data editing going on and the output and errors are logged to editData.log
and
editData.err
. (2) It fetches
GNSS products for you (unless you specify yourself), it tells you that it got them from
https://sideshow.jpl.nasa.gov/pub/JPL_GNSS_Products/Final
(stored in the
directory GNSSinitValues
. The output / errors from the position processing run
are saved to rtgx_ppp_0.tree.log0_0
and rtgx_ppp_0.tree.err0_0
,
respectively. It's always good to check the error files to ensure all went well
(It didn't if anything FATAL occurred).
After this run, your directory will now contain many files! A listing should look like this:
$ ls allStations.xyz GNSSinitValues rtgx_ppp_0.tree.err0_0 ambres.stats gnssList.txt rtgx_ppp_0.tree.log0_0 constraints.txt P595.gde.debug.tree runAgain dataRecordFile.gz P595.gde.stats smooth0_0.tdp debug.tree postfitResiduals.out smoothFinal.cov editData.err preamb_final.pos smoothFinal.tdp editData.log preamb_finalResiduals.out stations.txt filter.tdp preamb_smooth0_0.tdp Summary final.pos preamb_smoothFinal.cov Trees finalResiduals.out prefitResiduals.out treesUsed gde.tree rinexStaDb y2019
The file we're all interested in is smoothFinal.tdp
which contains
the final position solution. We'll get into that in a bit. First, have a look at
Trees/ppp_0.tree
(use less
or something like it).
This is the full configuration file for the PPP run! At the beginning
you will see that the satellite orbit files are set (GNSSORB_POS
,
which antenna phase center correction file is to be used (GNSSXYZ_FILE
),
that we're not applying ocean load corrections in this run (OCEANLOAD
)
(more about that next time) and so forth. Even if the syntax may be a bit odd, if you
read a bit you could get a sense for what's going on here.
Another important file that was created is rinexStaDb
. Have a look!
It's a database that contains five lines, KEYWORDS, ID, STATE, ANT, RX.
The first line defines in which order the keyword lines are to be expected.
The following 4 lines identify the database record in the first column (P595),
tell us its ID is P595 (duh!), its STATE refers to its position in XYZ (ECEF)
coordinates. Copy the three numbers that end in e+06
and convert it
to latitude, longitude, height by adding that text to the end of this command:
$> xyz2llh.py -xyz
Copy the command as you ran it and its output into the lab report. Add a map (maybe Google Earth) showing the location of the station. Where is it?
The remaining two lines in the rinexStaDb
file contain information about the
antenna type and the receiver. Remember that different antennas have different
phase centers, so we better know which antenna was used.
Let's get at that precise position that came out of the processing:
$>tdp2llh.py -tdp smoothFinal.tdp | head -1
Include this result in your report. This should be different from the
result you got from xyz2llh.py
before as this is the precise position.
How different?
$>tdp2envDiff.py < smoothFinal.tdp
This shows us that the final position is 27.5 cm to the west, 7.02 cm to the west and 1.6 cm down from the original position. Let's use the more precise estimate as our starting position and add this to StaDB:
$> grep P595.State.Pos smoothFinal.tdp | head -3
These are the last three lines of the final position solution that give us
a position estimate for station P595. One line for X,Y,Z component each. The first
column is the time in seconds past Jan. 1, 2000 UTC (you can use
the program sec2date
to convert this value to something that's
more easily digestible, e.g. $> sec2date 615384000
yields 2019-07-03 00:00:00
). The second column
is the nominal value, this should be the same as in rinexStaDb
. The third
column is the actual value, i.e. the new position. Fourth and fifth column
are sigma and id, respectively. Copy the values from the third column into the
correct positions in rinexStaDB of the STATE variable (there X is 5, Y is 6 and Z is 7th columns).
To make sure you're not losing your edits, rename the station database:
$> cp rinexStaDb myStaDb
We will make use of this next time (when processing rinex files directly,
gd2e.py
will always create its own new rinexStaDb
file
from the rinex headers and ignore the one you may give to it.)
We'll now try to generate a time series. Download rinex files for the
station we just processed one day for from June 20th to July 20, 2019. Adjust the
call we used above accordingly. This will take a while. Add the call you use
to download the data to your report. Eventually, the output should show:
29 daily files downloaded of 60 found on unavco
and
$ ls -d y2019/* y2019/d171 y2019/d176 y2019/d181 y2019/d186 y2019/d191 y2019/d196 y2019/d172 y2019/d177 y2019/d182 y2019/d187 y2019/d192 y2019/d197 y2019/d173 y2019/d178 y2019/d183 y2019/d188 y2019/d193 y2019/d198 y2019/d174 y2019/d179 y2019/d184 y2019/d189 y2019/d194 y2019/d199 y2019/d175 y2019/d180 y2019/d185 y2019/d190 y2019/d195 y2019/d200
Now process all these days with gd2e.py
(see above).
After the run, make sure you save the smoothFinal.tdp
contents. Do this:
$> cp smoothFinal.tdp DATE.tdp
where you replace DATE with the date of the file you just processed.
This processing can be easily automated using a for loop in the bash shell environment. If you know how to do that, go ahead. Else, process each file individually and work on your report in between processing calls. Don't forget to save the daily result!
Once this is done, plot it as a time series in east north and vertical components (3 panels).
For this, you will first need to extract a solution for each day and put this in
a new file. Let's use the change in position from nominal. Use the second to
last line from the tdp2envDiff.py
output (again, this can be automated in a loop).
and add this to a new file. Note that the values are in cm. You may use sec2date
or another tool to convert the first column in this file to a more easily legible
date.
Is the station moving? To where and by how much and on what day(s)? What happened at these time(s) in this region? (research!)
rg <at> nmt <dot> edu | Last modified: October 03 2019 00:16.