1차 정리된 파일 불러와서 자료 가공하기

기본 설정

setwd("~/Rclass/")

설정된 작업 디렉토리를 확인하는 방법은 ‘getwd()’ 입니다.

getwd()

## [1] "D:/OneDrive - hufs.ac.kr/Workspace_hufs/Rclass/Ehwa_240227"

위와 같이 “D:/Rclass/”가 작업 디렉토리로 설정된 것을 볼 수 있습니다.

다음으로는 분석에 필요한 라이브러리를 로드해줍니다.

library("lubridate")

## 
## 다음의 패키지를 부착합니다: 'lubridate'

## The following objects are masked from 'package:base':
## 
##     date, intersect, setdiff, union

library("plyr")
library("openair")

aero.meteo.CSV 불러오기 * aero.meteo.CSV 자료는 sample_data.csv

aero.meteo = read.csv(file = "sample_data.csv", header = TRUE)
aero.meteo$date = as.POSIXct(aero.meteo$Time, tz = "Asia/Seoul")
aero.meteo$date[1:10]

##  [1] "2020-12-15 00:00:00 KST" "2020-12-15 01:00:00 KST"
##  [3] "2020-12-15 02:00:00 KST" "2020-12-15 03:00:00 KST"
##  [5] "2020-12-15 04:00:00 KST" "2020-12-15 05:00:00 KST"
##  [7] "2020-12-15 06:00:00 KST" "2020-12-15 07:00:00 KST"
##  [9] "2020-12-15 08:00:00 KST" "2020-12-15 09:00:00 KST"

역궤적 자료 처리 및 기존 자료와 합치기

역궤적 자료는 trajectories.zip
기본적으로 linux shell scipt 기반으로 합치는 것이 빠르로 편함

gawk '$1 == 1 && $2 == 1 {print $0}' 2020*.txt > ../bulgwang.traj

윈도우에서 리눅스 쓰는 방법!

1. Windows 10에서 PowerShell(x64)을 관리자 권한으로 실행하고 아래 명령어를 실행한다.

-    Enable-WindowsOptionalFeature -Online -FeatureName Microsoft-Windows-Subsystem-Linux

1. WSL이 설치가 완료되면 재부팅 하라는 메세지가 뜬다. 간단히 재부팅…
1. 마이크로소프트 스토어에 들어가 Ubuntu를 찾아 설치한다.
1. 설치가 완료되면 ‘실행’ 버튼을 눌러 설치를 마무리.
1. PowerShell이 뜨고 설치가 마무리되게 된다. 바로 사용을 위한 계정 생성 프롬프트를 볼 수 있는데 원하는 user name과 passwd를 설정하면 바로 사용이 가능하다.
하지만, r로도 합칠 수 있음! (수십배는 시간걸림)

flist = list.files(path = "trajectories/", pattern = "txt")
all.dat = data.frame()

for (fn in 1:length(flist)){
  DAT = read.table(file = paste0("trajectories/", flist[fn]), skip = 11)
  all.dat = rbind(all.dat, DAT)
}

역궤적 자료 읽기 - 합쳐진 자료는 bulgwang.traj

traj = read.table("bulgwang.traj", header = F)
traj$V3 = traj$V3 + 2000
traj$date2 = paste0(traj$V3,traj$V4,traj$V5) # 간단한 날짜 형식으로 변환 (아직 날짜로 인식은 X)
traj$date2 = as.POSIXct(traj$date2, format = "%Y%m%d") + hours(traj$V6 + 9) # 날짜로 인식
traj$date2[1:10] #end point의 시간

##  [1] "2020-12-15 00:00:00 KST" "2020-12-14 23:00:00 KST"
##  [3] "2020-12-14 22:00:00 KST" "2020-12-14 21:00:00 KST"
##  [5] "2020-12-14 20:00:00 KST" "2020-12-14 19:00:00 KST"
##  [7] "2020-12-14 18:00:00 KST" "2020-12-14 17:00:00 KST"
##  [9] "2020-12-14 16:00:00 KST" "2020-12-14 15:00:00 KST"

역궤적 출발시간을 date 열로 붙이기!

traj$date = traj$date2 - hours(traj$V9)
traj$date[1:10]

##  [1] "2020-12-15 KST" "2020-12-15 KST" "2020-12-15 KST" "2020-12-15 KST"
##  [5] "2020-12-15 KST" "2020-12-15 KST" "2020-12-15 KST" "2020-12-15 KST"
##  [9] "2020-12-15 KST" "2020-12-15 KST"

