> ## 07-1 어느 지역이 제일 비쌀까?
>
> # Step 1 : 지역별 평균 가격 구하기
>
> setwd(dirname(rstudioapi::getSourceEditorContext()$path))
> load("./06_geodataframe/06_apt_price.rdata")
> library(sf)
> grid <- st_read("./01_code/sigun_grid/seoul.shp") # 서울시 1km 그리드 불러오기
Reading layer `seoul' from data source `D:\1_Study\1_BigData\01_R\02_Doit_R_Shiny\01_code\sigun_grid\seoul.shp' using driver `ESRI Shapefile'
Simple feature collection with 694 features and 2 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 126.7645 ymin: 37.42899 xmax: 127.1835 ymax: 37.70146
Geodetic CRS: GCS_unknown
> apt_price <- st_join(apt_price, grid, join = st_intersects) # 실거래+그리드 결합
> head(apt_price, 2)
Simple feature collection with 2 features and 15 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: 126.9674 ymin: 37.57442 xmax: 126.9689 ymax: 37.58621
CRS: +proj=longlat +datum=WGS84 +no_defs
ymd ym year code addr_1 apt_nm juso_jibun price
1 2021-01-14 2021-01-01 2021 11110 서울_종로 청운현대 서울특별시 종로구 청운동 56-45 청운현대 130000
2 2021-01-07 2021-01-01 2021 11110 서울_종로 광화문스페이스본 서울특별시 종로구 사직동 9-1 광화문스페이스본 150000
con_year area floor py cnt ID SIG_CD geometry
1 2000 130 2 3300 1 83556 1100 POINT (126.9674 37.58621)
2 2008 145 6 3414 1 83095 1100 POINT (126.9689 37.57442)
>
> kde_high <- aggregate(apt_price$py, by=list(apt_price$ID), mean)
> colnames(kde_high) <- c("ID", "avg_price")
> head(kde_high, 2)
ID avg_price
1 79520 1506.167
2 79765 3600.258
>
> # Step 2 : 평균 가격 정보 표시하기
>
> kde_high <- merge(grid, kde_high, by="ID")
> library(ggplot2)
> library(dplyr)
> kde_high %>% ggplot(aes(fill = avg_price)) +
+ geom_sf() +
+ scale_fill_gradient(low = "white", high = "red")
>
>
> # Step 3 : 지도 경계 그리기
>
> library(sp)
> kde_high_sp <- as(st_geometry(kde_high), "Spatial")
> x <- coordinates(kde_high_sp)[,1]
> y <- coordinates(kde_high_sp)[,2]
>
> l1 <- bbox(kde_high_sp)[1,1] - (bbox(kde_high_sp)[1,1] * 0.0001)
> l2 <- bbox(kde_high_sp)[1,2] + (bbox(kde_high_sp)[1,2] * 0.0001)
> l3 <- bbox(kde_high_sp)[2,1] - (bbox(kde_high_sp)[2,1] * 0.0001)
> l4 <- bbox(kde_high_sp)[2,2] + (bbox(kde_high_sp)[1,1] * 0.0001)
>
> install.packages("spatstat")
> library(spatstat)
> win <- owin(xrange = c(l1,l2), yrange = c(l3,l4))
> plot(win)
> rm(list = c("kde_high_sp", "apt_price", "l1", "l2", "l3", "l4"))
>
> # Step 4 : 밀도 그래프 표시하기
>
> # 포인트 패턴 분석
> p <- ppp(x, y, window = win) # 경계선 위에 좌표값 포인트 생성
> d <- density.ppp(p, weights = kde_high$avg_price, # 커널 밀도 함수로 변환
+ sigma = bw.diggle(p),
+ kernel = 'gaussian')
> plot(d) # 밀도 그래프 확인
> rm(list = c("x", "y", "win", "p"))
>
> # Step 5 : 래스터 이미지로 변환하기
>
> # 래스터 이미지 : 일반적으로 직사각형 격자의 화소, 색의 점을 모니터, 종이 등의 매체에 표시하는 자료 구조
> # 커널 밀도 추정에서 의미 있는 데이터는 보통 상위 20~30% 범위 내에 있다.
> # 여기에서는 전체 데이터의 25% 지점만 의미 있는 신호라고 가정.
> d[d < quantile(d)[4] + (quantile(d)[4] * 0.1)] <- NA
> library(raster)
> raster_high <- raster(d) # 래스터 변환
> plot(raster_high)
>
> # Step 6 : 불필요한 부분 자르기
>
> bnd <- st_read("./01_code/sigun_bnd/seoul.shp") # 서울시 경계선 불러오기
Reading layer `seoul' from data source `D:\1_Study\1_BigData\01_R\02_Doit_R_Shiny\01_code\sigun_bnd\seoul.shp' using driver `ESRI Shapefile'
Simple feature collection with 1 feature and 1 field
Geometry type: POLYGON
Dimension: XY
Bounding box: xmin: 126.7645 ymin: 37.42899 xmax: 127.1835 ymax: 37.70146
Geodetic CRS: GCS_unknown
> raster_high <- crop(raster_high, extent(bnd)) # 외곽선 자르기
> crs(raster_high) <- sp::CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0") # 좌표계 정의
> plot(raster_high)
> plot(bnd, col=NA, border="red", add=TRUE)
>
> # Step 7 : 지도 그리기
>
> install.packages("rgdal")
> library(rgdal)
> library(leaflet)
> leaflet() %>%
+ #---# 기본 지도 불러오기
+ addProviderTiles(providers$CartoDB.Positron) %>%
+ #---# 서울시 경계선 불러오기
+ addPolygons(data = bnd, weight = 3, color = "red", fill = NA) %>%
+ #---# 래스터 이미지 불러오기
+ addRasterImage(raster_high,
+ colors = colorNumeric(c("blue", "green", "yellow", "red"),
+ values(raster_high), na.color = "transparent"),
+ opacity = 0.4)
Warning messages:
1: In wkt(projfrom) : CRS object has no comment
2: In showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj = prefer_proj) :
Discarded ellps WGS 84 in Proj4 definition: +proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defs +type=crs
3: In showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj = prefer_proj) :
Discarded datum World Geodetic System 1984 in Proj4 definition
4: In wkt(pfrom) : CRS object has no comment
5: In rgdal::rawTransform(projfrom, projto, nrow(xy), xy[, 1], xy[, :
Using PROJ not WKT2 strings
6: In wkt(pfrom) : CRS object has no comment
7: In rgdal::rawTransform(projfrom, projto, nrow(xy), xy[, 1], xy[, :
Using PROJ not WKT2 strings
8: In rgdal::rawTransform(projto_int, projfrom, nrow(xy), xy[, 1], :
Using PROJ not WKT2 strings
9: In wkt(pfrom) : CRS object has no comment
10: In rgdal::rawTransform(projfrom, projto, nrow(xy), xy[, 1], xy[, :
Using PROJ not WKT2 strings
>
> # Step 8 : 평균 가격 정보 저장하기
>
> dir.create("07_map")
> save(raster_high, file = "./07_map/07_kde_high.rdata")
>
>
> ## 07-2 요즘 뜨는 지역이 어디일까?
>
> # Step 1 : 데이터 준비하기
>
> setwd(dirname(rstudioapi::getSourceEditorContext()$path))
> load("./06_geodataframe/06_apt_price.rdata")
> grid <- st_read("./01_code/sigun_grid/seoul.shp") # 서울시 1km 그리드 불러오기
Reading layer `seoul' from data source `D:\1_Study\1_BigData\01_R\02_Doit_R_Shiny\01_code\sigun_grid\seoul.shp' using driver `ESRI Shapefile'
Simple feature collection with 694 features and 2 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 126.7645 ymin: 37.42899 xmax: 127.1835 ymax: 37.70146
Geodetic CRS: GCS_unknown
> apt_price <- st_join(apt_price, grid, join = st_intersects) # 실거래+그리드 결합
> head(apt_price, 2)
Simple feature collection with 2 features and 15 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: 126.9674 ymin: 37.57442 xmax: 126.9689 ymax: 37.58621
CRS: +proj=longlat +datum=WGS84 +no_defs
ymd ym year code addr_1 apt_nm juso_jibun price
1 2021-01-14 2021-01-01 2021 11110 서울_종로 청운현대 서울특별시 종로구 청운동 56-45 청운현대 130000
2 2021-01-07 2021-01-01 2021 11110 서울_종로 광화문스페이스본 서울특별시 종로구 사직동 9-1 광화문스페이스본 150000
con_year area floor py cnt ID SIG_CD geometry
1 2000 130 2 3300 1 83556 1100 POINT (126.9674 37.58621)
2 2008 145 6 3414 1 83095 1100 POINT (126.9689 37.57442)
>
> # Step 2 : 이전/이후 데이터 세트 만들기
>
> kde_before <- subset(apt_price, ymd < '2021-07-01')
> kde_before <- aggregate(kde_before$py, by=list(kde_before$ID), mean)
> colnames(kde_before) <- c("ID", "before")
>
> kde_after <- subset(apt_price, ymd >= '2021-07-01')
> kde_after <- aggregate(kde_after$py, by=list(kde_after$ID), mean)
> colnames(kde_after) <- c("ID", "after")
>
> kde_diff <- merge(kde_before, kde_after, by = "ID")
> kde_diff$diff <- round( ((kde_diff$after - kde_diff$before) / kde_diff$before) * 100, 0) # 변화율 계산
> head(kde_diff, 2)
ID before after diff
1 79520 1451.000 1517.20 5
2 79765 3545.216 3681.72 4
>
> # Step 3 : 가격이 오른 지역 찾기
>
> library(sf)
> kde_diff <- kde_diff[kde_diff$diff > 0, ]
> kde_hot <- merge(grid, kde_diff, by = "ID") # 그리드에 상승 지역 결합
> library(ggplot2)
> library(dplyr)
> kde_hot %>%
+ ggplot(aes(fill = diff)) +
+ geom_sf() +
+ scale_fill_gradient(low = "white", high = "red")
>
> # Step 4 : 지도 경계 그리기
>
> library(sp)
> kde_hot_sp <- as(st_geometry(kde_hot), "Spatial")
> x <- coordinates(kde_hot_sp)[,1]
> y <- coordinates(kde_hot_sp)[,2]
>
> l1 <- bbox(kde_hot_sp)[1,1] - (bbox(kde_hot_sp)[1,1] * 0.0001)
> l2 <- bbox(kde_hot_sp)[1,2] + (bbox(kde_hot_sp)[1,2] * 0.0001)
> l3 <- bbox(kde_hot_sp)[2,1] - (bbox(kde_hot_sp)[2,1] * 0.0001)
> l4 <- bbox(kde_hot_sp)[2,2] + (bbox(kde_hot_sp)[1,1] * 0.0001)
>
> library(spatstat)
> win <- owin(xrange = c(l1,l2), yrange = c(l3,l4))
> plot(win)
> rm(list = c("kde_hot_sp", "apt_price", "l1", "l2", "l3", "l4"))
>
> # Step 5 : 밀도 그래프 표시하기
>
> # 포인트 패턴 분석
> p <- ppp(x, y, window = win) # 경계선 위에 좌표값 포인트 생성
> d <- density.ppp(p, weights = kde_hot$diff, # 커널 밀도 함수로 변환
+ sigma = bw.diggle(p),
+ kernel = 'gaussian')
> plot(d) # 밀도 그래프 확인
> rm(list = c("x", "y", "win", "p"))
>
> # Step 6 : 래스터 이미지로 변환하기
>
> d[d < quantile(d)[4] + (quantile(d)[4] * 0.1)] <- NA
> library(raster)
> raster_hot <- raster(d) # 래스터 변환
> plot(raster_hot)
>
> # Step 7 : 불필요한 부분 자르기
>
> bnd <- st_read("./01_code/sigun_bnd/seoul.shp") # 서울시 경계선 불러오기
Reading layer `seoul' from data source `D:\1_Study\1_BigData\01_R\02_Doit_R_Shiny\01_code\sigun_bnd\seoul.shp' using driver `ESRI Shapefile'
Simple feature collection with 1 feature and 1 field
Geometry type: POLYGON
Dimension: XY
Bounding box: xmin: 126.7645 ymin: 37.42899 xmax: 127.1835 ymax: 37.70146
Geodetic CRS: GCS_unknown
> raster_hot <- crop(raster_hot, extent(bnd)) # 외곽선 자르기
> crs(raster_hot) <- sp::CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0") # 좌표계 정의
> plot(raster_hot)
> plot(bnd, col=NA, border="red", add=TRUE)
>
> # Step 8 : 지도 그리기
>
> library(leaflet)
> leaflet() %>%
+ #---# 기본 지도 불러오기
+ addProviderTiles(providers$CartoDB.Positron) %>%
+ #---# 서울시 경계선 불러오기
+ addPolygons(data = bnd, weight = 3, color = "red", fill = NA) %>%
+ #---# 래스터 이미지 불러오기
+ addRasterImage(raster_hot,
+ colors = colorNumeric(c("blue", "green", "yellow", "red"),
+ values(raster_hot), na.color = "transparent"),
+ opacity = 0.4)
Warning messages:
1: In showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj = prefer_proj) :
Discarded ellps WGS 84 in Proj4 definition: +proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defs +type=crs
2: In showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj = prefer_proj) :
Discarded datum World Geodetic System 1984 in Proj4 definition
>
> # Step 9 : 평균 가격 변화율 정보 저장하기
>
> save(raster_hot, file = "./07_map/07_kde_hot.rdata")
>
>
> ## 07-3 우리 동네가 옆 동네보다 비쌀까?
>
> # Step 1 : 데이터 준비하기
>
> setwd(dirname(rstudioapi::getSourceEditorContext()$path))
> load("./06_geodataframe/06_apt_price.rdata")
> load("./07_map/07_kde_high.rdata")
> load("./07_map/07_kde_hot.rdata")
>
> library(sf)
> bnd <- st_read("./01_code/sigun_bnd/seoul.shp")
Reading layer `seoul' from data source `D:\1_Study\1_BigData\01_R\02_Doit_R_Shiny\01_code\sigun_bnd\seoul.shp' using driver `ESRI Shapefile'
Simple feature collection with 1 feature and 1 field
Geometry type: POLYGON
Dimension: XY
Bounding box: xmin: 126.7645 ymin: 37.42899 xmax: 127.1835 ymax: 37.70146
Geodetic CRS: GCS_unknown
> grid <- st_read("./01_code/sigun_grid/seoul.shp")
Reading layer `seoul' from data source `D:\1_Study\1_BigData\01_R\02_Doit_R_Shiny\01_code\sigun_grid\seoul.shp' using driver `ESRI Shapefile'
Simple feature collection with 694 features and 2 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 126.7645 ymin: 37.42899 xmax: 127.1835 ymax: 37.70146
Geodetic CRS: GCS_unknown
>
> # Step 2 : 마커 클러스터링 옵션 설정하기
>
> #---# 마커 클러스터링 함수(javascript avg.fomula) 등록
> load("./01_code/circle_marker/circle_marker.rdata")
> #---# 마커 클러스터링 색상 설정 : 상, 중, 하
> circle.colors <- sample(x=c("red","green","blue"), size = 1000, replace = TRUE)
>
> # Step 3 : 마커 클러스터링 시각화하기
>
> # install.packages("purrr")
> library(purrr)
> library(leaflet)
> leaflet() %>%
+ #---# 오픈스트리트맵 불러오기
+ addTiles() %>%
+ #---# 서울시 경계선 불러오기
+ addPolygons(data = bnd, weight = 3, color = "red", fill = NA) %>%
+ #---# 최고가 래스터 이미지 불러오기
+ addRasterImage(raster_high,
+ colors = colorNumeric(c("blue", "green", "yellow", "red"), values(raster_high)
+ , na.color = "transparent"), opacity = 0.4, group = "2021 최고가") %>%
+ #---# 급등지 래스터 이미지 불러오기
+ addRasterImage(raster_hot,
+ colors = colorNumeric(c("blue", "green", "yellow", "red"), values(raster_hot)
+ , na.color = "transparent"), opacity = 0.4, group = "2021 급등지") %>%
+ #---# 최고가/급등지 선택 옵션 추가하기
+ addLayersControl(baseGroups = c("2021 최고가", "2121 급등지"),
+ options = layersControlOptions(collapsed = FALSE)) %>%
+ #---# 마커 클러스터링 불러오기
+ addCircleMarkers(data = apt_price, lng = unlist(map(apt_price$geometry,1)),
+ lat = unlist(map(apt_price$geometry,2)), radius = 10, stroke = FALSE,
+ fillOpacity = 0.6, fillColor = circle.colors, weight = apt_price$py,
+ clusterOptions = markerClusterOptions(iconCreateFunction=JS(avg.formula)))
Warning messages:
1: In showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj = prefer_proj) :
Discarded ellps WGS 84 in Proj4 definition: +proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defs +type=crs
2: In showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj = prefer_proj) :
Discarded datum World Geodetic System 1984 in Proj4 definition
3: In wkt(projfrom) : CRS object has no comment
4: In wkt(pfrom) : CRS object has no comment
5: In rgdal::rawTransform(projfrom, projto, nrow(xy), xy[, 1], xy[, :
Using PROJ not WKT2 strings
6: In wkt(pfrom) : CRS object has no comment
7: In rgdal::rawTransform(projfrom, projto, nrow(xy), xy[, 1], xy[, :
Using PROJ not WKT2 strings
8: In rgdal::rawTransform(projto_int, projfrom, nrow(xy), xy[, 1], :
Using PROJ not WKT2 strings
9: In wkt(pfrom) : CRS object has no comment
10: In rgdal::rawTransform(projfrom, projto, nrow(xy), xy[, 1], xy[, :
Using PROJ not WKT2 strings
출처 : 김철민, ⌜공공데이터로 배우는 R 데이터분석 with 샤이니⌟, 이지스퍼블리싱, 2022
'데이터분석 > R' 카테고리의 다른 글
| [R 데이터분석 with 샤이니] 샤이니 입문하기 (0) | 2022.07.19 |
|---|---|
| [R 데이터분석 with 샤이니] 통계 분석과 시각화 (0) | 2022.07.12 |
| [R 데이터분석 with 샤이니] 지오 데이터프레임 만들기 (0) | 2022.07.03 |
| [R 데이터분석 with 샤이니] 카카오맵 API로 지오 코딩하기 (0) | 2022.07.02 |
| [R 데이터분석 with 샤이니] - 전처리 : 데이터를 알맞게 다듬기 (0) | 2022.07.02 |
