Note that you DO have to work on redoubt today!
Today we shift our focus from static positioning and looking at phase residuals for the individual processing parameters to kinematic, i.e. subdaily, processing. Following the detailed processing instructions in class today, you are well-equipped to take on a more complex problem. You will process 5Hz data recorded at 2 GPS stations near the 25-April-2015 Mw 7.8 Gorkha earthquake in Nepal.
Below I show the sequence of commands that should get us all to the same starting point.
Note that anything between redoubt and > is part of the prompt
and not to be typed. The commands you will use start at the first space and range until the
'#', which I use to comment what that line does.
redoubt:~> cd $GEOP555 #go to working directory
redoubt:/data/GEOP555/gps> #'gps' part should be your user name!
#I am working as user 'gps'
redoubt:/data/GEOP555/gps> mkdir lab06 #create folder to work in this week
redoubt:/data/GEOP555/gps> cd lab06 #go there
redoubt:/data/GEOP555/gps/lab06>setenv WORKDIR `pwd` #creates the variable WORKDIR to refer
#to this place
redoubt:/data/GEOP555/gps/lab06>echo $WORKDIR #print the contents of this variable
/data/GEOP/gps/lab06/ #this should be your output I will use
#this variable below!
Whenever you log out and log back in, make sure to have $WORKDIR set to the value above. Obviously you don't
have to create the lab06 directory every single time, as it will remain there.
If your $WORKDIR contains gps, you're doing it wrong!
Make a new directory for this run:
$> mkdir $WORKDIR/gorkha
$> cd $WORKDIR/gorkha
Now you'd need to get rinex files for the 2 stations KKN4 (Kakani 4) and NAST (site in Kathmandu). You can get the data from UNAVCO's high-rate archive. It was quite the effort to get the data off the receivers before it was overwritten (5Hz data on most receivers are stored in a ringbuffer that overwrites data older than usually 2 weeks). Remember, though, that the data are organized by YEAR and DAY OF YEAR directories. First, get the day of year value for the day of the earthquake (search for it, if you don't know it). We have a program that does just that conversion:
$> yymmmdd2doy YYMMMDD
The argument should be the 2-digit year, first 3 letters of the month, day of month. All in one string.
Armed with the correct day of year (doy) you can now ahead and download the data for the 2 stations (make sure to replace DOY with the value you've just calculated):
$> wget ftp://data-out.unavco.org/pub/highrate/5-Hz/rinex/2015/DOY/kkn4/kkn4\*.15d.Z
$> wget ftp://data-out.unavco.org/pub/highrate/5-Hz/rinex/2015/DOY/nast/nast\*.15d.Z
Feel free to explore the ftp-sever in a browser! There's lots of good stuff available and it's good to know how the data are organized once you need them.
You should now have 2 files in your directory (the stars indicate the DOY - I am not quite giving this away :) ...
$> ls
kkn4***0.15d.Z nast***0.15d.Z
Unzip them and convert them from compressed rinex (*.15d) to uncompressed rinex (*.15o)!
$> gunzip *.Z
$> crz2rnx nast*
$> crz2rnx kkn4*
You should now have 4 files in your directory (the stars indicate the DOY - I am not quite giving this away :) ...
$> ls
kkn4***0.15d kkn4***0.15o nast***0.15d nast***0.15o
Add necessary information for these two new stations to your sta_info database!
$> grep APPROX *o #gives you the apriori positions for sta_pos
$> grep ANT *o #gives you the antenna information for sta_svec
#(remember, 9 characters!,
# and the height of the antenna - DELTA H/E/N )
Since you have all the necessary information on the screen to create the antenna-phase-center correction files, lets go ahead and do that. Here is the example for NAST:
$> antex2xyz.py -antexfile $GIPSY/pcenter/latest.atx -xyzfile NAST.xyz \
-radcode SCIT -anttype TRM57971.00 -fel 0 -daz 5 -del 5 \
-extrap -recname NAST
To do the same for KKN4 you, obviously, need to replace all occurrences of NAST with KKN4, and replace the antenna type with the correct value
from KKN4's rinex file. The radome code SCIT remains the same for both! At this point you should have the following files in your
directory (still not giving the doy away!):
$> ls
kkn4***0.15d kkn4***0.15o KKN4.xyz nast***0.15d nast***0.15o NAST.xyz
Now add the 2 sites NAST, KKN4 with their respective information to your sta_id, sta_pos, sta_svec files!
If you followed the instructions in previous labs, you should find your sta_info database at $GEOP555/sta_info.
Yes, this is slightly tedious, and you may have to go back to LAB04 where we first did this (remember how important the formatting is!).
OK. Now for the orbits. This event was a few months ago now and we can expect the final orbits to be available. Use goa_prod_ftp.pl
to download the flinnR product for that day. Make sure to use the -hr flag! Your directory contents should be looking like this:
$> ls
2015-04-25.ant.gz 2015-04-25.shad.gz KKN4.xyz
2015-04-25.eo.gz 2015-04-25.tdp.gz nast***0.15d
2015-04-25.frame.gz 2015-04-25.wlpb.gz nast***0.15o
2015-04-25_hr.tdp.gz kkn41150.15d NAST.xyz
2015-04-25.pos.gz kkn41150.15o
I pre-calculated the ocean loading coefficients for you. Just as in previous cases, they are deposited at
/usr/local/GIPSY/ocnload/tpxo72atlas_CM/gipsy/SSSS.ocnld (where SSSS is to be
replaced with the 4-char station id!). If you wanted to do that yourself on my system for stations that don't have
ocean loading coefficients yet, you could invoke:
$> calc_oceanload.tpxo72atlas_CM
Usage: calc_oceanload.tpxo72atlas_CM 4charid lat(deg) long(deg) height-msl (m)
This will deposit a file SSSS.ocnld in the directory I gave above.
Well, now you've got everything lined up to process your data. Note that the organization of the files is suboptimal! I had you dump everything in one directory, because it's easier to show. However, for your own future processing, I recommend to organize everything a bit better (orbits, antenna calibrations, rnx files in separate directories etc.)!
Here's a is a template of the gd2p.pl command that you should be using:
#!/bin/csh -f
#The following environment variables were active:
#setenv GOA /usr/local/GIPSY/goa-6.3
#You may need the following source command to set your PATH if you have changed it
#source
gd2p.pl \
-i RINEXFILE \
-n SSSS \
-r 0.2 \
-trop_map GPT2 \
-w_elmin 10 \
-amb_res 1 \
-type k \
-d 2015-04-25 \
-orb_clk "flinnR /data/GEOP555/USER/lab06/gorkha"\
-sta_info /data/GEOP555/USER/sta_info \
-e " -a 1 -PC -LC -F -t1 483214260.0 -t2 483214380.0"\
-pb_min_slip 1.0e-3 \
-pb_min_elev 30 \
-tides WahrK1 FreqDepLove OctTid PolTid \
-add_ocnld " -c /usr/local/GIPSY/ocnload/tpxo72atlas_CM/gipsy/SSSS.ocnld" \
-OcnldCpn \
-add_ocnldpoltid \
-dwght 1.0e-5 1.0e-3 \
-kin_sta_xyz 1.0E-3 5.7E-7 0.2 RANDOMWALK\
-post_wind 5.0e-3 5.0e-5 \
-trop_z_rw 5.0e-8 \
-eldepwght SQRTSIN \
-wetzgrad 5.0e-9 \
-arp \
-AntCal SSSS.xyz
At this point this shouldn't be much of a surprise as it's very much in line with previous labs and last week's lecture. You
can take this and create 2 files from it: run_nast and run_kkn4. If you're a skilled
shell hacker, you'd create 1 run file that takes the site id and maybe a data as parameters and works out
things that need to get done internally. I will assume you have 2 files!
Here are some things that need to be replaced:
RINEXFILE: with the nast or kkn4 rinexfiles, as they are in the same directory, you can omit the pathSSSS: 4char-site ID, capitalized!USER: your user id as it shows up in /data/GEOP555Once the replacements are done you may need to make your files executable:
$> chmod +x run_*
I should explain some additions to the template
-e " -a 1 -PC -LC -F -t1 483214260.0 -t2 483214380.0": This string is passed on to the editor (ninja, you can call
ninja -h to get the help and more options). t1, t2 indicate the start and end times in seconds since J2000 that we're
interested in. If you were to take the difference of these 2 values, you'd find that we're looking at 240 seconds worth of data. -a
determines the minimum arc length we expect for each satellite in minutes. The default, which we used in previous labs, is 20. That causes issues
when we're looking at just 2 minutes of data. Therefore, we reduce this to 1 minute.-r 0.2: We are using 5Hz data, that's 1/5 seconds which is 0.2!-kin_sta_xyz 1.0E-3 5.7E-7 0.2 RANDOMWALK: Stochastic properties for the STAX, STAY, STAZ parameters. The format is
-kin_sta_xyz apriori-sigma (km) process-noise (km/sqrt(sec)) rate (secs) RANDOMWALK|WHITENOISE, we will play with the
process noise parameter later.At last! You can start with the processing!!!
Let's get started:
$> run_nast
The errors should go to stdout now. Fix any errors that come up - it should be straightforward. Then plot the results:
$> tdp2neu NAST > nast.neu
$> plot_neu nast.neu
tdp2neu takes the XYZ solution and rotates them into a local coordinates and gives north, east, and
vertical motion for the site with respect to the nominal position (in cm). nast.neu now contains
4 columns of data: time in minutes, north, east and vertical with respect to nominal position in sta_pos.
The reason that your lines don't all start at zero, is because we used a poor initial position. As I said earlier, the
proper way would be to run a static solution first (or you can average the data before the earthquake and subtract this
from the respective column).
Now, that plot doesn't look too exiting, does it? A bunch of lines - where are the earthquake wiggles? First you'll
see that there are tremendous offsets. Those are due to poor apriori position in the sta_pos database. If we were
doing this right, I'd have you do a static run for this day, but let's just take the first position in
tdp_final, and update NAST's entry sta_pos with it.
You'll likely get lines again. This time the scale is a lot smaller though. Ideally, you'd iterate on this for a while (or run the static solution). Certainly, at this resolution though, we should be seeing some deformation due to the earthquake, shouldn't we?
From lecture you should remember a crucial lesson about the process-noise, which determines how much a site can move in between epochs.
Here, we're using a RANDOMWALK model and the process noise is (see above) 5.7e-7 * sqrt(0.2) = 0.25 mm of motion
that we allow. We need to open this up! Let's allow 2 meters of motion each epoch: (2/1000) km / sqrt(0.2 sec) = 4.5e-3.
Replace the 5.7E-7 in the -kin_sta_xyz line with this new value and rerun:
$> run_nast
Whoa! Now do the same for KKN4! (the line)
You should get 2 very different plots for this earthquake. What causes this difference in motion at each respective location? Feel free to explore in Google earth. Here are the coordinates from the rinex files and the USGS:
| LAT | LON
-----+--------+------------
KKN4 | 27.8007 85.278804
NAST | 27.6567 85.327729
EQ | 28.230 84.731
gorkha directory - don't send it! I will look at it myself, but don't alter any files after the deadline! (I'll confirm file information)rg <at> nmt <dot> edu | Last modified: August 13 2018 21:30.