企业与高校空间分布

R
sf
ggspatial
spatstat
Author

Rui

Published

May 6, 2023

仅以绍兴市越城区的企业与高校数据为例。

导入相关的包

library(tidyverse)
library(sf)
library(geojsonio)
library(ggplot2)
library(ggspatial)
library(spatstat)
library(purrr)

读取数据

qiye <- read.csv("越城区内的企业.csv")
gaoxiao <- read.csv("绍兴高校.csv")
qiye %>% head(8) %>% knitr::kable()
company_type name clean_address longitude latitude
共建研究院 浙江大学绍兴研究院 浙江省绍兴市越城区迪荡湖隧道科学园3号楼26-27层 120.6082 30.01323
共建研究院 天津大学浙江国际创新设计与智造研究院 浙江省绍兴市越城区洋泾湖科创园1号楼 120.6354 30.05344
共建研究院 上海大学绍兴研究院 浙江省绍兴市越城区三江路78号 120.6241 30.10052
共建研究院 浙江工业大学绍兴研究院 浙江省绍兴市越城区洋泾湖科创园2号楼 120.6354 30.05344
共建研究院 江南大学(绍兴)产业技术研究院 浙江省绍兴市越城区洋江东路19号 120.6075 30.04988
共建研究院 中国纺织科学研究院江南分院 浙江省绍兴市越城区双堰路30号 120.6113 30.09832
共建研究院 浙江数字内容研究院 浙江省绍兴市越城区马欢路398号科创园综合楼1-3楼 120.7605 30.12733
共建研究院 浙江理工大学绍兴生物医药研究院 浙江省绍兴市越城区马欢路398号科创园C幢5楼 120.7492 30.13325 
gaoxiao %>% head(8) %>% knitr::kable()
名称 经度 纬度 类型
绍兴文理学院 120.5678 29.98551 本科
浙江越秀外国语学院 120.6024 29.96868 本科
浙江农林大学暨阳学院 120.2545 29.74984 本科
绍兴文理学院元培学院 120.5412 30.06775 本科
浙江工业职业技术学院 120.5791 30.03196 专科
浙江邮电职业技术学院 120.7543 30.14240 专科
浙江农业商贸职业学院 120.6452 30.07056 专科
绍兴职业技术学院 120.5611 29.98086 专科

数据合并

qiye <- qiye %>%
  select(-c(company_type, clean_address)) %>%
  rename("名称" = "name", "经度" = "longitude", "纬度" = "latitude") %>%
  mutate(数据来源 = "企业及研究院")

gaoxiao <- gaoxiao %>%
  select(-类型) %>%
  mutate(数据来源 = "高校")

all_data <- bind_rows(qiye, gaoxiao)

空间分布是否具有随机性

接下来使用 spatstat 包检验企业及研究院的空间分布是否具有随机性。

创建 ppp 对象，绘制点图：

xwindow <- owin(c(range(qiye$经度)), c(range(qiye$纬度)))
points <- ppp(qiye$经度, qiye$纬度, window = xwindow)
Note

owin(c(xmin, xmax), c(ymin, ymax))

plot(points)

density(points) %>% plot(main = "density")

样方分析

teststats <- quadrat.test(points, nx = 6, ny = 6)
teststats
## 
##  Chi-squared test of CSR using quadrat counts
## 
## data:  points
## X2 = 3150.8, df = 35, p-value < 2.2e-16
## alternative hypothesis: two.sided
## 
## Quadrats: 6 by 6 grid of tiles
plot(points, pch = 16, cex = 0.5, main = "样方分析")
plot(teststats, add = T, col = "red")

发现 p 值小于 0.05，拒绝原假设，即原始数据的空间分布非随机。

K 检验

K 函数

