IRIS_fdsnws_dataselect.sh

#!/bin/bash

wkdir=$(pwd)

# FDSN web services
#fdsnws_event="http://service.iris.edu/fdsnws/event/1/query"
fdsnws_station="http://service.iris.edu/fdsnws/station/1/query"
fdsnws_dataselect="http://service.iris.edu/fdsnws/dataselect/1/query"
irisws_sacpz="http://service.iris.edu/irisws/sacpz/1/query"

#======
usage(){
cat<<EOF
NAME

  IRIS_fdsnws_dataselect - a shell wrapper of IRIS web service to download waveform data.

SYNOPSIS

  IRIS_fdsnws_dataselect [-s <station_specs>] [-t <twin_specs>] [-o <out_dir>] <gcmt_catalog>

DESCRIPTION

  <gcmt_catalog>  a catalog list generated by fdsnws-event with catalog=GCMT

  Options:

  -s <station_specs> [b/20/60/90/150]
    box region: b/lat0/lat1/lon0/lon1
    radial region: r/lat0/lon0/radius0/radius1

  -t <twin_specs> [p,P/-500/1500]
    time window phase(s)/begin/end. If multiple phases specified,
    use the smallest traveltime as the reference time.

  -o <out_dir> [work directory]

NOTE
  1) TauP/bin/taup_time is called to calculate traveltime
  2) channel code is restricted to [FCHB]H?, i.e. only broadband high gain seismic data
  3) completeness of 3 components are checked, i.e. must have all ?H[ENZ] or ?H[12Z].

EOF
}

#====== parameters

# default option args
station_specs=b/20/60/90/150
twin_specs=p,P/-500/1500
out_dir=$wkdir

# parse options
while getopts ":s:t:o:h" opt
do
  case $opt in
  s)  
    station_specs="$OPTARG"
    ;;
  t)  
    twin_specs="$OPTARG"
    ;;
  o)
    out_dir=$(readlink -f "$OPTARG")
    ;;
  :)
    echo "Option -$OPTARG requires an argument." >&2
    usage
    exit 1
    ;;
  h)
    usage
    exit 1
    ;;
  \?)
    echo "Invalid option: -$OPTARG" >&2
    usage
    exit 1
    ;;
  esac
done

# get all mass-arguments
shift $((OPTIND-1))

event_list=${1:?[arg] need gcmt_catalog}
event_list=$(readlink -f $event_list)

# process arguments
geo_range=$(echo $station_specs | awk -F"/" '\
  $1=="b"{printf \
    "minlatitude=%s&maxlatitude=%s&minlongitude=%s&maxlongitude=%s",\
    $2,$3,$4,$5}; \
  $1=="r"{printf \
    "latitude=%s&longitude=%s&minradius=%s&maxradius=%s",$2,$3,$4,$5}')

phase_name=$(echo ${twin_specs} | cut -d / -f 1)
twin_begin=$(echo ${twin_specs} | cut -d / -f 2)
twin_end=$(echo ${twin_specs} | cut -d / -f 3)

if [ -z "$phase_name" -o -z "$twin_begin" -o -z "$twin_end" ]
then
  echo "[ERROR] time window: phase=$phase_name b=$twin_begin e=$twin_end"
  exit 1
fi

#====== get data for each event 

