카테고리 보관물: Programming

책: 실용주의 프로그래머 – 앤드류 헌트, 데이비드 토머스

‘실용주의 프로그래머’는 전문 프로그래머들을 위한 조언들이 담긴 책입니다. 제가 대학원생 때 읽고 책의 조언들 중 몇 가지를 직접 실천해보면서 많은 도움을 얻었습니다. 총 70가지의 조언들이 담겨 있는데 여기에서 내용을 확인할 수 있습니다. 제가 실천했던 조언들 중 몇 가지만 살펴보겠습니다.

  • 11. DRY (Don’t Repeat Yourself) – 반복하지 마라.
  • 12. 재사용하기 쉽게 만들라.
  • 21. 명령어 셸의 힘을 사용하라.

연구하면서 반복적으로 하는 작업들이 있습니다. 대학원생 때 반복적으로 하던 작업들 중 이진 파일로 된 탄성파 자료나 속도모델 등에서 trace나 profile을 추출하는 작업, 지하 영상에 미분이나 라플라시안 필터를 적용하는 작업, 깊이에 따라 속도가 선형으로 증가하는 속도모델을 만드는 작업 등이 있었습니다. 그때 그때 코드를 새로 작성하거나 예전에 사용했던 코드를 수정해서 작업했었는데 이러한 반복 작업을 없애기 위해 명령어 셸에서 옵션만 바꿔가며 쉽게 재사용할 수 있는 프로그램들을 작성했었습니다. 그리고 그런 프로그램들을 모아 gpl 라이브러리라고 이름을 붙였죠. 주로 당시 사용하던 포트란 프로그래밍 언어로 프로그램들을 작성했는데, 명령어 셸에서 사용할 수 있도록 option parser도 만들고 파일 입출력도 많이 다루면서 프로그래밍 연습을 할 수 있었습니다.

  • 22. 하나의 에디터를 잘 사용하라.

리눅스 Command line 상에서 프로그래밍을 했기 때문에 당시 사용할 수 있는 에디터로 VimEmacs가 있었습니다. 둘 중 어느 것을 사용할까 비교해보다가 결국 어디에나 설치되어 있는 Vi(Vim)를 사용하기로 했고, 각종 명령어와 plugin들을 이용해 Vim 고급 사용법을 익혔습니다.

  • 23. 언제나 소스코드 관리 시스템을 사용하라.

당시에는 Mercurial을 선택해서 사용하다가 몇 년 전부터는 Git만 사용하고 있습니다.

  • 28. 텍스트 처리 언어를 하나 익혀라.

당시 포트란과 C 언어만 사용하고 있었는데, 이 책을 읽은 후 Ruby, Perl, Python 중 프로그래머를 행복하게 하자는 철학을 내세웠던 Ruby를 선택해서 공부하고 사용했었습니다. Ruby를 배워놓았던 것이 도움이 돼서 나중에 Ruby와 Seismic Unix를 이용해 Muting/Interpolation 스크립트도 작성할 수 있었습니다. 참고로 지금은 Ruby는 사용하지 않고, 과학/공학 라이브러리가 많은 Python을 주로 사용합니다.

위의 조언들 외에도

  • 20. 지식을 일반 텍스트로 저장하라.
  • 29. 코드를 작성하는 코드를 작성하라.
  • 36. 모듈간의 결합도를 최소화하라.
  • 49. 소프트웨어를 테스트하라.
  • 61. 수작업 절차를 사용하지 말라.
  • 68. 문서가 애초부터 전체의 일부가 되게 하고, 나중에 집어넣으려고 하지 말라.

등 다양한 조언들을 실천해보면서 프로그래밍 실력을 향상시킬 수 있었습니다. 우리가 전문 프로그래머는 아니지만 주로 프로그래밍을 이용해 연구를 수행하고 있으니 프로그래밍에 대해 공부하는 것도 많은 도움이 됩니다. 단, 이것 자체가 연구는 아니니 연구하는 중에 틈틈이 연습해보는 것이 좋겠죠.

Fortran config parser

포트란 언어에서 사용할 수 있는 configuration file parser 모듈을 github에 공개했습니다.