接下来，我们可以计算 Ripley’s K(r) 函数（以下简称 K 函数），这是点模式分析中常用的函数之一。 它可以用于检验点的空间分布是否随机。 K 函数描述了在距离 r 内的点对之间的平均距离与随机分布的期望距离之间的差异。 如果点是随机分布的，则 K 函数应该是一个对角线（y = x）。 如果 K 函数在某些范围内高于对角线，则表示点分布在该范围内更密集；如果 K 函数在某些范围内低于对角线，则表示点分布在该范围内更稀疏。

# 计算K函数
K <- Kest(points)

# 绘制K函数图像
plot(K, main="K函数图")

蒙特卡罗模拟

可以使用蒙特卡罗模拟来检验点分布是否随机。蒙特卡罗模拟会生成大量的随机点集，然后计算它们的 K 函数，从而得到 K 函数的置信区间。如果原始数据的 K 函数值在置信区间内，则可以认为点分布是随机的；否则，就有证据表明点分布是非随机的。

# 进行蒙特卡罗模拟，生成999个随机点集
sim <- envelope(points, Kest, nsim = 999)
## Generating 999 simulations of CSR  ...
## 1, 2, 3, ......10.........20.........30.........40.........50.........60........
## .70.........80.........90.........100.........110.........120.........130......
## ...140.........150.........160.........170.........180.........190.........200....
## .....210.........220.........230.........240.........250.........260.........270..
## .......280.........290.........300.........310.........320.........330.........340
## .........350.........360.........370.........380.........390.........400........
## .410.........420.........430.........440.........450.........460.........470......
## ...480.........490.........500.........510.........520.........530.........540....
## .....550.........560.........570.........580.........590.........600.........610..
## .......620.........630.........640.........650.........660.........670.........680
## .........690.........700.........710.........720.........730.........740........
## .750.........760.........770.........780.........790.........800.........810......
## ...820.........830.........840.........850.........860.........870.........880....
## .....890.........900.........910.........920.........930.........940.........950..
## .......960.........970.........980.........990........ 999.
## 
## Done.

# 绘制蒙特卡罗模拟结果图
plot(sim, main="蒙特卡罗模拟结果图")

在绘制的图形中，原始数据的 K 函数曲线应该落在蒙特卡罗模拟生成的置信区间内，如果超出置信区间则说明点的分布是非随机的。

Note

仔细看上图粉色虚线（K 函数模拟值）上下有非常狭窄的灰色部分（置信区间）。上图中原始数据拟合的 K 函数曲线超出置信区间，说明企业及研究院的空间分布不是随机的。

L 函数与 G 函数

# 计算L函数
L <- Lest(points)

# 绘制L函数图像
plot(L, main = "L函数图")

如果点是随机分布的，则 L 函数应该是一条水平线。 如果 L 函数在某些范围内高于水平线，则表示点分布在该范围内更密集；如果 L 函数在某些范围内低于水平线，则表示点分布在该范围内更稀疏。

# 计算G函数
G <- Gest(points)

# 绘制G函数图像
plot(G, main="G函数图")

如果点是随机分布的，则 G 函数应该是一条对角线。 如果 G 函数在某些范围内高于对角线，则表示点分布在该范围内更聚集；如果 G 函数在某些范围内低于对角线，则表示点分布在该范围内更分散。

总体空间可视化

# 将经纬度转换为sf格式的点数据
all_data_sf <- st_as_sf(all_data, coords = c("经度", "纬度"), crs = 4326)
enterprise_sf <- all_data_sf[all_data_sf$数据来源 == "企业及研究院", ]
school_sf <- all_data_sf[all_data_sf$数据来源 == "高校", ]

使用 geojsonio 包读取 json 文件中的地理信息，再使用 sf 包创建空间矢量数据。

