R 시작
> 작업공간 (디렉토리)을 설정해줌.
setwd("D:/HUFS/Workspace/Rclass/CC_2023//")
> 오늘은 안면도 (태안) 기후변화감시소에서 장기간 측정한 CO2와 CH4
자료 분석 (다른 이름으로 저장)
> R의 기본은 자료 읽기와 쓰기
텍스트 데이터
- 텍스트 데이터는 데이터와 데이터를 구분하기 위해서 공백(blank),
콤마(comma), 탭(tab)와 같은 구분자(separator)를 사용한다. 텍스트
데이터를 R에서 읽어오기 위해서는 read.table() 함수를 사용하며, 그 사용
방법은 다음과 같다.
file |
R에서 읽어올 외부 데이터가 있는 파일의 위치(경로)와
파일명(확장자 포함)을 문자형으로 지정한다. |
header |
논리형으로 TRUE를 지정하면 외부 데이터에 있는 변수명을
R 데이터에서도 동일하게 사용하며, FALSE를 지정하면 외부 데이터의
변수명을 사용하지 않고 R에서 자체적으로 변수명을 지정한다. 그 이름은
“V1”, “V2” 식으로 된다. |
sep |
구분자로 문자형 형태로 지정한다. 공백이면 ” “,
콤마이면”,“, 탭이면”\t” |
stringsAsFactors |
문자형 데이터를 요인(factor)으로 자동으로 변경할 지를
설정한다. FALSE를 지정하면 문자형 데이터를 요인으로 변경하지 않고
문자형으로 읽어들인다. |
na.strings |
데이터에 결측치가 있거나 결측치로 지정할 내용을
문자형으로 지정한다 |
nrow |
읽어들일 행의 갯수를 지정 |
skip |
읽지 않고 넘길 행의 갯수 |
fileEncoding |
encoding 타입 윈도우는 WINDOWS-1252, 리눅스/맥은
UTF-8 |
co2 = read.table(file = "co2_tap_surface-flask_2_3001-9999_monthly.txt", skip = 209, header = F)
head(co2)
## V1 V2 V3 V4 V5 V6 V7 V8 V9 V10 V11 V12 V13 V14 V15 V16 V17
## 1 TAP 1990 11 1 0 0 0 -999 -9 -9 -9 -9 -9 360.79 -999.999 -9 36.73
## 2 TAP 1990 12 1 0 0 0 -999 -9 -9 -9 -9 -9 361.32 -999.999 -9 36.73
## 3 TAP 1991 1 1 0 0 0 -999 -9 -9 -9 -9 -9 362.87 -999.999 -9 36.73
## 4 TAP 1991 2 1 0 0 0 -999 -9 -9 -9 -9 -9 364.76 -999.999 -9 36.73
## 5 TAP 1991 3 1 0 0 0 -999 -9 -9 -9 -9 -9 364.85 -999.999 -9 36.73
## 6 TAP 1991 4 1 0 0 0 -999 -9 -9 -9 -9 -9 364.25 -999.999 -9 36.73
## V18 V19 V20 V21 V22 V23 V24 V25 V26 V27
## 1 126.13 -999.999 20 -999.999 -999.999 -999.999 1 -9 -9 158
## 2 126.13 -999.999 20 -999.999 -999.999 -999.999 1 -9 -9 158
## 3 126.13 -999.999 20 -999.999 -999.999 -999.999 1 -9 -9 158
## 4 126.13 -999.999 20 -999.999 -999.999 -999.999 1 -9 -9 158
## 5 126.13 -999.999 20 -999.999 -999.999 -999.999 1 -9 -9 158
## 6 126.13 -999.999 20 -999.999 -999.999 -999.999 1 -9 -9 158
> 컬럼 이름 변경 colnames()
- site_gaw_id, year, month, day, hour, minute, second, year, month,
day, hour, minute, second, value, value_unc, nvalue, latitude,
longitude, altitude, elevation, intake_height, flask_no, ORG_QCflag,
QCflag, instrument, measurement_method scale
colnames(co2) = c("site_gaw_id", "year", "month", "day", "hour", "minute", "second", "year1", "month1", "day1", "hour1", "minute1", "second1", "value", "value_unc", "nvalue", "latitude", "longitude", "altitude", "elevation", "intake_height", "flask_no", "ORG_QCflag", "QCflag", "instrument", "measurement_method", "scale")
head(co2)
## site_gaw_id year month day hour minute second year1 month1 day1 hour1 minute1
## 1 TAP 1990 11 1 0 0 0 -999 -9 -9 -9 -9
## 2 TAP 1990 12 1 0 0 0 -999 -9 -9 -9 -9
## 3 TAP 1991 1 1 0 0 0 -999 -9 -9 -9 -9
## 4 TAP 1991 2 1 0 0 0 -999 -9 -9 -9 -9
## 5 TAP 1991 3 1 0 0 0 -999 -9 -9 -9 -9
## 6 TAP 1991 4 1 0 0 0 -999 -9 -9 -9 -9
## second1 value value_unc nvalue latitude longitude altitude elevation
## 1 -9 360.79 -999.999 -9 36.73 126.13 -999.999 20
## 2 -9 361.32 -999.999 -9 36.73 126.13 -999.999 20
## 3 -9 362.87 -999.999 -9 36.73 126.13 -999.999 20
## 4 -9 364.76 -999.999 -9 36.73 126.13 -999.999 20
## 5 -9 364.85 -999.999 -9 36.73 126.13 -999.999 20
## 6 -9 364.25 -999.999 -9 36.73 126.13 -999.999 20
## intake_height flask_no ORG_QCflag QCflag instrument measurement_method scale
## 1 -999.999 -999.999 -999.999 1 -9 -9 158
## 2 -999.999 -999.999 -999.999 1 -9 -9 158
## 3 -999.999 -999.999 -999.999 1 -9 -9 158
## 4 -999.999 -999.999 -999.999 1 -9 -9 158
## 5 -999.999 -999.999 -999.999 1 -9 -9 158
## 6 -999.999 -999.999 -999.999 1 -9 -9 158
> 필요한 행만 추출하기 (str 함수로 전체적인 자료 형태를
보기)
str(co2)
## 'data.frame': 372 obs. of 27 variables:
## $ site_gaw_id : chr "TAP" "TAP" "TAP" "TAP" ...
## $ year : int 1990 1990 1991 1991 1991 1991 1991 1991 1991 1991 ...
## $ month : int 11 12 1 2 3 4 5 6 7 8 ...
## $ day : int 1 1 1 1 1 1 1 1 1 1 ...
## $ hour : int 0 0 0 0 0 0 0 0 0 0 ...
## $ minute : int 0 0 0 0 0 0 0 0 0 0 ...
## $ second : int 0 0 0 0 0 0 0 0 0 0 ...
## $ year1 : int -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 ...
## $ month1 : int -9 -9 -9 -9 -9 -9 -9 -9 -9 -9 ...
## $ day1 : int -9 -9 -9 -9 -9 -9 -9 -9 -9 -9 ...
## $ hour1 : int -9 -9 -9 -9 -9 -9 -9 -9 -9 -9 ...
## $ minute1 : int -9 -9 -9 -9 -9 -9 -9 -9 -9 -9 ...
## $ second1 : int -9 -9 -9 -9 -9 -9 -9 -9 -9 -9 ...
## $ value : num 361 361 363 365 365 ...
## $ value_unc : num -1000 -1000 -1000 -1000 -1000 ...
## $ nvalue : int -9 -9 -9 -9 -9 -9 -9 -9 -9 -9 ...
## $ latitude : num 36.7 36.7 36.7 36.7 36.7 ...
## $ longitude : num 126 126 126 126 126 ...
## $ altitude : num -1000 -1000 -1000 -1000 -1000 ...
## $ elevation : int 20 20 20 20 20 20 20 20 20 20 ...
## $ intake_height : num -1000 -1000 -1000 -1000 -1000 ...
## $ flask_no : num -1000 -1000 -1000 -1000 -1000 ...
## $ ORG_QCflag : num -1000 -1000 -1000 -1000 -1000 ...
## $ QCflag : int 1 1 1 1 1 1 1 1 1 1 ...
## $ instrument : int -9 -9 -9 -9 -9 -9 -9 -9 -9 -9 ...
## $ measurement_method: int -9 -9 -9 -9 -9 -9 -9 -9 -9 -9 ...
## $ scale : int 158 158 158 158 158 158 158 158 158 158 ...
co2 = co2[ , c("site_gaw_id", "year", "month", "day", "hour", "minute", "second","value", "QCflag", "scale")]
head(co2)
## site_gaw_id year month day hour minute second value QCflag scale
## 1 TAP 1990 11 1 0 0 0 360.79 1 158
## 2 TAP 1990 12 1 0 0 0 361.32 1 158
## 3 TAP 1991 1 1 0 0 0 362.87 1 158
## 4 TAP 1991 2 1 0 0 0 364.76 1 158
## 5 TAP 1991 3 1 0 0 0 364.85 1 158
## 6 TAP 1991 4 1 0 0 0 364.25 1 158
> CO2 농도 변화 그림 그려보기
plot(co2$value)
plot(co2$value, type = "l")
plot(co2$value, type = "b")
plot(co2$value, type = "l", col = "red")
plot(co2$value, type = "l", col = "red", lty = 2)
plot(co2$value, type = "l", col = "red", lty = 2, lwd = 2)
> 내부 데이터 저장
- 텍스트 파일은 주로 txt (write.table) 또는 csv로 저장함 (write.csv
함수 사용)
x |
저장할 데이터 프레임 |
file |
저장할 파일명 |
append |
기존 파일 뒤에 병합하여 저장하는 옵션, 논리형 (T 또는
F) |
sep |
각 열의 구분자 (seperator) |
row.names |
행 이름 (주로 NA로 지정) |
col.names |
열 이름 (주로 지정하지 않음) |
fileEncoding |
파일 encoding 형식 (주로 지정하지 않음) |
write.csv(file = "tap_co2.csv", co2, row.names = F) # 행번호를 없애기 위해 F를 씀.
CO2와 기상변수 시간자료 처리해보기
> 작업공간 (디렉토리)을 설정해줌.
setwd("D:/HUFS/Workspace/Rclass/CC_2023/")
> 안면도에서 측정한 CO2와 기상변수 다운로드 하기
> CO2 자료 읽기
co2 = read.table(file = "co2_amy_surface-insitu_39_9999-9999_hourly.txt", skip = 214, header = F, na.strings = "-999.999")
head(co2)
## V1 V2 V3 V4 V5 V6 V7 V8 V9 V10 V11 V12 V13 V14 V15 V16 V17 V18
## 1 AMY 1999 1 1 0 0 0 -999 -9 -9 -9 -9 -9 NA NA -9 36.5386 126.33
## 2 AMY 1999 1 1 1 0 0 -999 -9 -9 -9 -9 -9 NA NA -9 36.5386 126.33
## 3 AMY 1999 1 1 2 0 0 -999 -9 -9 -9 -9 -9 NA NA -9 36.5386 126.33
## 4 AMY 1999 1 1 3 0 0 -999 -9 -9 -9 -9 -9 NA NA -9 36.5386 126.33
## 5 AMY 1999 1 1 4 0 0 -999 -9 -9 -9 -9 -9 NA NA -9 36.5386 126.33
## 6 AMY 1999 1 1 5 0 0 -999 -9 -9 -9 -9 -9 NA NA -9 36.5386 126.33
## V19 V20 V21 V22 V23 V24 V25 V26 V27
## 1 108 42 66 NA NA 3 1 9 90
## 2 108 42 66 NA NA 3 1 9 90
## 3 108 42 66 NA NA 3 1 9 90
## 4 108 42 66 NA NA 3 1 9 90
## 5 108 42 66 NA NA 3 1 9 90
## 6 108 42 66 NA NA 3 1 9 90
> 컬럼 이름 변경 colnames()
- site_gaw_id, year, month, day, hour, minute, second, year, month,
day, hour, minute, second, value, value_unc, nvalue, latitude,
longitude, altitude, elevation, intake_height, flask_no, ORG_QCflag,
QCflag, instrument, measurement_method scale
colnames(co2) = c("site_gaw_id", "year", "month", "day", "hour", "minute", "second", "year1", "month1", "day1", "hour1", "minute1", "second1", "value", "value_unc", "nvalue", "latitude", "longitude", "altitude", "elevation", "intake_height", "flask_no", "ORG_QCflag", "QCflag", "instrument", "measurement_method", "scale")
head(co2)
## site_gaw_id year month day hour minute second year1 month1 day1 hour1 minute1
## 1 AMY 1999 1 1 0 0 0 -999 -9 -9 -9 -9
## 2 AMY 1999 1 1 1 0 0 -999 -9 -9 -9 -9
## 3 AMY 1999 1 1 2 0 0 -999 -9 -9 -9 -9
## 4 AMY 1999 1 1 3 0 0 -999 -9 -9 -9 -9
## 5 AMY 1999 1 1 4 0 0 -999 -9 -9 -9 -9
## 6 AMY 1999 1 1 5 0 0 -999 -9 -9 -9 -9
## second1 value value_unc nvalue latitude longitude altitude elevation
## 1 -9 NA NA -9 36.5386 126.33 108 42
## 2 -9 NA NA -9 36.5386 126.33 108 42
## 3 -9 NA NA -9 36.5386 126.33 108 42
## 4 -9 NA NA -9 36.5386 126.33 108 42
## 5 -9 NA NA -9 36.5386 126.33 108 42
## 6 -9 NA NA -9 36.5386 126.33 108 42
## intake_height flask_no ORG_QCflag QCflag instrument measurement_method scale
## 1 66 NA NA 3 1 9 90
## 2 66 NA NA 3 1 9 90
## 3 66 NA NA 3 1 9 90
## 4 66 NA NA 3 1 9 90
## 5 66 NA NA 3 1 9 90
## 6 66 NA NA 3 1 9 90
> 필요한 행만 추출하기 (str 함수로 전체적인 자료 형태를
보기)
str(co2)
## 'data.frame': 204984 obs. of 27 variables:
## $ site_gaw_id : chr "AMY" "AMY" "AMY" "AMY" ...
## $ year : int 1999 1999 1999 1999 1999 1999 1999 1999 1999 1999 ...
## $ month : int 1 1 1 1 1 1 1 1 1 1 ...
## $ day : int 1 1 1 1 1 1 1 1 1 1 ...
## $ hour : int 0 1 2 3 4 5 6 7 8 9 ...
## $ minute : int 0 0 0 0 0 0 0 0 0 0 ...
## $ second : int 0 0 0 0 0 0 0 0 0 0 ...
## $ year1 : int -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 ...
## $ month1 : int -9 -9 -9 -9 -9 -9 -9 -9 -9 -9 ...
## $ day1 : int -9 -9 -9 -9 -9 -9 -9 -9 -9 -9 ...
## $ hour1 : int -9 -9 -9 -9 -9 -9 -9 -9 -9 -9 ...
## $ minute1 : int -9 -9 -9 -9 -9 -9 -9 -9 -9 -9 ...
## $ second1 : int -9 -9 -9 -9 -9 -9 -9 -9 -9 -9 ...
## $ value : num NA NA NA NA NA NA NA NA NA NA ...
## $ value_unc : num NA NA NA NA NA NA NA NA NA NA ...
## $ nvalue : int -9 -9 -9 -9 -9 -9 -9 -9 -9 -9 ...
## $ latitude : num 36.5 36.5 36.5 36.5 36.5 ...
## $ longitude : num 126 126 126 126 126 ...
## $ altitude : int 108 108 108 108 108 108 108 108 108 108 ...
## $ elevation : int 42 42 42 42 42 42 42 42 42 42 ...
## $ intake_height : int 66 66 66 66 66 66 66 66 66 66 ...
## $ flask_no : logi NA NA NA NA NA NA ...
## $ ORG_QCflag : logi NA NA NA NA NA NA ...
## $ QCflag : int 3 3 3 3 3 3 3 3 3 3 ...
## $ instrument : int 1 1 1 1 1 1 1 1 1 1 ...
## $ measurement_method: int 9 9 9 9 9 9 9 9 9 9 ...
## $ scale : int 90 90 90 90 90 90 90 90 90 90 ...
co2 = co2[ , c("site_gaw_id", "year", "month", "day", "hour", "minute", "second","value", "QCflag", "scale")]
head(co2)
## site_gaw_id year month day hour minute second value QCflag scale
## 1 AMY 1999 1 1 0 0 0 NA 3 90
## 2 AMY 1999 1 1 1 0 0 NA 3 90
## 3 AMY 1999 1 1 2 0 0 NA 3 90
## 4 AMY 1999 1 1 3 0 0 NA 3 90
## 5 AMY 1999 1 1 4 0 0 NA 3 90
## 6 AMY 1999 1 1 5 0 0 NA 3 90
paste 함수에 구분자 ’-’를 사용하여 날짜 형식을 만들어줌. 그리고
’Time’라는 변수로 저장
co2$Time = paste(paste(co2$year, co2$month, co2$day, sep = "-"), paste(co2$hour, co2$minute, co2$second, sep = ":"))
’Time’를 1행부터 10행까지 찍어보면,
co2$Time[1:10]
## [1] "1999-1-1 0:0:0" "1999-1-1 1:0:0" "1999-1-1 2:0:0" "1999-1-1 3:0:0"
## [5] "1999-1-1 4:0:0" "1999-1-1 5:0:0" "1999-1-1 6:0:0" "1999-1-1 7:0:0"
## [9] "1999-1-1 8:0:0" "1999-1-1 9:0:0"
우리가 원하는 형식과 다른 것을 확인할 수 있음. 따라서 날짜 형식을
’날짜’로 변경해줘야함. 명령어는 as.POSIXct. ’tz’는 timezone
co2$Time = as.POSIXct(co2$Time, tz = "Asia/Seoul")
co2$Time[1:10]
## [1] "1999-01-01 00:00:00 KST" "1999-01-01 01:00:00 KST"
## [3] "1999-01-01 02:00:00 KST" "1999-01-01 03:00:00 KST"
## [5] "1999-01-01 04:00:00 KST" "1999-01-01 05:00:00 KST"
## [7] "1999-01-01 06:00:00 KST" "1999-01-01 07:00:00 KST"
## [9] "1999-01-01 08:00:00 KST" "1999-01-01 09:00:00 KST"
> 기상자료도 똑같이 처리
met = read.table(file = "co2_amy_surface-insitu_39_9999-9999_hourly_met.txt", skip = 140, header = F, na.strings = c("-9999.9", "-99.9", "-999.9", "-999.999"))
colnames(met) = c("site_gaw_id", "year", "month", "day", "hour", "minute", "second", "wind_direction", "wind_speed", "relative_humidity", "precipitation_amount", "air_pressure", "air_temperature", "dew_point_temperature", "sea_water_temperature", "sea_surface_water_temperature", "sea_water_salinity", "sea_surface_water_salinity", "latitude", "longitude", "altitude", "elevation")
str(met)
## 'data.frame': 181810 obs. of 22 variables:
## $ site_gaw_id : chr "AMY" "AMY" "AMY" "AMY" ...
## $ year : int 1999 1999 1999 1999 1999 1999 1999 1999 1999 1999 ...
## $ month : int 3 3 3 3 3 3 3 3 3 3 ...
## $ day : int 29 29 29 29 29 29 29 29 30 30 ...
## $ hour : int 16 17 18 19 20 21 22 23 0 1 ...
## $ minute : int 0 0 0 0 0 0 0 0 0 0 ...
## $ second : int 0 0 0 0 0 0 0 0 0 0 ...
## $ wind_direction : num NA NA 25.8 25.7 41.7 34.8 36.5 37 47.1 70.3 ...
## $ wind_speed : num NA NA 1.6 0.5 0.3 1.9 1.1 1.4 2 1.1 ...
## $ relative_humidity : logi NA NA NA NA NA NA ...
## $ precipitation_amount : logi NA NA NA NA NA NA ...
## $ air_pressure : logi NA NA NA NA NA NA ...
## $ air_temperature : logi NA NA NA NA NA NA ...
## $ dew_point_temperature : logi NA NA NA NA NA NA ...
## $ sea_water_temperature : logi NA NA NA NA NA NA ...
## $ sea_surface_water_temperature: logi NA NA NA NA NA NA ...
## $ sea_water_salinity : logi NA NA NA NA NA NA ...
## $ sea_surface_water_salinity : logi NA NA NA NA NA NA ...
## $ latitude : num 36.5 36.5 36.5 36.5 36.5 ...
## $ longitude : num 126 126 126 126 126 ...
## $ altitude : logi NA NA NA NA NA NA ...
## $ elevation : int 42 42 42 42 42 42 42 42 42 42 ...
met = met[ , c("site_gaw_id", "year", "month", "day", "hour", "minute", "second","wind_direction", "wind_speed")]
head(co2)
## site_gaw_id year month day hour minute second value QCflag scale
## 1 AMY 1999 1 1 0 0 0 NA 3 90
## 2 AMY 1999 1 1 1 0 0 NA 3 90
## 3 AMY 1999 1 1 2 0 0 NA 3 90
## 4 AMY 1999 1 1 3 0 0 NA 3 90
## 5 AMY 1999 1 1 4 0 0 NA 3 90
## 6 AMY 1999 1 1 5 0 0 NA 3 90
## Time
## 1 1999-01-01 00:00:00
## 2 1999-01-01 01:00:00
## 3 1999-01-01 02:00:00
## 4 1999-01-01 03:00:00
## 5 1999-01-01 04:00:00
## 6 1999-01-01 05:00:00
met$Time = paste(paste(met$year, met$month, met$day, sep = "-"), paste(met$hour, met$minute, met$second, sep = ":"))
met$Time = as.POSIXct(met$Time, tz = "Asia/Seoul")
met$Time[1:10]
## [1] "1999-03-29 16:00:00 KST" "1999-03-29 17:00:00 KST"
## [3] "1999-03-29 18:00:00 KST" "1999-03-29 19:00:00 KST"
## [5] "1999-03-29 20:00:00 KST" "1999-03-29 21:00:00 KST"
## [7] "1999-03-29 22:00:00 KST" "1999-03-29 23:00:00 KST"
## [9] "1999-03-30 00:00:00 KST" "1999-03-30 01:00:00 KST"
> 두 데이터 셋의 병합. (시간을 기준으로): co2와 met는 행의 갯수가
서로 다름
nrow(co2); nrow(met)
## [1] 204984
## [1] 181810
merge() 함수를 이용해서 두 데이터 셋을 합치기.
x, y |
합칠 데이터 프레임 두 개 |
by |
합칠 기준이되는 컬럼 이름. 지금이라면 “Time |
all |
모든 자료 다 합침 |
all.x |
x 데이터 기준 모든 자료 다 합침 |
all.y |
y 데이터 기준 모든 자료 다 합침 |
co2.met = merge(co2, met, by = "Time")
head(co2.met)
## Time site_gaw_id.x year.x month.x day.x hour.x minute.x
## 1 1999-03-29 16:00:00 AMY 1999 3 29 16 0
## 2 1999-03-29 17:00:00 AMY 1999 3 29 17 0
## 3 1999-03-29 18:00:00 AMY 1999 3 29 18 0
## 4 1999-03-29 19:00:00 AMY 1999 3 29 19 0
## 5 1999-03-29 20:00:00 AMY 1999 3 29 20 0
## 6 1999-03-29 21:00:00 AMY 1999 3 29 21 0
## second.x value QCflag scale site_gaw_id.y year.y month.y day.y hour.y
## 1 0 382.01 2 90 AMY 1999 3 29 16
## 2 0 383.47 2 90 AMY 1999 3 29 17
## 3 0 386.09 2 90 AMY 1999 3 29 18
## 4 0 NA 3 90 AMY 1999 3 29 19
## 5 0 384.64 2 90 AMY 1999 3 29 20
## 6 0 385.10 2 90 AMY 1999 3 29 21
## minute.y second.y wind_direction wind_speed
## 1 0 0 NA NA
## 2 0 0 NA NA
## 3 0 0 25.8 1.6
## 4 0 0 25.7 0.5
## 5 0 0 41.7 0.3
## 6 0 0 34.8 1.9
과제: co2.met 자료로 부터 매 시간평균 CO2를 구하고 (0시부터 23시),
그래프로 그려오시오.
필요한 열만 선택해서 붙여보기!
co2.met = merge(co2[,c("Time", "value")], met[,c("Time", "wind_direction", "wind_speed")], by = "Time")
co2.met를 csv 파일로 저장하기!
write.csv(file = "AMY_co2_met.csv", co2.met, row.names = F)
Openair
openair는 대기 자료 분석에 가장 많이 사용되는 패키지로써,
대기환경에서 가장 많이 사용되는 패키지중 하나이다. 특히 복잡한 통계 및
그래프를 짧은 명령어로 충분히 실행되기에 초보자에게도 진입장벽이 높지
않다.
openair 패키지 설치: openair가 처음이라면 아래 명령어를 사용하여
설치할 수 있다.
install.packages("openair")
패키지 불러오기: 기존에 설치가 되었거나, 설치가 완료되었다면, 아래
명령어로 lubridate 패키지와 함께 불러온다.
library("lubridate")
library("openair")
자료 불러오기: csv 파일이기에 read.csv 명령어를 이용하여 앞서 저장한
자료 (AMY_co2_met.csv)를 불러들인다. 우선 working directory
설정부터…
setwd("D:/HUFS/Workspace/Rclass/CC_2023/")
co2.met = read.csv(file = "AMY_co2_met.csv", header = T)
시간형태로 읽어들이기! as.POSIXct 명령어를 사용 (!! 주의할 점,
openair는 date를 시간으로 인식함!)
co2.met$date = as.POSIXct(co2.met$Time, tz = "Asia/Seoul") #연, 월, 일을 인식시켜줌으로써 손쉽게 시간형태로 변경
summary plot
- 시계열 자료의 전체적인 경향을 파악하기 좋은 그래프. 자료 및 결측치의
분포, 평균, 중간값 등에 대한 정보를 제공함함
summaryPlot(co2.met) # 안면의 summary plot
summaryPlot(co2.met, percentile = 0.95) # 상위 5% 값은 제거
summaryPlot(co2.met, na.len = 10) # 적어도 10개이상 연속적인 결측치를 보여줌
summaryPlot(co2.met, col.data = "green") # 자료 색깔을 녹색으로
summaryPlot(co2.met, col.mis = "yellow") # 결측치를 노란색으로
summaryPlot(co2.met, col.dens = "black") # density plot선을 검은색으로
summaryPlot(co2.met[, c(5,2,4)]) # 컬럼 2번째와 5-7번째만 보여주기
summaryPlot(subset(co2.met, select = c(date, value, wind_speed))) # 앞선 명령어의 대안
windrose
- 풍향, 풍속에 대한 정보를 제공 (!! 주의할 점, openair는 풍향, 풍속을
wd, ws로 인식함)
colnames(co2.met)[c(2, 3,4)] = c("co2","wd", "ws") # 풍향, 풍속의 컬럼 이름을 wd, ws로 교체. (교체 안할 시 일일이 지정해야하는 번거로움..)
windRose(co2.met) # 기본 바랑장미
windRose(co2.met, type = "month", layout = c(4, 3)) # 월별 바랑장미
windRose(co2.met[is.na(co2.met$co2) == F,], type = "co2", layout = c(4, 1)) # PM2.5 구간에 따른 바람장비
### pollution rose
pollutionRose(co2.met, pollutant = "co2") # PM2.5 농도 구간별 풍향의 빈도 수
pollutionRose(co2.met, pollutant = "co2", normalise = TRUE) # normalized된 PM2.5 농도 구간별 풍향의 빈도 수
Polar plot
- 풍향/풍속에 따른 오염물질의 평균 농도를 보여줌
polarPlot(co2.met, pollutant = "co2")
polarPlot(co2.met, pollutant = "co2", statistic = "cpf", percentile = 75) # Conditional Probability Function
polarPlot(co2.met, pollutant = "co2", statistic = "cpf", percentile = c(0,10)) # Conditional Probability Function
polarPlot(co2.met, pollutant = "co2", statistic = "cpf", percentile = c(50,60)) # Conditional Probability Function
polarPlot(co2.met, pollutant = "co2", type = "daylight") # ”year”, ”hour”, ”month”, ”season”, ”weekday”, ”weekend”, ”monthyear”, ”daylight”
Calendar plot
- 달력에 원하는 오염물질의 농도등을 표시하는 그림. 보다 직관적이기에
발표자료에 많이 사용됨.
st = as.POSIXct("2018-12-01") # 분석 시작 날짜
et = as.POSIXct("2019-03-31") # 분석 종료 날짜
calendarPlot(co2.met[co2.met$date >= st & co2.met$date <= et,], pollutant = "co2") # 용산 PM2.5 농도 기준, 표시
calendarPlot(co2.met[co2.met$date >= st + years(1) & co2.met$date <= et + years(1),], pollutant = "co2") # 분석 기간 기준을 1년 후로 미룸. 용산 PM2.5 농도 기준, 표시
calendarPlot(co2.met[co2.met$date >= st & co2.met$date <= et,], pollutant = "co2", annotate = "ws") # 관악 PM2.5 기준 농도가 아닌, 풍향/풍속 표시
Trend 분석 (smoothTrend, TheilSen)
smoothTrend(co2.met, pollutant = "co2", deseason = TRUE, simulate = F, ylab = "concentration (ppm)",
main = "monthly mean deseasonalised CO2 (bootstrap uncertainties)") # smooth trend를 계산. GAM 모델을 사용하여 smooth를 수행
## Warning in StructTS(myts): possible convergence problem: 'optim' gave code = 52
## and message 'ERROR: ABNORMAL_TERMINATION_IN_LNSRCH'
smoothTrend(co2.met, pollutant = "co2", deseason = TRUE, simulate = F, ylab = "concentration (ppm)",
main = "monthly mean deseasonalised CO2 (bootstrap uncertainties)", type = "wd") # 풍향을 기준으로 분류
## Warning: 5744 missing wind direction line(s) removed
## Warning: possible convergence problem: 'optim' gave code = 52 and message
## 'ERROR: ABNORMAL_TERMINATION_IN_LNSRCH'
## Warning: possible convergence problem: 'optim' gave code = 52 and message
## 'ERROR: ABNORMAL_TERMINATION_IN_LNSRCH'
TheilSen(co2.met, pollutant = "co2", deseason = TRUE, ylab = "CO2", main = "Theil-Sen Trend") # Theil-sen 기울기 계산산
## Warning in StructTS(myts): possible convergence problem: 'optim' gave code = 52
## and message 'ERROR: ABNORMAL_TERMINATION_IN_LNSRCH'
## [1] "Taking bootstrap samples. Please wait."
TheilSen(co2.met, pollutant = "co2", deseason = TRUE, ylab = "CO2", main = "Theil-Sen Trend", type = "season") # 계절절을 기준으로 분류
## Warning in StructTS(myts): possible convergence problem: 'optim' gave code = 52
## and message 'ERROR: ABNORMAL_TERMINATION_IN_LNSRCH'
## [1] "Taking bootstrap samples. Please wait."
## [1] "Taking bootstrap samples. Please wait."
## Warning in StructTS(myts): possible convergence problem: 'optim' gave code = 52
## and message 'ERROR: ABNORMAL_TERMINATION_IN_LNSRCH'
## [1] "Taking bootstrap samples. Please wait."
## [1] "Taking bootstrap samples. Please wait."
### timeVariation 분석
# 기본 그림
timeVariation(co2.met, pollutant = "co2")
# 여러 항목들의 시간변화 비교 (절대값이 다르기에 nomarlized 적용)
timeVariation(co2.met, pollutant = c("co2", "wd", "ws"), normalise = TRUE)
# factor (예: 기준)을 적용하여 전,후 차이를 비교. 아래는 코로나 전, 후의 NO2의 변화
co2.met[co2.met$date < as.POSIXct("2010-01-01"), "millenium"] = "Before 2010" # 코로나 전 지정
co2.met[co2.met$date >= as.POSIXct("2010-01-01"), "millenium"] = "After 2010" # 코로나 후 지정
timeVariation(co2.met, pollutant = "co2", group = "millenium", difference = TRUE)
timeVariation(co2.met, pollutant = "co2", group = "millenium", difference = TRUE, ylim = c(350,450))
그 외 유용한 함수
co2.met = rollingMean(co2.met, pollutant = "co2", hours = 8, new.name = "co2.8", data.thresh = 75) # 8시간 이동 평균
trendLevel(co2.met, pollutant = "co2", cols = "jet", y = "hour")