다음과 같은 configuration 파일에서 변수들을 읽어들일 수 있는 모듈입니다. 사용법은 github에 올렸습니다.

[DEFAULTS]
path = ../include
use_abs = True

[Section 1]
nmax = 30
# comment 1
vmin = 1.0
freqs = 5.0, 10.0, 30.0, 50.0
amps = 0.0, 1.0, 1.0, 0.0
path = ../text

[Section 2]
use_abs = no
; comment 2
my file = $(Section 1:path)/file.txt

속도 프로파일과 탄성파 트레이스 추출하여 그리기

속도모델에서 프로파일을 추출하여 깊이에 따라 속도 그림을 그려보겠습니다. 이진 형식의 속도파일에서 텍스트 파일로 프로파일을 추출한 후 그리는 방법과 이진 속도파일을 직접 읽어서 그리는 방법을 살펴보겠습니다. 참고로, 탄성파 공통송신원모음 등에서 트레이스를 추출하여 그리는 과정 또한 동일합니다.

텍스트 파일로 추출하여 그리기

바이너리 파일에서 프로파일 또는 트레이스를 추출하기 위해 gpl 라이브러리의 gplTracePick 프로그램을 사용하겠습니다. 이차원 단면(속도모델, 공통송신원모음 등)에서 세로 방향 트레이스를 추출할 때 사용하는 프로그램입니다. (가로방향 트레이스는 gplHTracePick 프로그램을 이용하면 됩니다.) 이 프로그램을 그냥 실행하면 아래와 같은 도움말이 나옵니다.


%%sh # 이 글을 쓰고 있는 jupyter notebook에서 shell 명령을 실행하기 위한 magic command입니다.
gplTracePick # 실제 터미널상에서 실행하는 명령어

 Gpl trace picker
 Required parameters:
     [i] n1=            : # of grids in fast dimension
     [s] fin=           : input binary file
     [s] fout=          : output binary file
     [i] pick=          : (=first), first pick (1~n2)
 Optional parameters:
     [i] last=first     : last pick (pick~n2)
     [i] step=1         : pick step
     [f] d1=1.0         : grid size
     [i] n2=calc        : # of grids in slow dimension
     [s] type=f         : data type [ifdcz]
     [s] otype=a        : output type [ab] (ascii/binary)

위에서 n1finfoutpick은 프로그램 실행시 필수적으로 넣어줘야 하는 값입니다.

  • n1은 세로 방향(fast dimension) 격자수
  • fin은 입력 파일 이름
  • fout은 출력 파일 이름
  • pick은 추출하고자하는 가로 방향(slow dimension) 격자 번호입니다. 격자 번호는 1번부터 시작합니다.

Marmousi 속도모델(nx=576, ny=188, h=0.016 km)에 대해 1.6 km 지점(격자번호 101)에서 시작하여 3.2 km 간격(200개 격자 간격)으로 3개의 속도 프로파일을 추출한다면 아래와 같이 실행할 수 있습니다.


%%sh
gplTracePick n1=188 d1=0.016 fin=marm16km.bin fout=vel_profile.txt pick=101 step=200 last=501

 n2=         576


     n1=188
     d1=0.016
     fin=marm16km.bin
     fout=vel_profile.txt
     pick=101
     step=200
     last=501

그 때 결과물은 아래와 같습니다. 첫 번째 열은 깊이 정보, 두 번째부터 네 번째 열까지는 추출한 속도 프로파일 정보입니다(1.6 km, 4.8 km, 8.0 km).

%%sh
head vel_profile.txt
   0.00000000       1.50000012       1.50000012       1.50000012    
   1.60000008E-02   1.50000012       1.50000012       1.50000012    
   3.20000015E-02   1.50000012       1.65800011       1.59800005    
   4.80000004E-02   1.66200006       1.66200006       1.60200012    
   6.40000030E-02   1.66600013       1.66600013       1.60600019    
   8.00000057E-02   1.67000008       1.73999715       1.69000006    
   9.60000008E-02   1.67400002       1.74399781       1.69400012    
  0.112000003       1.67800009       1.61800003       1.69800007    
  0.128000006       1.78200006       1.70200002       1.63200009    
  0.144000009       1.78600013       1.70600009       1.63600004    

