未记录事实核验 · 状态：现行内容

记录仅说明核验时间与维护状态，不保证适用于所有环境；使用前请留意正文中的版本、适用条件与风险提示。

R语言地理空间数据分析

---

一、空间数据的基础概念

1.1 三种基本几何类型

空间数据类型：

点（Point）     → 门店位置、用户地址、POI兴趣点
线（LineString）→ 道路、配送路线、河流
面（Polygon）   → 行政区划、商圈范围、配送区域

1.2 sf包：简单要素的R实现

sf（Simple Features）是处理空间数据的现代R包，遵循OGC标准，与tidyverse完全兼容。

library(sf)
library(tidyverse)

# 从数据框创建sf对象（点数据）
stores <- tibble(
  store_id   = 1:10,
  store_name = paste0("门店", LETTERS[1:10]),
  lng = c(116.40, 116.45, 116.38, 116.42, 116.35,
          116.48, 116.41, 116.44, 116.37, 116.50),  # 经度
  lat = c(39.90,  39.92,  39.88,  39.91,  39.93,
          39.89,  39.94,  39.87,  39.90,  39.92),   # 纬度
  monthly_sales = c(120, 85, 200, 150, 60, 180, 95, 220, 75, 130)
)

# 转换为sf对象
stores_sf <- st_as_sf(stores,
  coords = c("lng", "lat"),  # 先经度后纬度
  crs = 4326                 # WGS84坐标系
)

# sf对象的结构
print(stores_sf)
# 注意：geometry列存储了空间信息，其他列仍可正常用dplyr操作

---

二、坐标系：中国数据的特殊陷阱

坐标系对应关系

坐标系	用途	典型来源
WGS84（EPSG:4326）	全球标准GPS坐标	GPS设备、Google地图（境外）、OpenStreetMap
GCJ-02（火星坐标）	部分互联网地图服务采用的偏移坐标体系，并非国家大地坐标系	高德地图、腾讯地图等，具体以服务接口说明为准
BD-09（百度坐标）	百度自研坐标系	百度地图
CGCS2000（EPSG:4490）	中国国家大地坐标系	政府GIS数据、天地图

# 以下安装示例使用amapGeocode包进行高德地址查询
# remotes::install_github("womeimingzi11/amapGeocode")
library(amapGeocode)

# 地址地理编码示例：查询地址对应坐标，不等同于GCJ-02与WGS84之间的坐标转换
# key需要在高德开放平台申请
# converted <- getCoord("北京市朝阳区望京街道", key = "your_key")

# GCJ-02转WGS84的近似算法示例；“精度约1-5米”仅能作为待验证假设，本示例不提供精度保证
# 以下是GCJ-02 -> WGS84转换的核心算法
gcj02_to_wgs84 <- function(lng, lat) {
  a <- 6378245.0
  ee <- 0.00669342162296594323

  transform_lat <- function(x, y) {
    -100.0 + 2.0*x + 3.0*y + 0.2*y*y + 0.1*x*y + 0.2*sqrt(abs(x)) +
    (20.0*sin(6.0*x*pi) + 20.0*sin(2.0*x*pi)) * 2.0/3.0 +
    (20.0*sin(y*pi) + 40.0*sin(y/3.0*pi)) * 2.0/3.0 +
    (160.0*sin(y/12.0*pi) + 320*sin(y*pi/30.0)) * 2.0/3.0
  }

  transform_lng <- function(x, y) {
    300.0 + x + 2.0*y + 0.1*x*x + 0.1*x*y + 0.1*sqrt(abs(x)) +
    (20.0*sin(6.0*x*pi) + 20.0*sin(2.0*x*pi)) * 2.0/3.0 +
    (20.0*sin(x*pi) + 40.0*sin(x/3.0*pi)) * 2.0/3.0 +
    (150.0*sin(x/12.0*pi) + 300.0*sin(x/30.0*pi)) * 2.0/3.0
  }

  dlat <- transform_lat(lng - 105.0, lat - 35.0)
  dlng <- transform_lng(lng - 105.0, lat - 35.0)

  radlat <- lat / 180.0 * pi
  magic <- sin(radlat)
  magic <- 1 - ee * magic * magic
  sqrtmagic <- sqrt(magic)

  dlat <- (dlat * 180.0) / ((a * (1 - ee)) / (magic * sqrtmagic) * pi)
  dlng <- (dlng * 180.0) / (a / sqrtmagic * cos(radlat) * pi)

  list(lng = lng - dlng, lat = lat - dlat)
}