# 读取越城区边界geojson文件
geojson <- readLines("越城区.json", warn = FALSE)
yuecheng_boundary <- st_as_sf(st_read(geojson, quiet = TRUE))
# 绘制散点图和面图
ggplot() +
  # 区域
  geom_sf(data = yuecheng_boundary, fill = NA) +
  # 学校点
  geom_sf(data = school_sf, color = "blue",  size = 2) +
  # 企业点
  geom_sf(data = enterprise_sf, color = "red", size = 1, alpha = 0.3) +
  # 设置刻度范围
  coord_sf(xlim = c(120.45, 120.88), ylim = c(29.85, 30.25), expand = FALSE) +
  theme_minimal()

高校缓冲区

# 创建缓冲区
school_buffer_sf <- st_buffer(school_sf, dist = 5000)
# 绘制所有企业点和高校缓冲区
ggplot() +
  # 区域
  geom_sf(data = yuecheng_boundary, fill = NA) +
  # 高校点
  geom_sf(data = school_sf, color = "blue",  size = 2) +
  # 企业点
  geom_sf(data = enterprise_sf, color = "red", size = 1, alpha = 0.3) + 
  # 高校缓冲区
  geom_sf(data = school_buffer_sf, fill = NA, color = "blue", size = 0.5) +
  # 设置地图范围
  coord_sf(xlim = c(120.45, 120.88), ylim = c(29.85, 30.25), expand = FALSE) +
  theme_minimal()

计算每家企业在多少个高校缓冲区内