텍스트 파일로 추출한 결과는 gnuplot과 같은 프로그램을 이용해 빠르게 확인해볼 수 있습니다. 여기서는 파이썬의 Matplotlib을 이용하여 위의 속도 프로파일을 그려보겠습니다.


%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np

trc=np.loadtxt("vel_profile.txt")

h=0.016
fs='large'

plt.figure(figsize=[15,5])
for i,ix in enumerate([100,300,500]):
    plt.plot(trc[:,0],trc[:,i+1],label="{0} km".format(ix*h))

plt.legend(loc="upper left",fontsize=fs)
plt.xlabel("Depth (km)",fontsize=fs)
plt.ylabel("Velocity (km/s)",fontsize=fs)

<matplotlib.text.Text at 0x10cc66c88>

png

이진 파일을 직접 읽어서 그리기

이번에는 파이썬에서 이진 형식의 속도모델 파일을 직접 읽어서 그려보겠습니다.


nx=576
ny=188
vel=np.fromfile("marm16km.bin",dtype=np.float32)
vel.shape=(nx,ny)

h=0.016
fs='large'
depth=np.arange(ny)*h

plt.figure(figsize=[15,5])
for ix in [100,300,500]:
    plt.plot(depth,vel[ix,:],label="{0} km".format(ix*h))

plt.legend(loc="upper left",fontsize=fs)
plt.xlabel("Depth (km)",fontsize=fs)
plt.ylabel("Velocity (km/s)",fontsize=fs)

<matplotlib.text.Text at 0x10d13dd30>

png

참고로, 파이썬은 배열 인덱스가 0번부터 시작하기 때문에 가로방향 100, 300, 500번 속도 프로파일을 가져다가 그렸습니다(gplTracePick을 이용하는 앞의 예제에서는 101, 301, 501번 격자 위치에서 추출했죠).

탄성파 트레이스 그리기

공통송신원모음에서 탄성파 트레이스를 추출하여 그리는 과정은 속도모델에서 프로파일을 추출하여 그리는 경우와 동일합니다. 아래는 샘플 개수가 723개, 샘플링 간격 4 ms, 트레이스가 96개인 공통송신원모음 파일(marm3000.bin)에서 31번째와 61번째 트레이스를 그리는 예제입니다.


ntr=96
ns=723
dt=0.004
trc=np.fromfile("marm3000.bin",dtype=np.float32)
trc.shape=(ntr,ns)

fs='large'
time=np.arange(ns)*dt

plt.figure(figsize=[15,5])
for itr in [30,60]:
    plt.plot(time,trc[itr,:],label="trace {0}".format(itr+1))
plt.legend(loc="upper left",fontsize=fs)
plt.xlabel("Time (s)",fontsize=fs)
plt.ylabel("Amplitude",fontsize=fs)
plt.xlim([0,ns*dt])

(0, 2.892)

png

속도모델 그림 그리기

두 가지 방법으로 2차원 속도모델을 그려보겠습니다. 첫 번째 방법은 SU의 psimage를 이용하는 방법, 두 번째는 python의 matplotlib을 이용하는 방법입니다.

psimage로 그리기

첫 번째 방법부터 보겠습니다. psimage는 쉘에서 사용하는 명령어이지만, gpl 라이브러리의 psplot 모듈을 이용하면 python 명령을 통해 간편하게 속도모델을 그릴 수 있습니다. Marmousi 속도모델을 그림으로 그려보겠습니다.


from gpl.psplot import plot

nx=576
ny=188
h=0.016
fin="marm16km.drt"

opt = "n1={0} d1={1} d2={1} d1num=1 lbeg=1.5 lend=5.5".format(ny,h,h)
plot.velocity("marm16km.png", fin, opt)

psimage label1="Depth (km)" legend=1 d2s=0.5 lheight=1.0 lstyle="vertright" label2="Distance (km)" height=1.0 labelsize=8 lwidth=0.1 d1s=0.5 width=2.65  n1=188 d1=0.016 d2=0.016 d1num=1 lbeg=1.5 lend=5.5 < marm16km.drt > marm16km.eps

// adding velocity unit (km/s)

// fixing bounding box

// converting .eps to .png ..

vel(marm16km.png)

