1. 프로젝트 개요

본 보고서는 성남시의 지형적 특성(경사도)과 인구 동태(거주 및 유동 인구)를 결합하여, 응급 차량인 구급차가 실제 주행 시 겪는 물리적/사회적 제약 사항을 수치화하고 골든아워 확보를 위한 취약 구간을 도출하는 과정을 담고 있습니다.

2. 데이터 로드 및 전처리

2.1. 라이브러리 및 데이터 호출

KTDB 도로망 네트워크, 수치표고모델(DEM), 격자 인구 데이터를 호출하고 EPSG:5179 좌표계로 정밀 동기화합니다.

library(sf)
library(terra)
library(tidyverse)
library(scales)
library(leaflet)
library(viridis)
library(data.table)

# 1. 데이터 로드
grid_pop <- st_read("Lattice/vl_blk.shp", quiet = TRUE)
dem_raw <- rast("2025SUWON37709.img")

# 텍스트 파일 파싱 (KTDB 노드/링크)
nodes <- read.table("network_2025/01. Node_Link data/2025node.txt", skip = 1, fill = TRUE, stringsAsFactors = FALSE) %>% 
  filter(V1 %in% c("a", "a*")) %>% 
  mutate(node_id = as.character(V2), x = as.numeric(V3), y = as.numeric(V4)) %>% 
  select(node_id, x, y) %>% drop_na()

links <- read.table("network_2025/01. Node_Link data/2025link.txt", skip = 1, fill = TRUE, stringsAsFactors = FALSE) %>% 
  filter(V1 == "a") %>% 
  mutate(from_node = as.character(V2), to_node = as.character(V3)) %>% 
  select(from_node, to_node)

links_geom <- links %>%
  inner_join(nodes, by = c("from_node" = "node_id")) %>% rename(from_x = x, from_y = y) %>%
  inner_join(nodes, by = c("to_node" = "node_id")) %>% rename(to_x = x, to_y = y) %>%
  filter(!is.na(from_x), !is.na(from_y), !is.na(to_x), !is.na(to_y))

# KATEC PROJ4 문자열 정의
katec_proj4 <- "+proj=tmerc +lat_0=38 +lon_0=128 +k=0.9999 +x_0=400000 +y_0=600000 +ellps=bessel +towgs84=-115.80,474.99,674.11,1.16,-2.31,-1.63,6.43 +units=m +no_defs"

wkt_strings <- sprintf("LINESTRING (%.3f %.3f, %.3f %.3f)", 
                       links_geom$from_x, links_geom$from_y, links_geom$to_x, links_geom$to_y)

road_link <- st_sf(
  from_node = links_geom$from_node,
  to_node = links_geom$to_node,
  geometry = st_as_sfc(wkt_strings, crs = katec_proj4)
)

# EPSG:5179(도로명주소 표준)로 모든 데이터 강제 통일
target_crs <- 5179
grid_pop <- st_set_crs(grid_pop, target_crs) %>% st_transform(target_crs)
road_link <- st_transform(road_link, target_crs)
dem_raw <- project(dem_raw, paste0("EPSG:", target_crs))

boundary_sn <- st_union(grid_pop)
road_clipped <- st_intersection(road_link, boundary_sn)

2.2. 데이터 정밀 정제 (공간/물리 필터링)

분석의 신뢰도를 높이기 위해 구급차 진입이 불가능한 서울공항(군사시설) 부지를 마스킹 처리하고, 데이터 노이즈로 발생하는 1,200m 이상의 비정상 장대선을 제거합니다.

# 1. 서울공항 대략적인 바운더리 설정 
airport_mask <- st_polygon(list(matrix(c(
  964000, 1939000,
  967500, 1939000,
  967500, 1944000,
  964000, 1944000,
  964000, 1939000
), ncol = 2, byrow = TRUE))) %>% 
  st_sfc(crs = 5179)

# 2. 원본 뼈대인 road_clipped에서 미리 공항 구역 cut
road_clipped <- st_filter(road_clipped, airport_mask, .predicate = st_disjoint)

# 3. 길이 기반 필터링 (1.2km 이상 장대선 제거)
road_cleaned <- road_clipped %>%
  mutate(link_len = as.numeric(st_length(geometry))) %>%
  filter(link_len < 1200) %>%  
  st_filter(airport_mask, .predicate = st_disjoint) 