in_buffer <- st_within(enterprise_sf, school_buffer_sf)
in_buffer
## Sparse geometry binary predicate list of length 2390, where the
## predicate was `within'
## first 10 elements:
##  1: 1, 2, 5
##  2: 7
##  3: 7
##  4: 7
##  5: 5, 7
##  6: 7
##  7: 6
##  8: 6
##  9: 6
##  10: (empty)
num <- map_dbl(in_buffer, length)
num
##    [1] 3 1 1 1 2 1 1 1 1 0 4 3 2 2 2 0 2 2 0 1 0 1 1 1 2 3 0 2 1 1 1 1 1 1 1 1 1
##   [38] 2 1 0 3 1 0 4 0 2 0 1 1 0 1 1 1 2 3 0 2 1 1 1 1 1 1 1 3 2 1 1 0 1 1 1 0 4
##   [75] 0 1 1 1 3 0 1 2 1 2 1 0 1 0 0 2 2 1 2 3 3 0 1 3 0 3 1 0 1 1 0 2 0 3 1 1 0
##  [112] 0 1 2 0 1 1 1 1 1 1 1 1 3 3 1 1 1 3 0 3 1 0 4 0 4 1 0 1 1 1 1 1 3 1 2 1 2
##  [149] 0 1 0 1 0 4 1 1 0 0 3 0 1 0 1 0 2 2 0 2 2 2 2 2 3 1 1 1 1 0 3 3 0 0 0 2 1
##  [186] 2 1 4 1 0 0 3 0 0 0 2 0 1 4 1 0 2 1 2 2 3 0 1 1 1 2 1 0 0 0 0 0 1 1 1 1 1
##  [223] 2 2 1 0 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 4 2 1 0 4 3 3 3 2 1 1 1 3 0 0
##  [260] 0 0 0 0 2 1 1 1 2 3 1 1 0 2 2 3 2 1 0 3 3 1 2 1 4 1 3 4 3 1 3 3 0 4 3 4 0
##  [297] 3 3 3 1 1 2 1 2 3 2 4 3 2 3 1 3 4 0 2 1 3 3 1 3 3 3 3 3 3 3 3 3 2 0 3 3 3
##  [334] 3 3 3 3 3 1 3 3 3 3 3 3 3 3 3 2 3 3 3 4 4 3 3 2 3 3 0 3 1 2 3 2 2 2 3 4 3
##  [371] 3 2 2 1 2 1 2 3 4 3 3 2 3 3 2 3 3 4 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
##  [408] 3 3 3 3 3 2 3 3 3 3 3 3 3 1 4 2 2 3 3 3 3 3 2 3 3 3 3 3 3 3 3 3 0 3 3 3 3
##  [445] 3 3 3 3 3 3 3 3 1 4 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 2 3 3 3 3 3 3 3
##  [482] 3 3 3 3 3 3 3 4 3 4 1 3 4 4 3 4 4 4 4 4 4 4 3 3 1 4 1 4 4 4 2 1 1 3 4 1 1
##  [519] 4 3 3 4 4 3 4 3 3 4 3 1 4 3 4 4 3 4 2 4 3 3 2 4 2 1 2 4 3 2 4 2 1 2 2 3 3
##  [556] 4 4 4 4 4 3 2 2 3 2 3 3 3 3 4 2 4 3 3 1 4 2 1 2 2 1 4 1 2 2 3 4 4 4 1 3 3
##  [593] 3 3 3 3 3 4 3 2 3 4 3 1 3 3 3 4 3 4 3 1 0 0 0 0 1 2 0 1 1 1 1 1 1 3 0 0 0
##  [630] 1 1 1 0 0 0 0 1 0 1 0 0 1 0 0 0 0 0 0 1 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0
##  [667] 0 0 0 0 0 0 0 2 0 1 0 0 0 0 1 2 1 2 1 1 1 1 2 1 1 2 2 0 1 1 1 1 1 1 2 1 1
##  [704] 1 2 1 1 2 1 1 0 1 1 1 1 3 1 1 2 1 1 1 1 0 1 0 1 1 1 0 1 1 0 2 2 1 1 1 1 1
##  [741] 0 1 2 4 3 1 0 1 2 0 1 2 1 2 1 2 1 1 1 1 1 1 1 0 1 1 1 1 2 1 1 2 1 0 2 1 2
##  [778] 2 0 2 1 1 1 0 2 1 1 1 1 2 1 1 0 1 2 1 1 0 1 0 0 2 1 2 2 1 2 1 1 2 2 1 2 1
##  [815] 1 1 2 2 0 2 1 1 2 2 1 1 1 1 1 0 0 2 0 2 1 0 1 0 2 1 1 1 2 0 2 2 1 0 1 1 0
##  [852] 1 1 1 1 2 1 1 1 2 1 3 0 2 2 1 2 0 0 1 1 0 2 1 2 0 1 0 1 0 2 1 2 1 0 2 0 2
##  [889] 1 0 2 1 1 1 1 2 4 1 0 2 2 2 0 2 1 2 1 1 1 1 2 1 1 2 0 1 0 0 2 1 2 1 1 2 1
##  [926] 2 1 1 1 1 2 1 2 2 1 1 1 2 1 1 2 2 1 1 2 1 2 1 0 2 1 1 2 0 1 3 4 4 4 4 4 4
##  [963] 4 4 4 2 4 4 2 3 4 4 4 4 3 2 4 4 4 3 3 3 4 4 3 4 0 0 0 0 1 0 0 0 0 0 0 0 0
## [1000] 0 0 0 0 0 0 0 0 0 0 3 0 0 0 0 2 0 0 2 1 1 0 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0
## [1037] 0 0 0 0 4 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0
## [1074] 1 0 0 0 1 0 0 0 0 0 2 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 2 4 0 0 0 0 0 0 0 0 0
## [1111] 0 0 0 0 0 0 4 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## [1148] 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 4 2 3 2 2 2 4 3 2 1 3 3
## [1185] 2 3 2 2 2 2 2 3 2 3 2 2 2 2 2 2 2 2 3 3 2 2 2 2 3 1 3 3 1 2 1 3 3 3 1 2 2
## [1222] 1 3 2 2 3 3 3 3 1 2 3 2 2 3 3 2 3 3 3 3 2 2 3 3 4 2 1 3 3 2 2 2 1 3 1 3 1
## [1259] 0 2 2 2 2 1 4 2 3 2 3 3 2 3 3 3 3 2 2 2 2 2 3 2 3 2 2 3 0 4 2 2 2 2 2 2 3
## [1296] 3 2 2 1 2 2 2 3 3 2 2 3 3 2 2 1 1 2 2 3 2 2 2 1 2 2 2 3 2 2 1 1 0 0 3 0 0
## [1333] 2 0 1 1 0 0 0 0 1 1 0 1 1 0 0 1 1 1 1 1 0 0 1 1 3 1 1 1 1 1 3 2 1 0 3 2 0
## [1370] 3 3 3 0 0 3 0 1 1 0 1 1 1 1 0 1 1 1 1 1 1 1 0 1 1 1 0 1 3 1 1 1 1 1 1 0 1
## [1407] 1 1 1 1 1 1 0 0 1 0 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 0 1 1 0 0 1 1 1 0 0 1
## [1444] 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 0 1 0 0 0 1 0 1 1 1 0 1 0 0 1 1 1 1
## [1481] 1 0 1 1 0 1 1 1 1 0 0 1 1 1 1 1 1 1 1 1 1 1 0 1 2 2 2 2 1 2 2 2 2 2 1 2 2
## [1518] 2 1 1 2 2 2 2 2 1 1 2 2 2 1 2 3 1 1 2 1 1 2 2 1 2 2 2 2 2 1 4 1 2 1 1 1 2
## [1555] 1 1 1 2 1 1 3 1 1 2 1 1 2 1 2 1 1 1 1 1 1 1 1 2 2 1 2 2 2 2 1 1 1 1 1 2 1
## [1592] 1 1 1 1 1 1 1 2 1 1 2 4 1 1 1 0 1 1 1 2 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 2
## [1629] 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1
## [1666] 1 1 1 1 1 1 1 1 1 4 1 1 1 1 1 1 1 1 1 1 0 1 1 2 1 1 1 1 1 1 1 1 1 1 1 1 1
## [1703] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 0 1 1 1 1 1 1 1 1
## [1740] 1 1 1 1 1 2 1 1 1 1 0 1 1 0 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 0 0 1 1 1 1 1 1
## [1777] 0 0 1 0 1 0 0 0 0 0 1 1 0 2 4 4 3 3 3 4 4 1 4 3 2 4 3 3 4 3 3 4 4 4 4 4 4
## [1814] 4 3 4 4 4 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 4 0 0 0 0 0 0 0
## [1851] 0 3 4 3 3 3 3 3 3 3 4 3 3 3 3 2 2 3 3 3 3 3 3 3 3 3 4 2 1 4 3 1 4 1 2 4 4
## [1888] 3 3 3 3 4 1 4 4 4 0 0 3 4 3 3 4 4 4 4 3 0 0 0 0 0 0 0 1 1 1 0 1 0 1 1 0 0
## [1925] 1 1 1 1 0 1 2 1 1 1 1 4 1 2 1 2 1 1 1 1 2 2 1 2 1 1 1 1 1 2 0 0 1 0 2 1 2
## [1962] 1 2 2 1 1 1 1 0 0 2 0 1 1 1 1 1 1 2 2 1 2 2 2 2 1 1 1 1 1 1 0 1 1 1 1 1 2
## [1999] 1 0 0 0 1 1 1 2 1 1 1 2 1 1 2 0 2 2 2 1 0 2 2 1 1 1 0 1 1 0 2 0 1 2 0 2 2
## [2036] 1 1 4 4 3 3 4 4 4 3 0 1 0 0 0 4 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 4 0
## [2073] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 3 2 3 2 2 2 2 2 1 3 1
## [2110] 2 1 3 2 2 2 3 2 2 3 3 3 2 2 3 2 3 2 3 3 3 3 2 2 2 3 2 3 2 2 3 2 2 2 3 2 2
## [2147] 2 1 2 1 3 3 4 3 3 2 2 2 2 2 2 3 3 0 1 0 0 3 0 0 3 1 1 3 0 1 0 0 0 1 1 0 1
## [2184] 1 0 1 0 1 1 1 0 1 1 1 1 1 1 1 1 1 1 0 0 1 1 0 0 1 1 0 0 1 1 1 1 0 0 0 1 1
## [2221] 1 1 0 1 0 1 1 1 1 1 1 1 1 1 1 2 2 2 1 2 1 1 2 1 2 2 1 1 2 2 1 2 1 1 2 1 2
## [2258] 2 1 2 2 1 1 1 1 2 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [2295] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 1 1 1 1 1 4 1 1 1 1 1 1 1 1 1 1 0 1 1 1
## [2332] 2 0 1 1 1 0 0 2 4 4 4 3 4 3 0 0 4 0 0 0 0 0 0 0 0 0 1 1 1 1 1 0 0 1 1 1 1
## [2369] 1 1 0 3 2 3 4 3 1 1 0 1 1 1 3 1 2 1 3 0 3 1
ent_data <- qiye %>%
  mutate(num = num)
