R단계구분도: 2021 서울시 미세먼지 지도 그리기
완벽한 논문을 위한 R 프로그래밍 분석
앞서 2021 서울시 미세먼지 관련 분석을 총 4차례에 거쳐 글을 올렸다.
주제가 흥미롭기도 하고 활용할 수 있는 아이디어도 많았던 것도 있지만, 더 내가 열심히 하고 싶었던 이유는 논문에 들어갈 단계구분도 기능 때문이였다.
비록 논문에 내가 넣어야 할 것은 중국 지도지만... 어찌됐든 코딩 구현은 대동소이할 것이라 믿는다 ㅠ
교재에서는 단계구분도 까지는 설명하지 않았지만, 활용했던 분석 데이터셋과 구글링으로 독학한 지도 그리기로 서울시 행정구별 미세먼지 수준을 제작할 수 있었다...!
첫 지도 제작이라 많은 착오가 있었고, 분명 그대로 따라했는데 각자의 작업 환경이 달라서 에러가 나는 등 여러가지 장애요인이 많았지만... 역시 R 프로그래밍은 다들 문제에 대해 서로 커뮤니케이션 자료가 많아서 금방금방 나만의 코드를 짤 수 있었다.
서울시 단계구분도를 그리는 나같은 초보자들에게 도움이 될 수 있길 희망하며, process를 차근차근 적어내려가도록 하겠다 :)
* 참고한 글들, 다들 감사합니다!
https://hwangknock.tistory.com/6
# 사용한 패키지(중간중간 추가한 것은 본문에 있음)
install.packages("ggmap")
install.packages("ggplot2")
install.packages("raster")
install.packages("rgeos")
install.packages("maptools")
install.packages("rgdal")
총 두 가지 데이터 셋을 합쳐서 단계 구분도를 제작한다고 한다.
: 시각화 대상이 되는 내 데이터 셋 + 지도 그림 데이터 셋 = merge_data -> ggplot에 대입
우선, 내가 입력하고자 하는 데이터셋 가공 부터 설명하겠다.
1. 시각화 대상이 되는 데이터셋(내 자료)
먼저 내가 사용할 자료는, 앞서 2021 서울 미세먼지 분석(최종) 중 "구별 미세먼지 등급 분석" 항목이다.
처음에는 전체 미세먼지 데이터로 단계구분도를 만들었었는데, 아무래도 같은 서울시이고, 극단치(평균 30~40, 극단치 500~)가 꽤 존재해서 전체 지도의 색의 구분이 없어지는 문제가 발생했었다....
그래서, 고민 끝에 구별 미세먼지 등급을 이용해 지도를 만들었으니 이 점 참고바람!
시각화할, 구별 good등급 빈도 데이터 셋구별 미세먼지 등급 비교
: 구별로 미세먼지 등급을 분류한 후에 good의 빈도와 백분율을 나타냄. 이 글에서는 빈도만 이용했음.
seoulair %>%
filter(!is.na(pm10)) %>%
mutate(pm_grade = ifelse(pm10<=30, 'good'
ifelse(pm10<=81, 'normal',
ifelse(pm10<=150, 'bad','worse')))) %>%
group_by(district, pm_grade) %>% # district별로 pm_grade 등급 분류summarise(cnt = n()) %>% # district 별로 빈도 총계 구하기
mutate(total = sum(cnt),
pct = round(cnt/total*100,1)) %>% # district별로 pm_grade등급별 백분율
filter(pm_grade == 'good') %>%
select(district, cnt, pct) %>%
arrange(desc(pct)) %>%
head(5)
# 시각화할 데이터셋에, 행정구별 id 데이터셋 합치기
앞서 말했듯, 우리는 시각화할 내 데이터와 지도 그림 데이터를 합쳐서(merge) 하나의 데이터셋으로 만든 후, ggplot에 대입해 표현하게 된다.
left_join이나 bind_rows를 해봤으면 알겠지만, 하나의 기준을 부여해서 합치게 된다. 이번 서울 지도에서는 뒤의 지도 데이터셋에서 확인할 수 있겠지만 각 행정구별 id를 기준으로 합치게 된다.
그런데, 위의 내가 사용할 데이터셋을 보면알 수 있듯, id 컬럼이 없다.
그래서 따로 id 컬럼을 합쳐줘야하는데, 행정구별 id 데이터는 인터넷에서 금방 찾을 수 있었다.
25개 행정구별 id 데이터우선 너무 복잡하기 때문에, 내 데이터는 p2(첫 시도에 p로 실패해서..)라는 새로운 객체에 데이터셋을 넣어준다.
p2<- seoulair %>%
filter(!is.na(pm10)) %>%
mutate(pm_grade = ifelse(pm10<=30, 'good',
ifelse(pm10<=81, 'normal',
ifelse(pm10<=150, 'bad','worse')))) %>%
group_by(district, pm_grade) %>%
summarise(cnt = n()) %>%
mutate(total = sum(cnt),
pct = round(cnt/total*100,1)) %>%
filter(pm_grade == 'good') %>%
select(district, cnt, pct) %>%
arrange(desc(pct))
그리고, district라는 공통 컬럼을 기준으로 p데이터셋에 구별 id컬럼을 추가시켜줬다.
그런데, id가 integer라서 error가 뜨니, numeric으로 변경 후 join해 준다.
p2$id <- as.numeric(p2$id)
p2<- left_join(p2, id, by = "district")
p2<- p2 %>% relocate(province,district,id,cnt,pct)
id 컬럼이 추가된 p2, 변수 순서는 상관 없다!추가로, 나는 참고했던 글에서 province열을 추가했길래 따라서 추가한건데, 사실 필요없으니 자신의 데이터셋에 따라 활용하면 되겠다.
나처럼 일괄적으로 seoul 이라는 컬럼을 넣고 싶은 분들은 아래 코드를 참고하시길! str_detect라는 함수는 앞으로도 굉장히 유용할듯 하다 :)
install.packages('tidyverse')
library('tidyverse')
p2 <- p2%>%
select(district, cnt, pct) %>%
mutate(province = case_when(
str_detect(district,'구')~'seoul'
))
이제 내가 시각화할 데이터셋 p2까지 완성했으니, 다음은 짝꿍인 지도 데이터 셋을 가공해야한다. 그냥 있는 지도 데이터셋을 합치는 줄 알았더니, 세상은 역시 쉽지 않다...후... :)
2. 지도 그림 데이터 셋
2-1 밑바탕 그림 데이터 셋 만들기
본문 초장에 적은 패키지들을 구동시켰다면 아래 순서대로 따라해주기만 하면 된다. 그런데 나는 뭔가 error가 많이 나서 이것저것 또 다 시도해봤다... 이상해
아래 사이트에서, 2017년 3월 지도그림을 다운받는다고 생각하면 된다. zip을 풀면 하나의 폴더 안에 여러 형식의 파일들이 있는데 서로가 보충해주는 역할을 하기 때문에 고대로~ 두고 하나씩 뽑아서 사용해 줄 것이다.
http://www.gisdeveloper.co.kr/?p=2332
아래 두 방법 중 하나로 map이라는 그림 데이터가 들어간 새로운 객체를 만들어 준다.
왜냐하면 작업환경마다 되는 함수가 있고 안되는 듯하기 때문이다. 구글링을 하다보니 이런 사람들이 꽤 많음.... 참고로 나는 후자로 성공했다! 이때 사용되는 파일은 shp라는 지도를 그려줄 파일, 그러지만 이 파일만 빼서 다른 폴더로 이동시키면 error가 뜨니 유의하자.
map <- readShapePoly('C:/R/Start R/book/OpenData/SIG_201703/TL_SCCO_SIG.shp') #error: readShapePoly is deprecated; use rgdal::readOGR or sf::st_read
map <- readOGR("C:/R/Start R/book/P1_2021서울미세먼지/SIG_201703/TL_SCCO_SIG.shp") # 성공
다 안된다면.. 'sf'라는 패키지 도전 추천...
잘 들어갔는지 확인해보자.
입력한 후에 뭐가 막 주저리주저리 아래와 같이 나와야 성공이다. error는 NA라고 뜸 ㅠㅠ
이제 지도 밑바탕 데이터까지는 만든 셈이다!
따로 어떤 항목이 있나 궁금하신 분들은 아래로 확인해보라~
class(map)
slotNames(map)
map_info <- map@data
head(map_info)
그림 그리는 지리정보 데이터셋 중 data 항목이다2-2 좌표변환하기
처음에는 이게 왜 필요한지 알지 못했다.
하다가 에러도 많이 나서, 위의 map 만 이용해서 p2랑 합쳐 ggplot을 그렸는데도 서울 지도가 만들어져서, 없어도 되나? 했는데... 이 과정을 거치지 않으면 지도 그림만 그리고 그 지도가 위도와 경도를 포함하지 않아서, 추가로 지도에 넣어야하는 행정구명이라던지 단위 같은거를 입력할 수가 없다!
그러니 꼭 필요한 단계이다.
우선 아래를 통해, 현재 지도의 default로 설정된 좌표를 확인한다.
map@proj4string
흠... 사실 봐도 뭔 말인지 모르겠다. 그냥 위에 처럼 확인하고 넘어가자
이제 변환을 해줄건데, 국제표준기준(WGS84)으로 변환하는거라, 아래처럼 변환코드 넣으면 된다.
map <- spTransform(map, CRSobj = CRS('+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs'))
그리고, 아래를 입력하면 좌표가 적절히 들어간 것을 확인할 수 있다.
map@polygons[[1]]@Polygons[[1]]@coords %>% head(n = 10L)
2-3 map의 지도정보 데이터를 data frame으로 바꾸기
계속 말하지만, 우리는 map데이터와 앞서 만든 나의 p2데이터를 합쳐서, 그 하나의 데이터셋을 ggplot에 넣는 것이다.
그런데 우리가 만든 map의 현 상태를 말하자면 좌표가 부여된 뭉텅이 지도 그림이다. 따라서 이걸 세밀한 data frame으로 만들어야, 우리가 이후 자유자재로 값을 부여해서 이쁜 그림을 만들수 있다.
이때 사용하는 함수가, ggplot2의 fortity라는 함수이다. 지도그리기에서 필수함수!
This function turns a map into a data frame that can more easily be plotted with ggplot2.
new_map <- fortify(map, region = 'SIG_CD') # map의 SIG_ID를 ID변수로 지정해준다
View를 통해서 확인해 보면, 아래와 같이 우리의 p2데이터셋과 비슷한 형식으로 바뀌었음을 알 수 있다. 그리고 좌표변환을 해줬기에 long&lat컬럼도 확인 가능하다. 원래 여는데 시간이 오래걸리니 error가 아니다!(내가 이랬음)
data frame으로 바뀐 map데이터셋!2-4 new_map 데이터 셋을 서울로 한정하기
위의 지리정보에는 전국의 모든 시군구가 포함되어 있다. 그러니 id를 중심으로 서울시로만 한정해서 따로 객체를 만든다.
id가 또 integer이니 numeric으로 우선 변환하고 서울시 데이터만 뽑는다. 11740이하까지가 서울시의 행정구별 id이다.
new_map$id <- as.numeric(new_map$id)
seoul_map <- new_map[new_map$id <= 11740,]
3. 시각화 대상이 되는 내 데이터 셋 + 지도 그림 데이터 셋 = merge_data
자 이제 이 두 데이터셋을 합쳐보자!
P2_merge <- merge(seoul_map, p2, by='id')
str(P2_merge)
지도와 내데이터가 합쳐진 p2_merge 지도를 그려보면,
아래와 같이 2021 서울시 행정구 미세먼지 'good'등급 빈도수를 나타낸 단계구분도가 나온다.
map_graph2<- ggplot() + geom_polygon(data = P2_merge, aes(x=long, y=lat, group=group, fill = cnt))
그런데 위의 범례를 보면 알 수 있듯, 값이 커질 수록 색이 연해진다.
내 데이터는 값이 클 수록 색이 진해지는게 더 적절하기 때문에, 높낮이에 따른 색을 부여할 수 있는 코드를 이용해서 다시 그림을 그려본다. 또 대제목과, 축 제목들도 달아주자!
#ggplot 코드: plot에 기본 틀을 넣고 ggplot의 여러 함수로 세부 옵션 조정. scale_fill_gradient()를 통해 숫자가 높으면 진한 색, 숫자가 낮으면 연한 색으로 구분할 수 있도록 변경, theme_bw()로 뒷배경 제거, labs()로 제목 추가, theme()으로 plot의 기본 격자 제거
map_graph2 +
scale_fill_gradient(low = "#C3D7A4", high = "#52854C", space = "Lab", guide = "colourbar") +
theme_bw() +
labs(title = "서울시 미세먼지 '좋음'수준") +
theme(panel.grid.major.x = element_blank(), panel.grid.minor.x = element_blank(), panel.grid.major.y = element_blank(), panel.grid.minor.y = element_blank(), plot.title = element_text(face = "bold", size = 18, hjust = 0.5))
색도 컨셉에 맞게, 구글에서 green 컬러 rgb를 찾아서 원하는 높낮이로 넣어줬다.
여기까지만 해도 되지만, 더 완벽한 지도 작성을 위해 구별 이름을 지정하도록 하면 좋을 것 같다.
3. 서울시 단계구분도에 행정구명 표시하기
행정구명은 geom_text()함수를 위의 그림 코드에 추가하기만 하면 된다. 그러나 이름이 지도에 표시될 때 정확히 구별 중앙에 들어가야 보기가 좋다. 아무 좌표로 넣으면 중구난방으로 옆 행정구 위치까지 침범하는 사태가 발생한다...
따라서 제일 먼저 구별 중심점 좌표 데이터를 구해야 한다.
열심히 구글에서 서울시 행정구 중심점 좌표 csv파일을 검색하던 중, 우연히 아래의 정보를 찾았다!! ㅠㅠ
사실 처음에 25개별로 제일위에 있는 좌표들을 뽑아서 그림그렸었는데, 망해서 중심점 좌표의 중요성을 알았다... 다들 나같은 실수 하지 말고 아래 데이터를 테이블로 끌어와 새 객체로 만들자!
나는 이렇게 위의 정보를 복붙하여 csv파일을 만든후 import하여 새 객체로 만들었다. 역시 첫 번째 시도가 실패해서 객체명에 2가 붙은 점은.. 참고바람 :)
seoul_district2<- read.csv('서울_행정구_좌표.csv', header = T)
이 데이터셋을 geom_text()로 우리가 그린 그림 코드에 추가만 하면 된다.
map_graph2 +
scale_fill_gradient(low = "#C3D7A4", high = "#52854C", space = "Lab", guide = "colourbar") +
theme_bw() +
labs(title = "서울시 미세먼지 '좋음'수준") +
theme(panel.grid.major.x = element_blank(), panel.grid.minor.x = element_blank(), panel.grid.major.y = element_blank(), panel.grid.minor.y = element_blank(), plot.title = element_text(face = "bold", size = 18, hjust = 0.5)) +
geom_text(data = seoul_district2, aes(x = long, y = lat, label= paste(district)))
짜잔! 아래와 같이 한눈에봐도 어느 구가 공기가 가장 좋은 편인지 알 수 있게 되었다~
정확한 수치를 같이 기입하는 방법은 위의 참고한 사이트 중 두 번째 사이트를 통해 시도해보는 것을 추천한다 :)
이로써, 나의 첫번째 지도 프로젝트는 끝이다!
초보자 여러분들 같이 화이팅~
