※ 본 문서는 PC 환경에 최적화되어있습니다.

1 문서를 둘로 나눈 이유..

추후에 추가로 설명하겠지만, 가시선 분석에 사용한 ‘ArcPy’ 모듈이 있는 Python이 32bit밖에 지원을 안 하기 때문에, 어쩔 수 없이 R 또한 32bit버전으로 사용하여 고해상도 격자형 자료에 필수적인 메모리 할당을 약 4Gb밖에 하지 못 하였다. 따라서 합본시 메모리 부족으로 인한 마크다운 에러가 발생하므로, 어쩔 수 없이 두 부분으로 나누게 되었다.

2 서울세계불꽃축제: 명당은 어디?

서울에서 해마다 열리는 대표적인 불꽃 축제인 서울세계불꽃축제는 코로나19의 영향으로 올해는 개최하지는 않지만, 해마다 열리므로 추후 구경하는 자리를 선정할 때 참고가 될겸 명당자리를 분석해보았다.

2.1 분석 방법

원래는 관측지점이 A이지만, 이번 분석에서는 B가 관측가능지점(다수, 면)이고 A에서(80m 오프셋) 폭죽이 터지는 가정이다.

위 그림과 같이 가시권 분석(viewshed)을 역이용하는 방식을 사용하였고, 분석에 사용한 자료는 eTL 에서 구득한 서울시 30m급 수치표고모델(DEM)과 국토교통부 산하 국가공간정보포털 GIS건물통합정보 에서 건물 높이 데이터를 합하여 간이 수치표면모델(DSM)을 만들어 사용하였다.

DEM과 DSM의 차이

가시권 분석이란?

2.2 사용한 라이브러리

library(tidyverse)  # 이제는 기본 베이스
library(sf)  # 공간 분석
library(rgdal)  #    ''
library(raster)  #    ''
library(rgeos)  #    ''
library(windfarmGA)  # r 코드로 가시권 분석 가능
library(tmap)  # 지도 도식
library(leafem)  #    ''

2.3 좌표정보 획득

TA 세션에서 구득한 서울시 구 행정구역 셰이프를 도식 및 좌표계 획득을 목적으로 불러왔다.

# seoul shp
sigungu_shp <- readOGR(dsn = "./data/shp/seoul", layer = "tl_emd._seoul_4326", encoding = "UTF-8", 
    stringsAsFactors = FALSE)

OGR data source with driver: ESRI Shapefile 
Source: "C:\Users\dydgn\OneDrive - 공주대학교\바탕 화면\Analytical_Methods_for_Spatial_Information_Assignment\Gini-09\data\shp\seoul", layer: "tl_emd._seoul_4326"
with 467 features
It has 3 fields

crs <- crs(sigungu_shp)
# CRS.new <- CRS('+proj=tmerc +lat_0=38 +lon_0=127 +k=1 +x_0=200000 +y_0=500000
# +ellps=bessel +units=m')
CRS.new <- CRS("+proj=utm +zone=52 +ellps=WGS84 +towgs84=0,0,0,-0,-0,-0,0 +units=m +no_defs")

2.4 DSM 제작

레스터 자료의 Cell size가 분석 단위의 기준이 되므로 모든 자료의 좌표투영 기준을 PCS로 진행하였다(dd -> meter).

2.4.1 건물 셰이프 불러오기

국가공간정보포털에서 구득한 GIS건물통합정보 셰이프를 R에다가 직접 불러와서 작업하려고 하였으나,

위와 같이 메모리 리소스를 셰이프 하나에 많이 할당하여, 보유중인 작업환경에서 분석이 불가능하다 판단, 상용 GIS 프로그램인 ArcGIS에서 백터형의 자료를 래스터로 변환한 후에 불러왔다.

ArcGIS에서 행한 과정은 위와 같다.
격자 크기는 5m by 5m로 지정하여 분석의 정확도와 소요시간이 적절히 균형을 이루게 하였다.

# buildings

building <- raster("./data/ras/buildings/buildings.img")
building <- reclassify(building, cbind(NA, 0))
crs(building)

CRS arguments:
 +proj=utm +zone=52 +ellps=WGS84 +towgs84=0,0,0,-0,-0,-0,0 +units=m
+no_defs

# tmap_mode('view') tm_shape(dem)+tm_raster(palette = 'RdBu')

래스터 연산을 수행하기 위해, 건물 래스터 자료에 존재하는 NULL 값을 0으로 설정하였다.

2.4.2 DEM 불러오기

eTL에서 구득한 서울시 30m DEM을 불러왔다.

# DEM
dem <- raster("./data/ras/dem/seoul_dem.img")
dem <- projectRaster(dem, crs = CRS.new)
# ext <- extent(dem) ncl <- ncol(dem) nrw <- nrow(dem)

dem5m <- disaggregate(dem, fact = 6)
res(dem5m)

[1] 5 5

crs(dem5m)

CRS arguments:
 +proj=utm +zone=52 +ellps=WGS84 +towgs84=0,0,0,-0,-0,-0,0 +units=m
+no_defs

origin(dem5m)