역궤적 열 이름을 붙여주기

colnames(traj)

##  [1] "V1"    "V2"    "V3"    "V4"    "V5"    "V6"    "V7"    "V8"    "V9"   
## [10] "V10"   "V11"   "V12"   "V13"   "date2" "date"

traj = traj[,c(-1,-7,-8)]
colnames(traj)[1:10] = c("receptor", "year", "month", "day", "hour", "hour.inc", "lat", "lon", "height", "pressure")

역궤적 (traj)와 측정자료 (aero) 합치기

aero.traj = merge(traj, aero.meteo[,c(6:16,26)], by = "date", all = T)

역궤적 군집분석 및 PSCF, CWT

trajCluster 함수를 이용 | argument | 설명| |:-| :- | |traj| 앞에서 만든 역궤적 파일| |method | 군집분석 방법 “Euclid” 또는 “Angle”| |n.cluster| 원하는 군집의 수| |type | 구분하고 싶은 자료의 특성 “season”, “year”, “weekday” | |orientation| 지도 회전 반경, c(45, 0, 30)과 같은 형식으로 입력|

clust = trajCluster(traj, n.cluster = 5)

## Warning in trajCluster(traj, n.cluster = 5): Trajectory lengths differ, using
## most common length.

clust = trajCluster(traj, n.cluster = 5, orientation = c(45, 0, 90))

## Warning in trajCluster(traj, n.cluster = 5, orientation = c(45, 0, 90)):
## Trajectory lengths differ, using most common length.

head(clust$data) ## note new variable 'cluster'

## $traj
## # A tibble: 33,142 × 16
## # Groups:   default [1]
##    receptor  year month   day  hour hour.inc   lat   lon height pressure
##       <int> <dbl> <int> <int> <int>    <dbl> <dbl> <dbl>  <dbl>    <dbl>
##  1        1  2020    12    11    15      -72  63.4  112.  4163.     543.
##  2        1  2020    12    11    16      -71  63.2  111.  4173.     544.
##  3        1  2020    12    11    17      -70  63.1  111.  4175      545.
##  4        1  2020    12    11    18      -69  63.0  111.  4170.     548.
##  5        1  2020    12    11    19      -68  62.8  111.  4159.     548.
##  6        1  2020    12    11    20      -67  62.7  111.  4150.     547.
##  7        1  2020    12    11    21      -66  62.5  110.  4145.     546.
##  8        1  2020    12    11    22      -65  62.4  110.  4132      546 
##  9        1  2020    12    11    23      -64  62.2  110.  4100.     548.
## 10        1  2020    12    12     0      -63  62.0  110.  4046.     555.
## # ℹ 33,132 more rows
## # ℹ 6 more variables: date2 <dttm>, date <dttm>, traj_len <int>, default <fct>,
## #   len <int>, cluster <chr>
## 
## $results
## # A tibble: 365 × 8
## # Groups:   cluster, hour.inc [365]
##    cluster hour.inc default            lat   lon date                    n  freq
##    <chr>      <dbl> <fct>            <dbl> <dbl> <dttm>              <int> <dbl>
##  1 C1           -72 15 12월 2020 to…  65.9  108. 2020-12-22 23:56:47  8176  24.7
##  2 C1           -71 15 12월 2020 to…  65.6  108. 2020-12-22 23:56:47  8176  24.7
##  3 C1           -70 15 12월 2020 to…  65.3  108. 2020-12-22 23:56:47  8176  24.7
##  4 C1           -69 15 12월 2020 to…  65.0  109. 2020-12-22 23:56:47  8176  24.7
##  5 C1           -68 15 12월 2020 to…  64.6  109. 2020-12-22 23:56:47  8176  24.7
##  6 C1           -67 15 12월 2020 to…  64.3  109. 2020-12-22 23:56:47  8176  24.7
##  7 C1           -66 15 12월 2020 to…  64.0  109. 2020-12-22 23:56:47  8176  24.7
##  8 C1           -65 15 12월 2020 to…  63.6  110. 2020-12-22 23:56:47  8176  24.7
##  9 C1           -64 15 12월 2020 to…  63.2  110. 2020-12-22 23:56:47  8176  24.7
## 10 C1           -63 15 12월 2020 to…  62.9  110. 2020-12-22 23:56:47  8176  24.7
## # ℹ 355 more rows
## 
## $subsets
## [1] "traj"    "results"

