Daily residuals, shocks, and cleaned low-rate data

Usage at the end of this page

Having to reproduce a series of  G1_garb_*-1s.mc  and  g*-OPNEND-1h.mc  files?

With the advent of the daily-resid procedure, this section of the documentation has become quite central

First, we create 1-second MC-files for every day, containing the following components:
GRAV|VAL
BARO|VAL
GRAV|PRED
GRAV|RES
from  RAW_o054/G1$date.054
to   d/G1_gar_$date-1s.mc
or   d/G1_garb_$date-1s.mc  (reduced for barometric effect)  
or   d/G1_garp_$date-1s.mc  (reduced for polar motion)

with the script
 daily-resid -pm -T garp RAW_o054/G1*.054
alt.
 daily-resid RAW_o054/G11011??.054
Note that -pm  requests to reduce polar motion, and  -T garp marks the file names accordingly.
This takes an awful lot of time (i.e. computing with run_urtip)

With the option to generate a barometric residual,
and with the recently added Tamura potential, the
following was obtained:

daily-resid -B F -T garb RAW_o054/G1100515.054
tsd -n 86400,1,c -L 'G|R' d/G1_garb_100515-1s.mc
tsd -n 86400,1,c -L 'G|B' d/G1_garb_100515-1s.mc





To subtract the static effect of pressure 1015 hPa, use
daily-resid -BS F -T garb RAW_o054/G1100515.054

 At this point we can make a graphic view of the micro-seismic spectrum (a time series of 2048s duration segments)

   useisp [-f G1_garp_] $year $m $d  (e.g. useisp -f G1_garp_ 2011 03 11 )
(The default is G1_gar_ )