velocity_color를 이용해 컬러로 그릴 수도 있습니다.


plot.velocity_color("marm16km_color.png",fin,opt)

psimage label1="Depth (km)" ghls="0.33,0.5,1" bps=24 bhls="0.67,0.5,1" d1s=0.5 lwidth=0.1 whls="0,0.5,1" legend=1 d2s=0.5 lheight=1.0 lstyle="vertright" label2="Distance (km)" height=1.0 labelsize=8 width=2.65  n1=188 d1=0.016 d2=0.016 d1num=1 lbeg=1.5 lend=5.5 < marm16km.drt > marm16km_color.eps

// adding velocity unit (km/s)

// fixing bounding box

// converting .eps to .png ..

vel(marm16km_color.png)

Matplotlib으로 그리기

두 번째 방법은 python의 matplotlib 라이브러리를 이용하는 방법입니다. 이를 위해서는 코드에서 numpy를 이용해 속도모델을 읽어들인 후에 matplotlib으로 그립니다. 속도모델을 그리는 부분은 함수로 작성하였는데, 필요에 따라 수정해서 사용하면 되겠습니다.


%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np

def plot_vel(vel, h, figsize=[15,4], unit='km/s', xticks=None, yticks=None, cticks=None, cmap='gray_r', fontsize=20):
    xmax=(vel.shape[0]-1)*h
    ymax=(vel.shape[1]-1)*h

    plt.figure(figsize=figsize)
    plt.imshow(vel.transpose(),extent=(0,xmax,ymax,0),cmap=cmap)

    # x,y labels
    plt.xlabel('Distance (km)',fontsize=fontsize)
    plt.ylabel('Depth (km)',fontsize=fontsize)

    # x,y ticks, tick labels
    plt.tick_params(labelsize=fontsize)
    plt.gca().xaxis.tick_top()
    plt.gca().xaxis.set_label_position("top")
    xticks and plt.xticks(xticks)
    yticks and plt.yticks(yticks)

    # colorbar
    cb=plt.colorbar(shrink=1.0,pad=0.01,aspect=10,ticks=cticks)
    plt.clim([vel.min(),vel.max()])
    cb.set_label(unit,fontsize=fontsize)
    ct=plt.getp(cb.ax,'ymajorticklabels')
    plt.setp(ct,fontsize=fontsize)

# 속도모델 읽기
vel=np.fromfile(fin,dtype=np.float32)
vel.shape=(nx,ny)

yticks=[0,1,2] # y축 ticks
cticks=[2,3,4,5] # colorbar ticks
plot_vel(vel,h,yticks=yticks,cticks=cticks)

png

# 컬러로 그리고(cmap='jet') 파일로 저장하기
plot_vel(vel,h,xticks=[0,3,6,9],cmap='jet')
plt.savefig("vel.png",bbox_inches='tight')

png

결과물로 저장한 vel.png 파일은 다음과 같습니다.

vel(vel.png)

Matplotlib을 이용한 탄성파 자료처리 그림 그리기

이전 글에서는 SU 명령어들을 이용해 탄성파 자료처리 결과 확인용 그림을 그리는 방법을 살펴보았습니다. 이번에는 Python의 Matplotlib을 이용하여 그린 그림 예제들을 보겠습니다. 그림은 IPython Processing 모듈을 이용해 그렸으며, 그릴 때 사용한 코드는 github에서 볼 수 있습니다.

속도모델, 구조보정 영상

우선, 다음과 같이 이진파일로부터 2차원 속도모델과 구조보정 결과를 그릴 수 있습니다. 기본적으로 속도모델은 컬러, 구조보정 영상은 흑백으로 그리도록 했지만, 필요에 따라 코드를 수정해서 색상을 바꿀 수 있습니다. 색상을 바꾸고 싶을 경우 imshow 함수의 cmap 인자를 이용하면 됩니다.

%matplotlib inline
from pkprocess import *
import numpy as np

vel = np.fromfile("marm16km.drt", dtype=np.float32)

nx = 576
nz = 188
h = 0.016
vel.shape = (nx, nz)

plot_vel(vel, h)

png


plot_mig(vel,h)

png

공통송신원 모음, 스펙트럼

