# 측정소-인근 AWS 매핑표 (필요 컬럼: 1~3, 5)
airkorea_station <- read_xlsx("./data/첨부2_측정소_인근_대기관측소(AWS).xlsx")[, c(1:3, 5)]
# 2025년 PM2.5 시간자료 (에어코리아 도시대기, CP949 인코딩)
pm25_korea <- read.csv(
"./data/첨부1_2025_PM2.5 시간자료.txt",
fileEncoding = "CP949", sep = "\t"
)① 동작구 PM2.5 심층분석
2025년 연평균·단기 환경기준 달성 검토 및 3월 황사 유입 사례 분석
1 개요
서울 동작구(측정소 코드 111241)의 2025년 PM2.5 시간자료로 다음을 확인한다.
- 연평균 환경기준(15 ㎍/㎥) 초과 여부
- 단기(일평균 35 ㎍/㎥, 99백분위) 환경기준 초과 여부
- 2025년 3월 황사 유입기의 시간별 농도 ↔︎ 인근 AWS 풍향·풍속 연동
2 데이터 로드
데이터는 2026년 봄학기 대기과학 전공 과정 ‘대기오염’ 과목에서 과제 용도로 제공된 데이터를 활용한다.
측정망의 종류에 따라 도시대기 측정소, 배경대기 측정소 등의 유형이 존재하는데 우리는 도시대기 관측이 목적이므로 이를 대상으로 한다.
airkorea_station %>%
group_by(mang_name) %>%
summarize(count = n(), .groups = "drop")# A tibble: 5 × 2
mang_name count
<chr> <int>
1 교외대기 27
2 국가배경농도(도서) 11
3 도로변대기 64
4 도시대기 530
5 항만 23# 도시대기 관측망으로 예상되는 측정소 개수
length(unique(pm25_korea$AREA_INDEX))[1] 5323 최근 설치 측정소 탐색
풍향/풍속 등 기상데이터와 연계된 분석을 하려면 도시대기 측정소에서 가장 가까운 방재기상 관측소(AWS)를 찾아야 하기 때문에 혹시 최근에 추가된 측정소가 있는지를 확인할 필요가 있다.
간단히 최근에 설치된 곳이라면 데이터 수가 적을 것이라는 점을 이용했다.
pm25_korea %>%
group_by(AREA_INDEX) %>%
summarize(count = n(), .groups = "drop") %>%
arrange(count) %>%
head(5)# A tibble: 5 × 2
AREA_INDEX count
<int> <int>
1 131543 1008
2 221251 2832
3 111121 8760
4 111123 8760
5 111131 87604 동작구 추출 및 일평균 계산
과제가 자기가 사는 곳 근처의 도시대기 측정소를 분석해보라는 것이어서 동작구(111241)를 추출하고, 기상데이터와 조인 시에 문제가 안 되도록 1~24시 기준을 기상데이터와 같이 0~23시 체계에 맞춰 시각을 변환해준다.
pm25_korea <- pm25_korea %>%
mutate(DateTime = ymd_h(YYYYMMDDHH),
hour = format(DateTime, "%Y%m%d%H"))
pm25_dongjak <- pm25_korea %>%
filter(AREA_INDEX == 111241)
# 일평균 (하루 18시간 이상 측정된 날만 유효 일평균으로 인정)
pm25_dongjak$date <- substr(pm25_dongjak$YYYYMMDDHH, 1, 8)
pm25_dongjak_daily_mean <- daily_valid_mean(pm25_dongjak)
head(pm25_dongjak_daily_mean)# A tibble: 6 × 4
date n_hour daily_mean count
<chr> <int> <dbl> <int>
1 20250101 24 21.8 24
2 20250102 24 20.1 24
3 20250103 22 15.4 24
4 20250104 24 19.5 24
5 20250105 24 24.2 24
6 20250106 24 45.3 24환경기준 판정에 쓰는 일평균은 하루 24시간 중 18시간(75%) 이상 측정된 날만 유효한 것으로 인정한다(
setup.R의daily_valid_mean). 여기서도 같은 기준을 적용해 일평균을 산출했다.
5 연평균 환경기준(15 ㎍/㎥) 판정
연평균 기준은 일별 평균치의 연간 평균으로 산출하도록 규정되어 있다.
year_critical <- mean(pm25_dongjak_daily_mean$daily_mean, na.rm = TRUE)
year_critical[1] 18.46442if (year_critical > PM25_ANNUAL_STD) "연평균 환경기준 초과" else "달성"[1] "연평균 환경기준 초과"2025년 1년간의 일평균 농도의 추이도 아래와 같은 코드로 시각화해 볼 수 있다.
pm25_dongjak_daily_mean %>%
mutate(DateTime = ymd(date)) %>%
ggplot() +
geom_line(aes(x = DateTime, y = daily_mean), lwd = 1.0) +
labs(title = "동작구 측정소 일평균 PM2.5 농도 변화 추이",
y = "PM2.5 농도(㎍/㎥)", x = "날짜") +
geom_hline(yintercept = PM25_DAILY_STD, col = "red")
6 단기 환경기준(일평균 35 ㎍/㎥, 99백분위) 초과 여부 판정
우리나라의 단기 환경기준은 단순히 일평균이 기준치보다 낮다고 달성되는 것이 아니다. 추가적으로 일평균의 99백분위수를 기준치(35)와 비교하고 그것이 초과되지 않아야 달성되는 것이라는 점이 환경기준에 명시되어 있다.
# 유효 일평균 개수 기준 99백분위 위치
n_valid <- nrow(pm25_dongjak_daily_mean)
p99_index(n_valid)99%
357 # 상위 농도일 확인
pm25_dongjak_daily_mean %>%
arrange(desc(daily_mean)) %>%
head(4)# A tibble: 4 × 4
date n_hour daily_mean count
<chr> <int> <dbl> <int>
1 20250121 24 89.3 24
2 20250120 24 75.3 24
3 20250122 24 66.5 24
4 20250310 24 63.9 24ascending_order <- order(pm25_dongjak_daily_mean$daily_mean)
critical_value <- pm25_dongjak_daily_mean[ascending_order[p99_index(n_valid)], ]
critical_value# A tibble: 1 × 4
date n_hour daily_mean count
<chr> <int> <dbl> <int>
1 20250310 24 63.9 24if (critical_value$daily_mean > PM25_DAILY_STD) "단기 환경기준 초과" else "달성"[1] "단기 환경기준 초과"숫자로만 보면 감이 안오고 복잡한 것 같아서 시각화를 통해 눈으로 확인할 수 있게 만들어 보았다. (규정에 따르면, 아래 그래프에서 파란색 수평선과 적색 수직선이 만드는 사분면 중 2사분면에 아무것도 없어야 단기 환경기준 달성이 된다.)
pm25_dongjak_daily_mean %>%
mutate(DateTime = ymd(as.character(date))) %>%
ggplot() +
geom_bar(aes(x = reorder(DateTime, daily_mean), y = daily_mean),
lwd = 1.0, stat = "identity") +
scale_y_continuous(limits = c(0, 100), breaks = seq(0, 100, 20)) +
labs(title = "동작구 일평균 PM2.5 농도 (오름차순 정렬)",
y = "PM2.5 농도(㎍/㎥)", x = "날짜") +
theme(axis.text.x = element_blank()) +
geom_hline(yintercept = PM25_DAILY_STD, col = "blue") +
geom_vline(xintercept = p99_index(n_valid), col = "red")
7 2025년 3월 황사 유입 사례
이것도 교수님께서 과제로 내주신 것인데, “최근의 황사일을 찾아서 초미세먼지 농도에 영향을 주는지 확인해 보라”는 것이었다.
어떻게 할까 고민하다, 도시대기 측정소 인근의 기상관측 데이터와 연계한 시각화 방안을 생각해 내었는데, 이 과제가 아니었으면, 만들지 못했을 것이라 생각한다. (결과적으로 기후부 AX 아이디어 공모전에서도 활용했다. 교수님께 감사드립니다!)
배경: 3월 11~15일경 중국 고비사막·네이멍구 발원 대규모 황사가 북서풍을 타고 한반도로 유입. 3월 13일 수도권·서쪽 지역 전국 황사 영향권, 환경부 황사 위기경보 ‘관심’ 단계 발령. 인천·세종·충청 일부는 ‘매우 나쁨’ 수준.
- 황사일 전후의 추세 파악에 필요한 3월 1~15일 구간과 동작구 도시대기 측정소에서 3km 이내에 위치한 AWS(425 남현) 기상데이터를 준비한다.
sand_data <- pm25_dongjak %>%
filter(date >= 20250301 & date <= 20250315)
# 동작구 인근 AWS(425) 기상데이터 — WD="-9"/WS="-9.0" 은 결측
aws_425 <- readRDS("./data/첨부4_AWS_425_동작구인근_기상데이터.RDS")
aws_425_event <- aws_425 %>%
select(STN, YYMMDDHHMI, date, hour, WD, WS) %>%
mutate(WS = ifelse(WD == "0.0", 0, WS),
WS = ifelse(WS == "-9.0", NA, WS),
WD = ifelse(WD == "-9", NA, WD)) %>%
filter(date >= "2025-03-01", date <= "2025-03-15")
# 결측 확인
c(WD_결측 = sum(is.na(aws_425_event$WD)),
WS_결측 = sum(is.na(aws_425_event$WS)))WD_결측 WS_결측
0 0 세 지표(PM2.5·풍향·풍속)를
patchwork라는 패키지로 시각화결과를 수직 결합해 황사 유입 시점(3/13)에서의 초미세먼지 농도와 풍향/풍속의 변화 관계를 한눈에 볼 수 있도록 만들어 보았다.
p1 <- sand_data %>%
mutate(DateTime = ymd_h(as.character(YYYYMMDDHH))) %>%
ggplot() +
geom_line(aes(x = DateTime, y = PM25), lwd = 1.0) +
labs(title = "동작구 PM2.5 시간별 농도", y = "PM2.5 농도(㎍/㎥)", x = NULL) +
geom_vline(xintercept = ymd_h("2025031301"), col = "red", linetype = "dashed") +
geom_vline(xintercept = ymd_h("2025031400"), col = "red", linetype = "dashed") +
annotate("text", x = ymd_h("2025031401"), y = 75,
label = "고농도 황사유입기(3/13)", col = "red", hjust = 0)
p2 <- aws_425_event %>%
mutate(DateTime = ymd_hms(as.character(YYMMDDHHMI))) %>%
ggplot() +
geom_line(aes(x = DateTime, y = as.numeric(WD)), lwd = 1.0) +
labs(title = "남현 AWS 풍향", y = "풍향(도)", x = NULL) +
geom_vline(xintercept = ymd_h("2025031301"), col = "red", linetype = "dashed") +
geom_vline(xintercept = ymd_h("2025031400"), col = "red", linetype = "dashed") +
geom_hline(yintercept = 270, col = "blue") +
annotate("text", x = ymd_h("2025031401"), y = 355,
label = "고농도 황사유입기", col = "red", hjust = 0)
p3 <- aws_425_event %>%
mutate(DateTime = ymd_hms(as.character(YYMMDDHHMI))) %>%
ggplot() +
geom_line(aes(x = DateTime, y = as.numeric(WS)), lwd = 1.0) +
labs(title = "남현 AWS 풍속", y = "풍속(m/s)", x = "측정일시") +
geom_vline(xintercept = ymd_h("2025031301"), col = "red", linetype = "dashed") +
geom_vline(xintercept = ymd_h("2025031400"), col = "red", linetype = "dashed") +
annotate("text", x = ymd_h("2025031401"), y = 4,
label = "고농도 황사유입기", col = "red", hjust = 0)
p1 / p2 / p3
8 일평균 추이 ↔︎ 환경기준 비교 (종합)
마무리로 아래와 같은 두 개의 그래프만 그려보면 전국의 모든 도시대기 측정소의 단기 환경기준 달성여부를 알 수 있게 된다. 그래서 이러한 방법을 활용하여 전국에서 단기 환경기준을 달성한 도시대기 측정소가 몇 개가 되는지를 다음 장에서 검토해보기로 한다.
# 동작구 유효 일평균 (앞서 만든 pm25_dongjak_daily_mean 과 동일 기준)
dongjak_daily <- pm25_dongjak_daily_mean %>%
mutate(DateTime = ymd(as.character(date)))
dongjak_n <- nrow(dongjak_daily)
plot1_1 <- dongjak_daily %>%
ggplot() +
geom_line(aes(x = DateTime, y = daily_mean), lwd = 1.0) +
geom_vline(xintercept = ymd("20250313"), col = "red") +
labs(title = "동작구 일평균 PM2.5 농도 변화 추이", y = "PM2.5 농도(㎍/㎥)", x = "날짜") +
geom_hline(yintercept = PM25_DAILY_STD, col = "blue") +
annotate("text", x = ymd("20250313"), y = 75,
label = "고농도 황사유입일(3/13)", col = "red", hjust = 0)
plot1_2 <- dongjak_daily %>%
ggplot() +
geom_bar(aes(x = reorder(DateTime, daily_mean), y = daily_mean),
lwd = 1.0, stat = "identity") +
scale_y_continuous(limits = c(0, 100), breaks = seq(0, 100, 20)) +
labs(title = "동작구 일평균 PM2.5 농도 (오름차순 정렬)", y = "PM2.5 농도(㎍/㎥)", x = "날짜") +
theme(axis.text.x = element_blank()) +
geom_hline(yintercept = PM25_DAILY_STD, col = "blue") +
geom_vline(xintercept = p99_index(dongjak_n), col = "red")
plot1_1 / plot1_2
9 분석 결과 요약
- 연평균: 동작구 일평균의 연간 산술평균이 15 ㎍/㎥ 기준을 초과 → 연평균 환경기준 미달성
- 단기: 99백분위 일평균이 35 ㎍/㎥를 초과 → 단기 환경기준 미달성
- 3월 13일 황사 유입 구간에서 PM2.5 급등이 인근 AWS 풍향(북서풍 계열)·풍속 변화와 시점상 연동됨