ent_data %>% head() %>% knitr::kable()
名称 经度 纬度 数据来源 num
浙江大学绍兴研究院 120.6082 30.01323 企业及研究院 3
天津大学浙江国际创新设计与智造研究院 120.6354 30.05344 企业及研究院 1
上海大学绍兴研究院 120.6241 30.10052 企业及研究院 1
浙江工业大学绍兴研究院 120.6354 30.05344 企业及研究院 1
江南大学(绍兴)产业技术研究院 120.6075 30.04988 企业及研究院 2
中国纺织科学研究院江南分院 120.6113 30.09832 企业及研究院 1 
ggplot(data = ent_data, aes(x=num)) +
  geom_bar(fill = "royalblue", alpha = 0.8) +
  theme_minimal()

提取高校缓冲区数据

in_buffer_list <- list(NA)
for (i in 1:length(in_buffer)) {
  l <- length(in_buffer[[i]])
  temp <- data.frame(名称 = rep(qiye[i, 1], l), 缓冲区 = in_buffer[[i]])
  in_buffer_list[[i]] <- temp
}

列表转换为数据框：

in_buffer_df <- do.call(rbind, in_buffer_list)

数据联结：

temp <- gaoxiao %>%
  rename("高校名称" = "名称") %>%
  select(高校名称) %>%
  mutate(缓冲区 = c(1:8))