[1] -2.407558 -1.540132

격자 크기를 건물 자료와 같게 해주려고 격자크기 30m를 6으로 나눠줬다.

2.5 사례지역 한정하기(버퍼 분석)

2.5.1 AOI 지정

tm1 <- tm_shape(sigungu_shp) + tm_polygons(alpha = 0)
lf <- tmap_leaflet(tm1)
addMouseCoordinates(lf)

마우스 커서의 위치에 따른 경위도 좌표가 지도의 좌상단에 나오는 위의 지도를 활용하여 Fireworks Area를 좌표로 획득하였다.

다음은 획득한 좌표를 가지고 Area of Interest(AOI) 중심 포인트를 만드는 과정이다.

### epicenter

pts <- data.frame(lon = 126.94861, lat = 37.5229)
pts$lon <- as.double(pts$lon)
pts$lat <- as.double(pts$lat)
epicenter <- SpatialPointsDataFrame(pts[, c(1:2)], pts, proj4string = crs)

pcsaoi <- spTransform(epicenter, CRS.new)
# plot(epicenter)
tmap_mode("view")
tm_shape(pcsaoi) + tm_dots()

buf_5km <- gBuffer(pcsaoi, width = 2500)

중심 포인트를 가지고 PCS로 변환한 다음, 2.5km 반경의 버퍼를 만들었고,

### echo buffer
tm_shape(buf_5km) + tm_polygons(alpha = 0.3, col = "green") + tm_shape(pcsaoi) + 
    tm_dots()

writeOGR(pcsaoi, paste0(getwd(), "/data/tmp"), "aoi", driver = "ESRI Shapefile", 
    layer_options = "encoding=UTF-8", overwrite_layer = TRUE)

추후 ArcPy에서 작업을 목적으로 셰이프를 저장하였다.

2.5.2 마스킹을 통한 추출

버퍼를 가지고 반경 2.5km 지역만 추출하는 과정이다.

DEM

### extract by mask
em1 <- crop(dem5m, extent(buf_5km), snap = "out")
fr1 <- rasterize(buf_5km, em1)
lr1 <- mask(x = em1, mask = fr1)

plot(lr1)

건물 고도

em2 <- crop(building, extent(buf_5km), snap = "out")
fr2 <- rasterize(buf_5km, em2)
lr2 <- mask(x = em2, mask = fr2)

plot(lr2)

2.6 Map Algebra

서로 다른 레스터의 중심점을 맞춰주고, 격자형 자료 연산을 해주었다.

- 건물 높이 자료 + DEM = DSM

# map algebra
origin(lr2) <- origin(lr1)
aoidsm <- lr1 + lr2

plot(aoidsm)

setwd("./data/tmp/")
writeRaster(aoidsm, filename = "aoi.img", overwrite = TRUE)
setwd("./../../")

만든 DSM을 저장하여, ArcPy로 불러올 수 있게 하였다.

2.7 Viewshed

2.7.1 ArcPy란?

ArcPy는 ESRI사의 ArcGIS/ArcPRO에 포함된 Python 모듈로, 유용하고, 생산적인 방식의 지리정보분석, 데이터 변환, 데이터 관리, 지도 자동화 등을 제공한다. 병렬 연산 등을 지원하며, 잘 활용하면 강력한 지오프로세싱 도구가 될 수 있다. 현재 32bit 버전의 Python을 사용중이라, 해당 코드를 사용할 때, R 또한 32bit 버전을 사용하여 제한사항이 많았다.

2.7.2 메모리 확보

32비트 프로그램의 고질적인 문제인, 시스템 메모리를 약 4Gb밖에 할당을 못하여, 최대한 공간을 확보하고자, 더 이상 필요없는 Env변수들은 삭제하였다.

rm(dem5m)
rm(building)

2.7.3 가시권 분석

라이브러리 로드

R에서 Python 명령을 사용할 수 있게 해주는 라이브러리를 불러오고, 작업 PC에 설치된 다른 버전의 Python(3.X)과 혼동되지 않도록 ArcGIS 내의 Python 실행파일을 명시하였다.

# Sys.setenv(RETICULATE_PYTHON = 'C:/Program Files/GRASS GIS 7.8/Python37')
library(reticulate)
use_python("C:/Python27/ArcGIS10.5/python.exe", required = TRUE)
# py_run_file('./data/tmp/viewshed.py')

파이선 스크립트

다음은 가시권 분석과, Tmap에 도식하기 위하여 무거운 고해상도 격자형 가시범위를 Raster to Polygon을 사용하여 셰이프로 만들어주는 작업이다.

ArcPy의 Spatial Analysist 익스텐션을 사용하였고, ArcGIS의 한글 배척을 고려하여, 경로를 D 드라이브로 잡고, 위에서 출력한 DSM과 중심점 셰입을 변수로 넣어주었다.

가시선 분석에서 사용된 변수로는: 관측지점 오프셋 = 80M(폭죽이 터지는 고도), 지표면 오프셋 = 1.5M(사람 눈높이) 를 사용하였다.

import arcpy
from arcpy import env
from arcpy.sa import *

