class 24 pptx

Earth Science Applications of Space Based Geodesy
DES-7355
Tu-Th 9:40-11:05
Seminar Room in 3892 Central Ave. (Long building)
Bob Smalley
Office: 3892 Central Ave, Room 103
678-4929
Office Hours – Wed 14:00-16:00 or if I’m in my office.
http://www.ceri.memphis.edu/people/smalley/ESCI7355/ESCI_7355_Applications_of_Space_Based_Geodesy.html
Class 24
1
First --- when you are going to be entering the same
command many times, even if it is a one liner:
Put it in a shell script
AUTOMATE, AUTOMATE, AUTOMATE
Don’t depend on the history feature and continuously be
searching for the same command over and over and over.
2
Try to give your shell script a name whose first few
characters are not the same as a system command or
another command you are running.
You can then use the “!x…” construct to rerun it while in an
edit, run, look at it loop.
eg.
>
>
>
>
>
>
vi my_gmt_script.sh
my_gmt_script.sh
gs my_gmt_script.sh.ps
!vi
!my
!gs
while developing a GMT script
3
If you are going to do a test run and then do the same
command for a larger data set – use variables.
Also - put in some comments.
#!/bin/sh
YR=2009
START=015
FINIS=365
YREXT=-yrext
#daily position estimation of your data plus some stations from global nets
sh_gamit -s $YR $START $FINIS -expt same -orbit IGSF -copt x k p -dopt c ao -remakex Y $YREXT > sh_gamit.log
4
Things that trip you up.
Metadata:
You need correct station information in
station.info (metadata about receiver, antenna, firmware,
offsets to reference point, etc.)
and the
a-priori file where you need a good starting location
(velocity can be zero)
5
You can get the good starting location from one of the rinex
processing/location services on the web page.
6
Getting the metadata for station info can be challenging.
First see if someone you know has it in a station.info file
already and ask (beg?) them for it.
If not you will have to construct it yourself by trying to get
copies of the station logs (usually available in various
formats from the GPS archive data servers – but accuracy
and currency sometimes questionable).
7
If you have a bad a priori estimate (or a rinex file has bad
observation data) you will probably get a bad rms in the
qfile.
----If you are missing metadata in station.info - it may just run,
ignoring the station that does not have the metadata – so
the only error will be no results for that station.
It may try to construct the metadata from the rinex header –
things work.
8
Now is also a good time to select and check your reference
stations.
Make an initial selection of at least 6 well distributed IGS
core sites (select from the map on the igs core sites web
page) within a few thousand km of your area of interest.
braz
chpi
conz
harb
ispa
lpgs
mbar
nklg
sant
unsa
areq
cfag
Copo
cslo
kour
lhcl
9
Getting tabular metadata on
data availability from sopac
10
Things looked “ok” for the sites we selected, so we
continue on -Process with sh_gamit and then sh_glred.
Look at repeatabilities.
11
Output of gl_red: problem – 2 regions with “jumps” common
to all sites, plus several small sets of jumps
This is most likely due to a
“bad” reference station.
“bad” in this case does not
necessarily mean the station
has bad data (you would catch
that in the gamit processing –
probably has big rms residual),
but in this case the “bad”
reference station probably has
two periods of missing data.
12
When you have a small number
of reference sites this is what
happens when one or more
reference sites are missing – the
whole network jumps.
If we used 100 reference sites
(it would have taken several
weeks for gamit to run), the
network would still jump when
there is missing data, but the
jump would be very small.
13
We have to fix this before we estimate velocities
14
Quick check on data (after sh_gamit processing).
capybara:rinex bob$ files=`ls *1500.09o | nawk '{print substr($1,1,4)}'`
capybara:rinex bob$ for file in $files; do echo $file; ls *$file* | wc; done
areq 357 357 4641
braz 338 338 4394
cfag 346 346 4498
chpi 363 363 4719
conz 365 365 4745
copo 361 361 4693
cslo 367 367 4765
harb 358 358 4654
ispa 351 351 4563
kour 361 361 4693
lhcl 340 340 4420
lpgs 363 363 4719
mbar 365 365 4745
nklg 328 328 4264
sant 365 365 4745
unsa 365 365 4745
A couple of sites are “perfect” (365 days, but don’t know if
100% data for each day).
Worst site, nklg, is missing total of 37 days.
15
Nothing particularly obvious here (no single bad site from
~125 to ~160.
16
Look for ref sites with missing data
17
Go to the database where you will obtain this data and try
to find the time series.
Lots of missing data, but not correlated with periods of
jumps (and we have data for 338 days – lots more than
shown here. This plot is from scripps, but we are using
orbits and h-files from mit).
18
Looks pretty good.
19
Looks pretty good.
20
Looks pretty good.
21
Gaps match up with problem areas on gl_red results.
Potential problem station identified.
22
Looks pretty good.
23
Looks pretty good.
24
Looks pretty good.
25
Looks pretty good.
26
So there seems to be a lone culprit -- IPSA.
27
Remove IPSA from the list of reference stations in
sites.defaults (just remove the “glreps” code from the line
for station ipsa) and remove it from the list of stabilization
sites in your stab site file (look in glorg_comb.cmd for a line
that looks like
source ../tables/stab_site.global
To identify to stabilization site file. Here is what is in my file
(after removing ispa)
*Global stabilization list for ITRF05
stab_site clear
stab_site kour braz chpi unsa lpgs conz nklg mbar harb
You will need to place a file like this in your tables directory
– list your stabilization sites.
28
Re-run sh_glred.
(this is fast, no need to re-run sh_gamit, which is slow)
29
Look at outputs.
Problem fixed.
Still a few single day
problems – ignore for
now.
30
The rest of the command file.
#!/bin/sh
YR=2009
START=015
FINIS=365
YREXT=-yrext
#daily position estimation of your data plus some stations from global nets
sh_gamit -s $YR $START $FINIS -expt same -orbit IGSF -copt x k p -dopt c ao -remakex Y $YREXT > sh_gamit.log
#now combine your daily results with those of processing center
#first clean up a little so you don’t use old cmd files, etc.
cd gsoln
\rm globk_comb.cmd
\rm glorg_comb.cmd
\rm *.org
\rm *.prt
\rm *.log
cd ..
sh_glred -expt same -s $YR $START $YR $FINIS -net MIT -opt H G E F $YREXT > sh_glred.log
#make list stations to process for velocities
cd gsoln
create_gdl.sh
##make *vel.cmd, make uselist
clean up
\rm same.prt
\rm same.log
\rm globk_vel.org
#estimate velocities
globk 6 same.prt same.log new.gdl globk_vel.cmd
sh_plotvel -f globk_vel.org
31
How to make the velocity command data file.
This is dependent on your processing tree and file naming.
#!/bin/csh -f
set a=1
ls -1 ../glbf/h$1*_same.glx >! newtemp1
set nn=`cat newtemp1 | wc -l`
#../glbf/h0901011200_same.glx
#000000000111111111122222222223
#123456789012345678901234567890
while($a <= $nn)
set mm=`head -n$a newtemp1 | tail -1 | cut -c12-13`
set dd=`head -n$a newtemp1 | tail -1 | cut -c14-15`
set yr=`head -n$a newtemp1 | tail -1 | cut -c10-11`
set jd1=`doy $yr $mm $dd | grep DOY | cut -c32-35`
set jdd=`echo $jd1 | awk '{printf "%03s",$1}'`
echo "../glbf/h"$yr$mm$dd"1200_same.glx +" >> new.gdl
echo "../glbf/H"$yr$jdd"_MIT.GLX 1.0 " >> new.gdl
set a = `expr $a + 1`
end
rm newtemp1
32