指定した一時間の中央値を計算

python

200523_N.MWDH_U_bp10-20_sqr.sの15時〜１６時の中央値を求める場合
$ python sac_median.py 200523_N.MWDH_U_bp10-20_sqr.s 15

sac_median.py

import obspy
from obspy import read
import numpy as np
import sys

#コマンドラインから引数
file = sys.argv[1]
hour = sys.argv[2]
#数値に変換
hour = int(hour)

#サンプリング周波数
samplerate = 100

#sacファイル読み込み
sac1 = read(file, debug_headers=True)#読み込みたいファイルをここに入力
#振幅データ（縦軸）
y = sac1[0].data


##時間（横軸）．
#time = np.arange(0, len(y))/samplerate

#指定した時刻〜一時間後の時刻までの中央値を計算
median = np.median(y[hour*360000:(hour+1)*360000])
median = median*1e-9 #m/sからnm/sへ変換
print(median)

fortran

200523_N.MWDH_U_bp10-20_sqr.sの15時〜１６時の中央値を求める場合
$ gfortran -o median sac_median.f95 #コンパイル
$ ./median 200523_N.MWDH_U_bp10-20_sqr.s 15 #実行

sac_median.f95

program main
 implicit none
 integer i, a(40), max, npts, hour
 real b(70),c(2), dt, med
 real, allocatable :: x(:), y(:)
 character :: file1*80, arg*10, khd(48)*4

  if(iargc()>0) then
    call getarg (1, file1) !コマンドラインで引数!
    call getarg (2, arg) !コマンドラインで引数所得
    read(arg,*) hour
 endif

 max=10000000
 allocate (x(max),y(max))
 call arsac(file1,b,a,khd,x,max)
  npts = a(10)
  dt = b(1)
  do i = hour*360000+1,(hour+1)*360000+1
     y(i-hour*360000) = x(i)
  enddo
  call median(y, med, 360000)
  med = med*1e-9 !m/sからnm/sへ変換
 write(*,*)med
end program main
!########################################################
subroutine arsac(fn,fhd,nhd,khd,y,maxy)
  real fhd(70),y(maxy)
  integer nhd(40),i,i1,i2,i3,npts
  character khd(48)*4,fn*80

  open(1,file=fn,access='direct',recl=632,status='old')
  read(1,rec=1)(fhd(i1),i1=1,70),(nhd(i2),i2=1,40),(khd(i3),i3=1,48)
  npts=nhd(10)
  close(1)
  open(1,file=fn,access='direct',recl=4,status='old')
  do i=1,npts
!     read(1,rec=i+158,end=100)y(i)     Mar.1 H.A
     read(1,rec=i+158)y(i)
  enddo
! 100 continue
  close(1)
  return
end subroutine arsac
!#####################################################
subroutine median(a, med, n)
  implicit none
  integer,intent(in):: n
  real,intent(inout):: a(n)
  real,intent(inout):: med
  integer i, m
  real am

!バブルソート（遅い）
!  do i = 1, n-1 !sort a
!     am = minval(a(i+1:n))
!     m  = minloc(a(i+1:n), 1) + i
!     if (a(i) > am) then
!        a(m) = a(i)
!        a(i) = am
!     endif
!  enddo

  call quicksort(a,1,360001)

  if(mod(n, 2) /= 0) then !median of a(:)
     med = a((n-1)/2+1)
  else
     med = 0.5*(a(n/2)+a(n/2+1))
  endif
end subroutine median

!##############################################
recursive subroutine quicksort(a, first, last)
  implicit none
  real  a(*)
  integer first, last
  integer i, j
  real  x, t

  x = a( (first+last) / 2 )
  i = first
  j = last
  do
     do while (a(i) < x)
        i=i+1
     end do
     do while (x < a(j))
        j=j-1
     end do
     if (i >= j) exit
     t = a(i);  a(i) = a(j);  a(j) = t
     i=i+1
     j=j-1
  end do
  if (first < i - 1) call quicksort(a, first, i - 1)
  if (j + 1 < last)  call quicksort(a, j + 1, last)
end subroutine quicksort