# 使用示例
result <- gcj02_to_wgs84(116.397128, 39.916527)  # 天安门GCJ-02坐标
cat("WGS84经度:", result$lng, "WGS84纬度:", result$lat, "\n")

---

三、用ggplot2绘制地图

3.1 使用行政区划底图

library(ggplot2)
library(sf)
library(rnaturalearth)  # 全球行政区划数据

# 获取中国省级行政区
china <- ne_states(country = "China", returnclass = "sf")

# 基础地图
ggplot(china) +
  geom_sf(fill = "#F5F5F5", color = "white", size = 0.3) +
  coord_sf(xlim = c(73, 135), ylim = c(18, 54)) +
  theme_void()

# 添加数据层（各省GDP可视化）
set.seed(42)
province_data <- tibble(
  name = unique(china$name_zh),
  gdp  = runif(length(unique(china$name_zh)), 5000, 120000)
) %>% filter(!is.na(name))

china_with_data <- china %>%
  left_join(province_data, by = c("name_zh" = "name"))

ggplot(china_with_data) +
  geom_sf(aes(fill = gdp), color = "white", size = 0.3) +
  scale_fill_viridis_c(
    name = "GDP（亿元）",
    labels = scales::comma,
    option = "plasma"
  ) +
  coord_sf(xlim = c(73, 135), ylim = c(18, 54)) +
  labs(title = "中国各省GDP分布示意图") +
  theme_void() +
  theme(legend.position = "right")

3.2 城市级别的门店分布地图

# 使用leaflet做交互式地图（更适合探索性分析）
library(leaflet)

leaflet(stores_sf) %>%
  addTiles() %>%  # 默认OpenStreetMap底图
  addCircleMarkers(
    radius = ~sqrt(monthly_sales) / 2,  # 大小映射到销售额
    color = ~colorQuantile("RdYlGn", monthly_sales)(monthly_sales),  # 颜色映射
    popup = ~paste0(
      "<b>", store_name, "</b><br>",
      "月销售额：", monthly_sales, "万元"
    ),
    label = ~store_name
  ) %>%
  addLegend(
    position = "bottomright",
    colors = colorQuantile("RdYlGn", stores$monthly_sales)(
      quantile(stores$monthly_sales, seq(0.1, 0.9, 0.2))
    ),
    labels = c("低", "", "中", "", "高"),
    title = "销售额水平"
  )

---

四、空间分析操作

4.1 空间距离计算

# 计算每个门店到市中心的距离
city_center <- st_as_sf(
  data.frame(lng = result$lng, lat = result$lat),  # 使用前文近似转换后的WGS84坐标
  coords = c("lng", "lat"), crs = 4326
)

# 转换为适用于北京附近的投影坐标系，便于以米进行局部距离和缓冲区计算；经纬度下的st_distance也可计算地表距离，并非必然误差大
stores_proj <- st_transform(stores_sf, crs = 32650)  # UTM Zone 50N，北京适用
center_proj <- st_transform(city_center, crs = 32650)

stores_sf <- stores_sf %>%
  mutate(
    dist_to_center = as.numeric(st_distance(stores_proj, center_proj)) / 1000  # 转为公里
  )

4.2 空间缓冲区（商圈范围分析）

# 创建每个门店3公里服务范围
stores_buffer <- stores_proj %>%
  st_buffer(dist = 3000)  # 3000米 = 3公里

# 检查哪些门店的服务范围有重叠（竞争分析）
overlaps <- st_intersects(stores_buffer, sparse = FALSE)
diag(overlaps) <- FALSE  # 排除自身

stores_sf <- stores_sf %>%
  mutate(
    competing_stores = rowSums(overlaps)  # 3公里服务圈与本店服务圈相交的其他门店数，仅作为潜在竞争的几何代理指标
  )

# 可视化
ggplot() +
  geom_sf(data = stores_buffer, aes(fill = stores_sf$monthly_sales),
          alpha = 0.3, color = NA) +
  geom_sf(data = stores_sf, aes(size = monthly_sales, color = competing_stores)) +
  scale_fill_viridis_c(guide = "none") +
  scale_color_gradient(low = "green", high = "red", name = "周边竞争门店数") +
  scale_size_continuous(name = "月销售额（万）") +
  labs(title = "门店服务范围与竞争分析") +
  theme_void()

4.3 空间聚类：门店选址辅助

library(dbscan)

# 使用DBSCAN进行空间聚类（识别客户密集区域）
set.seed(42)
customers <- tibble(
  lng = rnorm(2000, 116.41, 0.05),
  lat = rnorm(2000, 39.91, 0.03)
)