road_clipped <- road_cleaned
cat("정제 완료: 남은 유효 도로 링크 수 =", nrow(road_clipped), "개\n")
## 정제 완료: 남은 유효 도로 링크 수 = 8281 개

3. 하이브리드 페널티 모델링

구급차의 주행 한계를 모사하기 위해 물리적 요인(경사)과 사회적 요인(거주/유동 인구)을 결합한 다중 감속 모델을 구축합니다. 최종 체감 주행 속도(\(V_{adj}\))는 다음과 같이 산출됩니다.

\[V_{adj} = V_{base} \times e^{-0.035 \times \max(0, \text{Slope} - 3)} \times P_{resident} \times P_{floating}\]

# [지형 분석] 도로별 평균 경사도 매칭
boundary_vect <- vect(boundary_sn)
dem_masked <- mask(crop(dem_raw, boundary_vect), boundary_vect)
slope_sn <- terrain(dem_masked, v = "slope", unit = "degrees")

road_slope <- terra::extract(slope_sn, vect(road_clipped), fun = mean, na.rm = TRUE)
road_clipped$avg_slope <- ifelse(is.na(road_slope[, 2]), 0, road_slope[, 2])

# [인구/통신 결합] 
grid_pop_valid <- grid_pop %>% filter(!is.na(val))
road_final <- st_join(road_clipped, grid_pop_valid %>% select(val), join = st_nearest_feature)

# 통신 데이터 로드 및 전처리 (18~20시 평균)
telecom_raw <- fread("T11_성남시.csv", encoding = "UTF-8")
telecom_agg <- telecom_raw %>%
  filter(O_TIME_CD %in% c("18", "19", "20")) %>% 
  group_by(O_ADMI_NM, O_CENTER_X, O_CENTER_Y) %>%
  summarise(avg_floating_pop = sum(CNT, na.rm = TRUE), .groups = "drop")

telecom_sf <- st_as_sf(telecom_agg, coords = c("O_CENTER_X", "O_CENTER_Y"), crs = 5179)
road_hybrid <- st_join(road_final, telecom_sf %>% select(avg_floating_pop), join = st_nearest_feature)

# NA 처리
avg_pop <- mean(grid_pop$val, na.rm = TRUE)
global_avg_floating <- mean(telecom_agg$avg_floating_pop, na.rm = TRUE)

# [핵심 수식] 하이브리드 페널티 적용
road_hybrid <- road_hybrid %>%
  mutate(
    val = ifelse(is.na(val) | val <= 0, avg_pop, val),
    avg_floating_pop = ifelse(is.na(avg_floating_pop), global_avg_floating, avg_floating_pop),
    link_length_m = as.numeric(st_length(geometry)),
    
    # 1. 기존 페널티들
    base_speed_mpm = 833, 
    slope_penalty = exp(-0.035 * pmax(0, avg_slope - 3)),
    
    # 2. 격자 기반 거주 인구 페널티 
    pop_threshold = quantile(val, 0.9, na.rm = TRUE),
    resident_penalty = 1 - (0.15 * (pmin(val, pop_threshold) / pop_threshold)),
    
    # 3. 통신 데이터 기반 유동 인구 페널티 
    floating_threshold = quantile(avg_floating_pop, 0.9, na.rm = TRUE),
    floating_penalty = 1 - (0.15 * (pmin(avg_floating_pop, floating_threshold) / floating_threshold)),
    
    # 4. 최종 주행 속도 및 시간 산출
    adj_speed_mpm = base_speed_mpm * slope_penalty * resident_penalty * floating_penalty,
    travel_time_min = link_length_m / adj_speed_mpm,
    time_per_100m = (travel_time_min / link_length_m) * 100
  )

# 밤 시간대용 데이터도 동일하게 val 전처리 적용
road_final <- road_final %>%
  mutate(
    val = ifelse(is.na(val) | val <= 0, avg_pop, val),
    link_length_m = as.numeric(st_length(geometry)),
    base_speed_mpm = 833,
    slope_penalty = exp(-0.035 * pmax(0, avg_slope - 3)),
    threshold_val = quantile(val, 0.9, na.rm = TRUE),
    pop_penalty = 1 - (0.2 * (pmin(val, threshold_val) / threshold_val)),
    adj_speed_mpm = base_speed_mpm * slope_penalty * pop_penalty,
    travel_time_min = link_length_m / adj_speed_mpm,
    time_per_100m = (travel_time_min / link_length_m) * 100
  )

