서울시 월평균 지도위에 표출하기

setwd("D:/HUFS/Workspace/Rclass/CC_2023/")
Sys.getlocale()
## [1] "LC_COLLATE=Korean_Korea.utf8;LC_CTYPE=Korean_Korea.utf8;LC_MONETARY=Korean_Korea.utf8;LC_NUMERIC=C;LC_TIME=Korean_Korea.utf8"
Sys.setlocale("LC_ALL","korean")
## Warning in Sys.setlocale("LC_ALL", "korean"): using locale code page other than
## 65001 ("UTF-8") may cause problems
## [1] "LC_COLLATE=Korean_Korea.949;LC_CTYPE=Korean_Korea.949;LC_MONETARY=Korean_Korea.949;LC_NUMERIC=C;LC_TIME=Korean_Korea.949"

서울지역 측정소 월평균 불러오기 seoul_mm.csv 다운로드

m.seoul_latlon = read.csv(file = "seoul_mm.csv")

서울 shapefile 가져오기 서울 shp 다운로드

install.packages("rgdal")

측정소 위치를 지도위에 표출

필요한 패키지 부르기

library("maps")
library("rgdal")
## 필요한 패키지를 로딩중입니다: sp
## Please note that rgdal will be retired during 2023,
## plan transition to sf/stars/terra functions using GDAL and PROJ
## at your earliest convenience.
## See https://r-spatial.org/r/2022/04/12/evolution.html and https://github.com/r-spatial/evolution
## rgdal: version: 1.6-5, (SVN revision 1199)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 3.5.2, released 2022/09/02
## Path to GDAL shared files: C:/Users/YChoi/AppData/Local/R/win-library/4.2/rgdal/gdal
## GDAL binary built with GEOS: TRUE 
## Loaded PROJ runtime: Rel. 8.2.1, January 1st, 2022, [PJ_VERSION: 821]
## Path to PROJ shared files: C:/Users/YChoi/AppData/Local/R/win-library/4.2/rgdal/proj
## PROJ CDN enabled: FALSE
## Linking to sp version:1.6-0
## To mute warnings of possible GDAL/OSR exportToProj4() degradation,
## use options("rgdal_show_exportToProj4_warnings"="none") before loading sp or rgdal.

shape file 읽기

Seoul = readOGR("Seoul.shp")
## Warning: OGR support is provided by the sf and terra packages among others
## Warning: OGR support is provided by the sf and terra packages among others
## Warning: OGR support is provided by the sf and terra packages among others
## Warning: OGR support is provided by the sf and terra packages among others
## Warning: OGR support is provided by the sf and terra packages among others
## Warning: OGR support is provided by the sf and terra packages among others
## Warning: OGR support is provided by the sf and terra packages among others
## OGR data source with driver: ESRI Shapefile 
## Source: "D:\HUFS\Workspace\Rclass\CC_2023\Seoul.shp", layer: "Seoul"
## with 25 features
## It has 3 fields
Seoul = spTransform(Seoul, CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"))

측정소 위치 그리기

{plot(m.seoul_latlon$lat ~ m.seoul_latlon$lon, xlab = "", ylab = "", pch = 21,
      xlim = c(126.75, 127.2), ylim = c(37.42, 37.72)
      ) # 위경도에 따른 심볼 그림
mtext(text = expression("Latitude ("*degree*"E)"), line = 2, side = 2) # y축 이름
mtext(text = expression("Longtitude ("*degree*"E)"), line = 2, side = 1) # x축 이름
map(Seoul, add = T) # 지도 추가
}
## Warning in SpatialPolygons2map(database, namefield = namefield): database does
## not (uniquely) contain the field 'name'.

PM2.5 농도에 따른 변화를 색깔별로 지도위에 표출

# 라이브러리가 없다면
install.packages("squash")
library("squash")
### PM.25의 농도 범위 확인
range(m.seoul_latlon$PM25)
## [1] 12.16667 22.16667
### color table 만들기
pm.25_map = makecmap(c(12, 23), colFn = bluered, outlier = 'red')

### 농도에 따른 사이즈 정하기
x1 = c(12, 23)
y1 = c(1.5, 5)

sl = lm(y1 ~ x1)

{plot(m.seoul_latlon$lat ~ m.seoul_latlon$lon, xlab = "", ylab = "", pch = 21,
      xlim = c(126.75, 127.2), ylim = c(37.42, 37.72),
      bg = cmap(m.seoul_latlon$PM25, map = pm.25_map),
      cex = m.seoul_latlon$PM25*sl$coefficients[2] + sl$coefficients[1]
      ) # 위경도에 따른 심볼 그림
mtext(text = expression("Latitude ("*degree*"E)"), line = 2, side = 2) # y축 이름
mtext(text = expression("Longtitude ("*degree*"E)"), line = 2, side = 1) # x축 이름
map(Seoul, add = T) # 지도 추가
vkey(pm.25_map, expression("PM"[2.5]), y = 37.43, x = 127.2, stretch = 1) # 범례 추가
}
## Warning in SpatialPolygons2map(database, namefield = namefield): database does
## not (uniquely) contain the field 'name'.

기본 상호작용하는 지도 생성하기

패키지 설치

install.packages("leaflet")

패키지 불러오기

library("leaflet")
{
  leaflet() %>%
  addTiles() %>%
  setView( lng= mean(m.seoul_latlon$lon), lat = mean(m.seoul_latlon$lat), zoom = 11) %>%
  addProviderTiles("Esri.WorldImagery") %>%
  addCircleMarkers(lng = m.seoul_latlon$lon, 
             lat = m.seoul_latlon$lat, 
             popup = paste(m.seoul_latlon$Name, "</br> PM2.5 </br>", m.seoul_latlon$PM25, "ug/m^3"))
}