# DBSCAN：基于密度的聚类，不需要预先指定簇数
coords_matrix <- as.matrix(customers[, c("lng", "lat")])
db_result <- dbscan(coords_matrix, eps = 0.01, minPts = 30)

customers <- customers %>%
  mutate(cluster = db_result$cluster)

ggplot(customers, aes(x = lng, y = lat, color = factor(cluster))) +
  geom_point(alpha = 0.5, size = 1) +
  scale_color_brewer(palette = "Set1",
                     labels = c("0" = "噪音点", "1" = "簇1", "2" = "簇2")) +
  labs(title = "客户分布空间聚类",
       subtitle = "DBSCAN算法识别密集客户群",
       color = "聚类结果") +
  theme_minimal() +
  coord_fixed(ratio = 1.5)  # 北京附近的近似显示比例，不是通用投影或精确地面距离校正

---

五、业务案例全流程：外卖配送区域优化

library(sf)
library(tidyverse)
library(leaflet)

# ============================================
# 场景：优化外卖配送区域划分
# 目标：基于模拟历史订单密度识别高密度候选区域；未建模路网、时效和容量约束，不输出最优配送区域
# ============================================

set.seed(123)

# 1. 模拟历史订单数据
orders_geo <- tibble(
  order_id  = 1:5000,
  del_lng   = rnorm(5000, 116.42, 0.03),
  del_lat   = rnorm(5000, 39.91, 0.02),
  delivery_time = rexp(5000, rate = 0.08)  # 配送时间（分钟）
)

# 2. 计算订单热力值（每个0.005度格子的订单密度）
orders_grid <- orders_geo %>%
  mutate(
    lng_bin = round(del_lng / 0.005) * 0.005,
    lat_bin = round(del_lat / 0.005) * 0.005
  ) %>%
  group_by(lng_bin, lat_bin) %>%
  summarise(
    order_count = n(),
    avg_delivery_time = mean(delivery_time),
    .groups = "drop"
  )

# 3. 按订单数的80%分位数筛选高密度格子；“top 20%格子覆盖约80%订单”仅是假设，实际格子比例与订单覆盖率应由数据计算
threshold <- quantile(orders_grid$order_count, 0.8)
core_area <- orders_grid %>%
  filter(order_count >= threshold)

# 4. 可视化：热力图 + 核心区域标注
ggplot() +
  geom_tile(data = orders_grid,
            aes(x = lng_bin, y = lat_bin, fill = order_count),
            alpha = 0.8) +
  geom_point(data = core_area,
             aes(x = lng_bin, y = lat_bin),
             color = "white", size = 0.5, alpha = 0.5) +
  scale_fill_gradient2(low = "#FFFFFF", mid = "#FFC107", high = "#E91E63",
                       midpoint = median(orders_grid$order_count),
                       name = "订单密度") +
  labs(
    title = "外卖配送热力图",
    subtitle = paste0("高密度区（白点标注）覆盖 ",
                     scales::percent(sum(core_area$order_count) / sum(orders_grid$order_count)),
                     " 的订单量"),
    x = "经度", y = "纬度"
  ) +
  theme_minimal() +
  coord_fixed(ratio = 1.5)

---

六、R vs ArcGIS/QGIS 的对比与协作

维度	R（sf/terra）	QGIS（开源）	ArcGIS（商业）
成本	免费	免费	商业许可，费用随产品、许可类型与采购条件变化；“数万/年”仅为未验证的预算假设
可视化	ggplot2，需编码	交互式GUI	专业制图GUI
数据处理	强，可与tidyverse融合	中等	强
可重复性	代码化有助于复现，仍依赖输入数据、环境和随机种子	可通过处理模型、历史记录或脚本辅助复现	有ModelBuilder等复现工具
3D/高级空间分析	需要额外包	插件丰富	功能完整
适合场景	数据驱动的统计+空间分析	快速可视化探索	专业GIS项目

最佳实践：用QGIS做初始空间数据探索和可视化，用R做可重复的统计分析，对接高德/腾讯地图API获取中国地理数据。

---

相关文档

• 09-ggplot2数据可视化深度实战 — 地图可视化的ggplot2基础
• 08-tidyverse数据处理全流程 — 空间数据与业务数据的关联
• 待创建：Python地理空间分析 — GeoPandas/Folium与sf包的对比
• 待创建：商业选址数据分析框架 — 结合POI数据的选址方法论