New Mexico Tech
Earth and Environmental Science

ERTH 455 / GEOP 555 - Geodetic Methods

Lab 6: gd2p.pl: kinematic position estimation

"I like my crust deformed."
UNAVCO bumper sticker

Note that you DO have to work on redoubt today!

Introduction

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.

Setup

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!

Ghorka Earthquake - SETUP

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.

Ghorka Earthquake - PROCESSING

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:

FIND ALL OCCURRENCES OF THESE PLACEHOLDERS AND REPLACE THEM WITH ACTUAL VALUES!

Once the replacements are done you may need to make your files executable:

    $> chmod +x run_*
    

I should explain some additions to the template

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 
	

Deliverables: (submit via canvas!)

rg <at> nmt <dot> edu | Last modified: August 13 2018 21:30.