R 用于 GIS - 计算邻近密度

Question

我最近开始将 R 用于 GIS 应用程序，并需要计算邻近密度的帮助，这是我不太熟悉的概念。为了说明这一点，假设我获得了学校和操场位置数据，我想计算这些数据在 0.25 英里、0.50 英里、0.75 英里和 1 英里半径处的邻近密度。如何使用 R 程序实现这一点？这是我到目前为止所尝试过的。

# load required packages
library(sp)
library(sf)

# create standard rectangular non-spatial data frames
schools <- data.frame(Location = c("Albertville, AL", "Albertville, AL", "Albertville, AL", "Albertville, AL", "Albertville, AL", "Albertville, AL"), Place = c("Albertville Middle School", "Albertville High School", "Albertville Intermediate School", "Albertville Elementary School", "Albertville Kindergarten and PreK", "Albertville Primary School"), Lat = c(34.25996, 34.260509, 34.27279, 34.25182, 34.29161, 34.25321), Long = c(-86.2074, -86.205116, -86.220818, -86.22271, -86.19384, -86.22192))
playgrounds <- data.frame(Location = c("Albertville, AL", "Albertville, AL", "Albertville, AL"), Place = c("Graham Park", "Sand Mountain Park", "Pecan Grove Playground"), Lat = c(34.26358, 34.275471, 34.27294), Long = c(-86.20289, -86.229111, -86.2264))

# define the coordinates
schools_coords <- cbind(schools$Long, schools$Lat)
playgrounds_coords <- cbind(playgrounds$Long, playgrounds$Lat)

# create the SpatialPointsDataFrame
schools_sp <- SpatialPointsDataFrame(schools_coords, data = schools, proj4string = CRS("+proj= longlat"))
playgrounds_sp <- SpatialPointsDataFrame(playgrounds_coords, data = playgrounds, proj4string = CRS("+proj= longlat"))

# convert to sf
schools_sf <- st_as_sf(schools_sp)
playgrounds_sf <- st_as_sf(playgrounds_sp)

这也只是我的起点。我可能会根据需要收集其他公开数据。

Answer 1

工作流程：

• 确定适当的 UTM 区域 CRS（或使用您认为合适的）。这假设您的坐标是 WGS84/EPSG:4326。
• 为学校和操场创建 sf 对象，并项目到上一步确定的 CRS
• 在指定距离（以米（半径）为单位）处创建缓冲区，例如400 == ~0.25 英里以及每个学校缓冲区的操场计数总和
• 合并所有缓冲区 sf 文件并计算近端密度

有必要将数据投影到投影坐标系 (PCS)，以确保缓冲区计算在整个研究区域范围内保持相对恒定。如果您的数据位于国家级别，或跨越多个 UTM 区域，您可能需要为每个 UTM 区域运行单独的计算。您还可以预先计算面积值并将它们分配给 buff* 文件，但现在让我们保持简单。

library(sf)
library(dplyr)
library(tidyr)

# Your sample data
schools <- data.frame(
  Location = c("Albertville, AL", "Albertville, AL", "Albertville, AL", "Albertville, AL", "Albertville, AL", "Albertville, AL"),
  Place = c("Albertville Middle School", "Albertville High School", "Albertville Intermediate School", "Albertville Elementary School", "Albertville Kindergarten and PreK", "Albertville Primary School"),
  Lat = c(34.25996, 34.260509, 34.27279, 34.25182, 34.29161, 34.25321),
  Long = c(-86.2074, -86.205116, -86.220818, -86.22271, -86.19384, -86.22192)
  )

playgrounds <- data.frame(
  Location = c("Albertville, AL", "Albertville, AL", "Albertville, AL"),
  Place = c("Graham Park", "Sand Mountain Park", "Pecan Grove Playground"),
  Lat = c(34.26358, 34.275471, 34.27294),
  Long = c(-86.20289, -86.229111, -86.2264)
  )

# Get UTM zone 
floor((schools$Long + 180) / 6) + 1
# [1] 16 16 16 16 16 16
floor((playgrounds$Long + 180) / 6) + 1
# [1] 16 16 16

# Create sf objects and project to correct UTM zone (or whatever CRS is appropriate)
# EPSG codes can be found @ https://epsg.io/
sf_schools <- schools %>%
  st_as_sf(coords = c("Long", "Lat"), crs = 4326) %>%
  st_transform(32616)

sf_playgrounds <- playgrounds %>%
  st_as_sf(coords = c("Long", "Lat"), crs = 4326) %>%
  st_transform(32616)

# School buffers
buff25 <- st_buffer(sf_schools, 400) %>%
  st_join(., sf_playgrounds, join = st_contains) %>%
  group_by(Place.x) %>%
  summarise(pg_count = sum(!is.na(Place.y)), .groups = 'drop') %>%
  mutate(buff = "buff025")

buff50 <- st_buffer(sf_schools, 800) %>%
  st_join(., sf_playgrounds, join = st_contains) %>%
  group_by(Place.x) %>%
  summarise(pg_count = sum(!is.na(Place.y)), .groups = 'drop') %>%
  mutate(buff = "buff050")

buff75 <- st_buffer(sf_schools, 1200) %>%
  st_join(., sf_playgrounds, join = st_contains) %>%
  group_by(Place.x) %>%
  summarise(pg_count = sum(!is.na(Place.y)), .groups = 'drop') %>%
  mutate(buff = "buff075")

buff10 <- st_buffer(sf_schools, 1600) %>%
  st_join(., sf_playgrounds, join = st_contains) %>%
  group_by(Place.x) %>%
  summarise(pg_count = sum(!is.na(Place.y)), .groups = 'drop') %>%
  mutate(buff = "buff100")

# Combine buffers, calculate area as square miles, calculate proximal density,
# convert to dataframe, pivot to wide and reorder buffers from smallest to largest
buff_all <- do.call(rbind, mget(ls(pattern = "^buff\\d{2}$"))) %>%
  mutate(area_sq_mile = as.numeric(st_area(.)) * 3.861e-7,
         pd = pg_count / area_sq_mile) %>%
  rename(Place = Place.x) %>%
  data.frame() %>%
  pivot_wider(id_cols = Place,
              names_from = buff,
              values_from = pd) %>%
  relocate(buff100, .after = last_col())

buff_all
# # A tibble: 6 × 5
#   Place                             buff025 buff050 buff075 buff100
#   <chr>                               <dbl>   <dbl>   <dbl>   <dbl>
# 1 Albertville Elementary School        0       0      0       0    
# 2 Albertville High School              5.16    1.29   0.573   0.322
# 3 Albertville Intermediate School      0       1.29   1.15    0.644
# 4 Albertville Kindergarten and PreK    0       0      0       0    
# 5 Albertville Middle School            0       1.29   0.573   0.322
# 6 Albertville Primary School           0       0      0       0

buff* 列中的值代表每平方英里每所学校周围游乐场的最近密度。