grep -v "^[ \t]*#" $event_list |\
while read event_line
do

  #------ event info
  IFS='|' read -ra evt <<< "$event_line"
  evid=${evt[8]#*,} # use GCMT event ID
  evodate=${evt[1]/T/ } # replace T with a whitespace
  evla=${evt[2]//[[:space:]]/} # remove whitespaces
  evlo=${evt[3]//[[:space:]]/}
  evdp=${evt[4]//[[:space:]]/}

  evoepoch=$(date -ud "$evodate UTC" +%s)
  evojday=$(date -ud "$evodate UTC" +%Y,%j,%H:%M:%S.%N)

  #------ create event directories
  event_dir=$out_dir/$evid
  mkdir -p $event_dir
  cd $event_dir
  mkdir data mseed sacpz

  #------ make log file 
  log_file=$event_dir/IRIS_fdsnws_dataselect.log
  echo "[LOG] ====== $event_line" >> $log_file
  echo "[LOG] $0 $event_list -s $station_specs -t $twin_specs -o $out_dir" >> $log_file
  echo "[LOG] wkdir=$wkdir" >> $log_file

  # save event info
  echo "$event_line" > $event_dir/data/event.txt

  #------ get channel list
  echo >> $log_file
  echo "[LOG] ------ get channel list" >> $log_file
  echo >> $log_file

  # channel operating time range (start 1 day before and end 1 day after origin time)
  startbefore=$(date -ud "$evodate -1 days" +%Y-%m-%dT%H:%M:%S)
  endafter=$(date -ud "$evodate 1 days" +%Y-%m-%dT%H:%M:%S)
  str_twin="startbefore=${startbefore}&endafter=${endafter}"

  # retrive channel info from IRIS web service
  channel_list=$event_dir/data/channel.txt
  echo >> $log_file
  echo "[LOG] get channel list" >> $log_file
  str_link="${fdsnws_station}?${str_twin}&${geo_range}&channel=FH?,CH?,HH?,BH?&level=channel&format=text"
  echo >> $log_file
  echo "[LOG] wget $str_link -O $channel_list" >> $log_file
  wget "$str_link" -O $channel_list -a $log_file

  #------ check channel list
  echo >> $log_file
  echo "[LOG] ------ check channel list" >> $log_file
  echo >> $log_file

  #--- remove records with incomplete components,
  # e.g. not all ?H[ENZ] or ?H[12Z] exist.
  echo >> $log_file
  echo "[LOG] --- check completeness of components" >> $log_file
  echo >> $log_file

  awk -F"|" '$1!~/#/{c=substr($4,1,2); printf "%s|%s|%s|%s\n",$1,$2,$3,c}' \
    $channel_list | sort -u |\
  while read seed_id
  do
    comp=$(grep "^${seed_id}" $channel_list |\
      awk -F"|" '{print substr($4,3,1)}' | sort | awk '{printf "%s",$1}')
    ncomp=$(echo -n "$comp" | wc -c)

    # discard stations which have less than 3 components
    if [ ${ncomp} -lt 3 ]
    then
      echo >> $log_file
      echo "[WARN] incorrect components found for ${seed_id}" >> $log_file
      grep "^${seed_id}" $channel_list >> $log_file
      sed -i "/^${seed_id}/s/^/#/" $channel_list
      continue
    fi
    
    # keep either ENZ or 12Z when only 3 components found
    if [ ${ncomp} -eq 3 ] && [ "$comp" != "ENZ" ] && [ "$comp" != "12Z" ]
    then
      echo >> $log_file
      echo "[WARN] incorrect components found for ${seed_id}" >> $log_file
      grep "^${seed_id}" $channel_list >> $log_file
      sed -i "/^${seed_id}/s/^/#/" $channel_list
      continue
    fi

    # when more than 3 components found, keep ENZ or 12Z
    if [ ${ncomp} -gt 3 ]
    then
      echo ${comp} | grep 'Z' > /dev/null
      if [ "$?" -eq 0 ]
      then
        echo ${comp} | grep 'EN'  > /dev/null
        if [ "$?" -eq 0 ]
        then
          sed -i "/^${seed_id}[12]|/s/^/#/" $channel_list
          echo >> $log_file
          echo "[WARN] keep ${seed_id}[ENZ] only" >> $log_file
          continue
        else
          echo ${comp} | grep '12' > /dev/null
          if [ "$?" -eq 0 ]
          then
            sed -i "/^${seed_id}[EN]|/s/^/#/" $channel_list
            echo >> $log_file
            echo "[WARN] keep ${seed_id}[12Z] only" >> $log_file
            continue
          else
            echo >> $log_file
            echo "[WARN] incomplete components found for ${seed_id}" >> $log_file
            grep "^${seed_id}" $channel_list >> $log_file
            sed -i "/^${seed_id}/s/^/#/" $channel_list
            continue
          fi
        fi
      else # Z does not exits
        echo >> $log_file
        echo "[WARN] incomplete components found for ${seed_id}" >> $log_file
        grep "^${seed_id}" $channel_list >> $log_file
        sed -i "/^${seed_id}/s/^/#/" $channel_list
        continue
      fi
    fi

  done

  #--- remove records with multiple sampling rates (e.g. HHZ and BHZ)
  # keep the first one found in the order of B H C F
  echo >> $log_file
  echo "[LOG] --- check multiple sampling rates" >> $log_file
  echo >> $log_file

  awk -F"|" '$1!~/#/{c=substr($4,2,2); printf "%s|%s|%s|.%s|\n",$1,$2,$3,c}' \
    $channel_list | sort -u |\
  while read seed_id
  do
    read -a band <<< $(grep "^${seed_id}" $channel_list |\
      awk -F"|" '{print substr($4,1,1)}' | sort)
    if [ "${#band[@]}" -gt 1 ]
    then
      echo >> $log_file
      echo "[WARN] multiple band codes found for $seed_id" >> $log_file
      grep "^${seed_id}" $channel_list >> $log_file
      sed -i "/^${seed_id}/s/^/#/" $channel_list
      for code in B H C F
      do
        echo "${band[@]}" | grep "${code}" > /dev/null
        if [ "$?" -eq 0 ]
        then
          seed_id=$(echo "$seed_id" | awk -F"|" -v c="$code" \
            '{x=substr($4,2,2); printf "%s|%s|%s|%s%s|",$1,$2,$3,c,x}')
          echo "[WARN] only keep $seed_id" >> $log_file
          sed -i "/^#${seed_id}/s/^#//" $channel_list 
          break
        fi
      done
    fi
  done

  #--- remove multiple location codes 
  # keep the first one from the sorted list of location codes
  echo >> $log_file
  echo "[LOG] --- check multiple location codes" >> $log_file
  echo >> $log_file

  awk -F"|" '$1!~/#/{printf "%s|%s|[^|]*|%s|\n",$1,$2,$4}' \
    $channel_list | sort -u |\
  while read seed_id
  do
    read -a loc <<< $(grep "^${seed_id}" $channel_list |\
      awk -F"|" '{if($3=="") $3="--"; print $3}' | sort -u)
    if [ "${#loc[@]}" -gt 1 ]
    then
      echo >> $log_file
      echo "[WARN] multiple location codes found for $seed_id" >> $log_file
      grep "^${seed_id}" $channel_list >> $log_file
      sed -i "/^${seed_id}/s/^/#/" $channel_list
      code=${loc[0]}
      if [ "$code" == "--" ]
      then
        code=""
      fi
      seed_id=$(echo "$seed_id" | awk -F"|" -v c="$code" \
        '{printf "%s|%s|%s|%s|",$1,$2,c,$5}')
      sed -i "/^#${seed_id}/s/^#//" $channel_list 
      echo "[WARN] only keep ${seed_id}" >> $log_file
    fi
  done

  #------ generate data selection list
  echo >> $log_file
  echo "[LOG] ------ make data selection list" >> $log_file
  echo >> $log_file

  dataselect_list=$event_dir/data/dataselect.txt
  cat /dev/null > $dataselect_list

  awk -F"|" '$1!~/#/&&$4~/.HZ/{if($3=="") $3="--"; \
    printf "%s %s %s %s %s %s\n",$1,$2,$3,$4,$5,$6}' $channel_list |\
  while read net sta loc cha stla stlo
  do
    # get traveltime
    ttp=$(taup_time -evt $evla $evlo -h $evdp \
        -sta $stla $stlo --time -ph $phase_name | head -n1 | awk '{print $1}')
    if [ -z "$ttp" ]
    then
      echo "[WARN] taup_time failed to get $phase_name travaltime for $net $sta $loc $cha" >> $log_file
      echo "taup_time -evt $evla $evlo -h $evdp -sta $stla $stlo --time -ph $phase_name" >> $log_file
      continue #exit 1
    fi
    # get time window
    t0=$(echo "$evoepoch + $twin_begin + $ttp" | bc -l)
    t1=$(echo "$evoepoch + $twin_end + $ttp" | bc -l)
    starttime=$(date -u -d @$t0 +%Y-%m-%dT%H:%M:%S)
    endtime=$(date -u -d @$t1 +%Y-%m-%dT%H:%M:%S)
    # print data selection list
    printf "%-6s %-6s  %-6s  %s  %s  %s\n" \
      $net $sta $loc "${cha:0:2}?" $starttime $endtime \
      >> $dataselect_list
  done

  #------ get waveform data
  echo >> $log_file
  echo "[LOG] ------ get waveform data" >> $log_file
  echo >> $log_file

  mseed_file=$event_dir/mseed/${evid}.mseed
  # query waveforms 
  str_link="${fdsnws_dataselect}"
  echo >> $log_file
  echo "[LOG] wget --post-file=$dataselect_list $str_link -O $mseed_file" >> $log_file
  echo >> $log_file

  wget --post-file=$dataselect_list "$str_link" -O $mseed_file -a $log_file

  #------ get sacpz
  echo >> $log_file
  echo "[LOG] ------ get sacpz" >> $log_file
  echo >> $log_file

  sacpz_dir=$event_dir/sacpz

  awk -F"|" '$1!~/#/{if($3=="") $3="--"; \
    printf "%s %s %s %s\n",$1,$2,$3,$4}' $channel_list |\
  while read net sta loc cha
  do
    read -ra x <<< $(grep "${net}[ ]*${sta}[ ]*${loc}" $dataselect_list)
    if [ "$?" -ne 0 ]
    then
      echo >> $log_file
      echo "[WARN] ${net}.${sta} is not found in $dataselect_list! SKIP" >> $log_file
      continue
    fi
    starttime="${x[4]}"
    endtime="${x[5]}"

    echo >> $log_file
    echo "SACPZ ${net}.${sta}.${loc}.${cha} $starttime $endtime" >> $log_file
  
    url="${irisws_sacpz}?net=${net}&sta=${sta}&loc=${loc}&cha=${cha}&starttime=${starttime}&endtime=${endtime}"

    if [ "$loc" == "--" ]
    then
      loc=""
    fi
    outfile="${sacpz_dir}/${net}.${sta}.${loc}.${cha}"
  
    wget "$url" -O $outfile -a $log_file
  
  done

done # for event_line in $(grep -v "^[ \t]*#" $event_list)

#END