# Set environment settings
env.workspace = "D:/arcpy"

# Check out the ArcGIS Spatial Analyst extension license
arcpy.CheckOutExtension("Spatial")

#input var

u'CheckedOut'

inRaster = Raster("D:/arcpy/dsm.img")
inObservers = "D:/arcpy/aoi.shp"
outAGL = ""
analysisType = "OBSERVERS"
verticalError = ""
outAnalysisRelationTable = ""
refractCoeff = ""
surfaceOffset = "1.5"
observerElevation = ""
observerOffset = "80"
innerRadius = ""
innerIs3D = "False"
outerRadius = ""
outerIs3D = "True"
horizStartAngle = ""
horizEndAngle = ""
vertUpperAngle = ""
vertLowerAngle = ""
analysisMethod = "ALL_SIGHTLINES"

# Execute Viewshed2
outViewshed2 = Viewshed2(inRaster, inObservers, outAGL, analysisType,
                         verticalError, outAnalysisRelationTable, refractCoeff,
                         surfaceOffset, observerElevation, observerOffset,
                         innerRadius, innerIs3D, outerRadius, outerIs3D,
                         horizStartAngle, horizEndAngle, vertUpperAngle,
                         vertLowerAngle, analysisMethod)

# Save the output
outViewshed2.save("D:/arcpy/output.img")

# Set local variables
inRaster = Raster("D:/arcpy/output.img")
outPolygons = "D:/arcpy/viewshed.shp"
field = "VALUE"

# Execute RasterToPolygon
arcpy.RasterToPolygon_conversion(inRaster, outPolygons, "SIMPLIFY", field)

<Result 'D:\\arcpy\\viewshed.shp'>

2.7.4 결과물 도식

arcviewshed <- readOGR(dsn = "D:/arcpy", layer = "viewshed", encoding = "UTF-8", 
    stringsAsFactors = FALSE)

OGR data source with driver: ESRI Shapefile 
Source: "D:\arcpy", layer: "viewshed"
with 2170 features
It has 2 fields
Integer64 fields read as strings:  Id gridcode

tm2 <- tm_shape(arcviewshed) + tm_polygons(alpha = 0.3, col = "green") + tm_shape(pcsaoi) + 
    tm_dots() + tm_scale_bar()
lf2 <- tmap_leaflet(tm2)
addMouseCoordinates(lf2)

초록색이 가시가능지역이다.

인터넷에서 쉽게 구할 수 있는 명당 지도이다.

위의 명당 지도 상의 R&D건물 옥상이 분석 결과 가시가능지역으로 나타났고..

63빌딩 뒤가 전부 가시 불가지역인 것으로 나타난다.

2.8 R Standalone Viewshed

참고로, R의 windfarmGA 라이브러리는 가시선 분석을 제공한다. 이를 사용하여 가시선 분석을 진행하는데, 64bit R에서는 잘 돌아갔지만, 32bit R에서는 메모리 부족으로, 실행이 불가하였기에, 코드와 결과물만 참고.

res <- viewshed(r = aoidsm, shape = buf_5km, turbine_locs = pcsaoi@coords, h1 = 80, 
    h2 = 1.7)
plot_viewshed(res, legend = T, interpolate = T, colNA = "black", col = topo.colors(15))

비슷하게 나온 것을 확인할 수 있다.

공간정보분석

ASS-09-2
‘Let me show you something’

by: 김용길, 백종찬, 서용훈

공간정보분석

1 문서를 둘로 나눈 이유..

2 서울세계불꽃축제: 명당은 어디?

2.1 분석 방법

2.2 사용한 라이브러리

2.3 좌표정보 획득

2.4 DSM 제작

2.4.1 건물 셰이프 불러오기

2.4.2 DEM 불러오기

2.5 사례지역 한정하기(버퍼 분석)

2.5.1 AOI 지정

2.5.2 마스킹을 통한 추출

DEM

건물 고도

2.6 Map Algebra

2.7 Viewshed

2.7.1 ArcPy란?

2.7.2 메모리 확보

2.7.3 가시권 분석

라이브러리 로드

파이선 스크립트

2.7.4 결과물 도식

2.8 R Standalone Viewshed

공간정보분석

ASS-09-2 ‘Let me show you something’ by: 김용길, 백종찬, 서용훈

공간정보분석

1 문서를 둘로 나눈 이유..

2 서울세계불꽃축제: 명당은 어디?

2.1 분석 방법

2.2 사용한 라이브러리

2.3 좌표정보 획득

2.4 DSM 제작

2.4.1 건물 셰이프 불러오기

2.4.2 DEM 불러오기

2.5 사례지역 한정하기(버퍼 분석)

2.5.1 AOI 지정

2.5.2 마스킹을 통한 추출

DEM

건물 고도

2.6 Map Algebra

2.7 Viewshed

2.7.1 ArcPy란?

2.7.2 메모리 확보

2.7.3 가시권 분석

라이브러리 로드

파이선 스크립트

2.7.4 결과물 도식

2.8 R Standalone Viewshed

ASS-09-2
‘Let me show you something’

by: 김용길, 백종찬, 서용훈