[Fortran, Gnuplot, ShellScript] 기상 자료를 이용한 통계 분석 및 가시화

 정보

• 업무명    : 기상 자료를 이용한 통계 분석 및 가시화

• 작성자    : 이상호

• 작성일    : 2019-08-31

• 설   명    :

• 수정이력 :

 

 내용

[특징]

• 기상 자료의 특성을 파악하기 위해 통계 분석 및 가시화 도구가 필요하며 이 프로그램은 이러한 목적을 달성하기 위해 고안된 소프트웨어

 

[기능]

• 구분자 (,)로 된 데이터셋을 읽기

• 날짜 (년, 월)에 따라 자료 정제 및 통계값 계산

• 통계값, 도수분포, 누적도수 자료 저장

• 저장된 자료를 바탕으로 빈도분포 및 상자 그림 가시화

 

[활용 자료]

• 자료 : 전천 일사, 직달 일사, 지구 복사, 기온

• 기간 : 2009년 12월 19일 - 2018월 12월 31일

• 해상도 : 매 1분

• 제공처 : 강릉원주대 생명대학 2호관 옥상 (복사 관측소)

• 관련 자료

 

표. 강릉원주대 복사 관측소의 일사계 세부 정보

	직달일사계	전천일사계	지구복사계
모델 (제조사)	CHP1 (Kipp & Zonen Inc.)	CMP21 (Kipp & Zonen Inc.)	CGR3 (Kipp & Zonen Inc.)
ISO 9060 분류	Secondary Standard	Secondary Standard	Secondary Standard
시간 해상도	매 1분	매 1분	매 1분
측정 파장	200-4000 nm	285-2800 nm	4500-42000 nm
관측 기간	2009-12-19 - 현재	2009-12-19 - 현재	2013-08-26 - 현재

 

그림. 강릉원주대 복사 관측소의 관측 장비 및 현황

 

그림. 강릉원주대 복사 관측소의 데이터셋 정보

 

[자료 처리 방안 및 활용 분석 기법]

• 기상 자료를 이용하여 통계 분석 및 가시화를 위해서 각 단계별로 수행과정은 다음과 같다.

• 1 단계 : 입력 자료를 구축

 

• 2 단계 : 각 날짜 (년, 월)마다 자료를 정제하여 통계값, 도수분포, 누적도수 계산하여 자료 저장

• 2010년 1-12월 통계값 계산 (2010_Total_STA.dat)

 

• 2010년 1-12월 도수분포, 누적도수 계산 (2010_Total_FRE.dat)

 

• 저장된 자료를 통해 가시화 및 분석

 

[사용법]

• 입력 자료를 동일 디렉터리에 위치

• 소스 코드를 실행 (csh Visualization_Of_Gaussian_Distribution_Using_Weather_Data.csh)

• 빈도분포 및 상자 그림 가시화 결과를 확인

 

[사용 OS]

• Linux

 

[사용 언어]

• Fortran

• Gnuplot

• ShellScript (csh)

 

 소스 코드

[기상 자료에서 통계 분석을 위한 전처리 및 가시화]

