University of Alaska Fairbanks
Department of Geosciences

GEOS F493 / F693 - Geodetic Methods

Lab 4: gd2e.py: static position estimation

"I like my crust deformed."
UNAVCO bumper sticker

Note that you DO have to work on the Linux cluster today today! THIS WILL NOT WORK ON badger!

Introduction

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!

Basic Setup

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.py
 	
If 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 | less
	
gives 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!

A simple Example

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 unavco
    
The 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.)

A more involved example - time series

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!)

Deliverables: (submit via blackboard!)

rg <at> nmt <dot> edu | Last modified: October 03 2019 00:16.