in_buffer_df <- in_buffer_df %>%
  left_join(temp, by = "缓冲区")

in_buffer_df <- in_buffer_df %>% distinct()

长数据转宽数据

wide_df <- tidyr::pivot_wider(
  in_buffer_df, 
  id_cols = 名称, 
  names_from = 高校名称, 
  values_from = 缓冲区, 
  values_fill = NA
)

df <- wide_df %>%
  select_if(is.numeric) %>%
  mutate_all(~ if_else(is.na(.), 0, 1))  
# 对每一列进行操作，将非NA值替换为1，将NA值替换为0
df %>% head() %>% knitr::kable()
绍兴文理学院 浙江越秀外国语学院 浙江工业职业技术学院 浙江农业商贸职业学院 浙江邮电职业技术学院 绍兴职业技术学院 绍兴文理学院元培学院
1 1 1 0 0 0 0
0 0 0 1 0 0 0
0 0 0 1 0 0 0
0 0 0 1 0 0 0
0 0 1 1 0 0 0
0 0 0 1 0 0 0

交集关系

library(ComplexUpset)

upset(
  df, 
  intersect = colnames(df),
  name = "缓冲区", 
  min_size = 5, #至少出现5次
  width_ratio = 0.16, # 左右图比例
  stripes = c("lemonchiffon", "deepskyblue"), # 下方图的背景颜色
  base_annotations = list(
    "Intersection sizes" = intersection_size(counts = FALSE, mapping = aes(fill = "bars_colour")) +
    scale_fill_manual(values = c("bars_colour" = "salmon"), guide = "none") # 需要{ggplot2}，上方柱状图的颜色，并且去除数字
  ),
)