Then we detect shocks,  -a -nh  appending the existing table without writing a header
(to improve this, the data is low-pass filtered in order to avoid too many losses due to microseismic noise;
the procedure adds the next day and truncates to 86400 smp in order not to lose data round midnight
change made 2010-12-15; the level 400 is way too high! New level 0.2 ~compatible with old 400.

  detect-shocks 0.2 -a -nh -m 600,1200 -f shocks-400.tse -L `label GRAV,RES` d/G1_gar_*-1s.mc
  detect-shocks 0.1 -a -nh -m 600,1200 -f shocks-200.tse -L `label GRAV,RES` d/G1_gar_*-1s.mc

(Test this before you really append the shock database - output to shocks-tst.tse, run without the -a -wh options.)
Check if the weeding-out makes sense using

  tsd -n 86400,1,c -L GRAV_________RES d/G1_gar_101130-1s.mc

The big shock table can be sorted by shocksize:
  awk '{printf "%10.1f %s\n",$20,$0}' shocks-200.tse | sort -n -k1 -r > shocks-200.sorted
so we get an overview.

We select shocks above a certain threshold for use in the resampling process:
  awk -v t=8 '{if ($20 > t){print $0}}' shocks-200.tse > shocks4use.tse
so we can live with an initially low detection level.

Then we can apply the following TSF EDIT sequences to interpolate linearly through shocked samples:
(editshocks.tse)
TSF EDIT D600S
NOECHO
CONT 37 B shocks4use.tse
ECHO
REPAIR L ALL
FILTER D:WD KAIS 60 1.5 0.0 0.03 0.5 0.0
DECIMATE 60 $MOVE
FILTER D:WD KAIS 10 1.5 0.0 0.18 0.5 0.0
DECIMATE 10 0
CLOSE 37
NOECHO
CONT 37 B shocks4use.tse
ECHO
REPAIR L ALL
END

TSF EDIT D600O
FILTER D:WD KAIS 60 1.5 0.0 0.03 0.5 0.0
FILTHR  24 12 6 1 0.5 0.01 0.0083333333 0.004
DECIMATE 60 $MOVE
FILTER D:WD KAIS 10 1.5 0.0 0.18 0.5 0.0
FILTHR  24 12 6 1 0.5 0.1 0.083333333 0.04
DECIMATE 10 0
END

Short-duration shocks might disturb the micro-seismic spectrum, so you could try
   useisp -T "-Euseshocks.tse,U -Erepair.tse,R" $year $m $d

Now let's get a long 10min time series: ( the purpose of  make-garb-mc )
  cp d/gar-600s.tsf d/gar-600s.tsf-bup
  setenv MOVE 210
  foreach file ( `all-but-the-last d/G1_gar_*-1s.mc` )
  app2days -wa d/gar-600s.tsf -E editshocks.tse,D600S -E editshocks.tse,D600O -E editshocks.tse,D600S \
           -T -Un144 $file 'GRAV|VAL' 'BARO|VAL' 'GRAV|BRES'
  end

Above:
If you want to limit the file list, use e.g.
foreach file ( `mfn -L -u 120101 d/G1_gar_ 090701 -1s.mc` )
mfn is in ~/TD

 Or simply ignore that the last operation will fail as the day beyond the last on the list is requested.

You may also use of the following to generate file names for a specific period:
   rm -f d/ap100916-600s.tsf
   set mjds = (  `do_mjds 2010 09 16 2010 12 18` )
   foreach file ( `mjd2yymmdd -p d/G1_gar_ -a -1s.mc $mjds` )
   app2days -wa d/ap100916-600s.tsf -E editshocks.tse,D600O -T -Un144 $file 'BARO|VAL'
   end
Example is to obtain 600s air pressure

Example to obtain multi-column most of the extent:
   rm -f d/garb-600s.tsf
   touch d/garb-600s.tsf
   setenv MOVE 540
   set mjds = (  `do_mjds 2009 07 01 2011 10 30` )
   foreach file ( `mjd2yymmdd -p d/G1_garb_ -a -1s.mc $mjds` )
      app2days -wa d/garb-600s.tsf -E editshocks.tse,D600S -E editshocks.tse,D600O \
                                   -E editshocks.tse,D600S -E editshocks.tse,D600S   \
               -T -Un144 $file 'GRAV|VAL' 'BARO|VAL' 'GRAV|RES' 'GRAV|BRES'
   end

Make mc files (you may use -Erepair.tse,R ):
This is the "manual" case _gar_:
  rm -f d/gar-600s.mc
  tslist d/gar-600s.tsf -gi4,2i3,i4,2i3,t37,e15.0 -k3 -I -O:`label GRAV,VAL` d/gar-600s.mc
  tslist d/gar-600s.tsf -gi4,2i3,i4,2i3,t52,e15.0 -k3 -I -O:`label BARO,VAL` d/gar-600s.mc
  tslist d/gar-600s.tsf -gi4,2i3,i4,2i3,t67,e15.0 -k3 -I -O:`label GRAV,RES` d/gar-600s.mc 

Advantage: The output file can be expanded.

The case _garb_ is managed with the script make-garb-mc, including all of the above, resulting in a 10-min sampled mc-file d/garb-600s.mc

At this point we can generate (and subtract) a polar-motion series:
Check the starting date:
    tslq d/garb-600s.mc
    set bdate = `tslq d/garb-600s.mc | awk '/File begin/{print substr($4,3,2)$5$6}'`
    polmotm -l 11.9,57.6 -d 1.16 d/garb-600s.mc  |\
    tslist - -gI4,5i3,t66,f11.0 -k3 -C3 -I -o PM/PMG${bdate}-600s.ts
    tslq PM/PMG${bdate}-600s.ts
    tslist d/garb-600s.mc -I -O:`label GPM,PRD` d/garb-600s.mc
(Carefull!!! - rather make a backup before. Check for compatibility of dates and lengths!)

    setenv FNSUB PM/PMG${bdate}-600s.ts
    tslist d/garb-600s.mc -L'G|BRES' -I -Esubts.tse,SUBU -Eeditshocks.tse,DRBIGSH -O:`label GRAV,BPMR` d/garb-600s.mc
 (Careful!!! - the PM-series is not least-squares adjusted)
The section DRBIGSH in editshocks.tse removes the drift and deletes a few repercussions around the coldhead failure in Feb. 2011..
The following is more prudent:
    tslist d/garb-600s.mc -L'G|BRES' -I -Esubts.tse,SUBU -Eeditshocks.tse,DRBIGSH -D -o d/garbp-600s.ts

How do we find the interpolated records? How can we count them?

   tslist d/gar-600s.mc -I -Eeditshocks.tse,COUNT

with
TSF EDIT COUNT
NOECHO
CONT 37 B shocks4use.tse
ECHO
REPAIR L ALL
END
This is a harmless wrap to edit out the  shocks4use.tse (although we do repair).
A benign wrap is  useshocks.tse
TSF EDIT USE
NOECHO
CONT 37 B shocks4use.tse
ECHO
END

Drift determination on the basis of d/gar-600s.mc

  cp newexp.tse newexp.tse.backup
  ln -s useshocks.tse newshocks.tse   # if not yet done
  polmotm -l 11.9,57.6 -d 1.16 d/gar-600s.mc |\
     tslist - -gI4,5i3,t66,f11.0 -k3 -C3 -I -o PM/PMG090615-600s.ts
  setenv FNSUB PM/PMG090615-600s.ts
  tslist d/gar-600s.mc -L'G|R' -Esubts.tse,SUBU -C3 -Ff10.3 -I -D -o d/grp-600s.ts

Augment the mc-file:
  tslist d/grp-600s.ts -I -D -O:`label GRAV,PMRES` d/gar-600s.mc
  tslist d/gar-600s.mc -L'G|PMR' -ML-3.1559,1,s:'B|V]Rwd' -I -D -o d/garp-600s.ts  

Eventually augment the mc-file:
  tslist d/garp-600s.ts -I -D -O:`label GRAV,BPMRES` d/gar-600s.mc

Now, d/garp-600.ts is the polar-motion and baro reduced residual. This is a nice time series to plot.
Use the old drift parameters:
  tsd -n 82000,1,c -E shockdrift.tse SD d/garp-600s.ts

Determine new:
  tslist d/garp-600s.ts -Euseshocks.tse,U -I -o d/garps-600s.ts

  setenv EXPFIT_INPUT d/garps-600s.ts
  setenv EXPFIT_OUT 090615-600s
  setenv EXPFIT_TSE expfit-try.tse
  expfitbm \!expfitm.ins '>BISECT>'

Reproducing G1_garb_*.mc  and g*-OPNEND-1h.mc
cd ~/TD
set days = ( `jdc -D -L0,21 -m -A1 -fs 2018 04 18` )
foreach d ( $days )
   daily-resid -B F -T garb RAW_o054/G1$d.054
end
rm -f d/tmp/g090615-OPNEND-1h.mc
/bin/cp -f d/g090615-OPNEND-1h.mc d/safe/
tslist d/g090615-OPNEND-1h.mc -LG -LB -I -U2018,4,17,5 -O1:=+2:= d/tmp/g090615-OPNEND-1h.mc
mv d/tmp/g090615-OPNEND-1h.mc d/g090615-OPNEND-1h.mc
In a new xterm, bash:
/home/hgs/TD/ts4openend GBR 


USAGE:    
           daily-resid [options] file [file ...]
 
PURPOSE:
           create an MC-file from e.g.
           an A1-, A2 or G1-ascii file
           or a xxYYMMDD*.mc file
           The output file will be named
           d/${type}_${theme}_YYMMDD-1s.mc
           by default with the following components:
           `label GRAV,VAR`   `label BARO,VAR`
           `label GRAV,PRED`  `label GRAV,RES`
           This can be varied, but GRAV|VAR is compulsory,
           see options -x and -L
           $type is copied from the input file (G1, A1, A2),
           $theme is set by the -T option, else "gar".
 
           Files: ASCII: may be a bz2-compressed file. Will be
           left behind uncompressed.
           Files: MC: default labels Grav-1 Baro-1
 
           Usually you would simply issue
             daily-resid RAW_o054/G1*.054
           and go and have your supper.
           Consider also
             daily-resid -B F -T garb RAW_o054/G111*.054
           to include atmospheric pressure reduction on the fly.
 
             daily-resid -x 5,11,19 -L `label GRAV,VAL` -L `label XTILT,PWR`
                                    -L `label YTILT,PWR` RAW_o054/A2*.054
 
           Once you have GRAV|RES you can subtract the drift:
             tslist d/G1_${theme}_$date-1s.mc -L'G|R' -Eexp.tse,EXPM ...
 
OPTIONS:
 
 -C  urtip      - The program to use for tide prediction
                  urtipg or urtipgt                                [urtipgt]
 
 -V             - Don't exec, show the commands.
 -k             - Keep the temporary ascii file.
 -N             - Don't overwrite existing MC output files.
 -pm            - (also: -PM) add polar motion to predicted tide
 -NT            - don't add columns for tide and residual
 
 -T theme       - theme part of file name                              [gar]
 
 -c cal         - calibration factor                              [-777.421]
 
 -B { #c | F | file }
                - barometric correction, creates an additional
                  (a fifth) channel. Coefficient will be taken
                  from FILE                                    [t/urtap.prl]
                  A numerical coefficient #c can be
                  specified explicitly
 
 -x c1,c2       - extract columns GRAV from c1
                              and BARO from c2                         [1,3]
 
 -x ca[,cb...]  - more columns can be specified
                  but then, all labels must be
                  defined, and 'GRAV,VAL' is compulsory

 -O odir        - Output directory                                       [d]
 
 -L label       - output labels,                         [GRAV,VAL BARO,VAL]
                  to be repeated such that
                  all columns get a label.
                  Else a default label                            [COL#n|ft]
                  will be injected, where
                  ft is the two first letters of
                  the first file's name (e.g. G1)
                  and #n is taken from the -x option
 
 -LI label      - input labels in the case of binary MC input.
                  to be repeated such that they match the
                  output labels                              [Grav-1 Baro-1]
 
ENVIRONMENT:
 
  RESIDEDIT      - normally empty. setenv RESIDEDIT -TLS
                  when processing days near leap-seconds
 
LOG FILES:
            Since this routine is most often called for a sequence
            of day files, logging will just keep the last turn of the loop
 
            run_urtip      logs to  /home/hgs/tslist-logs/daily-resid.pt.ts.log
                           where you can check that polar motion is o.k.
 
            Final assembly logs to  /home/hgs/tslist-logs/daily-resid-1s.mc.log