## use different distance matrix calculation, and calculate by season
clust = trajCluster(traj, method = "Angle", type = "year", n.cluster = 4)

## Warning in trajCluster(traj, method = "Angle", type = "year", n.cluster = 4):
## Trajectory lengths differ, using most common length.

trajLevel 함수를 이용 | argument | 설명| |:-| :- | |mydata| 앞에서 만든 역궤적 파일과 농도가 결합된 파일| |pollutant | 적용하고 싶은 오염물질| |smooth| smooth한 그래프를 원하는가? 그러면 T| |statistic | pscf, cwt등 사용| |lon.inc, lat.inc| 위경도 격자 크기 (중요할 수도..)| |percentile| pscf 계산 시 고려될 사례일 기준|

trajLevel(subset(aero.traj, lat > 20 & lat < 50 & lon > 100 & lon < 140),
  pollutant = "PM2.5", statistic = "cwt", smooth = TRUE, orientation = c(45, 0, 90)
)

## Warning: Missing dates detected, removing 12164 lines

trajLevel(subset(aero.traj, lat > 20 & lat < 50 & lon > 100 & lon < 140),
  pollutant = "PM2.5", statistic = "cwt", smooth = TRUE, lon.inc = 0.5, lat.inc = 0.5, orientation = c(45, 0, 90)
)

## Warning: Missing dates detected, removing 12164 lines

trajLevel(subset(aero.traj, lat > 20 & lat < 50 & lon > 100 & lon < 140),
   statistic = "pscf", percentile = 90, smooth = TRUE, lon.inc = 0.5, lat.inc = 0.5, orientation = c(45, 0, 90)
)

## Warning: Missing dates detected, removing 12164 lines

For 승미 (아라온 궤적그리기)

궤적 자료 불러오기sample_data.xlsx

library("readxl")
library("maps")

## 
## 다음의 패키지를 부착합니다: 'maps'

## The following object is masked from 'package:plyr':
## 
##     ozone

track = read_excel(path =  "track.xlsx", sheet = 1, col_names = TRUE)#, fileEncoding = "WINDOWS-1252")

{
  plot(Latitude ~ Longitude, track, type = "l", col = "darkorange", lwd = 2)
  points(Latitude ~ Longitude, track[is.na(track$Benzene) == F,], pch = 15, cex = 0.4, col = "darkred")
  map("world", add = T) # 지도 추가
}

* squash 라이브러리 호출

library("squash") # color bar를 줄 수 있는 라이브러리 호출
(zlim = range(track$Benzene, na.rm = T)) # ORG의 농도 범위 파악

## [1] 0.00 0.67

zlim[2] = 0.3
col.map = makecmap(zlim,  n = 100, colFn = jet, outlier = 'red') # 팔레트 제작 jet, heat, bluered, blueorange, rainbow, greyscale, coolheat
track$spc = track$Benzene
track[track$Benzene >= zlim[2], "spc"] = zlim[2]

{plot(track$Latitude ~ track$Longitude, type = "p", xlab = "", ylab = "", # 위경도 그림
      col = cmap(track$spc, col.map), main = expression("Benzene (ppm)"), pch = 15) # ㅒORG 농도에 따른 심볼 색 추가
mtext(text = expression("Latitude ("*degree*"E)"), line = 2, side = 2) # y축 이름
mtext(text = expression("Longtitude ("*degree*"E)"), line = 2, side = 1) # x축 이름
map("world", add = T) # 지도 추가
vkey(col.map, paste("Conc."),y = 35.0, x = 175, stretch = 0.05) # 레전드 추가
}

상관계수 그림

corrplot 패키지의 corrplot.mixed 함수를 이용

library("corrplot")

## corrplot 0.92 loaded

cor.data = cor(aero.meteo[,6:16], use = "pairwise.complete.obs")
corrplot.mixed(cor.data)

library("PerformanceAnalytics")

## Warning: 패키지 'PerformanceAnalytics'는 R 버전 4.3.2에서 작성되었습니다

## 필요한 패키지를 로딩중입니다: xts

## Warning: 패키지 'xts'는 R 버전 4.3.2에서 작성되었습니다

## 필요한 패키지를 로딩중입니다: zoo