4. 시각화 결과: 시계열 연동 하이브리드 대시보드

도출된 위험도(time_per_100m)를 백분위수(percent_rank)로 변환하여 0.0부터 1.0까지의 정규화된 스코어를 부여합니다. 유동 인구가 결합된 Day 시나리오와 거주/지형 중심의 Night 시나리오를 토글형 지도로 구현했습니다.

# [데이터 정제] Day vs Night
road_day_leaf <- road_hybrid %>%
  select(val, avg_floating_pop, avg_slope, time_per_100m) %>%
  st_transform(4326) %>%
  mutate(
    rank_score = percent_rank(time_per_100m),
    dynamic_weight = (rank_score * 4) + 1.5,
    dynamic_opacity = (rank_score * 0.6) + 0.3
  )

road_night_leaf <- road_final %>%
  select(val, avg_slope, time_per_100m) %>%
  st_transform(4326) %>%
  mutate(
    rank_score = percent_rank(time_per_100m),
    dynamic_weight = (rank_score * 4) + 1.5,
    dynamic_opacity = (rank_score * 0.6) + 0.3
  )

pal_vivid <- colorNumeric(palette = "RdYlGn", domain = c(0, 1), reverse = TRUE)


# [Leaflet 대시보드 렌더링]
map_hybrid_system <- leaflet(width = "100%") %>%
  # 텍스트 라벨을 제일 위로 올리기 위한 Pane 설정만 남김
  addMapPane("text_pane", zIndex = 400) %>%
  
  # 배경 지도 (Pane 할당 제거로 렌더링 충돌 방지)
  addTiles(group = "☀️ 낮 시간대 분석 (유동인구 반영)") %>%
  addProviderTiles(providers$CartoDB.DarkMatterNoLabels, group = "🌙 밤 시간대 분석 (지형/거주 중심)") %>%
  
  # 낮 데이터
  addPolylines(
    data = road_day_leaf,
    group = "☀️ 낮 시간대 분석 (유동인구 반영)",
    color = ~pal_vivid(rank_score), 
    weight = ~dynamic_weight, opacity = ~dynamic_opacity,
    popup = ~sprintf("<div style='color:black;'><b>☀️ 낮 시간 주행 분석</b><hr><b>유동인구:</b> %d명<br><b>경사도:</b> %.2f도<br><b>위험도:</b> 상위 %.1f%%</div>", 
                     as.integer(avg_floating_pop), avg_slope, (1-rank_score)*100) %>% lapply(htmltools::HTML)
  ) %>%
  
  # 밤 데이터
  addPolylines(
    data = road_night_leaf,
    group = "🌙 밤 시간대 분석 (지형/거주 중심)",
    color = ~pal_vivid(rank_score), 
    weight = ~dynamic_weight, opacity = ~dynamic_opacity,
    popup = ~sprintf("<div style='color:black;'><b>🌙 밤 시간 주행 분석</b><hr><b>거주인구:</b> %d명<br><b>경사도:</b> %.2f도<br><b>위험도:</b> 상위 %.1f%%</div>", 
                     as.integer(val), avg_slope, (1-rank_score)*100) %>% lapply(htmltools::HTML)
  ) %>%
  
  # 야간 지명 라벨 (최상단 배치)
  addProviderTiles(providers$CartoDB.DarkMatterOnlyLabels, group = "🌙 밤 시간대 분석 (지형/거주 중심)", options = providerTileOptions(pane = "text_pane")) %>%
  
  # 성남시 윤곽선
  addPolygons(
    data = st_transform(boundary_sn, 4326),
    color = "#00BFFF", fill = FALSE, weight = 2, opacity = 0.5
  ) %>%
  
  # 컨트롤러
  addLayersControl(
    baseGroups = c("☀️ 낮 시간대 분석 (유동인구 반영)", "🌙 밤 시간대 분석 (지형/거주 중심)"),
    options = layersControlOptions(collapsed = FALSE),
    position = "topright"
  ) %>%
  
  addLegend(pal = pal_vivid, values = c(0, 1), title = "주행 위험도", position = "bottomright")

map_hybrid_system