• 기상 자료에서 통계 분석을 위한 전처리

	#!/bin/csh
	rm -f Total_*.dat FIG/* DATA/*
	cat >! Preprocessing_For_Statistical_Analysis.f90 << EOF
	implicit none
	real, allocatable, dimension( : , : , : , : , : , :) :: var
	integer, allocatable, dimension( : ) :: SS
	character(200) :: dat, fn1, fn2
	integer :: yy, mon, dd, tt, mm, cc, i, j, k
	integer, parameter :: Max_row=99999
	real, allocatable, dimension( : , : , : , :) :: var2
	real :: Temp
	integer :: row
	allocate ( var(1:20,1:12,1:31,0:23,0:59,1:30), SS(1:100) )
	! call system("rm -rf OUTPUT/*")
	open(10,file='INPUT/DATA.dat',action='read',status='old')
	2 read(10,'(a)',end=99) dat
	read(dat(2:5),'(i4)') yy
	read(dat(7:8),'(i2)') mon
	read(dat(10:11),'(i2)') dd
	read(dat(13:14),'(i2)') tt
	read(dat(16:17),'(i2)') mm
	cc=0
	do i=22,200
	if(dat(i:i) == ',' .or. dat(i:i) == ' ') then
	cc=cc+1 ; SS(cc)=i
	endif
	if(dat(i:i) == ' ') goto 1
	enddo
	1 do i=1,cc-1
	if(SS(i+1)-1 >= SS(i)+1) then
	if(dat(SS(i)+1:SS(i)+1) /= '"') read(dat(SS(i)+1:SS(i+1)-1),*) var(yy,mon,dd,tt,mm,i)
	if(dat(SS(i)+1:SS(i)+1) == '"') var(yy,mon,dd,tt,mm,i)=-999.
	else
	var(yy,mon,dd,tt,mm,i) = -999.
	endif
	enddo
	if(var(yy,mon,dd,tt,mm,2) > -20) then
	write(fn1,'(i4,a,i2.2)') yy,'_',mon
	open(11,file='OUTPUT/'//trim(fn1)//'.dat')
	write(11,'(i5,4(1x,i2.2),f10.4)') yy,mon,dd,tt,mm,var(yy,mon,dd,tt,mm,2:2)
	endif
	goto 2
	99 end
	EOF
	pgf90 Preprocessing_For_Statistical_Analysis.f90 ; a.out
	foreach yy (2010)
	foreach mon (11)
	set mon1 = `printf "%02d" $mon`
	set fn1 = `echo ${yy}_${mon1}`
	echo $yy $mon
	cat >! Statistical_Analysis_Using_Weather_Data.f90 << EOF
	implicit none
	integer, parameter :: Max_row = 99999
	real, parameter :: pi = 4.d0*atan(1.d0)
	real, allocatable, dimension( : , : ) :: var2
	integer, allocatable, dimension( : ) :: var3
	integer, dimension(12) :: mm
	character(200) :: fn2
	real :: Temp, Xran
	integer :: row, i, j, k, yy, mon, cc
	real :: ST_MAX, ST_MIN, ST_AVE, ST_SUM, ST_K, ST_C, ST_C2, ST_MOD, ST_CC
	real :: ii, STA_MID, wa, STA_RF, STA_CC, wb
	real :: ST_VAR, ST_SDE, ST_CV, ST_SK, ST_SKEW, ST_KN, ST_KNRT
	real :: low, mid, upp, ST_LOW, ST_MID, ST_UPP, IQR, ST_LOW_IQR, ST_UPP_IQR
	real :: ST_GAU
	mm = reshape( (/31,28,31,30,31,30,31,31,30,31,30,31/), (/12/) )
	yy = $yy
	mon = $mon
	if( (mod(yy,4) == 0 .and. mod(yy,100) /= 0) .or. mod(yy,400) == 0 ) mm(2)=29
	write(fn2,'(i4,a,i2.2)') yy,'_',mon
	open(10,file='OUTPUT/'//trim(fn2)//'.dat')
	row = 0
	do j=1,Max_row
	read(10,*,end=100)
	row=row+1
	enddo
	100 rewind(10)
	allocate( var2(1:row, 1:10) )
	allocate( var3(1:row) )
	do i=1,row
	read(10,*) var2(i,1:6)
	enddo
	!********************* Sort Value **********************
	do i=1,row
	if(var2(i,6) /= -999) then
	do j=i+1,row
	if(var2(i,6) >= var2(j,6)) then
	temp=var2(j,6) ; var2(j,6)=var2(i,6) ; var2(i,6)=temp
	endif
	enddo
	Xran=real(mon)+((var2(i,3)-1.)/real(mm(mon)))+(var2(i,4)/(24.*real(mm(mon))))+(var2(i,5)/(3600.*real(mm(mon))))
	write(11,'(i5,4(1x,i2.2),2f10.5)') int(var2(i,1:5)), var2(i,6), Xran
	endif
	enddo
	!******************** STAT_FUN **********************************
	!************* ST_MOD Algo **********************
	cc=0
	do i=1,row
	if(var2(i,6)==var2(i+1,6)) then
	cc=cc+1 ; var3(i)=cc
	else
	cc=1
	endif
	enddo
	!*******************************************************
	ST_CC = count(var2(1:row,6) /= -999.)
	ST_SUM = sum(var2(1:row,6), mask=var2(1:row,6) /= -999.)
	ST_AVE = ST_SUM/real(ST_CC)
	ST_MAX = maxval(var2(1:row,6))
	ST_MIN = minval(var2(1:row,6))
	ST_K = 1.+(3.32*log10(real(ST_CC)))
	ST_C = (ST_MAX-ST_MIN)/ST_K
	ST_VAR = sum((var2(1:row,6)-ST_AVE)**2., mask=var2(1:row,6) /= -999.)/(ST_CC-1.)
	ST_SDE = sqrt(ST_VAR)
	ST_CV = ST_SDE/ST_AVE
	ST_SK = sum((var2(1:row,6)-ST_AVE)**3., mask=var2(1:row,6) /= -999.)/(ST_CC-1.)
	ST_KN = sum((var2(1:row,6)-ST_AVE)**4., mask=var2(1:row,6) /= -999.)/(ST_CC-1.)
	ST_SKEW = ST_SK/(ST_SDE**3.)
	ST_KNRT = ST_KN/(ST_SDE**4.)
	low=real(ST_CC)*0.25
	mid=real(ST_CC)*0.5
	upp=real(ST_CC)*0.75
	ST_LOW=var2(int(low)+1,6)
	ST_MID=var2(int(mid)+1,6)
	ST_UPP=var2(int(upp)+1,6)
	if(low==int(low)) ST_LOW = (var2(int(low),6)+var2(int(low)+1,6))/2.
	if(mid==int(mid)) ST_MID = (var2(int(mid),6)+var2(int(mid)+1,6))/2.
	if(upp==int(upp)) ST_UPP = (var2(int(upp),6)+var2(int(upp)+1,6))/2.
	IQR = ST_UPP-ST_LOW
	ST_LOW_IQR = ST_LOW-(1.5*IQR)
	ST_UPP_IQR = ST_UPP+(1.5*IQR)
	write(12,'(/,a,i5,a,i3)') 'Year =',yy, ' Mon =',mon
	write(12,*) ' LOW_IQR MIN AVE &
	MID MOD MAX UPP_IQR'
	write(12,'(7f12.2)') ST_LOW, ST_MIN, ST_AVE, ST_MID, var2(maxloc(var3(1:row)),6), ST_MAX, ST_UPP
	write(12,*) ' VAR SDE CV &
	SKEW KNRT COUNT'
	write(12,'(5f12.2,i13)') ST_VAR, ST_SDE, ST_CV, ST_SKEW, ST_KNRT, int(ST_CC)
	!********** Frequency Diagram *******************************
	write(13,'(/,a,i5,a,i3,3f10.2)') 'Year =', yy, ' Mon =',mon, (ST_MIN-ST_C)+(ST_C/2.), ST_MAX, ST_C
	write(13,*) 'Count Range Mid_value Fre &
	Relative_Fre Cumulative_Fre Cumulative_R_Fre'
	wa=0. ; wb=0. ; k=0
	do ii=ST_MIN-ST_C, ST_MAX, ST_C
	STA_MID = (ii+(ii+ST_C))/2.
	STA_CC = count(var2(1:row,6) >= ii)-count(var2(1:row,6) >= ii+ST_C)
	STA_RF = STA_CC/ST_CC
	k=k+1 ; wa=wa+STA_RF ; wb=wb+STA_CC
	write(13,'(i4,f10.2,a3,f6.2,f10.2,i10,f10.2,i16,f18.2)') &
	k, ii,' ~ ',ii+ST_C, STA_MID, int(STA_CC), STA_RF, int(wb), wa
	enddo
	write(14,'(2i5,7f10.2)') yy, mon, ST_LOW, ST_MIN, ST_MAX, ST_UPP, ST_AVE, ST_MID, var2(maxloc(var3(1:row)),6)
	!************* Over Value / Gaussian Distribution **************************
	do i=1,row
	if(var2(i,6) /= -999.) then
	if(var2(i,6) < ST_LOW_IQR .or. var2(i,6) > ST_UPP_IQR) write(15,*) int(var2(i,1:2)), var2(i,6)
	ST_GAU = (1./(ST_SDE*sqrt(2.*pi)))*exp(-(((var2(i,6)-ST_AVE)**2.)/(2.*(ST_SDE**2.))))
	write(16,*) var2(i,6), ST_GAU*10
	endif
	enddo
	end
	EOF
	pgf90 Statistical_Analysis_Using_Weather_Data.f90 ; a.out
	csh Visualization_Of_Gaussian_Distribution_Using_Weather_Data.csh
	mv -f 1.gif FIG/${fn1}.gif
	cat fort.11 >> Total_Time.dat
	cat fort.12 >> Total_STA.dat
	cat fort.13 >> Total_FRE.dat
	cat fort.14 >> Total_BOX.dat
	end
	csh Visualization_Of_Gaussian_Distribution_Using_Weather_Data.csh
	mv -f 2.gif FIG/${yy}.gif
	mv -f Total_STA.dat DATA/${yy}_Total_STA.dat
	mv -f Total_FRE.dat DATA/${yy}_Total_FRE.dat
	mv -f Total_BOX.dat DATA/${yy}_Total_BOX.dat
	mv -f Total_Time.dat DATA/${yy}_Total_Time.dat
	end
	rm -f Preprocessing_For_Statistical_Analysis.f90 Statistical_Analysis_Using_Weather_Data.f90 1.gs 2.gs fort.*

view raw Preprocessing_And_Visualization_For_Statistical_Analysis_Of_Weather_Data.csh hosted with ❤ by GitHub

 

• 기상 자료를 활용한 가우시안 분포 가시화

	#!/bin/csh
	set yy = `sed -n '2p' fort.13 | awk '{print $3}'`
	set mon = `sed -n '2p' fort.13 | awk '{print $6}'`
	set mon1 = `printf "%02d" $mon`
	set Xran_max = `sed -n '2p' fort.13 | awk '{print $8}'`
	set Xran_min = `sed -n '2p' fort.13 | awk '{print $7}'`
	set dXran = `sed -n '2p' fort.13 | awk '{print $9}'`
	set Xran_MM = `echo $Xran_min +11 | bc`
	set Xran_mm = `echo $Xran_min +1 | awk '{printf "%.2f", $1 +$2}'`
	set LOW_IQR = `sed -n '4p' fort.12 | awk '{print $1}'`
	set MIN = `sed -n '4p' fort.12 | awk '{print $2}'`
	set AVE = `sed -n '4p' fort.12 | awk '{print $3}'`
	set MID = `sed -n '4p' fort.12 | awk '{print $4}'`
	set MOD = `sed -n '4p' fort.12 | awk '{print $5}'`
	set MAX = `sed -n '4p' fort.12 | awk '{print $6}'`
	set UPP_IQR = `sed -n '4p' fort.12 | awk '{print $7}'`
	set VAR = `sed -n '6p' fort.12 | awk '{print $1}'`
	set SDE = `sed -n '6p' fort.12 | awk '{print $2}'`
	set CV = `sed -n '6p' fort.12 | awk '{print $3}'`
	set SKEW = `sed -n '6p' fort.12 | awk '{print $4}'`
	set KNRT = `sed -n '6p' fort.12 | awk '{print $5}'`
	set COUNT = `sed -n '6p' fort.12 | awk '{print $6}'`
	cat >! 1.gs << EOF
	set terminal post enhanced color font "Time-Roman, 18" background rgb "white"
	set output "1.gif"
	set multiplot
	set size 1,1
	set origin 0,0
	set label "Max = $MAX" at $Xran_mm, 9600 font "bold, 16"
	set label "Upper = $UPP_IQR" at $Xran_mm, 9200 font "bold, 16"
	set label "Ave = $AVE" at $Xran_mm, 8800 font "bold, 16"
	set label "Mid = $MID" at $Xran_mm, 8400 font "bold, 16"
	set label "Mod = $MOD" at $Xran_mm, 8000 font "bold, 16"
	set label "Lower = $LOW_IQR" at $Xran_mm, 7600 font "bold, 16"
	set label "Min = $MIN" at $Xran_mm, 7200 font "bold, 16"
	set label "Var = $VAR" at $Xran_MM, 9600 font " bold, 16"
	set label "Sde = $SDE" at $Xran_MM, 9200 font " bold, 16"
	set label "CV = $CV" at $Xran_MM, 8800 font " bold, 16"
	set label "Skew = $SKEW" at $Xran_MM, 8400 font " bold, 16"
	set label "Knrt = $KNRT" at $Xran_MM, 8000 font " bold, 16"
	#set label "N = $COUNT" at $Xran_MM, 7600 font "bold, 16"
	set key noopaque
	set key top right
	set key font "bold, 16"
	set format x "%1.0f" numeric
	set title "${yy}.${mon1}"
	set title font "bold, 26" textcolor rgb "black"
	set tics nomirror
	set xlabel "Temperature({\312C})"
	set xlabel font "bold, 18" textcolor rgb "black" norotate
	set xrange [ $Xran_min : $Xran_max ] noreverse nowriteback
	set xtics out $Xran_min, $dXran, $Xran_max
	set ylabel "Relative Frequency [$COUNT]"
	set ylabel font "bold, 18" textcolor rgb "black" rotate by -270
	set yrange [ 0 : 10000 ] noreverse nowriteback
	set ytics out 0, 1000, 10000
	set mytics 2
	set y2label "Gaussian Distribution"
	set y2label font "bold, 18" textcolor rgb "black" rotate by -270
	set y2range [ 0 : 1.0 ] noreverse nowriteback
	set y2tics out 0, 0.1, 1.0
	set my2tics 2
	plot "<(sed -n '3,100p' fort.13)" u 5:6 axes x1y1 w boxes lc "red" lw 4 t 'FRE',\
	"fort.16" u 1:2 axes x1y2 w l lc "blue" lw 3 t 'GAU'
	plot "<(sed -n '3,100p' fort.13)" u 5:6 axes x1y1 w l lc "red" lw 4 t "FRE"
	unset multiplot
	quit
	EOF
	gnuplot 1.gs
	mogrify -rotate 90 1.gif
	#display 1.gif

view raw Visualization_Of_Gaussian_Distribution_Using_Weather_Data.csh hosted with ❤ by GitHub

 

• 기상 자료를 이용한 상자 그림 가시화

	#!/bin/gnuplot
	set yy = `sed -n '2p' fort.13 | awk '{print $3}'`
	cat >! 2.gs << EOF
	set terminal post enhanced color font "Time-Roman, 18" background rgb "white"
	set output "2.gif"
	set multiplot
	set size 1,1
	set origin 0,0
	set key noopaque
	set key top right
	set key font "bold, 16"
	set title "$yy"
	set title font "bold, 26" textcolor rgb "black"
	set tics nomirror
	set xlabel "Month"
	set xlabel font "bold, 18" textcolor rgb "black" norotate
	set xrange [ 0 : 13 ] noreverse nowriteback
	set xtics out 1, 1, 12
	set mxtics 2
	set ylabel "Temperature({\312C})"
	set ylabel font "bold, 18" textcolor rgb "black" rotate by -270
	set yrange [ -30 : 50 ] noreverse nowriteback
	set ytics out -30, 10, 50
	set mytics 2
	set bar 0.5
	set boxwidth 0.4
	set style fill
	plot "Total_BOX.dat" u 2:3:4:5:6 with candlesticks notitle whiskerbars ,\
	"" u 2:7:7:7:7 w candlesticks lc 3 t "AVE" ,\
	"" u 2:9:9:9:9 w candlesticks lc 4 t "MOD" ,\
	"" u 2:8:8:8:8 w candlesticks lc 2 t "MID"
	plot "Total_BOX.dat" u 2:7 w l lc 3 t "AVE" ,\
	"" u 2:9 w l lc 4 t "MOD"
	unset multiplot
	quit
	EOF
	gnuplot 2.gs
	mogrify -rotate 90 2.gif
	#display 2.gif

view raw Visualization_Of_Box_Plot_Using_Weather_Data.csh hosted with ❤ by GitHub

 

 

 결과

[빈도 분포 및 가우시안 분포]

• 2010년 11월 분포는 12월보다 자료 개수가 적기 때문에 정규 분포를 따르지 않음

그림. 2010년 11월 빈도 분포 및 가우시안 분포 가시화 (적색 박스: 빈도 분포, 청색 선: 가우시안 분포).

 

그림. 2010년 12월 빈도 분포 및 가우시안 분포 가시화 (적색 박스: 빈도 분포, 청색 선: 가우시안 분포).

 

[상자 그림 가시화]

• 2010년 상자 그림에서 온도 차이는 여름보다 겨울이 큼

그림. 2010년 상자 그림 가시화 (청색 선: 평균값, 황색 선: 최빈값, 녹색 선: 중앙값, 보라색 선: 상자 그림).

 

 참고문헌

[논문]

• 없음

[보고서]

• 없음

[URL]

• 없음

 

블로그에 대한 궁금하신 점을 문의하시면 자세히 답변드리겠습니다.

E. sangho.lee.1990@gmail.com & saimang0804@gmail.com