## 
## 다음의 패키지를 부착합니다: 'zoo'

## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric

## 
## 다음의 패키지를 부착합니다: 'PerformanceAnalytics'

## The following object is masked from 'package:graphics':
## 
##     legend

pm.data = aero.meteo[,6:16]
chart.Correlation(pm.data, histogram = TRUE, pch = 19)

## Warning in par(usr): argument 1 does not name a graphical parameter

## Warning in par(usr): argument 1 does not name a graphical parameter

## Warning in par(usr): argument 1 does not name a graphical parameter

## Warning in par(usr): argument 1 does not name a graphical parameter

## Warning in par(usr): argument 1 does not name a graphical parameter

## Warning in par(usr): argument 1 does not name a graphical parameter

## Warning in par(usr): argument 1 does not name a graphical parameter

## Warning in par(usr): argument 1 does not name a graphical parameter

## Warning in par(usr): argument 1 does not name a graphical parameter

## Warning in par(usr): argument 1 does not name a graphical parameter

## Warning in par(usr): argument 1 does not name a graphical parameter

## Warning in par(usr): argument 1 does not name a graphical parameter

## Warning in par(usr): argument 1 does not name a graphical parameter

## Warning in par(usr): argument 1 does not name a graphical parameter

## Warning in par(usr): argument 1 does not name a graphical parameter

## Warning in par(usr): argument 1 does not name a graphical parameter

## Warning in par(usr): argument 1 does not name a graphical parameter

## Warning in par(usr): argument 1 does not name a graphical parameter

## Warning in par(usr): argument 1 does not name a graphical parameter

## Warning in par(usr): argument 1 does not name a graphical parameter

## Warning in par(usr): argument 1 does not name a graphical parameter

## Warning in par(usr): argument 1 does not name a graphical parameter

## Warning in par(usr): argument 1 does not name a graphical parameter

## Warning in par(usr): argument 1 does not name a graphical parameter

## Warning in par(usr): argument 1 does not name a graphical parameter

## Warning in par(usr): argument 1 does not name a graphical parameter

## Warning in par(usr): argument 1 does not name a graphical parameter

## Warning in par(usr): argument 1 does not name a graphical parameter

## Warning in par(usr): argument 1 does not name a graphical parameter

## Warning in par(usr): argument 1 does not name a graphical parameter

## Warning in par(usr): argument 1 does not name a graphical parameter

## Warning in par(usr): argument 1 does not name a graphical parameter

## Warning in par(usr): argument 1 does not name a graphical parameter

## Warning in par(usr): argument 1 does not name a graphical parameter

## Warning in par(usr): argument 1 does not name a graphical parameter

## Warning in par(usr): argument 1 does not name a graphical parameter

## Warning in par(usr): argument 1 does not name a graphical parameter

## Warning in par(usr): argument 1 does not name a graphical parameter

## Warning in par(usr): argument 1 does not name a graphical parameter

## Warning in par(usr): argument 1 does not name a graphical parameter

## Warning in par(usr): argument 1 does not name a graphical parameter

## Warning in par(usr): argument 1 does not name a graphical parameter

## Warning in par(usr): argument 1 does not name a graphical parameter

## Warning in par(usr): argument 1 does not name a graphical parameter

## Warning in par(usr): argument 1 does not name a graphical parameter

## Warning in par(usr): argument 1 does not name a graphical parameter

## Warning in par(usr): argument 1 does not name a graphical parameter

## Warning in par(usr): argument 1 does not name a graphical parameter

## Warning in par(usr): argument 1 does not name a graphical parameter

## Warning in par(usr): argument 1 does not name a graphical parameter

## Warning in par(usr): argument 1 does not name a graphical parameter

## Warning in par(usr): argument 1 does not name a graphical parameter

## Warning in par(usr): argument 1 does not name a graphical parameter

## Warning in par(usr): argument 1 does not name a graphical parameter

## Warning in par(usr): argument 1 does not name a graphical parameter

대기환경학자를 위한 R (실전편 - 역궤적 등)

이화여자대학교 2024

한국외국어대학교 최용주

2024-02-27

1차 정리된 파일 불러와서 자료 가공하기

기본 설정

역궤적 자료 처리 및 기존 자료와 합치기

윈도우에서 리눅스 쓰는 방법!

역궤적 군집분석 및 PSCF, CWT

For 승미 (아라온 궤적그리기)

상관계수 그림