Time series 그래프
필요한 라이브러리 설치 (일회성)
install.packages("readxl") # 엑셀파일을 읽는 라이브러리 설치
install.packages("lubridate") # 엑셀파일을 읽는 라이브러리 설치
라이브러리 불러오기
library("readxl") # 엑셀파일을 읽는 라이브러리 설치
library("lubridate") # 엑셀파일을 읽는 라이브러리 설치
##
## 다음의 패키지를 부착합니다: 'lubridate'
## The following objects are masked from 'package:base':
##
## date, intersect, setdiff, union
library("openair") # 엑셀파일을 읽는 라이브러리 설치
작업 디렉토리 지정
setwd("D:/HUFS/Workspace/Rclass/NIER_220918/") # 설정
getwd() #확인
## [1] "D:/HUFS/Workspace/Rclass/NIER_220918"
측정자료 (excel 파일) 불러오기
SH = read_xlsx(path = "시화산단AQMS.xlsx", sheet = 1, skip = 5) # 시화산단 자료
## New names:
## * `` -> ...1
## * 등급 -> 등급...2
## * `1시간` -> `1시간...3`
## * 등급 -> 등급...4
## * `1시간` -> `1시간...5`
## * ...
WS = read_xlsx(path = "원시동AQMS.xlsx", sheet = 1, skip = 5) # 원시동 자료
## New names:
## * `` -> ...1
## * 등급 -> 등급...2
## * `1시간` -> `1시간...3`
## * 등급 -> 등급...4
## * `1시간` -> `1시간...5`
## * ...
AWS = read_xlsx(path = "안산(545) AWS.xlsx", sheet = 1, skip = 1) # AWS 자료
불러온 자료의 컬럼 이름 확인 및 필요한 컬럼만 추출
colnames(SH)
## [1] "...1" "등급...2" "1시간...3" "등급...4" "1시간...5"
## [6] "등급...6" "1시간...7" "등급...8" "1시간...9" "등급...10"
## [11] "1시간...11" "등급...12" "1시간...13"
colnames(SH[,c(seq(1,13,2))])
## [1] "...1" "1시간...3" "1시간...5" "1시간...7" "1시간...9"
## [6] "1시간...11" "1시간...13"
SH = SH[,c(seq(1,13,2))]
WS = WS[,c(seq(1,13,2))]
colnames(SH) = colnames(WS) = c("date", "PM10", "PM2.5", "O3", "NO2", "CO", "SO2")
자료의 유형 파악
str(SH)
## tibble [600 x 7] (S3: tbl_df/tbl/data.frame)
## $ date : POSIXct[1:600], format: "2022-07-26 01:00:00" "2022-07-26 02:00:00" ...
## $ PM10 : chr [1:600] "22" "16" "27" "25" ...
## $ PM2.5: chr [1:600] "10" "9" "16" "17" ...
## $ O3 : chr [1:600] "0.020" "0.021" "0.010" "0.002" ...
## $ NO2 : chr [1:600] "0.017" "0.017" "0.028" "0.033" ...
## $ CO : chr [1:600] "0.4" "0.4" "0.4" "0.5" ...
## $ SO2 : chr [1:600] "0.001" "0.001" "0.001" "0.002" ...
모두 ’문자형’이기에 숫자형으로 바꿔줌
for (i in 2:ncol(SH)){
SH[,i] = as.numeric(unlist(SH[,i]))
WS[,i] = as.numeric(unlist(WS[,i]))
}
AWS 자료도 컬럼 이름을 수정
colnames(AWS)
## [1] "지점" "지점명" "일시" "기온(°C)"
## [5] "풍향(deg)" "풍속(m/s)" "강수량(mm)" "현지기압(hPa)"
## [9] "해면기압(hPa)" "습도(%)"
colnames(AWS) = c("site", "name", "date", "temp", "wd", "ws", "prec", "pres", "s.pres", "rh")
AWS 자료도 컬럼 이름을 수정ㅜ
SH = merge(SH, AWS, by = "date", all.x = T)
WS = merge(WS, AWS, by = "date", all.x = T)
summaryPlot(SH[,1:7])
## date1 date2 PM10 PM2.5 O3 NO2 CO SO2
## "POSIXct" "POSIXt" "numeric" "numeric" "numeric" "numeric" "numeric" "numeric"
plot(SH$PM10 ~ SH$date) # (1): 일반적인 점 그래프
plot(SH$PM10 ~ SH$date, type = "l") # (2): 일반적인 선 그래프
plot(SH$PM10 ~ SH$date, type = "l", ylim = c(0,50)) # (3): (2) + y축 범위 변경
plot(SH$PM10 ~ SH$date, type = "l", ylim = c(0,50), col = "darkorange") # (4): (3) + 색깔 변경
plot(SH$PM10 ~ SH$date, type = "l", ylim = c(0,50), col = "darkorange", xlab = "Time", ylab = "PM10") # (5): (4) + 이름 변경
plot(SH$PM10 ~ SH$date, type = "l", ylim = c(0,50), col = "darkorange", xlab = "Time", ylab = expression("PM"[10]*" ("*mu*"g m"^-3*")")) # (6): (5) + 특수문자 및 아래첨자
{
plot(SH$PM10 ~ SH$date, type = "l", ylim = c(0, 70), col = "darkorange", xlab = "Time", ylab ="")
points(SH$PM10 ~ SH$date, col = "darkorange", pch = 22, cex = 0.7)
mtext(text = expression("PM"[10]*" ("*mu*"g m"^-3*")"), line = 2, side = 2) # (6): (5) + 특수문자 및 아래첨자 간격 줄이기기
points(WS$PM10 ~ WS$date, col = "blue", pch = 21, cex = 0.7) # (8): (7) + SO4 자료 추가
lines(WS$PM10 ~ WS$date, col = "blue")
legend("topright", legend = c("시화공단", "원시"), col = c("darkorange", "blue"), pch = 22, cex = 0.7) # (7): (6) + 심볼 추가 + 레전드 추가
}
기타 그래프
Bar 그래프
DAT = SH[SH$date >= as.POSIXct("2022-08-15", tz = "UTC") & SH$date < as.POSIXct("2022-08-18", tz = "UTC"),]
DAT$PM10_25 = DAT$PM10 - DAT$PM2.5
barplot(t(DAT[,c("PM10_25", "PM2.5")]), col = c("darkorange", "skyblue"), beside = F, names.arg = as.Date(DAT$date))
바람장미 그래프
polarplot
Scatterplot
기본 상호작용하는 지도 생성하기
패키지 설치
install.packages("leaflet")
패키지 불러오기
library("leaflet")
library("geosphere")
측정소 위경도 및 AWS와 거리 계산
lon_nier = c(126.7272389, 126.7885028 , 126.8418193)
lat_nier = c(37.3351039, 37.3052549, 37.2802273)
loc = data.frame(group = c("A", "B", "A", "B"),
long = c(lon_nier[c(1,3)], lon_nier[c(2,3)]),
lat = c(lat_nier[c(1,3)], lat_nier[c(2,3)]))
dist.sh = distm(x = c(lon_nier[1], lat_nier[1]), y = c(lon_nier[3], lat_nier[3]))/1000
dist.ws = distm(x = c(lon_nier[2], lat_nier[2]), y = c(lon_nier[3], lat_nier[3]))/1000
{
leaflet() %>%
addTiles() %>%
setView( lng= mean(lon_nier), lat = mean(lat_nier), zoom = 13) %>%
addProviderTiles("Esri.WorldImagery") %>%
addMarkers(lng = 126.7272389,
lat = 37.3351039,
popup = paste("시화공단 </br>도시대기측정망</br> AWS 기준", round(dist.sh, 1), "km")) %>%
addMarkers(lng = 126.7885028,
lat = 37.3052549,
popup = paste("원시동 </br>도시대기측정망 </br> AWS 기준", round(dist.ws, 1), "km")) %>%
addMarkers(lng = 126.8418193,
lat = 37.2802273,
popup = "AWS 측정소</br>") %>%
addPolylines(data = loc, lng = ~long[1:2], lat = ~lat[1:2], col = "red") %>%
addPolylines(data = loc, lng = ~long[3:4], lat = ~lat[3:4], col = "blue")
}
지리학적 분석
무인드론에 설치된 OPC (optical particle couneter) 자료를 분석해보자
필요한 라이브러리 설치 (일회성)
install.packages("readxl") # 엑셀파일을 읽는 라이브러리 설치
install.packages("maps") # google map 을 불러오는 라이브러리 설치
install.packages("squash") # color bar를 줄 수 있는 라이브러리 설치
install.packages("plyr") # 다소 복잡한 통계 함수를 위한 라이브러리 설치
필요한 라이브러리 호출
library("readxl")
library("ggmap")
## 필요한 패키지를 로딩중입니다: ggplot2
## Google's Terms of Service: https://cloud.google.com/maps-platform/terms/.
## Please cite ggmap if you use it! See citation("ggmap") for details.
library("squash")
library("plyr")
작업 디렉토리 지정
setwd("D:/HUFS/Workspace/Rclass/NIER_220918/") # 설정
getwd() #확인
## [1] "D:/HUFS/Workspace/Rclass/NIER_220918"
측정자료 (excel 파일) 불러오기
DAT = read_xlsx(path = "opc_data.xlsx", sheet = 1, skip = 1)
불러온 자료 확인
head(DAT)
## # A tibble: 6 x 14
## 일자_텍스트 시간_텍스트 코드 지역 OPC_자료사용 일자_오전오후_~ TSP PM10
## <chr> <chr> <chr> <chr> <chr> <chr> <dbl> <dbl>
## 1 2022-07-26 9:55:22 2022-~ 1구역 X 2022-07-26 오전~ 50.6 49.9
## 2 2022-07-26 9:55:28 2022-~ 1구역 X 2022-07-26 오전~ 68.5 61.3
## 3 2022-07-26 9:55:34 2022-~ 1구역 X 2022-07-26 오전~ 53.1 52.7
## 4 2022-07-26 9:55:40 2022-~ 1구역 X 2022-07-26 오전~ 66.6 48.7
## 5 2022-07-26 9:55:46 2022-~ 1구역 X 2022-07-26 오전~ 86.3 67.1
## 6 2022-07-26 9:55:52 2022-~ 1구역 X 2022-07-26 오전~ 65.3 56.7
## # ... with 6 more variables: PM4 <dbl>, PM25 <dbl>, PM1 <dbl>,
## # OPC_Latitude <dbl>, OPC_Longitude <dbl>, OPC_Height <dbl>
자료의 컬럼별 통계값 확인
summary(DAT)
## 일자_텍스트 시간_텍스트 코드 지역
## Length:45402 Length:45402 Length:45402 Length:45402
## Class :character Class :character Class :character Class :character
## Mode :character Mode :character Mode :character Mode :character
##
##
##
##
## OPC_자료사용 일자_오전오후_시간 TSP PM10
## Length:45402 Length:45402 Min. : 5.60 Min. : 5.60
## Class :character Class :character 1st Qu.: 19.20 1st Qu.: 19.20
## Mode :character Mode :character Median : 25.95 Median : 25.70
## Mean : 36.97 Mean : 28.38
## 3rd Qu.: 37.23 3rd Qu.: 36.10
## Max. :4477.70 Max. :1135.10
## NA's :36862 NA's :36862
## PM4 PM25 PM1 OPC_Latitude
## Min. : 5.60 Min. : 5.50 Min. : 4.60 Min. :37.29
## 1st Qu.: 18.80 1st Qu.: 17.70 1st Qu.: 15.00 1st Qu.:37.31
## Median : 25.10 Median : 24.00 Median : 21.90 Median :37.33
## Mean : 27.12 Mean : 25.36 Mean : 22.43 Mean :37.33
## 3rd Qu.: 34.40 3rd Qu.: 31.90 3rd Qu.: 28.40 3rd Qu.:37.34
## Max. :1123.20 Max. :963.60 Max. :308.70 Max. :37.35
## NA's :36862 NA's :36862 NA's :36862 NA's :37258
## OPC_Longitude OPC_Height
## Min. :126.7 Min. :-19.00
## 1st Qu.:126.7 1st Qu.: 37.00
## Median :126.7 Median :105.00
## Mean :126.7 Mean : 90.92
## 3rd Qu.:126.8 3rd Qu.:126.00
## Max. :126.8 Max. :217.00
## NA's :37258 NA's :37258
PM2.5 자료릐 NA값이 아닌 자료만을 취합함 (엑셀 필터 기능)
DAT = DAT[is.na(DAT$PM25) == F, ]
단순한 그래프화
plot(DAT$OPC_Latitude ~ DAT$OPC_Longitude) # (1) y ~ x 형태의 일반적인 점 그래프
plot(DAT$OPC_Latitude ~ DAT$OPC_Longitude, xlab = "Longitude", ylab = "Latitude") # (2): (1) + 축 추가
plot(DAT$OPC_Latitude ~ DAT$OPC_Longitude, xlab = "Longitude", ylab = "Latitude", ylim = c(37.29, 37.36)) # (3): (2) + y축 범위 변경
plot(DAT$OPC_Latitude ~ DAT$OPC_Longitude, xlab = "Longitude", ylab = "Latitude", ylim = c(37.29, 37.36), col = "darkorange") # (4): (3) + 색깔 변경
plot(DAT$OPC_Latitude ~ DAT$OPC_Longitude, ylim = c(37.29, 37.36), col = "darkorange", ylab = expression("Latitude ("*degree*"N)"), xlab = expression("Longitude ("*degree*"E)")) # (5): (4) + 특수부호
단순한 그래프화 (축 이름 별도 지정)
{
plot(DAT$OPC_Latitude ~ DAT$OPC_Longitude, ylim = c(37.29, 37.36), col = "darkorange", ylab = "", xlab = "")
mtext(text = expression("Latitude ("*degree*"N)"), line = 2, side = 2) # y축 추가
mtext(text = expression("Longitude ("*degree*"E)"), line = 2, side = 1) # x축 추가
}
PM2.5 농도에 따른 color bar
PM2.5의 농도 범위 파악
(zlim = range(DAT$PM25, na.rm = T)) # 결과를 zlim에 저장
## [1] 5.5 963.6
최고농도의 범위를 90퍼센타일 값으로 변경
(zlim[2] = round_any(quantile(DAT$PM25, 0.90, na.rm = T), accuracy = 1, f = trunc))
## 90%
## 39
col.map = makecmap(zlim, n = 100, colFn = jet, outlier = 'red') # 팔레트 제작 jet, heat, bluered, blueorange, rainbow, greyscale, coolheat
DAT$PM25.1 = DAT$PM25
DAT[DAT$PM25 > zlim[2] & is.na(DAT$PM25) == F, "PM25.1"] = zlim[2]
{
plot(DAT$OPC_Latitude ~ DAT$OPC_Longitude, type = "p", xlab = "", ylab = "", # 위경도 그림
col = cmap(DAT$PM25.1, col.map), main = expression("PM"[2.5]*" ("*mu*"g m"^-3*")")) # ㅒORG 농도에 따른 심볼 색 추가
mtext(text = expression("Latitude ("*degree*"E)"), line = 2, side = 2) # y축 이름
mtext(text = expression("Longtitude ("*degree*"E)"), line = 2, side = 1) # x축 이름
vkey(col.map, expression("PM"[2.5]),y = 37.3, x = 126.81, stretch = 0.2) # 레전드 추가
}
구글 지도 위에 표시
key = "AIzaSyDGUtMAIyCV7ejSian3e3M3NgFTK4rr8IU" #본인이 할당 받은 구글 키를 입력 (참고: https://duopix.co.kr/google-map-key/)
register_google(key)
loc = c(mean(DAT$OPC_Longitude, na.rm = T), mean(DAT$OPC_Latitude, na.rm = T)) # 비행 경로의 평균 위경도를 구글 지도의 중심점으로 잡음. MA가 있기에 na.rm = T라는 명령어를 넣어줌
qmap( loc , zoom = 12, maptype = "satellite") + # 중심점을 기준으로 zoom 값으로 지도를 확대 (큰 수) 혹은 축소 (작은 수). maptype은 satellite, terrain, hybrid, roadmap가 있음
geom_point(aes(x = OPC_Longitude, y = OPC_Latitude), data = DAT, col = cmap(DAT$PM25.1, col.map))
qmap( loc , zoom = 12, maptype = "terrain") +
geom_point(aes(x = OPC_Longitude, y = OPC_Latitude), data = DAT, col = cmap(DAT$PM25.1, col.map))
qmap( loc , zoom = 12, maptype = "hybrid") + # hybrid
geom_point(aes(x = OPC_Longitude, y = OPC_Latitude), data = DAT, col = cmap(DAT$PM25.1, col.map))
qmap( loc , zoom = 12, maptype = "roadmap") + # roadmap
geom_point(aes(x = OPC_Longitude, y = OPC_Latitude), data = DAT, col = cmap(DAT$PM25.1, col.map))
3차원 그래프 완성
library("rgl") # OPEN GL (open graphics library) 2D, 3D 렌더링을 해주는 라이브러리
par3d(windowRect = c(10, 10, 800, 800)) # 3D 윈도우 생성
plot3d(DAT$OPC_Longitude, DAT$OPC_Latitude, DAT$OPC_Height, col = cmap(DAT$PM25.1, col.map),
size = 2, type = 'p', xlab = "Lon", ylab = "Lat", zlab = "Alt", zlim = c(0,300))
movie3d(spin3d(axis = c(0, 0, 0.5), rpm = 5), duration = 20, dir = getwd())