그리고 SU 파일로부터 공통송신원 모음이나 F-X, F-K 스펙트럼을 그릴 수 있습니다. 공통송신원 모음은 Wiggle trace 또는 이미지로 그릴 수 있고, 이미지 색상은 cmap으로 조절 가능합니다.


su = read_su("marm3000.su")

plot_wiggle(su, perc=97)

min=-616.05078125 max=613.4453125

png


plot_image(su, perc=97)

min=-616.05078125 max=613.4453125

png


plot_image(su, perc=97, cmap='bwr')

min=-616.05078125 max=613.4453125

png


specfx(su)

dt=0.004, fmax=125.0

png


specfk(su)

dt=0.004, fmax=125.0

dx=0.025, kmax=20.0

png

위의 그림들 모두 Matplotlib으로 그렸으므로, 수정이 필요할 경우 Matplotlib 문서를 참고하여 수정해서 사용하시면 되겠습니다.

탄성파 자료처리 그림 그리기

탄성파 자료처리 결과 그림을 쉽게 그리는 방법을 살펴보겠습니다.

탄성파 자료처리를 하다 보면 결과물을 그림으로 확인해야 하는 경우가 많습니다. 특별히 노력해서 그려야 하는 그림도 있지만 속도모델, 공통송신원모음 등 대부분의 그림은 거의 비슷한 명령으로 그릴 수 있습니다. 개인적으로 논문이나 발표자료에 넣을 그림을 그릴 때 Seismic Un*x(SU)를 많이 이용하는데, 몇 가지 자주 그리는 그림들을 쉽게 그릴 수 있도록 파이썬 모듈을 만들었습니다. 모듈은 gpl라이브러리에 포함되어 있습니다. 최근 python 3 용으로 수정하였습니다.

먼저 속도모델을 예로 들어보겠습니다. 그림을 그리기 위한 코드는 다음과 같습니다.


from gpl.psplot import plot

vel="marm16km.drt"
opt="n1=188 d1=0.016 d2=0.016 d1num=1 d2num=2"

plot.velocity_color("vel_color.png",vel,opt)

위 코드는 gpl.psplot 모듈에서 plot을 가져오고, marm16km.drt 파일로부터 opt 문자열의 옵션을 이용하여 vel_color.png 파일을 생성하는데, 컬러로 된 속도모델 그림으로 만들라는 코드입니다.

velocity_color는 그림 종류를 지정하는 명령인데, 현재 다음과 같은 명령들을 지원합니다.

  • velocity(target, source, option, unit=”km/s”)
  • velocity_color(target, source, option, unit=”km/s”)
  • gradient(target, source, option)
  • gradient_color(target, source, option)
  • migration(target, source, option)
  • contour(target, source, option)
  • seismogram(target, source, option)
  • spectrum(target, source, option)

위의 명령들은 SU를 이용해 해당 그림을 그리라는 명령으로, contourpscontour를 사용하고 나머지는 psimage를 사용합니다. 입력 파일이 SU 파일이라면 supscontour 또는 supsimage를 사용합니다.

인자들 중 target은 출력 파일, source는 입력 파일, option은 그림 그릴 때 사용할 옵션입니다. 그림 종류에 따라 기본적으로 몇 가지 옵션이 들어가있는데, n1, d1, d2와 같이 입력 파일에 따라 달라지는 옵션을 option에 넣어주면 됩니다. 그리고 기본 옵션을 덮어쓰고 싶은 경우에도 option에 추가해줍니다.

속도모델의 단위는 기본적으로 km/s로 지정해 놓았는데, 필요에 따라 바꿔서 사용할 수 있습니다. g/cc로 바꾸면 밀도 모델을 그릴 수도 있겠죠. migration은 snapshot을 그릴 때 사용할 수도 있습니다.

SU 명령은 기본적으로 eps 파일을 생성합니다. target을 eps 외의 다른 파일(png, tiff, jpg 등)로 지정하면ImageMagickconvert 명령을 이용해 eps 파일을 변환합니다.

따라서 본 모듈의 모든 기능을 이용하려면 Python, SU, ImageMagick이 필요합니다.

터미널에서 위의 코드를 실행했을 때 나오는 메시지는 다음과 같습니다.

psimage height=1.0 width=2.65 d2s=0.5 lwidth=0.1 lstyle="vertright" lheight=1.0
label2="Distance (km)" ghls="0.33,0.5,1" bps=24 whls="0,0.5,1" legend=1
bhls="0.67,0.5,1" labelsize=8 label1="Depth (km)" d1s=0.5  n1=188 d1=0.016
d2=0.016 d1num=1 d2num=2 < marm16km.drt > vel_color.eps

psimage: bclip=5.5 wclip=1.5

// adding velocity unit (km/s)

// fixing bounding box
Original:  %%BoundingBox: 66 41 353 207
Updated:   %%BoundingBox: 85 104 324 202

// converting .eps to .png ..

내용을 살펴보면 다음 순서로 실행됩니다.

  1. SU의 psimage 명령을 이용해 속도모델 eps 파일을 생성합니다. 옵션은 컬러 속도모델에 맞춰서 들어갑니다. 참고로, 그림 크기는 Geophysics 논문 기준에 맞춘 것입니다.
  2. km/s라는 단위를 넣어줍니다(postscript 수정).
  3. 그림 여백을 조절합니다(bounding box 수정).
  4. eps 파일을 png 파일로 수정합니다.

그리고, 결과물인 vel_color.png은 다음과 같습니다. output(vel_color.png)

아래 코드와 다른 그림 예시를 올리니 필요한 그림에 해당하는 명령을 사용하시면 되겠습니다.


from gpl.psplot import plot

vel="marm16km.drt"
opt="n1=188 d1=0.016 d2=0.016 d1num=1 d2num=2"

plot.velocity("vel.png",vel,opt+"lbeg=1.5 lend=5.5 lfnum=1.5")
plot.velocity_color("vel_color.png",vel,opt)
plot.velocity_color("density_color.png",vel,opt,unit="g/cc")
plot.gradient("grad.png",vel,opt)
plot.gradient_color("grad_color.png",vel,opt)
plot.migration("mig.png",vel,opt)
plot.contour("contour.png",vel,opt)

seismo="marm3000.su"
opt2="f2=0 d2=0.025 d1s=0.5 d2s=0.5"
plot.seismogram("seismo.png",seismo,opt2)

spec="marm3000fx.su"
plot.spectrum("spec.png",spec,opt2)

output(vel.png)

output(vel_color.png)

output(density_color.png)

output(grad.png)

output(grad_color.png)

output(mig.png)

output(contour.png)

output(seismo.png)

output(spec.png)

Muting scripts

SU 파일 뮤팅 스크립트를 공개하였습니다. SU 파일을 손으로 뮤팅할 때 모든 shot gather를 뮤팅하는 것이 아니라 건너 뛰면서 뮤팅한 후 건너뛴 부분은 interpolation하도록 만들었습니다. 해상 스트리머 탐사 자료에 사용할 수 있습니다. Github에서 받을 수 있습니다.

 

참고 문헌

Kim, A., D. Ryu, W. Ha, and C. Shin, 2015, An efficient first arrival picking procedure for marine streamer data, Geosystem Engineering,  18(5), 251-258.

IPython Processing

학부 탄성파 자료 처리 수업에서 사용하기 위한 파이썬 패키지입니다. 기존의 공개된 Matlab 패키지를 파이썬으로 변환하였습니다. 간단한 2차원 자료 처리(time processing)가 가능합니다.

육상, 해상 자료 처리 예제와 함께 공개하였고, Github에서 받으실 수 있습니다.

참고 문헌

하완수, 2015, 대화식 탄성파 자료 처리 수업을 위한 파이썬 패키지 개발, 한국자원공학회지, 52(4), 414-421.

SU 파일 입출력을 위한 포트란 라이브러리

포트란 사용자들이 Seismic Un*x의 SU 파일을 읽고, 수정하고, 쓸 수 있도록 하는 입출력 라이브러리를 작성하여 공개하였습니다.

Github에서 받으실 수 있습니다.

매뉴얼은 ReadTheDocs를 참고하세요.

참고문헌

하완수, 2015, SU 파일 입출력을 위한 포트란 라이브러리 개발, 한국자원공학회지, 52(1), 81-90.