章节概述

本章基于 sf 与 ggplot2 构建空间可视化基础，串联矢量数据结构、坐标参考系（CRS）与核心几何运算（如空间连接与裁剪），打通从数据清洗到制图输出的完整工作流。

实践层面，本章聚焦计算社会科学中高频的分层设色图（Choropleth Map）与多图层综合制图。重点解析连续数据离散化的算法策略（如自然断点法），并引入网格聚合等进阶技巧，以科学化解空间要素的视觉拥挤。

规范层面，本章归纳了地图辅助要素的排版美学，并着重划定中国公开地图出图的合规红线与标准画法，旨在输出兼具科学严谨、政治合规与发表友好的专业出版级图件。

本章学习内容

学习导航

本章学习内容

欢迎开始 第 4 章 的学习。请点击 下方导航卡 进入相应小节：

空间数据导论

空间要素与属性

空间数据读写

sf 框架与 I/O

空间数据处理

投影与几何运算

空间数据制图

分层设色与聚合

制图规范

规范与排版

章节练习

复现与挑战

💡 提示：学习完一个小节后，请再次点击 屏幕右下角的章节主页按钮回到本导航页

1.0 空间数据

1.0 空间数据导论

基本概念

本章进入空间可视化（geo-visualisation），也常被称为空间制图（spatial mapping）；在更广义的语境中，它与制图学（cartography）、地理信息科学（GIScience）及地理信息系统（GIS）等领域紧密相关。

学科拓展空间数据与空间制图覆盖面很广，本章仅聚焦于计算社会科学研究中最常见的空间数据概念与基于 R 的基础制图流程；若你对上述方向有兴趣，可进一步系统学习 GIS/GIScience 与空间分析方法。

与第 3 章主要处理的非空间数据（non-spatial data）不同，本章关注的数据包含明确的空间位置或空间形状信息，即空间数据（spatial data）。Figure 1 展示了常见的空间数据的核心概念。

Figure 1. 空间数据概念

简单概括：空间数据 = 几何（在哪里/形态） + 属性（是什么）。

1.1 空间数据主要类型

常见的空间数据（spatial data）通常分为两大类：

矢量数据（vector）：以点、线、面等几何对象来表达空间实体，适合描述离散对象（discrete features）的形状、边界与相互位置关系；
栅格数据（raster）：以规则网格单元（像素）组织空间信息，更适合刻画连续表面（continuous surfaces）及其空间梯度变化。

Figure 2. 矢量数据与栅格数据

本章主要面向矢量数据（vector）的空间可视化与基础处理流程。这是因为在计算社会科学的常见研究对象中，分析单元多为离散空间实体（discrete units）（如行政区、社区、学校、POI、道路网络等），其数据结构天然更适合用点、线、面等矢量几何来表达，并便于与人口普查、调查数据等表格型属性进行连接与解释。

栅格数据（raster）同样是重要的空间数据类型，但它更常用于表达连续空间表面（continuous surfaces）（如遥感影像、温度、降水、夜间灯光、PM2.5 网格、房价插值图等），其数据结构与制图逻辑与矢量不同。由于本课程主要面向计算社会科学应用场景，栅格在本课程中不是核心重点，因此我们仅在必要处进行简要说明，不做系统展开。

需要强调的是：在计算社会科学研究中，栅格数据依然有重要用途，例如用夜间灯光或遥感土地覆盖刻画区域发展水平、用网格化环境暴露（如空气污染、热暴露）解释健康或行为差异等。

在 R 的生态中，目前通常使用 sf 包来处理矢量空间数据，并使用 ggplot2 + geom_sf() 来完成可视化。
（注：其他空间可视化包也会在后续课程中提及）

1.2 空间数据的核心结构

基础 · 概念

如果用“表格”的视角来理解数据结构，那么空间数据与非空间数据最直观的区别在于：
空间数据会多出一列专门的几何列（geometry column），用于记录对象的空间形状与位置。

属性（attributes）：
类似 csv 的普通列，用于描述对象“是什么”
（例如区域名称、人口、收入、病例数等）
几何（geometry）：
用于描述对象“在哪里/长什么样”
（例如点的位置、道路的线形、行政区的边界面）

在 sf 对象中，几何通常以一列 geometry 形式存在，并与属性表按行一一对应（即“一行 = 一个空间对象”）

因此它仍然可以像 dataframe 一样进行 filter()、mutate()、select() 等 tidyverse 操作。

用示例数据理解 `geometry` 列

基础 · 示例

下面使用 sf 包自带的示例数据 nc（North Carolina counties）来直观看空间数据的表格结构。

# install.packages("sf")
library(sf)

# 1) 读取 sf 自带示例数据：North Carolina counties
# st_read() 是sf中读取空间数据的核心函数

nc <- st_read(system.file("shape/nc.shp", package = "sf"), quiet = TRUE)

# 2) 展示 nc 空间数据
head(nc)

Simple feature collection with 6 features and 14 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -81.74107 ymin: 36.07282 xmax: -75.77316 ymax: 36.58965
Geodetic CRS:  NAD27
   AREA PERIMETER CNTY_ CNTY_ID        NAME  FIPS FIPSNO CRESS_ID BIR74 SID74
1 0.114     1.442  1825    1825        Ashe 37009  37009        5  1091     1
2 0.061     1.231  1827    1827   Alleghany 37005  37005        3   487     0
3 0.143     1.630  1828    1828       Surry 37171  37171       86  3188     5
4 0.070     2.968  1831    1831   Currituck 37053  37053       27   508     1
5 0.153     2.206  1832    1832 Northampton 37131  37131       66  1421     9
6 0.097     1.670  1833    1833    Hertford 37091  37091       46  1452     7
  NWBIR74 BIR79 SID79 NWBIR79                       geometry
1      10  1364     0      19 MULTIPOLYGON (((-81.47276 3...
2      10   542     3      12 MULTIPOLYGON (((-81.23989 3...
3     208  3616     6     260 MULTIPOLYGON (((-80.45634 3...
4     123   830     2     145 MULTIPOLYGON (((-76.00897 3...
5    1066  1606     3    1197 MULTIPOLYGON (((-77.21767 3...
6     954  1838     5    1237 MULTIPOLYGON (((-76.74506 3...

关于示例数据 `nc`

nc 是 sf 自带的示例空间数据，来源于美国 North Carolina（北卡罗来纳州）的县级行政区边界。

每一行对应一个县的面要素（areal feature）：
- POLYGON 表示该县由一个连续的面构成；
- MULTIPOLYGON 表示该县由多个彼此分离的面共同构成（例如包含小岛或飞地）。
普通列是该县的属性信息（例如名称、代码、统计指标等）。
geometry 列记录该县边界的几何形状与地理位置坐标，用于地图绘制与空间计算。

注：此数据主要用于教学演示，并不是本章必须使用的研究数据

逐行读 `head(nc)` 的输出

head(nc) 的输出通常包含两部分：对象概览信息与属性表（含 geometry）。

1. 对象概览信息

你会看到类似下面这些行（不同版本的 sf 文字可能略有差异，但含义一致）：

Simple feature collection with 6 features and 14 fields：
表示当前打印的是一个简单要素集合（simple feature collection），包含 6 个空间对象（因为使用了 head() 仅展示前 6 行），并带有 14 个属性字段（即除 geometry 之外的“普通列”）。
Geometry type: MULTIPOLYGON：
表示几何类型为多面（MULTIPOLYGON）。它仍属于“面要素”，一个对象可能由多个彼此分离的面共同构成（例如包含小岛或飞地）。
Bounding box: xmin ... xmax ... ymin ... ymax ...：
表示该数据（当前展示部分）的空间范围边界框。
术语 bounding box：范围边界框，通常是一个四边形，用四个数值概括数据的最小外接矩形范围：xmin/xmax 对应左右边界，ymin/ymax 对应下上边界。
直观理解：就像给一张照片画一个“外框”，外框刚好把所有对象包住
Geodetic CRS: NAD27：
表示该数据采用的坐标参考系（CRS）是 NAD27（North American Datum 1927），属于一种地理坐标系统（geodetic / geographic CRS）。后续章节会专门解释 CRS 。
术语 Geodetic：指“基于椭球体/大地基准”的坐标表达方式，通常以经纬度或类似的地理坐标体系为基础（与投影坐标的概念相关）。

2. 属性表

接下来这一段看起来就像普通 data.frame 或 tibble：每行对应一个县（一个空间对象），每列对应一个属性字段。
关键区别在于：它多了一列 geometry，从而成为空间数据（spatial data）。

geometry：
这是几何列，打印时通常会被简写为 MULTIPOLYGON (((...)))。
你不需要在这里读懂每个坐标点的数值，只需知道：geometry 保存了用于绘图与空间计算的边界坐标与几何结构信息（例如边界点序列、面与面之间的组合关系等）。

常见误解

在 head(nc) 中看到的 AREA、PERIMETER 等字段只是该示例数据自带的属性列，并不等同于 sf 的几何信息本身。
真正决定“画出来是什么形状”的，是 geometry 列。

1.3 矢量空间数据的几何类型

基础 · 概念

在矢量空间数据中，几何信息通常以三类基本对象表示：点、线、面。

它们也可以从 “几何维度” 来理解：
点是 0 维（仅表示位置），线是 1 维（沿某一方向延展的路径，具有长度），面是 2 维（具有封闭边界的区域，具有边界与面积）。

同时，它们在空间数据中也有典型的“存储方式”（GIS领域的概念）：

点：通常由一对坐标表示，如 (x, y)（如经纬度 lon, lat）；
线：由一串按顺序连接的点构成，即 P1 → P2 → ... → Pn（最简单情形是两点连线）；
面：由闭合的边界点序列构成，即起点与终点相同，形成封闭边界
注：在复杂情形下还可能包含“空洞”或多个分离部分。

在数据可视化的概念中，这些差异直接决定了地图表达时应选择的视觉编码方式：例如点常用大小（size）与形状（shape），线常用线宽（linewidth）与线型（linetype），面则常用填充色（fill）来表达区域指标。

点（Point）：位置型对象

点（point）用于表示可被视作“一个位置”的空间对象，例如：GPS 定位点、医院、地铁站、事故发生点、调查样本点等。
在数据结构上，点要素通常由一对坐标表示（x, y）（或经纬度 lon, lat）；在地图表达上，点图常用以下视觉编码：

size：表达数量或强度（例如门诊量、病例数、访问量）；
colour / shape：表达类别差异（例如医院等级、站点类型、事件类别）。

下面我们使用 spData 包中的示例数据 cycle_hire_osm 来展示一个典型的点数据（共享单车租赁站点），并进行一次最基本的点地图可视化。

library(sf)
library(spData) # 需install
library(ggplot2)

# ------------------------------------------------------------
# 示例点数据：cycle_hire_osm
# 说明：这是 OpenStreetMap (OSM) 提取的“共享单车租赁站点”点数据，
# 区域：英国伦敦（Transport for London / Santander Cycles 所在区域的站点）
# 用途：教学演示 sf 的 POINT 几何与基础点地图绘制
# ------------------------------------------------------------
data("cycle_hire_osm", package = "spData")

head(cycle_hire_osm)    # 简单展示

Simple feature collection with 6 features and 5 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: -0.1293092 ymin: 51.52583 xmax: -0.090836 ymax: 51.53402
Geodetic CRS:  WGS 84
     osm_id                       name capacity cyclestreets_id description
1    108539            Windsor Terrace       14            <NA>        <NA>
2 598093293 Pancras Road, King's Cross       NA            <NA>        <NA>
3 772536185 Clerkenwell, Ampton Street       11            <NA>        <NA>
4 772541878                       <NA>       NA            <NA>        <NA>
5 781506147                       <NA>       NA            <NA>        <NA>
6 783824668      Finsbury Library, EC1       NA            <NA>        <NA>
                     geometry
1 POINT (-0.0933878 51.52913)
2 POINT (-0.1293092 51.53402)
3 POINT (-0.1182352 51.52729)
4  POINT (-0.090836 51.52583)
5 POINT (-0.1210572 51.53001)
6 POINT (-0.1038272 51.52594)

ggplot(cycle_hire_osm) + # 简单可视化展示
  geom_sf(alpha = 0.7, size = 0.8) +
  labs(
    title = "Cycle Hire Stations (London, UK)",
    subtitle = "Example sf POINT data from spData (OSM-derived)",
    caption = "Data: spData::cycle_hire_osm"
  ) +
  theme_minimal()

线（Line）：连接与路径型对象

线（line / linestring）用于表示连接关系或路径，常见于道路网络、河流、公交线路、迁移路径、通勤线路等。
在线几何中，一个线对象通常由一串按顺序连接的点构成，可写作 P1 → P2 → ... → Pn（最简单情形是两点连线：线段）。

在地图表达上，线图常用以下视觉编码：

linewidth：表达流量或强度（例如交通流量、OD 规模、道路容量）；
linetype / colour：表达类别或等级（例如道路等级、线路类型、河流等级）。

下面我们使用 spData 包中的示例数据 seine 来展示一个典型的线数据：它刻画了法国塞纳河（Seine）及其主要支流（如 Marne、Yonne）的河网线要素，并以此进行一次最基本的线地图可视化。

library(sf)
library(spData)
library(ggplot2)

# ------------------------------------------------------------
# 示例线数据：seine
# 说明：法国塞纳河（Seine）及其支流（Marne, Yonne）的河网线数据
# 数据类型：MULTILINESTRING（线要素）
# 用途：教学演示 sf 的 LINE 几何与基础线地图绘制
# ------------------------------------------------------------
data("seine", package = "spData")

# 快速查看
head(seine)

Simple feature collection with 3 features and 1 field
Geometry type: MULTILINESTRING
Dimension:     XY
Bounding box:  xmin: 518344.7 ymin: 6660431 xmax: 879955.3 ymax: 6938864
Projected CRS: RGF93 / Lambert-93
   name                       geometry
1 Marne MULTILINESTRING ((879955.3 ...
2 Seine MULTILINESTRING ((828893.6 ...
3 Yonne MULTILINESTRING ((773482.1 ...

# 简单可视化
ggplot(seine) +
  geom_sf(colour = "#1F78B4", linewidth = 0.7, alpha = 0.9) +
  labs(
    title = "River Network Lines in France (Seine, Marne, Yonne)",
    subtitle = "Example sf LINE data from spData: seine",
    caption = "Data: spData::seine"
  ) +
  theme_minimal()

面（Polygon）：边界与区域型对象

面（polygon / multipolygon）用于表示具有边界与面积的空间对象，例如：行政区、社区边界、土地利用地块、建筑轮廓（footprint）、洪泛区等。
在数据结构上，面要素通常由闭合的边界点序列构成，即边界坐标按顺序连接并最终回到起点（起点与终点相同），从而形成一个封闭区域。

面要素是空间制图中最常见的“底图结构”，尤其适合制作分层设色图（choropleth map）：
- 用 fill 表达区域指标（例如失业率、贫困率、人口密度）；
- 用边界线强调行政区轮廓

下面我们使用 spData 包中的示例数据 world 绘制一幅最基本的世界面边界图，以直观理解“面数据”的结构与地图呈现方式。

library(sf)
library(spData)
library(ggplot2)

# 示例面数据：全球各国家/地区边界（POLYGON / MULTIPOLYGON）
data("world", package = "spData")

head(world)

## Simple feature collection with 6 features and 10 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -180 ymin: -18.28799 xmax: 180 ymax: 83.23324
## Geodetic CRS:  WGS 84
## # A tibble: 6 × 11
##   iso_a2 name_long  continent region_un subregion type  area_km2     pop lifeExp
##   <chr>  <chr>      <chr>     <chr>     <chr>     <chr>    <dbl>   <dbl>   <dbl>
## 1 FJ     Fiji       Oceania   Oceania   Melanesia Sove…   1.93e4  8.86e5    70.0
## 2 TZ     Tanzania   Africa    Africa    Eastern … Sove…   9.33e5  5.22e7    64.2
## 3 EH     Western S… Africa    Africa    Northern… Inde…   9.63e4 NA         NA  
## 4 CA     Canada     North Am… Americas  Northern… Sove…   1.00e7  3.55e7    82.0
## 5 US     United St… North Am… Americas  Northern… Coun…   9.51e6  3.19e8    78.8
## 6 KZ     Kazakhstan Asia      Asia      Central … Sove…   2.73e6  1.73e7    71.6
## # ℹ 2 more variables: gdpPercap <dbl>, geom <MULTIPOLYGON [°]>

ggplot(world) +
  geom_sf(linewidth = 0.15) +
  labs(
    title = "World Countries (Polygon / MultiPolygon)",
    subtitle = "Example sf polygon data from spData: world",
    caption = "Data: spData::world"
  ) +
  theme_minimal()

在城市形态与规划研究中，建筑轮廓与地块边界等面数据也常用于绘制“实—虚”图（figure–ground），用于突出建成环境与开放空间的结构关系。

1.4 坐标参考系（Coordinate Reference System, CRS）

在空间数据中，坐标参考系（Coordinate Reference System, CRS）是一个基础且关键的概念。
空间数据的 geometry 本质上是一组坐标；只有在明确“这些坐标属于哪一种坐标系”时，它们才具有可解释的空间意义（例如位置、方向、距离、面积，以及对象之间的相对空间关系/拓扑关系）。

现实中存在多种 CRS，它们基于不同的表达方式与应用需求来定义坐标。
因此在空间制图中，尤其当你需要同时使用多份空间数据（例如底图 + 点位 + 线网）时，必须检查它们是否处于同一 CRS。只有在统一坐标参考系的前提下，不同数据的叠加、比较与进一步计算（如距离、缓冲区、面积）才是可靠的。

CRS 的两种常见类型：地理坐标系 vs 投影坐标系

CRS 在实际研究中主要会遇到两类（参考 Figure 3）：

地理坐标系（Geographic CRS / GCS）：
以经度与纬度来表达位置，数值单位通常是度（可用十进制度或度分秒表示）。这种坐标便于全球定位，但在距离、面积等度量计算上需要谨慎，因为“度”并不是线性距离单位。
投影坐标系（Projected CRS / PCS）：
通过某种地图投影把地球表面转换到平面坐标系，坐标单位通常是米（有时也可能是千米或英尺/foot）。这类坐标更适合进行距离、面积、缓冲区等空间度量与制图表达。

Figure 3. 地理坐标系 vs 投影坐标系

在可视化层面，你会观察到一个直观现象：经纬网（graticule）在图上是“直的还是弯的”，与地图采用的 CRS/投影方式密切相关。
更具体地说，经纬线在不同投影下会发生不同的几何变换，因此其在图上的形态也会随之改变。

当地图采用经纬度的“直绘”展示（将经度近似作为 x、纬度近似作为 y）时，经纬网通常呈现为规则的水平/竖直线条；
当地图采用其他世界投影（如 Robinson 等）时，经纬网往往会出现弯曲或收敛的形态，这是投影变换的正常结果。

下面用一段示例进行对比：使用同一份世界底图，分别以 WGS84（EPSG:4326）与 Robinson 投影绘制，并将两幅图上下排列，直观看到经纬网的差异。

# install.packages("ggmapcn")
# install.packages("patchwork")

library(ggplot2)
library(sf)
library(ggmapcn)
library(patchwork)

# ============================================================
# 对比目的：同一份世界底图，在不同 CRS/投影下绘制，经纬网形态会不同
# - 图 1：WGS84 (EPSG:4326) —— 经纬度“直绘”，经纬网通常呈水平/竖直正交
# - 图 2：Robinson 投影      —— 世界制图常用投影，经纬网会发生弯曲（投影效应）
# ============================================================

# ------------------------------------------------------------
# CRS 1: WGS84（EPSG:4326）
# 说明：这是最常见的地理坐标系（经纬度）。这里用它作为“直绘对照组”。
# ------------------------------------------------------------
p_wgs84 <- ggplot() +
  geom_world() +
  annotation_graticule(
    lon_step     = 60,  # 经线间隔（单位：度）
    lat_step     = 30,  # 纬线间隔（单位：度）
    label_offset = 5    # 标签偏移
  ) +
  coord_sf(crs = 4326) +
  labs(title = "WGS84 (EPSG:4326)") +
  theme_void() +
  theme(panel.ontop = TRUE)

# ------------------------------------------------------------
# CRS 2: Robinson（Robinson projection）
# 说明：Robinson 是常见的世界地图投影，用于更“视觉均衡”的全球展示。
# 关键点：同一世界底图在该投影下会发生几何变换，因此经纬网会弯曲/收敛。
# ------------------------------------------------------------
crs_robin <- "+proj=robin +datum=WGS84"

p_robin <- ggplot() +
  geom_world(crs = crs_robin) +
  annotation_graticule(
    crs          = crs_robin,
    lon_step     = 30,   # 投影图下可用更密的经纬网
    lat_step     = 15,
    label_offset = 3e5   
  ) +
  coord_sf(crs = crs_robin) +
  labs(title = "Robinson projection") +
  theme_void() +
  theme(panel.ontop = TRUE)

# ------------------------------------------------------------
# 上下对比（两行一列）
# ------------------------------------------------------------
p_wgs84 / p_robin

为什么要统一 CRS？以及如何通过 `sf` 统一

重点

当你需要把两份空间数据叠加在同一张地图上（例如“底图面 + 点位/线网”），它们必须处在同一坐标参考系（CRS）下，才能实现正确的空间对齐。
更重要的是：许多空间操作（如空间连接、距离/面积计算、缓冲区）都隐含地依赖 CRS 的定义；若 CRS 不一致，结果往往不可比较，甚至没有意义。

在绘图或空间分析之前，先 st_crs() 检查两份数据的 CRS，然后用 st_transform() 统一到同一个 CRS。

方法1 以“目标 CRS”为标准，转换两份数据

library(sf)
library(spData)

# 示例：面数据（world） + 线数据（seine）
data("world", package = "spData")
data("seine", package = "spData")

# 1) 先检查 CRS
st_crs(world)
st_crs(seine)

# 2) 设定一个目标 CRS（例如 WGS84，经纬度 EPSG:4326）
crs_target <- st_crs(4326)

# 3) 同时转换
world_4326 <- st_transform(world, crs_target)
seine_4326 <- st_transform(seine, crs_target)

# 4) 再次确认：两份数据 CRS 应一致
st_crs(world_4326)
st_crs(seine_4326)

方法2 以其中一个数据为“标杆”，把另一个转换到它的 CRS

library(sf)
library(spData)

data("world", package = "spData")
data("seine", package = "spData")

# 1) 以 world 的 CRS 作为标杆
crs_ref <- st_crs(world)

# 2) 把 seine 转换到 world 的 CRS
seine_to_world <- st_transform(seine, crs_ref)

# 3) 检查：两者应一致
st_crs(world)
st_crs(seine_to_world)

常见 CRS 的例子

常用

不同国家/地区在制图与空间分析中常用的 CRS 并不相同（尤其是投影坐标系）。
现阶段的重点不是记参数，而是建立“编号—名称—用途”的基本对应关系，并能够用 sf 进行查看与转换。

全球通用（地理坐标系）：EPSG:4326 = WGS 84（经纬度表达）
英国常用（投影坐标系）：EPSG:27700 = OSGB36 / British National Grid (BNG)（英国国家网格）
中国常用（地理坐标系）：EPSG:4490 = CGCS2000（经纬度表达）
- 历史基准仍可能出现：EPSG:4610 = Xian 1980；EPSG:4214 = Beijing 1954
美国常见情况：常见基准包括 NAD83 / NAD27，并在实践中大量使用各类 State Plane 或 UTM 等投影系统（不同州/分区对应不同 EPSG 代码）。

如需其他国家或地区的常用坐标系，可通过询问AI或查阅资料

本节小结：课程定位与进一步学习

空间数据与空间制图本身具有完备的方法体系，能够支撑起独立的学科与课程体系。
从学术结构上看，它们与制图学、地理信息科学、地理信息系统、空间分析等方向紧密相连，也是典型的交叉学科知识框架之一。

许多更深层的内容（如投影原理、测度误差、拓扑规则、制图表达理论、空间统计与模型等）需要在后续学习中练习老师或进行系统补充。
若对空间方向有进一步兴趣，也建议自行查阅相关专业教材与资料，并结合实际问题逐步拓展。

2.0 空间数据读写

上一节已经明确：矢量空间数据本质上是“属性表 + 几何列（geometry）”。

本节聚焦于建立一个基于 sf 与 tidyverse 的空间数据读写工作流，目标是将不同来源的数据整理为可直接用于制图的 sf 输入，并支持后续复用与共享。该工作流主要包含以下四类操作：

读入：用 sf::st_read() 读取常见矢量格式（如 shp/geojson/gpkg 等），并保持为 sf 对象；
转换：将表格型数据（如 .csv 中的 x/y 或 lon/lat 坐标列）转换为空间点数据（sf::st_as_sf()），并显式设置 CRS；
连接：将外部表格数据与空间数据按键连接（attribute join），形成带有制图指标的“空间属性表”；
写出：将处理结果导出为常用空间格式（如 gpkg/geojson），以便复用、共享与后续制图。

2.1 常见矢量数据格式

基础 · 概念

矢量空间数据存在多种存储格式。结合 sf 的常见读写需求，以下格式需要建立基本认识：

shp（Shapefile）：传统且仍被广泛使用的矢量格式；通常由一组同名文件共同构成（如 .shp/.shx/.dbf/.prj 等）。

注意迁移该类型文件需要将其整体一并处理，否则可能出现属性缺失或 CRS 信息丢失。
geojson：基于 JSON 的开放矢量格式，便于跨平台交换与网页端展示；结构直观，但在要素规模较大时文件体积可能较大。
gpkg（GeoPackage）：现代的单文件矢量格式，可在一个文件中容纳多个图层与属性表，适合项目组织与可复现分析。
csv + 坐标列：严格来说并非空间文件格式，但在实践中常见（例如点数据以 lon/lat 或 x/y 两列记录）；读入后可用 st_as_sf() 转换为 sf 点对象。
json：并非典型空间文件格式，但在数据接口与网页数据中较常见；若包含可解析的坐标或几何字段（如 lon/lat 或坐标数组），可在整理后转换为 sf 对象。
xml（常见为 GML / KML）：基于 XML 的空间数据家族，常用于标准化交换与可视化共享。
- gml（Geography Markup Language）：OGC 标准之一，常见于机构数据发布；结构规范但相对“重”。
- kml/kmz：Google Earth 生态常用格式；kmz 为压缩形式。

Figure 4. 常见空间矢量数据格式

2.2 使用 `sf` 读写矢量空间数据

基础 · 操作

在 sf 体系中，矢量空间数据的读写主要由两类函数完成：
- st_read()：负责从文件系统读入多种常见矢量格式；
- st_write()：负责将处理后的 sf 对象写出为指定空间格式，便于复用与共享。

其中，st_read() 基于 GDAL 驱动，能够在同一套接口下读入多种格式（如 shp/geojson/gpkg/kml/gml 等），并返回标准 sf 对象（属性表 + geometry）。

读写 `.shp` 文件

常用

本节以 sf 包内置的 Shapefile 示例数据 nc 作为演示对象。该数据以 .shp 形式存储，可用于展示 Shapefile 的读入与写出流程：先用 st_read() 读入为 sf 对象，再用 st_write() 将其写出为新的 .shp 文件。

st_read() 的输入是文件路径；读入后返回标准 sf 对象（属性表 + geometry）。
st_write() 将 sf 对象写出到指定路径；若目标文件已存在，可用 delete_dsn = TRUE 覆盖写出（谨慎使用）。

library(sf)

# ------------------------------------------------------------
# 1) 读入：sf 内置示例 Shapefile（nc）
# ------------------------------------------------------------
path_nc <- system.file("shape/nc.shp", package = "sf")
path_nc  # 查看该 .shp 在本机的存储位置

nc <- st_read(path_nc, quiet = TRUE)

# 最基本检查
head(nc)

# ------------------------------------------------------------
# 2) 写出：将 sf 对象导出为新的 .shp 文件
#    说明：Shapefile 会生成一组同名文件（.shp/.shx/.dbf/.prj 等）
# ------------------------------------------------------------

# 这里使用了临时文件，实践阶段改为需要储存的文件地址
out_dir  <- tempdir() 
# 输出的文件名是 nc_out.shp；.shp是格式后缀
out_shp  <- file.path(out_dir, "nc_out.shp")

st_write(nc, out_shp, delete_dsn = TRUE, quiet = TRUE)

# 验证：再读入一次，确认写出成功
nc_out <- st_read(out_shp, quiet = TRUE)
head(nc_out)

读写 `.gpkg` 文件

本节以 spData 包内置的 GeoPackage 示例文件 world.gpkg 作为演示对象。gpkg 属于现代的单文件矢量格式，其特点是：一个文件可包含多个图层（layer）（可理解为同一文件中包含多个数据层）。因此在读入前通常先用 st_layers() 查看图层信息，再用 st_read(layer = ...) 指定目标图层。
写出时可用 st_write() 将处理后的 sf 对象导出为新的 .gpkg 文件，并可通过 layer = 指定写出的图层名称。

library(sf)
library(spData)

# ------------------------------------------------------------
# 1) 读入：spData 内置示例 GeoPackage（world.gpkg）
# ------------------------------------------------------------
path_gpkg <- system.file("shapes/world.gpkg", package = "spData")
path_gpkg  # 查看该 .gpkg 在本机的存储位置

# 查看图层信息（名称、要素数、几何类型等）
st_layers(path_gpkg)  # 此示例通常仅包含一个图层

# 读取指定图层（layer 名称以 st_layers() 的输出为准）
world_gpkg <- st_read(path_gpkg, layer = "world", quiet = TRUE)

# 最基本检查
head(world_gpkg)

# ------------------------------------------------------------
# 2) 写出：将 sf 对象导出为新的 .gpkg 文件（可指定 layer 名称）
# ------------------------------------------------------------
out_dir  <- tempdir()
out_gpkg <- file.path(out_dir, "world_out.gpkg")

st_write(world_gpkg, out_gpkg, layer = "world_out", delete_dsn = TRUE, quiet = TRUE)

# 验证：查看写出文件的图层，并读入确认
st_layers(out_gpkg)
world_out <- st_read(out_gpkg, layer = "world_out", quiet = TRUE)
head(world_out)

读写 `.geojson` 文件

重点 · 常用

本节以 spData 包内置的 GeoJSON 示例数据 cycle_hire_osm.geojson 作为演示对象。该数据以 .geojson 格式存储，可用于展示 GeoJSON 的读入与写出流程：先用 st_read() 读入为 sf 对象，再用 st_write() 写出为新的 .geojson 文件。

library(sf)
library(spData)

# ------------------------------------------------------------
# 1) 读入：spData 内置示例 GeoJSON（cycle_hire_osm.geojson）
# ------------------------------------------------------------
path_geojson <- system.file("shapes/cycle_hire_osm.geojson", package = "spData")
path_geojson  # 查看该 .geojson 在本机的存储位置

cycle_hire_osm_gj <- st_read(path_geojson, quiet = TRUE)

# 最基本检查
head(cycle_hire_osm_gj)

# ------------------------------------------------------------
# 2) 写出：将 sf 对象导出为新的 .geojson 文件
# ------------------------------------------------------------
out_dir     <- tempdir()
out_geojson <- file.path(out_dir, "cycle_hire_osm_out.geojson")

st_write(cycle_hire_osm_gj, out_geojson, delete_dsn = TRUE, quiet = TRUE)

# 验证：再读入一次，确认写出成功
cycle_out <- st_read(out_geojson, quiet = TRUE)
head(cycle_out)

从 `csv` 到 `sf`：`read_csv()` + `st_as_sf()` + 写出

重点 · 常用

点位数据在实践中常以表格形式（如 .csv）提供：每一行对应一个对象，同时包含两列坐标信息（如 lon/lat 或 x/y）。
此类数据的典型工作流包括三步：
1) 读入为 tibble（属性表）；
2) 用 sf::st_as_sf() 将坐标列转换为点几何，并显式指定 CRS；
3) 视需求进行 CRS 转换，并将结果写出为空间文件格式，便于复用与共享。

下例以 spData 的示例点数据 cycle_hire_osm 为参照构造一个“表格型点数据”用于演示：先提取坐标列形成 tibble（模拟 .csv 读入后的结构），再转换为 sf 点对象；随后将 WGS84（EPSG:4326）转换到英国常用的 BNG（EPSG:27700），并将结果写出为 GeoPackage 文件。

library(sf)
library(spData)
library(tidyverse)

# ------------------------------------------------------------
# 示例参照：spData::cycle_hire_osm（sf 点数据）
# 目标：演示 “表格坐标列 → sf 点对象 → CRS 转换 → 写出”
# ------------------------------------------------------------
data("cycle_hire_osm", package = "spData")

# 1) 将 sf 点数据“还原”为表格：提取 lon/lat 坐标列 + 属性列（模拟 csv 读入结果）
coords <- st_coordinates(cycle_hire_osm)   # 提取点坐标矩阵（X/Y）

cycle_tbl <- cycle_hire_osm %>%
  st_drop_geometry() %>%                  # 去掉 geometry，仅保留普通属性列
  mutate(
    lon = coords[, 1],                    # X = 经度（lon）
    lat = coords[, 2]                     # Y = 纬度（lat）
  )

head(cycle_tbl)                           # 表格型点数据：属性列 + lon/lat 坐标列

# 2) 表格 → sf 点对象：指定坐标列名 + 指定 CRS（WGS84，经纬度）
cycle_sf_wgs84 <- cycle_tbl %>%
  st_as_sf(coords = c("lon", "lat"), crs = 4326)

# 3) CRS 转换：WGS84 → British National Grid（BNG, EPSG:27700）
cycle_sf_bng <- st_transform(cycle_sf_wgs84, 27700)

# 4) 检查：对比两个 CRS 下 geometry 的坐标数值表达
head(cycle_sf_wgs84)
head(cycle_sf_bng)

# 5) 写出：将结果导出为 GeoPackage（推荐单文件格式）
out_dir  <- tempdir()
out_gpkg <- file.path(out_dir, "cycle_hire_bng.gpkg")

st_write(cycle_sf_bng, out_gpkg, layer = "cycle_hire_bng", delete_dsn = TRUE, quiet = TRUE)

# 验证：再读入确认
cycle_bng_out <- st_read(out_gpkg, layer = "cycle_hire_bng", quiet = TRUE)
head(cycle_bng_out)

格式说明：读入格式 ≠ 写出格式

示例采用“同格式读写”只是演示流程；在实际工作中，读入格式与写出格式并不需要一致。

sf::st_write() 支持将同一个 sf 对象写出为多种矢量格式（如 shp/geojson/gpkg 等），因此空间数据在不同存储格式之间可以相互转换。

同时，sf 对象也可以写出为纯表格格式（如 .csv），但需注意：

.csv 通常不会以标准空间结构保留几何信息；若仅保存属性列，则几何会丢失。
若需要“写出后还能恢复为空间对象”，通常应保留坐标列（点数据的 lon/lat），或将几何编码为 WKT/WKB 字符串再写出；后续可再用 st_as_sf() / st_as_sfc() 还原几何。

实践建议 gpkg 适合项目内复用与管理，geojson 适合跨平台交换；若需要进入专业 GIS 软件环境进行可视化、编辑与空间分析，也可将数据写出为 shp 以便在 ArcGIS / QGIS 等软件中直接使用。

2.3 按照地理键链接空间与非空间数据

重点 · 常用

在实践中，空间边界数据 (boundary data) 通常只包含地理几何信息（ geometry），其属性字段相对有限；而数据分析所需的统计指标往往以非空间表格形式单独存储（如 .csv 指标表）。

因此，在空间制图与空间数据分析中，一个常见的步骤是：依据一致的地理键（geographic key / geographic identifier），将外部表格中的指标字段并入空间边界数据，形成“边界 + 指标”的 sf 数据结构（如 Figure 4.5 所示）。

该操作等价于第二章的表格合并，通常使用 dplyr::left_join() 完成。

Figure 5. 通过地理键链接数据

地理键的常见形式

基础 · 概念

地理键（geographic key / identifier）用于为每一个空间单元提供唯一标识：在典型的边界数据中，每一行对应一个地理单元，因此地理键通常应当是唯一的（可视为该空间单元的“编号”）。

基于地理键的属性连接，本质上是在不同数据源之间建立“同一空间单元”的对应关系。

常见地理键形式主要包括两类：

地理编号/代码（code）：
以标准化编码唯一标识空间单元，通常更稳定、更适合作为连接键。
例如：英国统计地理体系中的 OA / LSOA / MSOA 代码、英国邮编（postcode）；美国的州/县/普查区代码体系；以及各类行政区划代码等。
地名字段（name）：
以名称描述空间单元，例如“上海市·浦东新区”“某某区/县”“某某州”等。
名称字段更直观，但在连接时更容易受到同名、别名、拼写差异、大小写、空格与符号差异等因素影响。

建议若同时提供 code 与 name，通常优先使用 code 作为连接键：其稳定性更高、歧义更少；name 更适合作为地图标注或图例展示字段。

常用连接方式：`left_join()` 将指标并入空间边界

常用

sf 对象的属性部分遵循数据框结构，因此可直接使用 dplyr::left_join() 进行属性连接。该操作的核心是：将外部指标表按地理键并入空间边界数据，同时保留 geometry 列不变。

基本写法

地理键同名（如，地理键所在的那一列都叫做 'geo_key'）：

boundary_sf %>% left_join(attr_tbl, by = "geo_key")
地理键不同名（如，地理键所在的那一列分别叫做 'geo_key_sf' 与 'geo_key_tbl'）：
boundary_sf %>% left_join(attr_tbl, by = c("geo_key_sf" = "geo_key_tbl"))

为何常用 `left_join()`？

空间制图中通常将空间边界数据作为左表（sf），将外部指标表作为右表（tibble/csv）。使用 left_join() 的主要原因在于：
- 边界完整保留：即使部分地理单元在指标表中缺失匹配，geometry 仍被保留，对应指标为 NA；
- 缺失可被表达：地图仍显示完整边界结构，缺失单元可通过 NA（或缺失色）明确呈现。

相对地，inner_join() 仅保留两表都能匹配到的空间单元：未匹配的边界将被剔除，地图上对应区域也会“消失”。该方式适用于明确需要排除缺失单元的场景，应结合研究目的选择。

示范案例

下面以 sf::nc 的边界数据为例演示“按地理键连接”的最基本流程。为使流程更贴近真实场景，这里将 nc 的属性部分抽取出来，构造一个“外部指标表”（等价于从 .csv 读入的指标表），再按地理键并回空间边界数据。

说明示例中使用 FIPSNO 作为地理键：该字段为县级单元的数值型标识，适合作为连接键。指标字段以“出生数 ÷ 面积”构造一个示意性密度指标（units 取决于原始字段定义；此处仅用于演示 join 流程）。

library(sf)
library(tidyverse)

# ------------------------------------------------------------
# 示例空间边界：sf 内置 nc（North Carolina counties）
# ------------------------------------------------------------
path_nc <- system.file("shape/nc.shp", package = "sf")
nc <- st_read(path_nc, quiet = TRUE)

# ------------------------------------------------------------
# 构造一个“外部指标表”（等价于从 csv 读入的指标表）
# 关键点：st_drop_geometry() 会移除 geometry 列
# - 移除后对象不再是 sf，而是普通表格（tibble/data.frame）
# - 此处仅用于模拟“外部属性表”的来源，便于演示 join 流程
# ------------------------------------------------------------
attr_tbl <- nc %>%
  st_drop_geometry() %>%          # 去掉 geometry：sf → 普通表格
  transmute(
    FIPSNO,                       # 地理键（county identifier）
    bir79_density = BIR79 / AREA  # 示例密度指标：1979 出生数 / 面积（示意）
  )

head(attr_tbl)

  FIPSNO bir79_density
1  37009     11964.912
2  37005      8885.246
3  37171     25286.713
4  37053     11857.143
5  37131     10496.732
6  37091     18948.454

# ------------------------------------------------------------
# 按地理键连接：将外部指标字段并回空间边界（sf 对象）
# 说明：以 nc（边界）为左表，确保所有边界保留；未匹配则指标为 NA
# ------------------------------------------------------------
nc_joined <- nc %>%
  left_join(attr_tbl, by = "FIPSNO")

# 检查：新指标字段已并入 (bir79_density)，geometry 保持不变
head(nc_joined)

Simple feature collection with 6 features and 15 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -81.74107 ymin: 36.07282 xmax: -75.77316 ymax: 36.58965
Geodetic CRS:  NAD27
   AREA PERIMETER CNTY_ CNTY_ID        NAME  FIPS FIPSNO CRESS_ID BIR74 SID74
1 0.114     1.442  1825    1825        Ashe 37009  37009        5  1091     1
2 0.061     1.231  1827    1827   Alleghany 37005  37005        3   487     0
3 0.143     1.630  1828    1828       Surry 37171  37171       86  3188     5
4 0.070     2.968  1831    1831   Currituck 37053  37053       27   508     1
5 0.153     2.206  1832    1832 Northampton 37131  37131       66  1421     9
6 0.097     1.670  1833    1833    Hertford 37091  37091       46  1452     7
  NWBIR74 BIR79 SID79 NWBIR79 bir79_density                       geometry
1      10  1364     0      19     11964.912 MULTIPOLYGON (((-81.47276 3...
2      10   542     3      12      8885.246 MULTIPOLYGON (((-81.23989 3...
3     208  3616     6     260     25286.713 MULTIPOLYGON (((-80.45634 3...
4     123   830     2     145     11857.143 MULTIPOLYGON (((-76.00897 3...
5    1066  1606     3    1197     10496.732 MULTIPOLYGON (((-77.21767 3...
6     954  1838     5    1237     18948.454 MULTIPOLYGON (((-76.74506 3...

3.0 空间数据处理

sf 对象在结构上可理解为“data.frame + geometry 列”，因此大多数 tidyverse 的数据操作可直接沿用，例如 select()、filter()、mutate()、arrange()、summarise() 等。

这一特性使得矢量空间数据能够在同一套语法体系下完成属性整理与制图所需的派生变量构造。

在属性层操作之外，sf 还提供面向 geometry 的空间处理能力，即常说的空间处理（geoprocessing）与空间分析操作：
- 叠置运算：相交、裁剪、并集/合并与差集等；
- 几何生成与变换：缓冲区、重心/代表点等；
- 空间匹配：将空间关系转化为属性并入的空间连接（spatial join）。

3.1 常规属性数据处理

参考第2章内容

在 sf 对象中，除 geometry 外的字段本质上仍是一张标准属性表，因此大多数属性层面的整理与构造可直接沿用第二章学习的 dplyr 的数据处理语法。

常用操作包括：

select()：筛选并保留制图/分析所需字段；
filter()：按条件筛选空间对象（例如限定空间范围或限定类别）；
mutate()：构造派生变量（如标准化、比率、密度、分组标签等）；常用
arrange()：按指标排序，用于质量检查与结果浏览；
rename()：统一字段命名，提高可读性与可复现性；
distinct()：检查或提取唯一键值（例如核对地理键是否唯一）；常用
summarise()：汇总统计（常与 group_by() 配合，用于形成上级尺度指标）。

在变量构造中，常见需求是基于一个或多个字段生成分类/标签变量，可通过 if_else() / case_when() 实现：
- if_else() 适用于两类条件分支；
- case_when() 适用于多条件分支（“分组规则表”写法）。常用

示例代码

下面使用 spData::world 作为示例数据。该数据包含国家边界与若干常用属性字段：
- continent：用于构造更粗粒度的洲别分组；
- pop 与 area_km2：用于派生人口密度指标。

本示例聚焦三个目标：
1) 用 case_when() 将 continent 重编码为 cont_group（将 North America 与 South America 归并为 Americas）；
2) 在 mutate() 中计算人口密度 pop_density = pop / area_km2；
3) 用 arrange() 对结果排序，以便快速检查高/低密度国家分布。

library(sf)
library(spData)
library(tidyverse)
options(sci = 9999)
data("world", package = "spData")

world2 <- world %>%
  # 1) 字段筛选：保留制图/检查常用字段（几何列也保留，用于后续空间制图）
  select(iso_a2, name_long, continent, area_km2, pop, lifeExp, gdpPercap, geom) %>%
  
  # 2) 派生变量：用 case_when() 进行“多条件分类/重编码”
  #    目标：把原始的 continent 进一步归并为更粗的分组标签 cont_group
  #    - North America + South America → Americas（美洲）
  #    - 其余保持为 Europe / Asia / Africa / Oceania
  #    - 未覆盖或缺失的情况统一标记为 Other
  mutate(
    cont_group = case_when(
      continent %in% c("Europe")                         ~ "Europe",
      continent %in% c("Asia")                           ~ "Asia",
      continent %in% c("Africa")                         ~ "Africa",
      continent %in% c("North America", "South America") ~ "Americas",
      continent %in% c("Oceania")                        ~ "Oceania",
      TRUE                                               ~ "Other"
    ),
    
    # 示例：构造人口密度（people per km^2）
    # 说明：area_km2 为面积（平方千米），pop 为人口数
    pop_density = pop / area_km2
  ) %>%
  
  # 3) 排序：便于检查
  arrange(desc(pop_density))


# 检查：查看分组结果与新变量
world2 %>%
  select(iso_a2, name_long, continent, cont_group, pop_density) %>%
  head()

Simple feature collection with 6 features and 5 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 3.314971 ymin: -2.917858 xmax: 129.4683 ymax: 53.5104
Geodetic CRS:  WGS 84
# A tibble: 6 × 6
  iso_a2 name_long    continent cont_group pop_density                      geom
  <chr>  <chr>        <chr>     <chr>            <dbl>        <MULTIPOLYGON [°]>
1 BD     Bangladesh   Asia      Asia             1192. (((92.67272 22.04124, 92…
2 PS     Palestine    Asia      Asia              853. (((35.39756 31.48909, 35…
3 LB     Lebanon      Asia      Asia              555. (((35.8211 33.27743, 36.…
4 KR     Republic of… Asia      Asia              512. (((126.1748 37.74969, 12…
5 RW     Rwanda       Africa    Africa            486. (((30.4191 -1.134659, 29…
6 NL     Netherlands  Europe    Europe            421. (((6.90514 53.48216, 6.0…

基于 `geometry` 的派生变量

重点

sf 的关键特征在于：几何列 geometry 可直接参与计算；因此 mutate() 不仅能生成普通属性变量，也能生成基于几何度量的空间指标。常见操作包括：

面要素：st_area() 计算面积；
边界/线要素：st_length() 计算长度（对面数据通常先取边界再计算周长）；
派生指标：将属性量与几何度量结合形成密度类指标（“属性 ÷ 面积/长度”）。

注意 st_area() / st_length() 的单位受 CRS 影响。若需要以“米/平方千米”等度量单位解释结果，通常应在投影坐标系下进行计算，并在计算前用 st_transform() 统一到合适的投影 CRS。

示例代码

下面使用 sf 自带的示例面数据 nc：读入后仅保留最基本的属性字段与 geometry，并在投影坐标系下计算面积与周长。这里使用 EPSG:3857 作为演示性投影坐标系（单位为米），以保证面积/长度的单位可直接换算。

library(sf)
library(tidyverse)

# ------------------------------------------------------------
# 示例面数据：sf::nc（North Carolina counties）
# 目标：基于 geometry 计算面积（km^2）与周长（km）
# ------------------------------------------------------------
path_nc <- system.file("shape/nc.shp", package = "sf")
nc <- st_read(path_nc, quiet = TRUE)

nc_geom <- nc %>%
  # 1) 字段整理：仅保留最基本字段 + geometry（便于聚焦“几何派生量”）
  select(NAME, geometry) %>%
  
  # 2) 投影转换：度量计算更适合在投影坐标系下进行
  #    这里用 EPSG:3857（单位：米）作为示例
  st_transform(3857) %>%
  
  # 3) 几何度量：面积 + 周长
  #    周长的计算：对面要素先取边界，再求边界长度
  mutate(
    area_m2       = st_area(geometry),                    # 面积（m^2，units 对象）
    area_km2      = as.numeric(area_m2) / 1e6,            # 换算为 km^2
    
    perimeter_m   = st_length(st_boundary(geometry)),     # 周长（m，units 对象）
    perimeter_km  = as.numeric(perimeter_m) / 1e3         # 换算为 km
  )

# 最小检查：去掉几何列便于打印
nc_geom %>%
  st_drop_geometry() %>%
  select(NAME, area_km2, perimeter_km) %>%
  head(10)

          NAME  area_km2 perimeter_km
1         Ashe 1760.2305     176.2149
2    Alleghany  946.4178     149.2654
3        Surry 2202.2231     199.5868
4    Currituck 1074.3327     375.1378
5  Northampton 2352.5281     263.5126
6     Hertford 1495.1612     199.9215
7       Camden  952.7121     187.0913
8        Gates 1398.7857     153.2864
9       Warren 1823.7884     175.4619
10      Stokes 1906.8947     174.8570

3.2 基础空间处理入门

重点

空间处理（geoprocessing）关注的是：如何基于 geometry 的空间关系与几何运算，生成新的空间对象或新的空间对应关系。

本节聚焦三类最常用操作：
- 空间关系判断：判断“是否包含/是否相交/是否接触”等空间关系；
- 空间连接：将一个图层的属性按空间关系写入另一个图层（如“点归属到面”, point in polygon）；
- 几何操作：通过缓冲、相交、裁剪、合并等操作改变几何形状或范围。

3.2.1 空间关系

空间关系（Spatial relationship）用于描述两个空间对象在位置上的“相对关系”（例如是否包含、是否相交、是否接触），是空间处理与空间连接的基础。

在 sf 中，这类判断通常通过一组函数完成（常被称为 spatial predicates）。重点掌握以下三类：

st_within()：判断 A 是否位于 B 的内部
- 典型理解：点是否落在某个行政区内。
st_intersects()：判断 A 与 B 是否有任何交集
- 交集可以是：共享面积、共享边界的一段，或仅共享一个点（只要“碰到”就算）。
st_touches()：判断 A 与 B 是否“只接触边界”
- 典型理解：两个面共享边界/顶点，但内部不重叠（没有共享面积）。

这类函数的输出本质上是在回答“哪些对象彼此满足某种空间关系”。它们通常被用来：
- 筛选：保留满足空间关系的要素；
- 匹配：为每个对象找到与之对应的对象；
- 连接：把匹配到的属性并入数据中（空间连接）。

Figure 6. 常见空间关系

3.2.2 空间连接

重点 · 常用

空间连接（spatial join）的核心是：
不依靠“相同列名/编号”进行匹配 (即，地理键，本章2.3节内容)，而是依靠对象之间的空间关系（是否落在内部、是否相交、是否接触）来完成“匹配”，并将匹配到的属性字段并入目标图层。

在 sf 中，空间连接最常用的形式是：

将 y 图层的属性，按指定空间关系，连接到 x 图层上：
st_join(x, y, join = st_within) 或 st_join(x, y, join = st_intersects) 也有相应的 st_within(), st_intersects()等

Figure 7. Spatial Join示意图

典型案例

点 → 面：
将点位归属到所在的面单元（经典 point-in-polygon 问题）。
常用于将 POI/事件点“分配”到行政区或规则网格，再进一步计算计数、密度或比率并制图。常用
线 → 面：
将道路、河流、轨迹等线要素与面单元匹配，识别其在面内的穿越/覆盖部分，并据此生成面尺度指标。
常见指标包括：线要素在面内的总长度、长度密度（length per area），或按类型分解后的网络供给强度。
面 ↔︎ 面：
处理两套面边界的叠置关系，用于边界体系转换、跨尺度对应与面尺度指标重分配。
常见于：将一套统计单元的指标转换到另一套规划/管理单元（例如“普查区 → 学区/服务区/规划分区”）；或计算两类区域之间的交叠面积比例，以支持后续的加权汇总(面积加权平均数计算)。较常用

示例：point in polygon问题

示例

下面构造一个矩形面作为“区域边界”，再随机生成一批点位。
目标是：用 st_join() 将点位归属到该面中，并进一步统计“面内点的数量”。
该流程与实践中“事件点/POI → 行政区/缓冲区”一致。

library(sf)
library(tidyverse)

# ============================================================
# 目标：演示“点 → 面”的空间连接（point in polygon）
# 核心函数：st_join(points, polygon, join = st_within)
# CRS：示例使用 EPSG:4326（WGS84，经纬度；单位为“度”）
# ============================================================

# ------------------------------------------------------------
# 1) 构造一个矩形面（POLYGON）
#    - st_polygon()：用“闭合坐标序列”定义一个面
#    - rbind()：按顺序给出边界点；最后一个点必须回到起点以闭合
# ------------------------------------------------------------
poly <- st_polygon(list(rbind(
  c(0, 0),
  c(10, 0),
  c(10, 10),
  c(0, 10),
  c(0, 0)   # 闭合：起点=终点，形成封闭边界
)))

# st_sfc()：将几何对象放入“几何列容器”（simple feature column）
# crs = 4326：声明该几何使用 WGS84 经纬度坐标系（EPSG:4326）
region <- st_sf(
  region_id = "A",
  geometry  = st_sfc(poly),
  crs       = 4326
)

# ------------------------------------------------------------
# 2) 随机生成点位（POINT）
#    - set.seed()：固定随机数种子，保证每次运行生成相同结果（可复现）
#    - runif(n, min, max)：从[min, max]区间生成 n 个均匀分布随机数
#    - 这里把范围设为 -2 到 12：故意让一部分点落在矩形外部
# ------------------------------------------------------------
set.seed(123)

pts_tbl <- tibble(
  id = 1:30,
  x  = runif(30, -2, 12),  # x 坐标：范围略大于矩形（0~10）
  y  = runif(30, -2, 12)   # y 坐标：范围略大于矩形（0~10）
)

# st_as_sf(coords = c("x","y"))：
# - 将表格中的 x/y 两列转换为 POINT 几何
# - crs = 4326：声明这些坐标也处于 WGS84（保证与 region 一致）
pts <- st_as_sf(pts_tbl, coords = c("x", "y"), crs = 4326)

# ------------------------------------------------------------
# 3) 空间连接：把面(region)的属性并入点(pts)
#    - st_join(x, y)：将 y 的属性连接到 x（输出行数通常以 x 为基准）
#    - join = st_within：判定“点是否位于面内部”
#      * 面内点：获得 region_id = "A"
#      * 面外点：region_id = NA（未匹配）
# ------------------------------------------------------------
pts_joined <- st_join(pts, region, join = st_within)

# 快速检查：统计匹配成功（A）与未匹配（NA）的点数量
pts_joined %>%
  st_drop_geometry() %>%     # 去掉几何列，仅看属性表（便于 count/打印）
  count(region_id)           # 查看 x 个点在 A内，x点不在

# A tibble: 2 × 2
  region_id     n
  <chr>     <int>
1 A            15
2 <NA>         15

# ------------------------------------------------------------
# 4) 点 → 面的汇总：统计每个面单元内部点的数量
#    - count(region_id)：按面ID计数（此处只有 "A" 与 NA 两类）
#    - right_join(region)：确保“面边界完整保留”
#      * 即使面内没有点，也仍保留该面，并把计数补为 0
# ------------------------------------------------------------
region_count <- pts_joined %>%
  st_drop_geometry() %>%
  count(region_id, name = "n_points") %>%
  right_join(st_drop_geometry(region), by = "region_id") %>%
  mutate(n_points = replace_na(n_points, 0))

region_count

# A tibble: 1 × 2
  region_id n_points
  <chr>        <int>
1 A               15

# ------------------------------------------------------------
# 5) 简单可视化 （后续小节会细讲）：面边界 + 点位
#    - 面使用浅色填充（仅用于示意范围）
#    - 点按 is.na(region_id) 着色：TRUE 表示在面外，FALSE 表示在面内
# ------------------------------------------------------------
ggplot() +
  geom_sf(data = region, fill = "darkgreen", alpha = 0.2, linewidth = 0.6) +
  geom_sf(
    data = pts_joined,
    aes(colour = is.na(region_id)),  # TRUE=面外（NA），FALSE=面内（匹配到 A）
    size = 2,
    alpha = 0.8
  ) +
  labs(
    title = "Spatial join (Point in Polygon) with st_join()",
    subtitle = "Points inside the polygon receive region_id; points outside remain NA",
    colour = "Outside polygon?"
  ) +
  theme_minimal()

3.2.3 简单空间处理

在空间处理（geoprocessing）中，常见任务可以有两类：
（1）对几何对象做构造/变形；（2）对多个几何对象做叠置运算（overlay）。

sf 通过一组以 st_ 开头的函数提供这些能力；入门阶段优先认识与了解以下类操作即可。

Figure 8. 常见空间处理示意图

裁剪（Clip）与相交（Intersection）

常用

参考 Figure 8 右上角 Intersect 与左下角 Clip。

设有两份矢量空间数据 x 与 y（均为 sf 对象，可为点/线/面图层）
在叠置运算中，st_intersection() 用于提取两者的共同空间部分，并据此生成新的几何结果。

Intersect（相交）：st_intersection(x, y) 返回 x ∩ y：

仅保留 x 与 y 在空间上同时覆盖的部分。
- 若 x/y 为面：输出为重叠区域的面片；
- 若涉及线/面：输出为位于面内的线段；
- 若涉及点/面：输出为落在面内的点。

Clip（裁剪）：可视为 Intersect 的“方向化”用法
将 y 视为裁剪边界（mask/boundary），仅保留 x 在 y 范围内的部分：
\[ \text{clip}(x \mid y) \equiv x \cap y \]
在 sf 中，裁剪通常还是用 st_intersection(x, y) 实现；差异主要体现在解释方式：
- Intersect 强调“两者共同部分”；
- Clip 强调“用 y 约束 x 的空间范围”。

典型范例
以行政区边界 y 作为裁剪范围：对道路、POI、轨迹或网格 x 进行裁剪，仅保留边界内部分，以保证后续制图与统计分析在同一空间范围内进行。研究区域的划定

并集（Union）与溶解（Dissolve）

参考 Figure 8 左上角 Union 与右下角 Split的反义词。

设有两份矢量空间数据 x 与 y（均为 sf 对象，可为点/线/面图层）。
Union（并集）的语义是合并两者的空间覆盖范围，保留 x ∪ y。在 sf 中主要通过 st_union() 实现。

Union（并集）：st_union(x, y)
返回 x ∪ y：输出几何覆盖范围包含 x 与 y 的全部区域（无论两者相交或分离）。
当边界存在交叠或穿越时，几何会发生叠置切割并重组，从而形成新的几何片段（overlay effect）。
溶解（dissolve）：st_union() 的常见应用方式
更常见的任务是在同一图层内部按分组将多个单元合并为更高层级的几何，并消除组内相邻单元的内部边界。
该操作通常称为 dissolve，可通过 group_by() + summarise() + st_union() 实现；同时需对组内属性字段进行一致的汇总规则（如求和、均值、计数）。

区分要点
- Union（并集）：面向两份数据/两个图层，强调空间范围的合并（x ∪ y）；
- Dissolve（溶解）：面向同一图层的分组整合，强调按组汇合几何并消除内部边界。
例如，将上海市各区边界按“上海市”合并为一个整体轮廓

差集（Difference）与对称差（Symmetric Difference）

设有两份矢量空间数据 x 与 y（均为 sf 对象，可为点/线/面图层）。这一组操作的核心在于：从 x 中“扣掉” y 的影响范围，或提取两者不重叠的部分。

参考 Figure 8 下中部 Erase 与左中部 Symmetrical Difference。

1）差集（Difference/Erase）

st_difference(x, y)：返回 x \ y
即：保留 x 中不被 y 覆盖/相交的部分；x 与 y 的重叠部分会被“抠除”。
若 x 与 y 相交：输出为 x 被切割后的“剩余部分” Figure 4.8 Erase部分；
若 x 与 y 不相交：输出通常与 x 相同（因为没有可扣除的部分）；
若 y 完全覆盖 x：输出可能为空几何（empty geometry）。

常见语义 Erase（擦除）：用 y 作为“遮罩”从 x 中擦除重叠区域，本质上就是 st_difference(x, y)。
若 y 由多个要素组成，实践中常先 st_union(y) 再做差集，以避免重复切割：
st_difference(x, st_union(y))。

2）对称差（Symmetric Difference）

st_sym_difference(x, y)：返回 (x \ y) ∪ (y \ x)
即：保留两者不重叠的部分，并将其合并；等价于“把交集挖空后，剩下的都要”。

区分
- Difference：只保留 x 的“剩余部分”（方向性：x \ y）；
- Symmetric Difference：保留 x 与 y 各自独有的部分（对称性：两边都保留）。

3.2.4 常用空间数据生成

缓冲区（Buffer）

重点 · 常用

st_buffer() 用于在点/线/面要素周围生成给定距离的缓冲区（buffer zone），常用于表达邻近范围与支撑基础的可达性/暴露分析。

点缓冲：
以每个点为圆心，通过设定半径生成圆形（或近似圆形）影响范围
（例如地铁站点的辐射区）；
线缓冲：
以线要素为“中心轴”，对线的两侧按给定距离进行外扩；从几何角度看，可理解为对线上的各位置做等距扩张，并在端点处形成圆帽或平帽（端点样式取决于实现参数/默认设置），从而得到走廊状区域；
面缓冲：
以面要素的边界为基准，对边界进行等距平移：
- 距离为正时，边界向外扩张（外扩）；
- 距离为负时，边界向内收缩（内缩）；
  结果是得到一个新的面几何，用于表示“边界扩大/缩小”的区域范围。

注意缓冲距离的单位由 CRS 决定。涉及距离的缓冲操作通常应在投影坐标系下进行（如，单位为米/千米或英尺等），以保证距离含义可解释。

Figure 9. 缓冲区分析示意图

重心（Centroid）

常用

在空间分析与制图中，常需要将面要素 (或线要素) 转换为一个代表性点位，用于标注、连线(网络分析里很常见)、简单直线距离计算或与点数据统一表达尺度。

Centroid 在中文语境中常译为重心或质心；在本章作为入门，可将其理解为：由面/线几何计算得到的“代表点”。

在 sf 中，可使用 st_centroid() 由线/面/多面要素生成质心点。需要注意两点：

st_centroid() 的输出仍是 sf 对象，并默认保留原有属性字段（即“属性表不变，几何从面变为点”）；
对于形状凹陷或带孔的多边形，几何质心可能落在面外（入门阶段了解即可）。若需要“保证点落在面内”，可使用 st_point_on_surface()。

Figure 10. 面生成重心示意图

下面以 sf 内置的 nc 面边界为例：生成各县的 centroid，并将“边界面 + centroid 点”叠加可视化。

library(sf)
library(ggplot2)

# ------------------------------------------------------------
# 1) 读取示例面数据：North Carolina counties（MULTIPOLYGON）
# ------------------------------------------------------------
path_nc <- system.file("shape/nc.shp", package = "sf")
nc <- st_read(path_nc, quiet = TRUE)

# ------------------------------------------------------------
# 2) 生成 centroid：面 → 点
#    - 输出仍为 sf
#    - 默认保留原始属性字段（如 NAME/FIPS 等）
# ------------------------------------------------------------
nc_cent <- st_centroid(nc)

# 检查：几何类型已变为 POINT，属性列仍在
st_geometry_type(nc)[1] # 注意观察，nc本来的空间类型是 面

[1] MULTIPOLYGON
18 Levels: GEOMETRY POINT LINESTRING POLYGON MULTIPOINT ... TRIANGLE

st_geometry_type(nc_cent)[1] # 注意观察，通过函数生成的空间数据类型是 点

[1] POINT
18 Levels: GEOMETRY POINT LINESTRING POLYGON MULTIPOINT ... TRIANGLE

head(nc_cent) # 注意观察，生成的重心也具有多边形的属性！

Simple feature collection with 6 features and 14 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: -81.49823 ymin: 36.36142 xmax: -76.02719 ymax: 36.49111
Geodetic CRS:  NAD27
   AREA PERIMETER CNTY_ CNTY_ID        NAME  FIPS FIPSNO CRESS_ID BIR74 SID74
1 0.114     1.442  1825    1825        Ashe 37009  37009        5  1091     1
2 0.061     1.231  1827    1827   Alleghany 37005  37005        3   487     0
3 0.143     1.630  1828    1828       Surry 37171  37171       86  3188     5
4 0.070     2.968  1831    1831   Currituck 37053  37053       27   508     1
5 0.153     2.206  1832    1832 Northampton 37131  37131       66  1421     9
6 0.097     1.670  1833    1833    Hertford 37091  37091       46  1452     7
  NWBIR74 BIR79 SID79 NWBIR79                   geometry
1      10  1364     0      19  POINT (-81.49823 36.4314)
2      10   542     3      12 POINT (-81.12513 36.49111)
3     208  3616     6     260 POINT (-80.68573 36.41252)
4     123   830     2     145 POINT (-76.02719 36.40714)
5    1066  1606     3    1197 POINT (-77.41046 36.42236)
6     954  1838     5    1237 POINT (-76.99472 36.36142)

# ------------------------------------------------------------
# 3) 叠加可视化：县级边界 + centroid 点
# ------------------------------------------------------------
ggplot() +
  geom_sf(data = nc, fill = "grey95", linewidth = 0.2) +
  geom_sf(data = nc_cent, size = 0.8, alpha = 0.8, color = 'blue') +
  labs(
    title = "Polygon boundaries and their centroids",
    subtitle = "Example: sf::nc (county polygons) + st_centroid() points",
    caption = "Data: sf built-in example (North Carolina counties)"
  ) +
  theme_minimal()

4.0 空间数据可视化

4.0 空间可视化指南

制图导航

空间数据可视化图谱

“无地图，不空间”

本节将正式梳理矢量空间数据的完整制图工作流：以 sf 为底层数据框架（负责几何特征与属性管理），重点解析研究中最为高频的分层设色图（Choropleth Map）及其静态渲染核心 ggplot2::geom_sf()；在此基础上，进一步剖析数据分级与色彩映射背后的算法逻辑（如自然断点法等）。随后，我们将引入网络底图（Basemap）丰富制图背景，并最终将图面表达从静态渲染延伸至动态交互式地图的构建。

点击下方卡片，快速跳转至相应模块：

可视化基础

sf 对象与要素

分层设色图

空间映射渲染

数据分级

算法与色彩断点

底图绘制

网络瓦片底图

交互地图

缩放与动态交互

💡 提示：学习完一个小节后，请再次点击 屏幕右下角的制图导航按钮回到本导航页

制图基础

4.1 静态制图基础：`sf` × `ggplot2::geom_sf()`

Figure 11 展示了 sf 与 ggplot2 配合进行矢量空间数据可视化的流程示意图。

sf 负责矢量几何（点/线/面）及其坐标参考系（CRS）的记录与管理；ggplot2 负责图层式绘图语法与图形呈现。在两者的衔接上，ggplot2::geom_sf() 提供面向 sf 对象的地图图层接口，coord_sf() 用于控制地图的显示范围与投影/坐标表达。

因此，静态地图的基本写法可沿用第三章已掌握的 ggplot2 结构：以 ggplot() 建立画布，通过 + 逐层叠加 geom_sf()、coord_sf()、labs() 与 theme_*() 等组件。整体流程可概括为：数据准备 → 图层叠加 → 坐标控制 → 主题布局 → 导出输出。

Figure 11. sf + ggplot2 空间数据出图示意图

基础制图框架

重点

数据准备：
读取为 sf 对象；检查几何类型与 CRS；
注意: 在多图层叠加前，必要时使用 st_transform() 统一 CRS。
图层叠加：
以 ggplot() 建立画布，使用多个 geom_sf() 分层绘制；常见顺序为“底图 → 主题图层 → 标注/辅助图层”，以保证信息层级清晰。
坐标控制：
用 coord_sf() 控制地图显示范围与坐标表达：
- xlim/ylim 设置绘图范围；
- crs 指定投影/坐标参考系（用于地图表达与图层对齐）。
视觉与布局：
用 labs() 与 theme_*() 设置标题、副标题、图例与版式；必要时补充制图要素：
- 图例：
  说明颜色/符号含义（常规由 ggplot2 自动生成，但可进一步微调）；
- 比例尺与指北针：
  用于距离与方向的基本参照（可通过扩展组件添加）。
导出输出：
使用 ggsave() 输出可复现的图形文件（PNG/PDF 等）；也可在 RStudio 的 Plots 面板中通过 Export 按钮进行手动导出（更适合快速预览与临时输出）。推荐

代码示例

下述代码使用 spData 包内置的示例数据，演示在 sf × ggplot2::geom_sf() 框架下绘制矢量图层，并展示基础的图层叠加方式。

lnd：大伦敦地区行政区边界（面要素 / polygon）；
cycle_hire_osm：伦敦共享单车租赁站点（点要素 / point，OSM 数据整理）；

代码 · 示例1

单层底图（面）

注：ggspatial 包提供地图制图的要素

library(sf)
library(spData)
library(ggplot2)
library(ggspatial)  # 为 ggplot2 地图提供“制图要素”注释层：比例尺、指北针、经纬网等

data("lnd", package = "spData")  # 示例面数据：大伦敦行政区边界（sf polygon）

ggplot(lnd) +
# ------------------------------------------------------------
# 1) 底图图层：绘制行政区边界线
#    linewidth 控制边界线粗细
# ------------------------------------------------------------
geom_sf(linewidth = 0.2) +
  
# ------------------------------------------------------------
# 2) 比例尺（scale bar）
#    annotation_scale() 会依据当前地图的 CRS/坐标单位推算距离
#    - location   : 放置位置（"bl"=左下角，"br"=右下角，"tl"/"tr" 类似）
#    - width_hint : 比例尺宽度相对绘图区宽度的比例（0~1），0.25 表示约占 1/4
#    - line_width : 比例尺线条粗细
#    - height     : 比例尺条的高度（使用 grid::unit() 指定物理单位）
# ------------------------------------------------------------
ggspatial::annotation_scale(
  location   = "bl",
  width_hint = 0.25,
  line_width = 0.5,
  height     = unit(0.15, "cm")
) +
  
# ------------------------------------------------------------
# 3) 指北针（north arrow）
#    annotation_north_arrow() 用于添加方位指示
#    - location    : 放置位置（这里 "br"=右下角）
#    - which_north : "true" 表示真北（一般地图默认即可）
#    - style       : 指北针样式（ggspatial 提供若干预设样式函数）
# ------------------------------------------------------------
ggspatial::annotation_north_arrow(
  location    = "br",
  which_north = "true",
  style       = ggspatial::north_arrow_fancy_orienteering
) +
  
# ------------------------------------------------------------
# 4) 最小文本要素：标题 + 数据来源（caption 常用于数据源/时间口径）
# ------------------------------------------------------------
labs(
  title   = "London borough boundaries (spData::lnd)",
  caption = "Data: spData"
) +
  
# ------------------------------------------------------------
# 5) 主题：最简洁的版式
# ------------------------------------------------------------
theme_minimal()

代码 · 示例2

图层叠加（面 + 点）—— 底图 + 点位

library(sf)
library(spData)
library(ggplot2)
library(ggspatial)

data("lnd", package = "spData")
data("cycle_hire_osm", package = "spData")

# 1) CRS 对齐
cycle_hire_osm2 <- st_transform(cycle_hire_osm, st_crs(lnd))

# 2) bbox: bounding box
bb <- st_bbox(lnd)

ggplot() +
# ----------------------------------------------------------
# 底图（面）
# ----------------------------------------------------------
geom_sf(data = lnd, fill = "grey95", color = "grey60", linewidth = 0.2) +
  
# ----------------------------------------------------------
# 叠加点位（点）：
# 1. 在 aes() 内部映射 color = "字符串名称"
# 2. size 和 alpha 依然可以在外面控制
# ----------------------------------------------------------
geom_sf(data = cycle_hire_osm2, 
        aes(color = "Cycle Hire Stations"), 
        size = 0.5, 
        alpha = 0.7) +
  
# ----------------------------------------------------------
# 手动设置颜色图例
# values: 指定实际颜色 (这里设为红色)
# name: 图例标题 (设为 NULL 去掉标题，或写 "Legend")
# ----------------------------------------------------------
scale_color_manual(name = "Legend", 
                   values = c("Cycle Hire Stations" = "#E41A1C")) +
  
# ----------------------------------------------------------
# 调整图例样式 (可选)：让图例里的点更大、不透明，方便看清
# ----------------------------------------------------------
guides(color = guide_legend(override.aes = list(size = 2, alpha = .7))) +
  
# ----------------------------------------------------------
# 绘图窗口范围控制
# ----------------------------------------------------------
coord_sf(
  xlim   = c(bb["xmin"], bb["xmax"]),
  ylim   = c(bb["ymin"], bb["ymax"]),
  expand = FALSE
) +
  
# ----------------------------------------------------------
# 比例尺
# ----------------------------------------------------------
ggspatial::annotation_scale(
  location   = "bl",
  width_hint = 0.25,
  style      = "ticks", # 换个样式
  line_width = 0.5,
  pad_x       = unit(0.45, "cm"),
  pad_y       = unit(0.5, "cm")
) +
  
# ----------------------------------------------------------
# 指北针
# ----------------------------------------------------------
ggspatial::annotation_north_arrow(
  location    = "bl", # 改到右上角，避免遮挡
  which_north = "true",
  style       = ggspatial::north_arrow_fancy_orienteering,
  height      = unit(1, "cm"),
  width       = unit(1, "cm"),
  pad_x       = unit(0.55, "cm"),
  pad_y       = unit(1.00, "cm")
) +
  
  labs(
    title    = "Cycle Hire Locations in London",
    subtitle = "Spatial distribution of hire stations within inner boroughs",
    caption  = "Source: spData package"
  ) +
  theme_minimal() +
  # 将图例放在合适的位置（例如底部或右侧）
  theme(legend.position = "bottom")

综合：基于 OpenStreetMap 的多图层空间可视化

在实际城市空间研究中，完整的专题地图往往需要融合多种类型的地理要素（如点、线、面）。本节将展示一个综合性的微观城市空间制图案例。

在获取真实世界的精细矢量数据时，除了读取本地存储的空间文件，还可以利用 osmdata 包。该工具能够通过从 OpenStreetMap (OSM) 开源数据库中按需提取指定范围的空间要素（如城市水系、绿地、道路网与建筑轮廓等）。

知识拓展：OpenStreetMap (OSM) 与 osmdata 包

OpenStreetMap (OSM) 是一个开源、协作式的全球地理空间数据库。它包含了由全球研究者与志愿者共同维护的极高精度的微观地理要素（如建筑轮廓、路网层级、兴趣点 POI、土地利用类型等），是计算社会科学、城市规划等领域不可或缺的开源数据库。

在 R 语言生态中，osmdata 包是连接底层数据库的核心桥梁。它封装了底层的 Overpass API，允许研究者通过简单的逻辑语句（如指定 Bounding Box 和特征标签）直接抓取所需的矢量要素，并将其无缝转化为标准的 sf 空间对象，直接用于 ggplot2 渲染与空间统计分析。

📚 推荐参考资料与进阶指南：

OSM Map Features ：https://wiki.openstreetmap.org/wiki/Map_features
osmdata 官方教程：https://docs.ropensci.org/osmdata/

OSM 的标签体系（Tags）极为庞大且精细。在实际研究中，若不确定某种具体的空间要素（例如“三甲医院”、“地铁站点”、“特定等级的公路”）在 OSM 中对应的准确键值对，建议直接查阅上述官方维基，或者直接向 AI 进行提问（例如：“我想用 R 语言的 osmdata 包获取某个城市的咖啡馆分布，应该使用哪个 key 和 value？”）

在完成要素提取后，依托 ggplot2 的图层叠加逻辑即可构建综合地图。

本案例在制图过程中结合了以下几个关键点：
1. 多图层叠加与色彩映射：
2. 标准制图要素引入：使用 ggspatial 包读取底图的坐标参考系，添加比例尺与指北针；
3. 网格置顶渲染：通过调整底层主题参数（如设定 panel.ontop = TRUE），强制将虚线经纬网格置于所有图层之上，以确保坐标参考的完整性与图面工程质感。

⚠️ 运行提示：
由于 osmdata 实时抓取大量矢量数据高度依赖外部 API 与网络连接的稳定性，为保证教程的流畅阅读，下述综合案例的代码不在本文档中进行实时编译演示。读者可自行复制并在本地 R 环境中运行，其生成的最终地图效果如下方图片所示。

代码 · 示例3

# 如果未安装，请先运行 install.packages("osmdata")
library(tidyverse)
library(sf)
library(osmdata) # 用于通过 API 实时获取 OpenStreetMap 矢量数据
library(ggspatial)

# ============================================================
# 1. 定义查询范围 
# ============================================================
# 获取曼哈顿下城的经纬度边界框
bbox_manhattan <- getbb("Lower Manhattan, New York")

# ============================================================
# 2. 构建查询并分别获取各个图层的数据
# ============================================================
# [提示] 这里会向 API 发送 3 次请求，可能需要 10-30 秒，请耐心等待

# 水系 (Water)
water_sf <- opq(bbox = bbox_manhattan) %>%
  add_osm_feature(key = "natural", value = "water") %>%
  osmdata_sf()

# 公园绿地 (Parks)
parks_sf <- opq(bbox = bbox_manhattan) %>%
  add_osm_feature(key = "leisure", value = "park") %>%
  osmdata_sf()


# 建筑轮廓 (Buildings)
bldg_data <- opq(bbox = bbox_manhattan) %>%
  add_osm_feature(key = "building") %>%
  osmdata_sf()

# ============================================================
# 3. 提取对应的矢量几何对象 
# ============================================================

water_poly <- water_sf$osm_polygons
water_multi <- water_sf$osm_multipolygons
parks_poly <- parks_sf$osm_polygons
buildings_poly <- bldg_data$osm_polygons

# ============================================================
# 4. 图层叠加与城市综合地图绘制
# ============================================================
ggplot() +
  
  # --- 图层 0：手动绘制底层陆地背景 ---
  # [技巧] 为了把经纬网放到最顶层，必须让默认背景透明。因此我们用 annotate 手动铺一层底色
  annotate("rect", xmin = -Inf, xmax = Inf, ymin = -Inf, ymax = Inf, fill = "#f2efe9") +
  
  # --- 图层 1：底层水系 ---
  # 
  geom_sf(data = water_poly, aes(fill = "Water"), color = NA) +
  geom_sf(data = water_multi, aes(fill = "Water"), color = NA) +
  
  # --- 图层 2：公园绿地 ---
  geom_sf(data = parks_poly, aes(fill = "Parks"), color = NA) +
  
  # --- 图层 3：表层建筑物 ---
  geom_sf(data = buildings_poly, aes(fill = "Buildings"), color = NA) +
  
  # --- 图层 4：指北针与比例尺 (依赖 ggspatial) ---
  # 指北针放在左上角 (tl: top-left)
  annotation_north_arrow(location = "tl", which_north = "true", 
                         pad_x = unit(0.1, "in"), pad_y = unit(0.1, "in"),
                         style = north_arrow_fancy_orienteering()) +
  # 比例尺放在左下角 (bl: bottom-right)
  annotation_scale(location = "br", width_hint = 0.2, style = "ticks",
                   text_col = "black", line_col = "black") +
  
  # --- 颜色映射与图例控制 ---
  scale_fill_manual(
    name = "Urban Features",      # 图例的大标题
    values = c(
      "Water" = "#caddf5",
      "Parks" = "#c8dfb3",
      "Buildings" = "#cfc8c0"
    ),
    breaks = c("Buildings", "Parks", "Water") # 强制图例显示顺序
  ) +
  
  # --- 视角控制：裁剪边界 ---
  coord_sf(xlim = c(bbox_manhattan[1,1], bbox_manhattan[1,2]),
           ylim = c(bbox_manhattan[2,1], bbox_manhattan[2,2]),
           expand = FALSE) +
  
  # --- 主题美化：网格置顶与细节排版 ---
  theme_bw() + 
  theme(
    # 1. 强制将经纬网格置于所有图层之上，并将面板背景彻底设为透明
    panel.ontop = TRUE, 
    panel.background = element_rect(fill = NA, color = NA),
    
    # 2. 经纬网格：加粗的黑色点线 (dotted)
    panel.grid.major = element_line(color = "black", linetype = "dotted", linewidth = 0.5),
    panel.grid.minor = element_blank(),
    
    # 3. 边框与坐标文字
    panel.border = element_rect(color = "grey20", fill = NA, linewidth = 1.8),
    axis.text = element_text(color = "grey40", size = 8),
    axis.ticks = element_line(color = "grey40"),
    
    # 4. 图例框线与位置微调
    legend.position = "right",
    legend.title = element_text(face = "bold", size = 10),
    legend.background = element_rect(fill = "white", color = "grey80", linewidth = 0.3),
    legend.margin = margin(5, 5, 5, 5),
    
    # 5. 标题排版
    plot.title = element_text(face = "bold", size = 16, hjust = 0.5, margin = margin(t = 15, b = 5)),
    plot.subtitle = element_text(color = "grey50", hjust = 0.5, margin = margin(b = 15)),
    plot.caption = element_text(color = "grey60", hjust = 1)
  ) +
  labs(
    title = "Multi-layer Urban Map: Lower Manhattan",
    subtitle = "Architectural Layout with Top-level Graticules & Map Elements",
    caption = "Data sourced from OpenStreetMap via {osmdata}",
    x = NULL, y = NULL 
  )

上述代码成果图

分层设色

4.2 分层设色图（Choropleth map）

重要 · 常用

分层设色图（choropleth map）通常以互不重叠的面单元（统计区/行政区等 enumeration units）作为制图单元，将“每个面单元对应的一个数值指标”编码为填充色（fill）的色阶强度，从而呈现指标在空间上的差异分布。其核心可概括为：区域单元 + 区域指标 + 颜色编码。

术语提示点要素或线要素同样可以用颜色深浅表达数值差异，但这类图通常归入“分级符号图/分级线（graduated symbols/lines）”等专题制图类型；为避免概念混用，本章将“choropleth map”专指面单元的分层设色表达。

1. 核心用途

空间差异表达：直观呈现区域间的高低差异与梯度变化；
空间格局识别：辅助识别集聚、分异与分区结构；

2. 基本绘制逻辑：geom_sf(aes(fill = ...))

在 ggplot2 中，分层设色图的核心结构通常由三部分构成：

面边界图层：geom_sf() 负责绘制边界；
设色映射：在 aes() 中将指标映射到 fill（连续/离散都可）；
色标与图例：用 scale_fill_*() 控制配色与图例显示，并通过 labs(fill = ...) 明确指标含义。

提示
分层设色图更适合映射率/密度/比率（rate / density / ratio）等可比指标；直接映射总量（count）往往会被面积或人口规模主导，解释需谨慎。

Figure 12. 常见的分层设色图出图流程示意图

代码示例 1：`spData::world` 绘制人口密度分层设色图

示例代码

下述示例以 spData::world 的国家边界数据为底图，利用 pop（人口）与 area_km2（面积，km²）构造人口密度指标 pop_density = pop / area_km2（单位：人/平方千米），并将该指标映射到 fill 绘制分层设色图。

library(sf)
library(spData)
library(tidyverse)
library(viridis)
library(units)
library(ggnewscale) # 新增：允许在同一图形中使用多个相同的标尺（如两个 fill 标尺）

options(scipen = 9999) # 禁用科学计数法显示

# ------------------------------------------------------------
# 1. 数据获取与预处理
# ------------------------------------------------------------
data("world", package = "spData")

# 指标计算：人口密度 (人/平方公里)
world_all <- world %>%
  mutate(pop_density = pop / area_km2)

# 数据分层 
# 将缺失值 (NA) 与有效观测值分离，以便为两类数据应用独立的视觉映射和图例
world_valid <- world_all %>% filter(!is.na(pop_density))
world_na    <- world_all %>% filter(is.na(pop_density))

# ------------------------------------------------------------
# 2. 可视化构建
# ------------------------------------------------------------
ggplot() +
  
  # --- 图层 1：缺失数据层 ---
  # 使用离散标尺绘制 NA 区域
  geom_sf(
    data = world_na, 
    aes(fill = "No Data"),  # 映射至固定字符串以生成离散图例项
    linewidth = 0.22,       # 地块边界线粗细
    colour = "white"        # 地块边界线颜色
  ) +
  
  # 标尺定义 1：缺失数据的离散颜色映射
  scale_fill_manual(
    name = NULL,            # 不单独命名图例 （因为不是重点）
    values = c("No Data" = "grey55"), 
    guide = guide_legend(order = 2) # 指定图例顺序：置于下方
  ) +
  
  # --- 标尺重置 ---
  # 关键操作：初始化新的 fill 标尺，解除后续图层对上一标尺的依赖
  new_scale_fill() +
  
  # --- 图层 2：有效数据层 ---
  # 使用连续标尺绘制人口密度
  geom_sf(
    data = world_valid, 
    aes(fill = pop_density), 
    linewidth = 0.22, 
    colour = "white"
  ) +
  
  # 标尺定义 2：有效数据的连续颜色映射
  # 采用对数变换 (Log transformation) 以增强偏态分布数据的可视化对比度
  scale_fill_viridis_c(
    option    = "magma",
    direction = -1,         # 颜色反转：深色表示高密度
    trans     = "log10",    # 对数变换
    labels    = scales::label_number(accuracy = 0.1), # 图例数值的位数，保留1位小数
    name      = "Population density\n(people / km²)", # 标注标题， 单位 （\n代表换行）
    guide     = guide_colourbar(order = 1) # 指定图例顺序：置于上方
  ) +
  
  # --- 图表布局与修饰 ---
  labs(
    title    = "Choropleth map: Population density", # 参考第三章设置中文标题，注释等
    subtitle = "Fill encodes pop_density = pop / area_km² (log10 scale); darker = higher",
    caption  = "Data source: spData::world"
  ) +
  theme_minimal() + # 可以自行更换主题
  theme(
    panel.grid.major  = element_line(linewidth = 0.25),
    legend.key.height = unit(0.7, "cm"),
    legend.spacing.y  = unit(0.2, "cm") # 调整图例项之间的垂直间距
  )

代码示例 2：基于 `tidycensus` 的纽约市贫困率空间可视化

拓展 · 专业 · 常用

本示例演示如何利用 tidycensus 包获取美国人口普查局（U.S. Census Bureau）发布的 ACS（American Community Survey） 5-year 数据，并以纽约市(New York City, NYC)为例，绘制 Census Tract 尺度下的人口贫困率分层设色图。

ACS 是美国最权威的官方社会经济统计体系之一，数据口径统一且公开，广泛应用于社区尺度的空间差异分析与制图。

Census Tract 是美国人口普查局定义的统计分区，人口通常维持在 4,000 人左右，其空间范围随人口密度变化（高密度区面积更小），是社会科学中表征邻里/社区（Neighborhood）最常用的精细地理尺度。

核心处理流程：

数据获取：
调用 tidycensus::get_acs() 拉取纽约市2023年份的贫困人口与总人口数据；
空间关联：
设置 geometry = TRUE 直接返回包含地理边界的 sf 对象（属性表 + 几何字段）；
指标计算：
构造贫困率指标 poverty_rate（贫困人口 / 总人口）；
离散分级：
利用 classInt 包的 Fisher-Jenks 算法计算自然断点，将连续的贫困率划分为若干离散等级（最优阈值划分）；
制图输出：
使用 ggplot2::geom_sf() 将等级映射至填充色 fill，完成可复现的专业地图绘制。

说明 tidycensus 接口支持从 Block Group (美国官方最小地理普查单元)到 State 等多种地理尺度（通过 geography 参数控制），极易实现跨尺度对比。本节重点在于全流程可视化实现，关于 API 配置及变量编码细节，建议查阅 ?tidycensus 或访问官方文档。

library(sf)          
library(tidyverse)   
library(tidycensus)  # 美国人口普查数据接口: 用于直接拉取 ACS/Decennial Census 数据及配套边界
library(tigris)      # 美国地理边界数据包: 提供 TIGER/Line 矢量边界文件
library(viridis)     # 色彩映射包
library(classInt)    # 单变量分级统计包: 用于计算 Fisher-Jenks 等断点
library(ggspatial)   # 地图制图扩展包: 添加比例尺、指北针等地图要素
library(units)       # 计量单位处理包: 用于控制比例尺单位换算

# 设置选项：启用 tigris 缓存以加速边界数据加载；禁用科学计数法
options(tigris_use_cache = TRUE)
options(scipen = 9999)

# ==============================================================================
# 1. 实验设计与参数设定 
# ==============================================================================

# 定义数据年份与地理范围
# 采用 5-year ACS 数据以确保小尺度 (Census Tract) 估计值的稳定性
acs_year <- 2023
# 定义纽约市 (NYC) 五大区对应的 County 名称
nyc_counties <- c("New York", "Kings", "Queens", "Bronx", "Richmond")

# ==============================================================================
# 2. 数据获取与指标构建 
# ==============================================================================

# 通过 tidycensus API 获取 ACS 贫困数据
# 地理尺度: Census Tract (普查小区)，适用于细粒度的城市内部空间差异分析
# 变量定义:
#   - B17001_001: 确定贫困状态的总人口 
#   - B17001_002: 收入低于过去12个月贫困线的人口 
nyc_tracts <- get_acs(
  geography = "tract",
  state     = "NY",
  county    = nyc_counties,
  variables = c(tot = "B17001_001", pov = "B17001_002"),
  year      = acs_year,
  survey    = "acs5",
  geometry  = TRUE,   # 同时获取空间几何信息 (sf对象)
  output    = "wide"  # 宽格式输出，便于直接进行列运算
)

# 指标计算：贫困率 (Poverty Rate)
# 数据清洗：对于总人口极少 (如机场、公园、工业区等 totE <= 100) 的区域，
# 其比率极不稳定，将其设为 NA 以避免误导性高值。
# NA_real_ 意思是：数值型的缺失值； NA本身是逻辑型

nyc_tracts2 <- nyc_tracts %>%
  mutate(poverty_rate = if_else(totE > 100, povE / totE, NA_real_))

# ==============================================================================
# 3. 空间数据预处理 
# ==============================================================================

# 获取 NYC 轮廓边界
# 美化用途：用于在地图最上层绘制深色轮廓，增强地图的结构感和区域辨识度
nyc_outline <- counties(state = "NY", cb = TRUE, year = acs_year, class = "sf") %>%
  filter(NAME %in% nyc_counties) %>%
  summarise(geometry = st_union(geometry), .groups = "drop") # 融合五个区的边界为一个整体

# 投影转换 
# 目标坐标系: EPSG:2263 (NAD83 / New York Long Island (ftUS))
# 说明: 这是一个局部投影坐标系，相比于经纬度 (EPSG:4326)，它能提供更准确的长度和面积量测，
# 且适用于绘制带线性单位 (ft/m) 的比例尺。
crs_nyc <- 2263
nyc_tracts2 <- st_transform(nyc_tracts2, crs_nyc)
nyc_outline <- st_transform(nyc_outline, crs_nyc)

# ==============================================================================
# 4. 统计分级：Fisher-Jenks 自然断点法 
# ==============================================================================

# 设定分级数量 建议不要超过7
n_class <- 6

# 提取有效数值用于计算断点 (剔除 NA 和 Inf)
# 剔除 Inf（无穷大）、-Inf（负无穷大）、NA（缺失值）和 NaN（非数字）
vals <- nyc_tracts2$poverty_rate
vals <- vals[is.finite(vals)]

# 计算 Fisher-Jenks 断点
# 算法逻辑: 迭代寻找使组内方差最小化、组间方差最大化的阈值
fi   <- classInt::classIntervals(vals, n = n_class, style = "fisher", warnSmall = FALSE)
brks <- fi$brks

# 构建图例标签 (百分比格式)
# 格式示例: "0% – 5%", "5% – 12%"
brk_labs <- paste0(
  scales::label_percent(accuracy = 1)(brks[-length(brks)]), " – ",
  scales::label_percent(accuracy = 1)(brks[-1])
)

# 数据离散化
nyc_tracts2 <- nyc_tracts2 %>%
  mutate(
    # 将连续变量 poverty_rate 切分为离散区间
    pov_class = cut(
      poverty_rate,
      breaks = brks,
      include.lowest = TRUE,
      right  = FALSE,
      labels = brk_labs
    ),
    # 处理缺失值: 将 NA 显式转换为 "No data" 因子水平
    pov_class = forcats::fct_na_value_to_level(pov_class, level = "No data")
  )

# ==============================================================================
# 5. 可视化参数设定 
# ==============================================================================

# 定义色彩映射 (Color Palette)
# 色板: Viridis "mako" (适合表达社会经济数据的冷色调方案)
# 方向: direction = -1 (深色表示高贫困率，浅色表示低贫困率)
pal_cols <- viridis::viridis(n_class, option = "mako", direction = -1)

# 将颜色与标签绑定，并追加 "No data" 的颜色定义
names(pal_cols) <- brk_labs
# 将颜色向量转换为命名向量，确保 scale_fill_manual 能正确匹配
pal_cols_final <- c(pal_cols, "No data" = "grey55")

# 定义图例排序
# 逻辑: 将高值区间放在图例顶部，低值在下，No data 置于最底
legend_breaks <- c(rev(brk_labs), "No data")

# ==============================================================================
# 6. 图层构建与制图 
# ==============================================================================

ggplot() +
  # --- 图层 1: 贫困率填色层  ---
  geom_sf(
    data = nyc_tracts2,
    aes(fill = pov_class),
    linewidth = 0, # 去除内部 Tract 边界线，避免线条过密干扰视觉
    colour = NA
  ) +
  
  # --- 图层 2: 区域轮廓层  ---
  geom_sf(
    data = nyc_outline,
    fill = NA,
    linewidth = 0.5,
    colour = "black" # 黑色轮廓勾勒 NYC 整体形状
  ) +
  
  # --- 标尺设置 ---
  scale_fill_manual(
    values = pal_cols_final,
    breaks = legend_breaks, # 应用自定义图例顺序
    drop   = FALSE          # 强制显示所有定义的图例项 (即使数据中不存在)
  ) +
  
  # --- 地图整饰要素 ---
  # 比例尺 (右下角)
  ggspatial::annotation_scale(
    location   = "br",
    width_hint = 0.25,
    line_width = 0.5,
    height     = unit(0.15, "cm"),
    pad_x      = unit(0.25, "cm"),
    pad_y      = unit(0.25, "cm")
  ) +
  
  # 指北针 (右下角，位于比例尺上方)
  ggspatial::annotation_north_arrow(
    location    = "br",
    which_north = "true",
    style       = ggspatial::north_arrow_nautical,
    height      = unit(1.1, "cm"),
    width       = unit(1.1, "cm"),
    pad_x       = unit(1.0, "cm"),
    pad_y       = unit(1.0, "cm")
  ) +
  
  # --- 标题与标注 ---
  labs(
    title    = "Choropleth map: Poverty rate by Census tract (NYC)",
    subtitle = paste0(
      "Data: ACS ", acs_year,
      " 5-year estimates; Classification: Fisher–Jenks (n = ", n_class, ")"
    ),
    fill     = "Poverty rate\n(Fisher breaks)",
    caption  = "Source: US Census Bureau via tidycensus"
  ) +
  
  # --- 坐标系与主题 ---
  coord_sf(expand = TRUE) +
  theme_bw() +
  theme(
    panel.grid.major  = element_line(linewidth = 0.3, colour = "grey90"),
    axis.title        = element_blank(),
    legend.key.height = unit(0.7, "cm")
  )

制图数据分级策略：为何选择 Fisher-Jenks 算法？

下一节详讲

分层设色图的核心不仅在于配色，更在于连续数值的离散化策略。面对存在极端值的偏态分布数据（如房价、人口），直接线性映射常导致色带被高值“拉伸”，从而掩盖了大部分区域的空间差异。

上述代码采用 Fisher-Jenks （自然断点法）进行分级。该算法通过迭代寻找类内差异最小、类间差异最大的阈值，能最大程度还原数据的自然聚类特征。

原理机制：计算断点 breaks （阈值序列），将连续数值划分为离散区间，同一区间映射为同一颜色；
实现路径：调用 classInt::classIntervals() 获取最优断点，再通过 cut() 函数将连续变量转化为等级变量。

建议 分级数量不宜过多，通常控制在 5–7 级，以在视觉区分度与认知负担之间取得平衡。下一小节将进一步介绍其他常用分类方法。

数据分类分色

4.3 数据分类分色

重要 · 常用

在计算社会科学的研究中，制作分层设色图（Choropleth Map）时，我们通常需要将连续的数值变量划分为离散的色彩等级（Data Classification），而非直接映射为连续的渐变色带。

思考：为什么要把连续数据变成离散等级？

虽然连续色带（Continuous Color Ramp）理论上保留了所有数据细节，但在空间格局的有效表达上往往效果不佳（特别是在社会科学研究中），主要基于以下三个方面：

降低认知负荷（Cognitive Load）
人眼对色彩明度的分辨能力是有限的。读者很难区分“深蓝”与“稍深一点的蓝”具体代表多少数值差异（例如 15.1% 与 15.9%）。将数据归纳为 5–7 个离散等级，能帮助读者快速锁定“高值区”与“低值区”，显著提升读图效率。
对抗数据偏态（Handling Skewness）
社会经济数据（如收入、贫困率）常呈现典型的长尾分布。若使用连续色带，极少数的极端高值会“吃掉”色带中最深的部分，导致剩余 90% 的区域颜色过浅且趋同，从而掩盖内部差异。离散分级（如 Fisher-Jenks）能强制拉开这些区域的色差；在极端情况下，研究者甚至需要手动干预分级，才能让主体数据的空间结构被肉眼识别。
平滑局部噪声（Generalisation）
若直接利用连续色带映射原始数值，邻近单元间微小的数值波动（例如 12.1% 与 12.3%）会导致地图上出现大量细碎的颜色差异，形成类似图像处理中的“椒盐噪声”（Salt-and-pepper noise，即地图看起来像撒了胡椒粉一样斑驳、破碎）。离散分级本质上是一种空间泛化策略，它过滤了这些无关紧要的微小波动，将相近数值合并为同一色块，从而突显出更宏观、连贯的空间集聚模式。

补充值得注意的是，在自然地理与环境科学领域（如高程 DEM、气温场、空气污染浓度），由于变量本身具有强烈的物理连续性且空间分布相对平滑，连续色带依然是主流且科学的表达方式。离散分级主要针对人文社会科学中分布不均、异质性强的社会经济人口指标。

R 语言中的 classInt 包提供了最全面的分级算法支持，直接决定了地图的空间格局表达效果。

核心函数：classIntervals(var, n, style = "...")

var：数值向量（待分级数据）；

n：期望的分级数量（通常 5–7 级）；

style：分级算法名称（如 "fisher", "quantile", "equal"）。

以下介绍三种最常用的分级方法及其适用场景：

1. Fisher-Jenks 自然断点法 (`style = "fisher"`)

2. 等量分位法 (`style = "quantile"`)

3. 等间距法 (`style = "equal"`)

原理机制：
基于数学上的绝对数值跨度（Max - Min），算法将整个数据的取值范围进行等距离切分，目标是确保每个色彩等级跨越的数值区间大小完全相同（类间距均等化）。例如，若数据范围是 0 到 100，分为 5 级，则断点严格固定为 20、40、60、80。

通俗理解：等间距法生成的断点就像是一个等差数列。与前面介绍的等量分位法（Quantile）截然不同，它完全不关心每个色彩区间内最终分配到了多少个观测单元（例如国家），而是极其机械地仅按数值本身的绝对大小进行物理等分。
适用场景：
适用于数据分布非常均匀（近似均匀分布或正态分布），或者指标本身具有绝对的线性物理意义的情况（如气温、降水、海拔高程，或标准的 0-100% 评分）。在大多数呈现典型长尾特征的计算社会科学研究中（如财富分布、人口密度），这种方法不常用。
主要优点：
1. 通俗易懂：分组逻辑最为直接，图例上的数值通常是整齐的等差数列，普通读者极易理解；
2. 便于横向对比：如果人为锁定一个全局的最大值和最小值，那么不同年份的断点将完全一致，非常适合制作时间序列地图（Time-series Maps）进行绝对数值的变化对比。
潜在缺点：
1. 对极端值/偏态分布极度敏感：这是它在社会科学中被弃用的主因。若数据集中存在极少数的高值，这几个高值会把整个跨度拉得极长。结果是 90% 的地理单元都被迫挤在最低的一个或两个颜色等级中，导致地图失去展现内部空间差异的能力（即地图看起来“糊成一片”）。

代码演示：
演示将再次使用 gapminder 数据集，通过等间距法对 2007 年全球预期寿命数据进行分级。

直方图部分

library(tidyverse)
library(gapminder)
library(classInt)
library(RColorBrewer) 

# ============================================================
# 1. 数据准备
# ============================================================
clean_data <- gapminder %>% 
  filter(year == 2007)

target_var <- clean_data$lifeExp

# ============================================================
# 2. 计算等间距断点 (Equal Interval Breaks)
# ============================================================
# 使用 Equal 算法将数据分为 5 类 (确保每个数值区间的跨度完全相等)
equal_res <- classIntervals(target_var, n = 5, style = "equal")

# 提取计算出的断点数值 (用于后续切分和绘图)
breaks_vals <- equal_res$brks 

# ============================================================
# 3. 数据离散化处理
# ============================================================
clean_data <- clean_data %>%
  mutate(life_class = cut(lifeExp, 
                          breaks = breaks_vals,       # 使用等间距断点
                          include.lowest = TRUE,      # 关键：闭合区间，确保包含最小值
                          dig.lab = 3))               # 标签保留3位有效数字

# ============================================================
# 4. 绘制精细控制的直方图
# ============================================================
ggplot(clean_data, aes(x = lifeExp, fill = life_class)) +
  
  # --- 图层1：直方图主体 ---
  geom_histogram(binwidth = 2,        # 组距设为2
                 color = "grey5",     # 柱子边框颜色
                 alpha = 0.9,         # 透明度
                 linewidth = 0.3) +   # 边框线条粗细
  
  # --- 图层2：限高分割线 ---
  # 使用 geom_segment 自定义线条起止点 (y=0 到 y=23)
  geom_segment(data = data.frame(x = breaks_vals[2:5]), # 排除首尾极值，只取中间断点
               aes(x = x, xend = x, y = 0, yend = 23), 
               linetype = "dashed",   # 线型：虚线
               color = "darkred",     # 颜色：深红
               linewidth = 1.2,       # 线条宽度
               inherit.aes = FALSE) + # 关键：不继承主图的 fill 属性，防止报错
  
  # --- 图层3：配色设置 ---
  # 使用 RColorBrewer 色盘 (注意 guide 上大下小)
  scale_fill_brewer(palette = "YlGnBu", 
                    name = "Equal Intervals \n(5 Classes)", # 更新为等间距图例名称
                    guide = guide_legend(reverse = TRUE)) +
  
  # --- 图层4：数值标签 ---
  # 将标签放置在分割线上方 (y=24)
  annotate("text", x = breaks_vals[2:5], y = 24, 
           label = round(breaks_vals[2:5], 1), # 保留1位小数
           fontface = "bold", color = "black", size = 4) +
  
  # --- 坐标轴控制 ---
  scale_y_continuous(
    limits = c(0, 25),            # 强制固定 Y 轴范围 0 - 25个
    expand = expansion(mult = 0)  # 移除 Y 轴底部的留白，让柱子紧贴 X 轴
  ) +
  
  # --- 主题美化 ---
  labs(title = "Equal Interval Breaks on Life Expectancy (2007)",
       subtitle = "Histogram colored by 5 equal-width classes",
       x = "Life Expectancy (Years)",
       y = "Count of Countries") +
  theme_bw() +
  theme(
    legend.position = "right", 
    plot.title = element_text(face = "bold")
  ) +
  
  # --- 坐标系比例 ---
  # 强制 x轴和 y轴显示比例为 1:1
  coord_equal()

地图部分

library(tidyverse)
library(gapminder)
library(classInt)
library(RColorBrewer) 
library(sf)
library(spData) 
library(countrycode) 
library(ggmapcn) # 用于生成贴合地图边界的曲面经纬网与标签

# ============================================================
# 1. 属性数据准备、断点计算与离散化
# ============================================================
clean_data <- gapminder %>% 
  filter(year == 2007) %>%
  # 将国家名称映射为 ISO-2 国际标准代码，以确保后续精准连接
  mutate(iso_a2 = countrycode(sourcevar = country, 
                              origin = "country.name", 
                              destination = "iso2c"))

# 计算等间距断点 (Equal Interval Breaks)
target_var <- clean_data$lifeExp
# style = "equal" 确保将整体数据跨度(最大值-最小值)机械地等分为 5 个物理区间
equal_res <- classIntervals(target_var, n = 5, style = "equal")
breaks_vals <- equal_res$brks 

# 将连续型预期寿命变量转化为离散的等级因子
clean_data <- clean_data %>%
  mutate(life_class = cut(lifeExp, 
                          breaks = breaks_vals, 
                          include.lowest = TRUE, # 确保包含最小值的左闭合区间
                          dig.lab = 3))

# ============================================================
# 2. 空间数据加载、属性挂载与投影转换
# ============================================================
data(world)

# 定义 Robinson 投影的 PROJ 字符串
crs_robin <- "+proj=robin +datum=WGS84"

map_data <- world %>%
  # 基于 ISO-2 代码连接空间与属性数据
  left_join(clean_data, by = "iso_a2") %>%
  # 提前将空间数据转换为 Robinson 投影，以提高渲染效率并确保辅助图层精准对齐
  st_transform(crs = crs_robin)

# ============================================================
# 3. 绘制分层设色地图 
# ============================================================
ggplot(data = map_data) +
  
  # --- 图层 1：绘制地图面要素 ---
  geom_sf(aes(fill = life_class), 
          color = "black",       # 设定国家边界线颜色
          linewidth = 0.3) +     
  
  # --- 图层 2：生成贴合弧线的经纬网与标签 ---
  annotation_graticule(
    crs = crs_robin, 
    lon_step = 30,         # 经度线间隔 (度)
    lat_step = 15,         # 纬度线间隔 (度)
    label_offset = 3e5,    # 标签与地图边缘的偏移量 (依据出图比例微调)
    color = "grey55",      
    linetype = "dashed",
    line_width = 0.2
  ) +
  
  # --- 图层 3：色彩映射 ---
  scale_fill_brewer(palette = "YlGnBu", 
                    name = "Life Expectancy\n(Equal Intervals)",  # 更新为等间距的图例说明
                    guide = guide_legend(reverse = TRUE),         # 图例顺序设为上大下小
                    na.value = "grey85") +                        # 缺失数据的国家填充为浅灰色
  
  # --- 图层 4：坐标参考系 (CRS) 声明 ---
  # datum = NA 用于关闭 ggplot2 默认的直角网格，避免与 annotation_graticule 冲突
  coord_sf(datum = NA) +
  
  # --- 图层 5：主题与排版美化 ---
  labs(title = "Global Life Expectancy (2007)",
       subtitle = "Choropleth map using Equal Interval classification (Robinson Projection)", # 更新副标题说明
       caption = "Data source: Gapminder & spData") +
  theme_void() +           # 移除默认背景、坐标轴线和刻度
  theme(
    legend.position = "bottom",          
    plot.title = element_text(face = "bold", size = 14, margin = margin(b = 5)),
    plot.subtitle = element_text(margin = margin(b = 5)),  # 增加副标题下方的留白
    plot.margin = margin(10, 10, 10, 10)                   # 设置图形整体边距，防止贴边标签被截断
  )

延伸：更多数据分级与自定义策略”

以上介绍的自然断点法 、等量分位法 和 等间距法 是制作分层设色地图时最经典、最常用的。

然而，在面对极其复杂的现实数据时，软件（如 ArcGIS、QGIS 或 R 的 classInt 包）还提供了许多其他的分级选项：

标准差法 (Standard Deviation)：
计算所有数据的平均值和标准差，并以平均值为中心，按 0.5 或 1 个标准差的步长向两端划分区间。它不强调数值的绝对大小，而是非常适合用来突出展示哪些地区“异常偏高”或“异常偏低”。
手动/自定义分级 (Manual / Fixed Breaks)：
在实际的学术研究或数据新闻中，比较常用的分类往往是“算法 + 人工”的结合。例如：
- 细化极值区间：在 Fisher-Jenks 跑出自然断点的基础上，如果最高数值区间内部依然存在极大的悬殊（比如包含了 80 分和 1000 分的异常值），研究者可以通过手动增加断点，将最高区间人为再切分出几个子区间，从而精准地将极高值单独剥离出来。同理，对极低数值区间也可以进行这种人为的细化操作。
- 对齐政策标准：如果研究贫困率，可以直接手动将断点设定为官方的“贫困线”数值(比如官方有明确规定的数值)，而不必理会算法自动切在哪。在 R 的 classInt 包中，可以通过设定 style = "fixed" 并配合具体的 fixedBreaks = c(...) 向量来实现。

📚 推荐参考资料：

ArcGIS Pro 官方文档：在搜索引擎中搜索 “ArcGIS Pro 对数值字段进行分类”，官方中文文档对上述所有算法的数学逻辑和适用场景有着极其权威且图文并茂的解释。
《Geocomputation with R》（R语言地理计算）。这是一本开源的在线教科书，详细讲解了 sf 与各种制图包的高阶联动。

底图绘制

4.4 底图：空间语境

常用

在空间可视化中，仅绘制几何形状（如多边形或散点）往往会使图面丧失地理参照。引入网络瓦片底图（Basemap Tiles）能够提供至关重要的空间语境（Spatial Context）。

高质量的底图不仅能辅助读者建立直观的地理定位与空间连接，还能显著增强制图的专业度与真实感。

制图思路

在叠加底图时，必须严格遵循“数据为主，底图为辅”的视觉层级原则，切忌让底图“喧宾夺主”。高水准的学术底图配置通常遵循以下设计规范：

色彩退隐（Color Muting）：
避免使用色彩丰富或饱和度过高的常规地图（如标准的 Google Map 或高德地图），以免与主图层（如分层设色的色阶）发生严重的视觉冲突。带有透明度处理的灰度底图（Grayscale）或低饱和度的浅色底图是学术制图的首选。
标签克制（Label Control）：
底图中过密的道路名称或地名标签会造成图面杂乱，严重干扰数据读取。建议优先选择无标签（No Labels）或少量标签的极简底图版本。
暗色模式（Dark Theme）：
当主图数据呈现高亮特征（如夜间灯光分布、明亮的点密度发光图）时，采用暗黑风格底图（如 Carto Dark Matter）能够形成强烈的明暗反差。这种基于特殊制图目的的暗调排版，能极大提升图面的现代感与美学张力。

工具包推荐

在基于 sf 与 tidyverse 的空间分析工作流中，加载在线底图通常依赖以下常用集成方案：

ggspatial 包：
作为 ggplot2 体系的生态扩展，其 annotation_map_tile() 函数可作为普通几何图层通过 + 号直接叠加至绘图结构中。该函数无需附加 API 配置，默认提供 CartoDB（如 Light/Dark Matter 灰度与暗色主题）及 OpenStreetMap 底图数据，并能与 sf 主图层的坐标参考系自动对齐。
tmap 包：
专注于专题地图绘制，语法逻辑与 ggplot2 相近并直接兼容 sf 对象。可通过 tm_basemap() 函数叠加多种开源底图。此外，tmap 支持在静态制图与动态交互模式间无缝切换，适用于探索性空间数据分析（ESDA）。
maptiles 包：
适用于静态出版级图件的底图处理。当需要对底图的下载缓存、缩放级别（Zoom Level）与分辨率进行精细控制时，该包可基于 sf 对象的边界范围（Bounding Box）精准裁剪，并生成高质量的本地栅格底图对象。

代码示例 1

本示例展示了基于开源数据获取与底图渲染的制图流程。首先，利用 osmdata 提取东京千代田区咖啡馆POI；针对商业品牌的长尾分布特征，借助 forcats::fct_lump_n() 函数执行数据聚合，保留 Top 7 连锁品牌。制图阶段的核心操作为调用 ggspatial::annotation_map_tile() 嵌入极简在线底图。

# ---- 加载核心工具包 ----
library(tidyverse) 
library(sf)        
library(ggspatial) # 提供底图与空间制图组件（比例尺/指北针）
library(osmdata)   # 提取 OpenStreetMap 开源矢量数据

# ---- 1. 数据获取与预处理 ----

# 构建查询：限定东京千代田区范围，并提取 cafe 设施要素
chiyoda_query <- opq(bbox = "Chiyoda, Tokyo, Japan") %>%
  add_osm_feature(key = "amenity", value = "cafe")

# 执行 API 查询并直接转化为 sf 空间点对象
chiyoda_cafes_raw <- osmdata_sf(chiyoda_query)$osm_points

# 空间属性清洗与分类规约
plot_data <- chiyoda_cafes_raw %>%
  filter(!is.na(name)) %>% # 剔除无效点位
  mutate(
    # 判别逻辑：依据英文品牌字段的缺失状态，划分独立门店与连锁品牌
    Brand_Initial = ifelse(is.na(`brand:en`), "Independent Cafe", `brand:en`),
    
    # 长尾聚合 (Lumping)：保留频次 Top 7 的头部品牌，其余归并为 "Other Chains"
    Cafe_Type = fct_lump_n(Brand_Initial, n = 7, other_level = "Other Chains")
  ) %>%
  # 因子重排：将长尾聚合项与独立门店置于图例底层，优化视觉层级逻辑
  mutate(Cafe_Type = fct_relevel(Cafe_Type, "Other Chains", "Independent Cafe", after = Inf))


# ---- 2. 空间可视化绘制 ----

ggplot() +
  # [1] 底图层：加载极简在线瓦片底图以提供空间语境（需置于主图层前）
  annotation_map_tile(type = "cartolight", zoom = 14, alpha = 0.9) +
  
  # [2] 主数据层：叠加几何要素，并通过颜色映射商业类型
  geom_sf(data = plot_data, aes(color = Cafe_Type), size = 1, alpha = 0.9) +
  
  # [3] 视觉映射：选用 ColorBrewer 离散色盘 (Set1) 增强多类别区分度
  scale_color_brewer(palette = "Set1", name = "Business Type") +
  
  # [4] 制图规范组件：挂载动态比例尺与花体指北针
  annotation_scale(
    location = "bl", width_hint = 0.25, style = "ticks", 
    text_cex = 0.8, text_col = "gray30"
  ) +
  annotation_north_arrow(
    location = "br", which_north = "true", 
    style = north_arrow_fancy_orienteering(
      fill = c("gray40", "white"), line_col = "gray20"
    )
  ) +
  
  # [5] 主题排版与网格视觉调和
  theme_bw() + 
  theme(
    legend.position = "right", 
    legend.title = element_text(face = "bold", size = 9),
    legend.text = element_text(size = 8),
    legend.key.size = unit(0.5, "cm"), 
    plot.title = element_text(face = "bold", hjust = 0.5, margin = margin(b = 10)),
    plot.subtitle = element_text(hjust = 0.5, color = "gray30", margin = margin(b = 15)),
    # 弱化经纬坐标网格，采用浅灰虚线以避免与底图真实路网产生视觉混淆
    panel.grid.major = element_line(color = "gray85", linetype = "dashed"),
    panel.grid.minor = element_blank()
  ) +
  labs(
    title = "Spatial Distribution & Brand Analysis of Cafes",
    subtitle = "Chiyoda Ward, Tokyo (Top Brands Listed)",
    x = "Longitude", y = "Latitude",
    caption = "Basemap: CartoDB Light | Data: OpenStreetMap"
  )

上述代码成果图

代码示例2

本示例展示了基于微观坐标边界与“暗图浅框”混合排版的空间制图流程。为优化获取效率，示例直接传入经纬度边界框，提取上海市中心约 3×3 公里核心区的便利店空间POI。属性重构阶段借助 forcats::fct_lump_n() 聚合长尾品牌，并将分类频次动态拼接入图例标签（如 Brand (n)）以提升信息密度。制图阶段调用 ggspatial::annotation_map_tile() 嵌入极简暗黑底图，配合亮色离散色盘赋予商业热点发光质感；同时保持图表外围的浅色设计，无缝兼容学术期刊的白底排版规范。

# ---- 加载核心工具包 ----
library(tidyverse) 
library(sf)        
library(ggspatial) # 集成底图、比例尺与指北针等制图组件
library(osmdata)   # 提供 OpenStreetMap 开源空间数据提取接口

# ---- 1. 数据获取与预处理 ----

# 构建查询：直接传入经纬度边界框 (Bounding Box) 锁定微观研究区，规避地名解析失败风险
lujiazui_bbox <- c(121.47, 31.22, 121.52, 31.26)
shanghai_query <- opq(bbox = lujiazui_bbox) %>%
  add_osm_feature(key = "shop", value = "convenience")

# 执行查询：提取 API 数据并强制转化为 sf 空间点图层
shanghai_stores_raw <- osmdata_sf(shanghai_query)$osm_points

# 空间属性清洗与高密度图例重构
plot_data <- shanghai_stores_raw %>%
  filter(!is.na(name)) %>% # 剔除无名无效锚点
  mutate(
    # [1] 属性重构：处理 brand 字段缺失值，将其统一归类为“独立门店”
    Brand_Initial = ifelse(is.na(brand), "Independent Store", brand),
    
    # [2] 长尾聚合 (Lumping)：提炼 Top 4 头部品牌，尾部数据规约为 "Other Chains"
    Store_Type = fct_lump_n(Brand_Initial, n = 4, other_level = "Other Chains")
  ) %>%
  # [3] 视觉层级优化：强制将长尾聚合项与独立门店置于因子底端，优化图例排版逻辑
  mutate(Store_Type = fct_relevel(Store_Type, "Other Chains", "Independent Store", after = Inf)) %>%
  
  # [4] 图例信息增容：计算各类别绝对频次 (n)，并动态拼接至图例标签 (例如 "FamilyMart (15)")
  add_count(Store_Type) %>%
  mutate(
    Store_Label = paste0(Store_Type, " (", n, ")"),
    # 因子锁序：利用 fct_reorder 确保拼接后的新标签严格继承原 Store_Type 的内在排序
    Store_Label = fct_reorder(Store_Label, as.integer(Store_Type))
  )


# ---- 2. 空间可视化绘制 ----

ggplot() +
  # [1] 语境基底：渲染 CartoDB 暗黑瓦片底图，zoom=15 适配微观街区尺度
  annotation_map_tile(type = "cartodark", zoom = 15, alpha = 0.9) +
  
  # [2] 主题数据层：映射带有频次信息的重构标签，适当放大要素(size=2)以防被底图吞噬
  geom_sf(data = plot_data, aes(color = Store_Label), size = 2, alpha = 0.9) +
  
  # [3] 色彩映射：采用 ColorBrewer Set2 离散色盘，利用高明度色彩在暗背景中激发荧光质感
  scale_color_brewer(palette = "Set2", name = "Store Brand (Count)") +
  
  # [4] 空间参照组件：比例尺与指北针采用高亮灰白配色，确保在暗调底图上的视觉可见性
  annotation_scale(
    location = "bl", width_hint = 0.25, style = "ticks", 
    text_cex = 0.8, text_col = "gray80", line_col = "gray80"
  ) +
  annotation_north_arrow(
    location = "br", which_north = "true", 
    style = north_arrow_fancy_orienteering(
      fill = c("gray50", "white"), line_col = "gray80", text_col = "gray80"
    )
  ) +
  
  # [5] 混合排版设计 (Dark Map on Light Canvas)
  theme_bw() + 
  theme(
    # 核心区填黑：强制将绘图面板底色涂黑，避免底图加载边界出现生硬白边
    panel.background = element_rect(fill = "#242424", color = NA),
    
    # 外围区留白：图例与标题保持常规浅色系，无缝适配标准学术期刊排版
    legend.position = "right", 
    legend.title = element_text(face = "bold", size = 9),
    legend.text = element_text(size = 8),
    legend.key.size = unit(0.5, "cm"), 
    
    plot.title = element_text(face = "bold", hjust = 0.5, margin = margin(b = 10)),
    plot.subtitle = element_text(hjust = 0.5, color = "gray50", margin = margin(b = 15)),
    
    # 辅助线调和：经纬网格采用低明度虚线，维持坐标参照的同时不破坏整体暗调氛围
    panel.grid.major = element_line(color = "gray40", linetype = "dashed"),
    panel.grid.minor = element_blank()
  ) +
  labs(
    title = "Spatial Distribution of Convenience Stores",
    subtitle = "Lujiazui & The Bund, Shanghai (Dark Map on Light Canvas)",
    x = "Longitude", y = "Latitude",
    caption = "Basemap: CartoDB Dark Matter | Data: OpenStreetMap"
  )

上述代码成果图

交互式地图

4.5 交互式地图（概览）

拓展 · 进阶

交互式地图可作为成果展示(特别是进行ppt汇报)的补充方案。其优势在于支持缩放与漫游、鼠标悬停标签、点击弹窗与高亮反馈，从而在浏览过程中即时查看空间单元的属性信息；同时可与在线底图叠加，更适合网页发布、Dashboard 展示与面向公众的可视化传播。

交互式分层设色图与静态地图的核心一致：仍需完成“指标 → 分级区间 → 颜色”的映射，以保证色阶层次清晰、读图更稳定（分级方法可参考上一小节的 breaks 思路）。

本节以 leaflet 作为示例框架：

leaflet 提供基于 HTML 的交互地图容器，可直接在 R Markdown knit 后输出为可交互页面；
sf 负责边界几何与 CRS 管理；交互地图通常统一到 EPSG:4326（经纬度）以确保底图正确对齐；
分层设色通过 addPolygons() 映射面要素颜色，并可进一步配置：
- label：鼠标悬停显示属性标签；
- popup：点击弹出更完整的属性信息；
- highlightOptions：悬停时高亮边界，增强交互可读性；
在线底图可通过 addProviderTiles() 调用公开底图（如 CartoDB DarkMatter）。

进阶提示 leaflet 的交互输出最终运行在浏览器环境中，常涉及 HTML/CSS/JavaScript 的渲染逻辑（例如标签与弹窗的 HTML 片段）。本章以“能用、能读、能解释”为目标，理解基本机制即可；更深度的样式定制与交互开发可在后续专题中展开。

# ============================================================
# 交互式分层设色图（ACS 5-year）：家庭收入中位数（纽约州 Census tract）
# - 地理尺度：Census tract（单州，下一级尺度更细，更能展示州内差异）
# - 分级方式：Fisher breaks（classInt）
# - 可视化  ：leaflet + 暗色底图（CartoDB.DarkMatter）
# ============================================================

library(sf)
library(tidyverse)
library(tidycensus)
library(tigris)
library(leaflet)

library(classInt)       # ★Fisher 分级断点
library(RColorBrewer)   # 调色板
library(htmltools)      # HTML()

options(tigris_use_cache = TRUE)
options(scipen = 9999)


# ------------------------------------------------------------
# 1) 参数：年份、州、分级数
# ------------------------------------------------------------
year_sel  <- 2023
state_sel <- "NY"
n_class   <- 7

# ------------------------------------------------------------
# 2) 获取 ACS：家庭收入中位数（B19013_001）
#    - geometry = TRUE：直接返回 sf（属性 + 边界）
# ------------------------------------------------------------
inc_tract <- get_acs(
  geography   = "tract",
  state       = state_sel,
  variables   = "B19013_001",
  survey      = "acs5",
  year        = year_sel,
  geometry    = TRUE,
  cache_table = TRUE
) %>%
  transmute(
    GEOID,
    NAME,
    med_hh_income = estimate,   # 家庭收入中位数（USD）
    moe           = moe,        # 误差范围（MOE）
    geometry
  ) %>%
  st_transform(4326)            # leaflet 底图：经纬度（EPSG:4326）

# ------------------------------------------------------------
# 3) Fisher 分级（classInt）
#    - 仅用有效值生成断点，避免 NA/Inf 干扰
# ------------------------------------------------------------
x_ok <- inc_tract$med_hh_income %>%
  as.numeric() %>%
  (\(z) z[is.finite(z) & !is.na(z)])()

brks <- classInt::classIntervals(x_ok, n = n_class, style = "fisher")$brks
brks <- unique(brks)  # 极端情况下避免重复断点
n_bins <- length(brks) - 1

# 若断点重复导致 n_bins 变小，至少保证有 3 级（否则颜色会过少）
if (n_bins < 3) {
  brks <- classInt::classIntervals(x_ok, n = 5, style = "quantile")$brks
  brks <- unique(brks)
  n_bins <- length(brks) - 1
}

# ------------------------------------------------------------
# 4) 颜色映射：按“断点区间”设色（低值浅，高值深）
#    - 生成“刚好 n_bins 个颜色”，避免长度不匹配
# ------------------------------------------------------------
base_cols <- RColorBrewer::brewer.pal(9, "Blues")
pal_vec   <- grDevices::colorRampPalette(base_cols)(n_bins)

pal <- leaflet::colorBin(
  palette  = pal_vec,
  bins     = brks,
  domain   = inc_tract$med_hh_income,
  na.color = "grey55"
)

# ------------------------------------------------------------
# 5) 预生成 label / popup（关键：不要在 formula 里 lapply）
#    - 先把每个要素的字符串准备好，再转为 HTML 列表
# ------------------------------------------------------------
inc_tract2 <- inc_tract %>%
  mutate(
    label_html = if_else(
      is.na(med_hh_income),
      paste0(NAME, "<br/>", "NA"),
      paste0(NAME, "<br/>$", formatC(med_hh_income, format = "f", digits = 0, big.mark = ","))
    ),
    popup_html = paste0(
      "<strong>", NAME, "</strong><br/>",
      "Median HH income: ",
      if_else(is.na(med_hh_income), "NA", paste0("$", formatC(med_hh_income, format = "f", digits = 0, big.mark = ","))),
      "<br/>",
      "MOE: ",
      if_else(is.na(moe), "NA", paste0("±$", formatC(moe, format = "f", digits = 0, big.mark = ","))),
      "<br/>",
      "GEOID: ", GEOID
    )
  )

lab_list   <- lapply(inc_tract2$label_html, htmltools::HTML)
popup_list <- lapply(inc_tract2$popup_html, htmltools::HTML)

# ------------------------------------------------------------
# 6) 图例（Legend）顺序：上大下小（符合“由下到上从少到多”的读图习惯）
#    - leaflet 图例是“从上到下依次列出”
#    - 反转 labels/colors：让高值在上、低值在下
# ------------------------------------------------------------
legend_labels <- paste0(
  "$", formatC(brks[-length(brks)], format = "f", digits = 0, big.mark = ","),
  " – $", formatC(brks[-1],         format = "f", digits = 0, big.mark = ",")
)

legend_labels_rev <- rev(legend_labels)
legend_colors_rev <- rev(pal_vec)

# ------------------------------------------------------------
# 7) leaflet 交互地图
# ------------------------------------------------------------
m <- leaflet(inc_tract2) %>%
  addProviderTiles(providers$CartoDB.DarkMatter) %>%
  addPolygons(
    fillColor    = pal(inc_tract2$med_hh_income),
    fillOpacity  = 0.82,
    color        = "#f2f2f2",
    weight       = 0.25,
    opacity      = 0.7,
    smoothFactor = 0.15,
    
    # 明确打开交互（保险）
    options = pathOptions(interactive = TRUE),
    
    # 悬停标签（hover）
    label = lab_list,
    labelOptions = labelOptions(
      direction = "auto",
      textsize  = "12px",
      opacity   = 0.95,
      sticky    = TRUE
    ),
    
    # 点击弹窗（click）
    popup = popup_list,
    
    # 高亮（hover 时描边变亮）
    highlightOptions = highlightOptions(
      weight = 1.0,
      color  = "#ffd54f",
      bringToFront = TRUE
    )
  ) %>%
  addLegend(
    colors   = legend_colors_rev,
    labels   = legend_labels_rev,
    title    = "Median HH income (USD)",
    opacity  = 0.9,
    position = "bottomright"
  )

# 在 R Script 里务必显式打印对象，确保 Viewer 渲染
m

交互式地图在学术论文的常规图件呈现中并非主流形式（多数期刊仍以静态图为主），因此本节不再深入展开。
若对交互式空间可视化有进一步兴趣，可检索 leaflet 、tmap 包的用法，并结合 AI 辅助探索交互式地图与动态可视化的常见实现技巧。

5.0 制图规范与美化

5.0 规范与美化指南

规范导航

地图出图规范与美化图谱

“科学严谨，专业呈现”
极为重要

本节将系统梳理地图的出图规范与版式美化原则。首先明确学术制图的基础要素配置，并重点解析中国公开地图出图的合规红线与标准表示方法；在此基础上，进一步探讨如何通过视觉层级控制、空间拥挤消解等高阶排版技巧，显著提升空间可视化的可读性、专业度与美学质感。

点击下方卡片，快速跳转至相应模块：

制图规范

通用要素与合规要求

地图美化

视觉层级与高级版式

💡 提示：学习完一个小节后，请再次点击 屏幕右下角的规范导航按钮回到本导航页

制图规范

5.1 出图规范

5.1.1 通用地图出图规范

重点

地图不仅是“展示”，也是一种证据表达。若制图要素缺失或表达不当，容易引发误读（如尺度误导、分类误导、符号含义不清），甚至在研究结论与政策解读中造成争议。因此，遵循通用出图规范的意义在于：确保地图表达准确、可核查、可比较，并降低读者因版式与编码差异产生的理解偏差。

实践经验 避免信息重复
如果文章中的题注(Figure caption) 已经完整交代地图的主题、数据来源与年份，那么地图画面内通常无需再重复“标题/来源/年份”等文字，以免分散注意力与造成版式拥挤。

建议按以下逻辑检查地图要素与表达的规范：

地图范围与定位：
明确研究区范围（study area），并确保与研究问题一致。必要时添加定位小图（locator map）帮助读者建立空间位置感与尺度层级。

在 SCI/SSCI 论文的地图中，常见做法是用“闵行区 → 上海市 → 中国”的多级定位组合图，清晰交代研究位置与空间背景。

图例：
图例（Legend）应清楚说明变量含义与单位（如 count/km²），并保证符号、色阶与图面表达一致。图例顺序建议符合常规读图习惯（低值在下，高值在上）。
若为分级设色（离散分组），建议在图例或题注中注明分级方法（如 Natural breaks / Quantile / Equal interval 等），以增强可解释性与可复现性。

缺失值区域应使用区别于主配色的颜色或纹理单独标识：在彩色设色图中常用不同深浅的灰色（如 grey55 / grey85）表示 NA，避免被误读为“低值”或“背景”。

比例尺、方向与经纬网：
比例尺（Scale）与方向（通常指指北针 North）并非每张图都必须，但在多数研究型地图中建议优先考虑加入。
- 局部尺度、需要距离感知时建议给出比例尺；
- 全球尺度或投影变形明显的场景，可弱化或不放比例尺，以避免误导。
- 比例尺/指北针/经纬网属于非主体元素，必须避免遮挡专题信息：如比例尺压住关键地块颜色，经纬网线过粗“盖住”主图层等（应控制图层顺序与视觉权重）。
指北针与经纬网有时可同时出现，但需注意避免过多的信息重复：经纬网已隐含南北方向时，通常保留其一即可；若指北针用于视觉引导或版式平衡，应采用简洁样式并放置在不干扰主体的位置。

数据来源、时间口径与底图署名：
数据来源与年份/时期应至少在题注或图面某处出现一次（择一即可，避免重复）。
若叠加外部底图服务（OSM、Carto、Google、Esri、高德等），应在图面或题注的合适位置标注底图来源与服务商署名（避免遗漏版权/署名要求）。

注记与标注：
注记与标注（Labels & Annotations）的原则是服务于读图而非装饰：控制数量与层级，避免密集重叠。尤其在分层设色图中，过多地名/编号可能遮挡色块（破坏主体信息表达）。必要时可仅标注少数关键区域，或使用外置标注/缩放图替代密集覆盖。

可视化一致性：

(参考第三章内容)

同一论文/报告内多张地图应尽量统一字体体系、字号层级、配色逻辑、图例位置、边界线宽与缺失值颜色，以保证跨图对比有效，并提升研究输出的整体性与专业感。

5.1.2 中国地图出图合规指南

(极为重要！)

在中国进行学术出版或公开展示时，地图内容的正确性是必须坚守的底线，而涉及国家版图的错误表达则是绝对不可触碰的红线。不规范的地图使用不仅会导致论文被拒稿、报告被驳回，严重时甚至涉及法律责任。

根据《地图管理条例》、《公开地图内容表示规范》等文件，以下是使用 R 语言绘制中国地图时的核心合规要点。

官方标准

严禁“自行改绘”国界线。 涉及中国全图、南海诸岛、台湾省、藏南地区、钓鱼岛等敏感区域时，必须严格遵循自然资源部发布的标准地图或界线画法。

官方合规参照体系：

核心要素的正确表示

在绘制中国全图时，以下要素必须完整呈现，不得遗漏或错误表达：

完整的领土范围：
- 必须包含：中国大陆、海南岛、台湾岛。
- 不可遗漏： 南海诸岛（包括十段线）、钓鱼岛及其附属岛屿、赤尾屿等重要岛屿。
南海诸岛的表示方式：
- 主图包含： 若主图范围足够大（如竖版地图），应直接在主图中表示南海诸岛。
- 附图表示（Inset）： 若使用横版地图，主图南部通常截取至海南岛最南端，并必须在右下角配置南海诸岛附图。附图必须包含“十段线”及主要岛礁。
台湾省的表示：
- 台湾省必须按省级行政单位表示。
- 在制作分省填色图（Choropleth Map）时，若台湾省数据缺失，严禁将其留白或与其他国家颜色一致。必须将其归类为国内省份颜色，或在图例/Caption中明确标注“台湾省资料暂缺”。

坐标系与投影

在 R (sf, terra) 中处理中国地图数据时，应遵循国家测绘标准：

坐标系（CRS）：
- 法定大地基准：2000国家大地坐标系 (CGCS2000)。
- EPSG 代码：推荐使用 EPSG:4490 (CGCS2000, Lat/Lon)。
- 注：虽然 WGS84 (EPSG:4326) 是国际通用标准，但在国内正式出版物中，标注为 4490 更为规范。两者在一般可视化精度下差异极小，但法理意义不同。
投影建议：
- 为了减少面积变形，中国全图通常采用 Albers 等面积圆锥投影。
- 常用参数：+proj=aea +lat_1=25 +lat_2=47 +lon_0=105 +x_0=0 +y_0=0 +ellps=GRS80 +units=m +no_defs (中心经线 105°E)。

辅助要素与图面规范

依据地图学规范及出版习惯，建议如下配置：

比例尺 (Scale Bar)：
- 必选要素。 一般使用线段比例尺。
- 单位：千米 (km)。
- 长度：通常控制在图幅宽度的 1/5 至 1/10 (约 1-3cm)，避免喧宾夺主。
经纬网与指北针 (Graticules & North Arrow)：
- 经纬网： 若主图绘制了经纬网，附图（南海诸岛）也必须绘制。单位应为“度、分、秒”或十进制度。为了图面整洁，刻度值通常只需标注在图框的左侧和下侧（或单边）。
- 指北针： 在拥有经纬网的地图中，经线指示南北，纬线指示东西，因此原则上可以省略指北针。若地图进行了非正向旋转，则必须加注指北针。
审图号 (Map Approval Number)：
- 若地图用于公开发表（期刊论文、书籍、公开报告），必须通过自然资源主管部门的地图审核，并标注审图号（如：GS(2023)xxxx号）。
- 学术类非正式图表建议： 在图注中声明“底图基于自然资源部标准地图服务网站下载的审图号为 GS(20xx)xxxx 号的标准地图制作，底图边界无修改”。

敏感信息避让

公开地图不得表示以下内容：

军事设施、军事禁区、国家安全要害部门。
涉密的水利、交通、电力、通信等设施的属性信息。
精度限制：地面遥感影像分辨率不得优于0.5米等。

自查清单

全图范围： 是否包含了南海诸岛（附图）、钓鱼岛、赤尾屿？
界线画法： 国界线是否使用了标准矢量数据？有无随意平滑或简化导致边界变形？
台湾省： 是否正确着色？若无数据，是否标注了“资料暂缺”？
文字注记： 是否避让了国界线？涉外地名是否符合官方译名？
引用声明： 是否在 Caption 中注明了数据来源及底图审图号？

具体注意细节请参考：2023年第1号自然资源部关于印发《公开地图内容表示规范》的通知 (以最新版本为准)

中国矢量地图数据资源库

在使用 R 进行空间数据可视化时，可从以下平台获取基础矢量边界（GeoJSON/Shapefile）：

注意： 必须核实数据来源的权威性，确保国界与争议领土表达符合自然资源部标准。

代码示例

本节代码将演示如何基于 sf, ggplot2 绘制一幅完整的中国省级人口分层设色地图。我们通过 API 调用 geojson.cn 的省级行政边界数据，并结合第七次全国人口普查的部分截取数据（仅包含部分省份作为演示）进行可视化。

本示例重点展示的技术点：

子图拼接： 针对中国地图绘制中必须包含的“南海诸岛等”要素，演示如何利用 cowplot 包的 draw_plot() 函数，通过图层叠加的方式，将南海诸岛以“画中画”形式嵌入主图右下角。

合规性免责声明

本案例仅作为代码技术演示

所使用的底图边界来自开源接口，未经自然资源部审核，未标注审图号

# ==============================================================================
# 系统环境配置与空间制图依赖库加载
# ==============================================================================
library(sf)          
library(tidyverse)   
library(cowplot)     # 复合图形排版 (用于实现主副图的空间嵌套/Inset Map)
library(ggspatial)   # 制图规范组件库 (提供标准比例尺与指北针)
library(units)       # 空间测度单位换算支持
library(scales)      # 视觉通道标度与数值格式化管理
library(classInt)    # 空间数据单变量离散化与自然断点计算 (Fisher-Jenks)

# ==============================================================================
# 1. 空间对象实例化与属性挂载 (Spatial Instantiation & Attribute Join)
# ==============================================================================

# 1.1 获取基底多边形并执行投影对齐
china_prov <- read_sf("https://geojson.cn/api/china/100000.topo.json") %>%
  st_set_crs(4326) %>%
  # [制图规范] 强制转换为中国国家大地坐标系 (CGCS2000, EPSG:4490)，保障本土几何测算准确性
  st_transform(4490)

# 1.2 构建截面统计数据表 (2020年人口规模，单位：万人)
province_pop <- tibble::tribble(
  ~province, ~pop_2020_wan, ~pop_2019_wan,
  "河北省", 7461, 7447, "山西省", 3492, 3497, "辽宁省", 4259, 4277,
  "吉林省", 2407, 2448, "江苏省", 8475, 8469, "浙江省", 6457, 6375,
  "安徽省", 6103, 6092, "福建省", 4154, 4137, "江西省", 4519, 4516,
  "山东省", 10153, 10106, "河南省", 9937, 9901, "湖北省", 5775, 5927,
  "湖南省", 6644, 6640, "广东省", 12601, 12489, "海南省", 1008, 995,
  "四川省", 8367, 8351, "贵州省", 3856, 3848, "云南省", 4721, 4714,
  "陕西省", 3953, 3944, "甘肃省", 2502, 2509, "青海省", 592, 590,
  "黑龙江省", 3185, 3255
)

# 1.3 属性表连接 (Attribute Join)
china_prov_pop <- china_prov %>%
  left_join(province_pop, by = c("fullname" = "province"))

# ==============================================================================
# 2. 空间统计分级与数据离散化 (Fisher-Jenks Classification)
# ==============================================================================

# 提取有效属性向量 (剔除 NA 避免算法失效)
pop_values <- na.omit(china_prov_pop$pop_2020_wan)

# 2.1 计算自然断点
# 依托 Fisher-Jenks 算法寻找数据内在断层，最大化组间异质性与组内同质性
fisher_breaks <- classIntervals(pop_values, n = 5, style = "fisher")$brks
brks <- round(fisher_breaks, 0) # 取整以保障图例标量界限的整洁度

# 2.2 构造严谨的数学区间标签 (Mathematical Interval Labels)
# 首区间采用闭区间 [a, b]，后续区间采用左开右闭 (b, c]
labels_vec <- paste0(
  ifelse(seq_along(brks[-1]) == 1, "[", "("),
  brks[-length(brks)], ", ", brks[-1], "]"
)

# 2.3 连续变量因子化 (Discretization)
china_prov_pop <- china_prov_pop %>%
  mutate(pop_class = cut(pop_2020_wan, 
                         breaks = brks, 
                         labels = labels_vec, 
                         include.lowest = TRUE))

# ==============================================================================
# 3. 视觉通道映射与视窗约束 (Visual Scale & Bounding Box)
# ==============================================================================

# 3.1 预设离散型分层设色比例尺 (Choropleth Color Scale)
fill_scale_discrete <- scale_fill_brewer(
  palette  = "Reds",
  na.value = "grey85", # 设缺失值(如港澳台区域)为中性灰
  name     = "2020 Population\n(10k, Fisher-Jenks)",
  # [核心修复] 强制倒置图例排列，确保高值在上、低值在下，顺应空间读图直觉
  guide    = guide_legend(reverse = TRUE)
)

# 3.2 定义空间视窗边界 (Spatial Extents)
bb_main <- c(xmin = 73,  xmax = 135.5, ymin = 18,  ymax = 54) # 中国大陆主图范围
bb_scs  <- c(xmin = 105, xmax = 125,   ymin = 3,   ymax = 26) # 南海诸岛附图范围

# ==============================================================================
# 4. 空间图层构建 (Map Layer Construction)
# ==============================================================================

# --- 4.1 渲染核心主图 (Main Extent) ---
p_main <- ggplot() +
  geom_sf(data = china_prov_pop, aes(fill = pop_class), colour = "grey35", linewidth = 0.25) +
  coord_sf(xlim = c(bb_main["xmin"], bb_main["xmax"]), 
           ylim = c(bb_main["ymin"], bb_main["ymax"]), 
           expand = TRUE) +
  fill_scale_discrete +
  theme_bw() +
  theme(
    plot.background   = element_rect(fill = "white", colour = NA),
    plot.margin       = margin(5, 5, 5, 5),
    legend.key.height = unit(0.7, "cm")
  ) +
  # 植入空间测度制图要素
  ggspatial::annotation_scale(location = "bl", width_hint = 0.25, style = "ticks") +
  ggspatial::annotation_north_arrow(location = "tl", which_north = "true", 
                                    style = ggspatial::north_arrow_fancy_orienteering,
                                    height = unit(1, "cm"), width = unit(1, "cm"))

# --- 4.2 渲染附属插图 (Inset Map: South China Sea) ---
p_scs <- ggplot() +
  geom_sf(data = china_prov_pop, aes(fill = pop_class), colour = "grey35", linewidth = 0.15) +
  # 严格锁定视窗并剥离拓展空白 (expand = FALSE)，保障附图紧凑性
  coord_sf(xlim = c(bb_scs["xmin"], bb_scs["xmax"]), 
           ylim = c(bb_scs["ymin"], bb_scs["ymax"]), 
           expand = FALSE) +
  fill_scale_discrete +
  theme_bw() +
  theme(
    panel.grid.major = element_line(colour = "grey85", linewidth = 0.25),
    panel.grid.minor = element_blank(),
    axis.title = element_blank(), axis.text = element_blank(), axis.ticks = element_blank(),
    panel.border    = element_rect(fill = NA, colour = "grey25", linewidth = 0.5),
    plot.margin     = margin(0, 0, 0, 0),
    legend.position = "none" # 附图隐去图例以避免冗余
  )

# ==============================================================================
# 5. 主副图空间嵌套排版 (Layout & Composition)
# ==============================================================================

# 利用 cowplot 基于相对坐标体系 (0~1) 将南海诸岛附图悬浮嵌入主图右下侧
map_with_inset <- ggdraw() +
  draw_plot(p_main) +
  draw_plot(p_scs, x = 0.61, y = 0.22, width = 0.17, height = 0.17)

# 渲染最终复合出版图
map_with_inset

拓展阅读与实战指南：R语言中国地图制图教程

在实际的研究与发表中，获取合规的地理边界数据并使用 R 语言实现完美的主副图嵌套（如南海诸岛附图的规范排版）往往存在一定的技术门槛。

为了更好地掌握中国地图的数据获取、投影坐标处理以及高阶出图技巧，强烈建议参考以下优秀的中文社区教程与实战记录：

1. 官方空间数据获取与解析
重点讲解如何利用 R 语言（如 sf 包）直接读取阿里云 DataV 或民政部等权威接口的 GeoJSON 格式标准边界数据：

2. 南海诸岛及九段线附图（Inset Map）嵌套排版技巧
重点讲解如何利用 sf 进行边界裁剪（Bounding Box），并结合 ggplot2 与 cowplot 包实现大陆主图与南海小地图的空间画中画（Inset）排版，这是学术论文出图的刚需技术：

3. sf 包绘制中国地图综合实战

以下知乎与微信公众号文章提供了更为丰富的代码范例，涵盖了分层设色、标签添加、数据离散化及各种视觉美化细节，适合作为本章练习的代码参考：

制图规范

5.2 地图美化

可参考第三章可视化美化内容

在第三章中，本教程已系统探讨了数据可视化的通用美学原则。然而，地图包含复杂的地理基底与多维度的空间拓扑关系，其美化逻辑具有极强的特殊性。

必须明确的核心前提是：所有地图的美化，都必须建立在绝对严谨的“制图规范”基础之上。

精心打磨的地图“颜值”并非仅为视觉愉悦。在学术报告、政策简报或商业竞标中，高品质的空间可视化能够显著提升作品的专业度与可信度，一定程度上可提高研究的录用倾向或项目的中标率。

针对空间制图，核心的美化与版式原则主要包括以下三个维度：

1. 视觉层级与底图退隐

与普通统计图表不同，地图的画板通常是广袤且复杂的地理空间。

在常见的空间制图中，必须确保核心数据图层处于视觉最前端。底图（如行政边界、地形或路网）仅提供空间参考，建议使用低饱和度、低对比度的中性色调（如浅灰色 #F5F5F5 或暗黑色 #1A1A1A）予以视觉弱化。同时，面要素的边界线（Borders）宜采用极细且贴近底色的线条（或直接取消边界），以有效消解画面的割裂感。

在高精度微观地理尺度下（如英国 OA、美国 Census Tract 或 Block Group 级别），建议直接取消面要素的边界绘制，这能有效避免密集边界线对分层设色视觉表达的严重干扰。

2. 空间拥挤的视觉消解

当进行空间大数据制图时（如海量的微观定位点或 OD 轨迹线），若直接将原始要素映射到图面上，在局部高度聚集的区域极易产生空间拥挤与遮挡（Overplotting）。这种缺乏修饰的粗暴映射不仅会导致底图被严重遮挡、画面糊成一团丧失美感，更会掩盖数据本身的绝对密度差异，使得图面完全无法有效传达出空间分布的“冷热点”信息。

为了精准展现真实的空间聚集特征并兼顾制图美学，推荐采用以下策略：

透明度叠加与暗场发光：
在暗色背景基底上，为数据要素设置极高的透明度与高明度色彩。依托海量点位或轨迹线的物理重叠与色彩累加，自然透出高亮光晕（即“伪热力发光”效果）。此方法能够在完全保留原始空间结构的同时，通过亮度差异直观反映出热点区域。

$Figure 13 透明度叠加与暗场发光示意图$

Figure 13 透明度叠加与暗场发光示意图

既有空间单元聚合（Spatial Aggregation）：
将微观离散点位直接归集至既有行政或统计单元（如英国的 OA、LSOA 或美国的 Census Tract）。通过执行空间连接（Spatial Join）或拓扑相交运算（st_intersects）汇总地块内的要素频数。为消除各单元因面积或人口基数悬殊造成的统计偏差，强烈建议将基础频数转化为空间密度或人均指标，再进行分层设色图的渲染。

Figure 14 空间单元聚合示意图

几何网格聚合（Tessellation & Binning）：
果断放弃散点映射，将地理空间划分为标准化的聚合单元（如六边形或正方形密铺网格），并统计网格内包含的要素频数或属性指标。辅以连续或离散的渐变色带进行分层设色，能有效规避相互遮挡，极大提升分布密度的精确可读性与学术严谨度。

Figure 15 几何网格聚合示意图

核密度估计（Kernel Density Estimation, KDE）：
运用空间统计学算法，将离散的点集直接转化为平滑、连续的密度表面分布图。此方法能极其精准且柔和地刻画出空间“热点”的渐变轮廓与核心区范围（注：核密度图通常涉及 raster 栅格数据格式的处理，本节仅作概念提及，不作代码展开）。

Figure 16 核密度估计示意图

3. 制图辅助要素的主题融合与版式一致性

地图的三大辅助要素（比例尺、指北针、图例）不仅是测绘规范组件，更是排版设计的有机组成部分。首先，需根据地图底色动态适配其视觉对比度。其次，对于长条形的连续型或离散化统计图例，推荐采用底部水平居中的排版方式，这能与地图整体的横竖向比例形成视觉配平，避免侧边垂直放置挤占核心图幅的有效空间。

此外，尽管 ggspatial 等扩展包提供了种类繁多的指北针与比例尺样式，但在同一项研究或同一份报告中，务必保持出图风格的高度一致性。这意味着不同地图间的辅助要素样式、所处位置、颜色搭配以及经纬线风格应严格统一，这也精准契合了数据可视化中对颜色与字体一致性的核心要求（参考第三章内容）。

扩展学习与灵感资源

空间制图美学需要持续的阅图积累与代码实践。以下精选了涵盖技术指南、制图工具与前沿灵感的高质量资源：

1. R 语言空间制图权威指南
- Geocomputation with R（第 9 章：Making maps）：R 语言空间分析开源教科书，系统讲解了静态与交互式高级地图的绘制规范。
- ggplot2: Elegant Graphics for Data Analysis（第 6 章：Maps）：ggplot2 官方指南地图专章，深度剖析空间矢量渲染与投影几何逻辑。

2. 经典工具与学术范例
- ColorBrewer 2.0 (colorbrewer2.org)：专为地图设计的经典配色基准工具，支持“色盲友好”与“打印友好”等科学筛选。
- DataShine (datashine.org.uk)：英国 UCL 团队开发的人口普查可视化平台，是微观尺度（如 OA 级别）高精度分层设色地图的顶级参考。

3. 制图美学与艺术灵感
- London: The Information Capital（James Cheshire & Oliver Uberti）：城市空间可视化的标杆级专著，展示了如何将海量城市数据转化为极致的制图艺术。
- David Rumsey Map Collection (davidrumsey.com)：全球最大历史古地图库，可从中汲取传统制图学的文本排版与古典色彩美学沉淀。
- #30DayMapChallenge（GitHub/X 检索同名标签）：每年 11 月由全球空间数据社区发起的开源制图挑战，汇集大量前卫设计灵感与复现代码。

6.0 本章练习

6.0 动手实战

实战目标

本节练习使用 Inside Airbnb 提供的开源真实房源数据与 spData 包中的大伦敦行政区数据，进行多图层与多尺度的空间可视化分析：

（1）以网络读取方式获取房源观测值，清洗并转化为空间点要素（st_as_sf），与面要素执行地理坐标对齐（st_transform）；
（2）出图 1 绘制暗夜发光风格的点分布地图，对高密度空间点要素进行层次化渲染（点径与透明度控制）；
（3）出图 2 构建六边形空间网格（Hexagonal Tessellation），执行空间拓扑关联（Point-in-Polygon）与均值过滤，并应用 Fisher-Jenks 自然断点法绘制学术规范的分层设色图（含几何计算与视觉呈现的坐标系逆向对齐逻辑）。

将巩固以下能力：
- 空间数据流处理：远程数据抓取、二维坐标映射、CRS 投影系统转换与逆向对齐。
- 空间统计与聚合：标准几何网格生成、基于空间相交特征的频数计算与阈值过滤。
- 制图美学与规范：高密度数据的过度绘制规避策略、离散数据分级表达、多主题（暗场/亮场）图面排版定制。

6.1 数据与研究问题设定

数据来源：
1. Inside Airbnb 全球短租房源开源数据（本练习选取大伦敦地区，含经纬度坐标的 .gz 完整版房源切片数据）。
2. spData::lnd 大伦敦区行政边界矢量面数据。

数据说明：Inside Airbnb 与 `lnd` 面数据

Inside Airbnb (https://insideairbnb.com/get-the-data/) 是一个独立的非商业开源项目，旨在探究并揭示短期租赁平台对全球城市住宅社区的真实影响。该平台定期抓取并发布全球主要城市的详尽短租数据，字段涵盖完整的房源属性、空间坐标、日历可用性及用户历史评论等（本练习仅调用伦敦地区最新一期的全量房源属性与坐标数据）。此类高精度微观数据目前已被广泛应用于城市空间结构、住房政策评估等领域的量化研究中。
lnd 是 spData 教学包中内置的经典空间数据集，包含了伦敦各行政区（Boroughs）的边界多边形。

研究任务：

制图任务1｜微观空间形态与分类映射（点分布图）：不同租赁类型（room_type）在城市空间中呈现何种微观分布特征？

本任务要求清洗网络读取的源数据并生成 sf 对象。在暗色基底上，应用“深空星云”配色方案直接映射房源分布。通过极小点径与高透明度配置，自然生成类似热力聚集的“发光”图面。
制图任务2｜空间网格离散化与聚合统计（六边形密铺分层设色图）：房源在标准物理空间单元（约 0.5 km² 的六边形网格）内的聚集强度如何？重要

本任务需将底层数据转换至投影坐标系（EPSG:27700）以生成精确面积的几何网格。计算网格内的包含要素数量，过滤低于均值的单元，并采用 Fisher-Jenks 法进行 7 级离散化。最终需将网格逆向转换回地理坐标系（WGS 84），配合纯净明亮的学术主题输出空间密度图。

这是一个“照着做”的练习：先复现图形，再逐行理解每一层空间映射、拓扑运算与主题定制代码的底层逻辑。

请在 R Project 中新建一个 R Markdown 文档：

文件名建议：Lab4_spatial_vis.Rmd
核心要求：
1. 加载 tidyverse、sf、spData、ggspatial 与 classInt 包（未安装请先安装）。
2. 使用 read_csv() 直接读取 Airbnb 的 .gz 网络数据，去除坐标缺失值，完成初始空间要素实例化与底图坐标对齐。
3. 完成图 1 绘制：叠加纯黑底色、暗灰边界与荧光色点要素，配置暗场主题（theme_void 深度定制）及相应制图标尺。
4. 完成图 2 绘制：
  - 执行投影转换生成六边形网格（st_make_grid）。
  - 利用 st_intersects 计算聚合数量，并保留高于均值的数据。
  - 使用 classIntervals 提取自然断点并分箱。
  - 逆向投影至 WGS 84 后，采用明亮学术主题与热力色偏映射离散密度，确保横向图例排版。
5. 生成 html 输出，并检查两幅地图的学术规范性。

6.2 可复现模板

注意

从第三章开始，练习不再提供完整的 .Rmd 成品模板。
练习要求是在自建的 R Markdown 文档中自行插入代码块并运行（可参考第2章的格式配置）。
请参考本章前述小节的核心代码逻辑，逐步构建并完成两类不同风格与分析维度的城市空间制图。

几何拓展：如何推算六边形的 `cellsize`？

在执行空间离散化时，我们往往需要指定具有明确物理意义的网格面积（例如本案例的 $0.5\text{ km}^2$）。在 sf::st_make_grid(square = FALSE) 函数中，生成六边形网格所用的 cellsize 参数，指的是正六边形对边之间的距离（即内切圆直径 $d$）。

根据正六边形的面积公式，面积 $A$ 与对边距离 $d$ 的关系为： \[A = \frac{\sqrt{3}}{2} d^2\]

假设我们需要生成一个面积精确为 0.5 平方千米（即 $500,000\text{ m}^2$）的六边形网格，通过逆向推导： \[500,000 = \frac{\sqrt{3}}{2} d^2\] \[d^2 = \frac{1,000,000}{\sqrt{3}} \approx 577,350.27\] \[d = \sqrt{577,350.27} \approx 759.835\text{ (m)}\]

因此，在代码中设定 cellsize = 759.83，即可精确生成面积约为 0.5 $\text{km}^2$ 的空间统计单元。

# ============================================================
# 系统环境配置与核心依赖库加载
# ============================================================
library(tidyverse)
library(sf)
library(ggspatial)
library(spData)   # 提供大伦敦行政区划 (lnd) 边界数据集
library(classInt) # 执行 Fisher-Jenks 自然断点分类

# ============================================================
# 1. 远程获取与解析全量微观空间数据 (Remote Data Acquisition)
# ============================================================
# 指定由 Inside Airbnb 提供的 GZIP 压缩格式原始数据源
airbnb_gz_url <- "https://data.insideairbnb.com/united-kingdom/england/london/2025-09-14/data/listings.csv.gz"

# 依托 readr 引擎执行内存级解压与读取，直接将网络流解析为数据框 (DataFrame)
# 提示：由于样本量及属性维度较大，网络 I/O 耗时取决于当前连接带宽
london_airbnb_full <- read_csv(airbnb_gz_url)

# ============================================================
# 2. 空间对象实例化与数据清洗 (Spatial Object Instantiation)
# ============================================================
airbnb_sf <- london_airbnb_full %>%
  # 剔除坐标字段缺失的无效观测值，保障底层几何转换的稳健性
  filter(!is.na(longitude) & !is.na(latitude)) %>%
  # 将二维数值坐标映射为标准的几何点要素 (POINT)，并显式声明 WGS 84 (EPSG:4326) 坐标系
  st_as_sf(coords = c("longitude", "latitude"), crs = 4326)

# 检验实例化结果：确认 Simple feature collection 结构及其包含的属性字段
print(airbnb_sf)

# ============================================================
# 3. 载入基底要素与地理坐标系对齐 (Geographic CRS Alignment)
# ============================================================
# 载入基底多边形要素 (Polygon)：大伦敦各区边界 (默认使用 WGS 84 经纬度系统)
data("lnd", package = "spData")

# [空间拓扑原则] 强制投影对齐
# 确保点数据与基底面数据的坐标系完全一致，避免后续可视化产生空间漂移
airbnb_sf_transformed <- airbnb_sf %>%
  st_transform(crs = st_crs(lnd))

# ============================================================
# 4. 预设高对比度分类色谱 (Categorical Color Palette)
# ============================================================
# 针对高密度点群特征，构建基于暗场 (Dark Mode) 环境发光特性的定性色盘
# 旨在通过高饱和度冷暖对比，剥离不同房源类型在微观尺度下的空间异质性
nebula_colors <- c(
  "Entire home/apt" = "#00E5FF", # 荧光青 
  "Private room"    = "#9D00FF", # 深空紫 
  "Shared room"     = "#FF1493", # 炽热粉 
  "Hotel room"      = "#FFD700"  # 星尘金 
)

# ============================================================
# 5. [出图 1] 绘制暗夜发光主题的综合点面图 (Point Distribution)
# ============================================================
ggplot() +
  
  # --- 5.1 基础空间框架渲染 ---
  # 设定深色行政边界作为基底，降低边界亮度 (grey30) 以抑制对主体数据的视觉干扰
  geom_sf(data = lnd, fill = "black", color = "grey30", linewidth = 0.5) +
  
  # --- 5.2 核心要素空间映射 ---
  # 采用极小点径 (size = 0.01) 与低透明度 (alpha = 0.1) 策略控制过度绘制 (Overplotting)
  # 依托高密度区域的色彩物理叠加，自然生成伪热力聚集的视觉特征
  geom_sf(data = airbnb_sf_transformed, 
          aes(color = room_type), 
          size = 0.01,    
          alpha = 0.1) +  
  
  # --- 5.3 空间测量参考配置 ---
  annotation_north_arrow(location = "tl", which_north = "true", 
                         pad_x = unit(.2, "in"), pad_y = unit(.2, "in"),
                         style = north_arrow_fancy_orienteering(text_col = "white", 
                                                                line_col = "white", 
                                                                fill = "white")) +
  annotation_scale(location = "bl", width_hint = 0.25, 
                   text_col = "white", line_col = "white", 
                   pad_x = unit(.2, "in"), pad_y = unit(.2, "in")) +
  
  # --- 5.4 视觉通道与图例配置 ---
  scale_color_manual(values = nebula_colors, name = "Room Type") +
  # 解耦图面与图例的尺寸映射，独立放大图例标示点以保障学术可读性
  guides(color = guide_legend(override.aes = list(size = 3, alpha = 1))) +
  
  # --- 5.5 空间范围裁切 (Bounding Box) ---
  coord_sf(xlim = st_bbox(lnd)[c("xmin", "xmax")], 
           ylim = st_bbox(lnd)[c("ymin", "ymax")],
           expand = TRUE) +
  
  # --- 5.6 绘图引擎主题重构 ---
  theme_void() +
  theme(
    plot.background  = element_rect(fill = "black", color = NA),
    panel.background = element_rect(fill = "black", color = NA),
    panel.grid.major = element_line(color = "grey15", linetype = "dotted", linewidth = 0.4),
    
    axis.text  = element_text(color = "grey40", size = 7), 
    axis.title = element_blank(), 
    axis.ticks = element_blank(), 
    
    legend.position   = "bottom",
    legend.background = element_rect(fill = "black", color = NA),
    legend.text       = element_text(color = "white"),
    legend.title      = element_text(color = "white", face = "bold"),
    
    plot.title    = element_text(color = "white", face = "bold", size = 16, hjust = 0.5, margin = margin(t = 15, b = 5)),
    plot.subtitle = element_text(color = "grey70", hjust = 0.5, margin = margin(b = 15)),
    plot.caption  = element_text(color = "grey50", hjust = 1, size = 8, margin = margin(t = 10)),
    plot.margin   = margin(10, 15, 10, 15) 
  ) +
  
  labs(
    title = "Airbnb Spatial Distribution in London",
    subtitle = paste0("Total Listings: ", scales::comma(nrow(airbnb_sf)), " | Data Updated: 14 Sep 2025"),
    caption = "Data Source: Inside Airbnb & spData"
  )

# ============================================================
# 6. 空间网格离散化处理 (Hexagonal Tessellation & Classification)
# ============================================================

# 6.1 投影平面坐标系转换 (Projected CRS)
# 必须将数据转入英国国家格网系统 (BNG, EPSG:27700)，以便使用物理单位(米)进行几何面积分割
lnd_bng <- st_transform(lnd, crs = 27700)
airbnb_sf_bng <- st_transform(airbnb_sf_transformed, crs = 27700)

# 6.2 生成等面积六边形网格
# 在 BNG 投影下生成 cellsize = 759.83m 的网格，确保单体面积约为 0.5 km²
hex_grid <- st_make_grid(lnd_bng, cellsize = 759.83, square = FALSE) %>% 
  st_sf(hex_id = 1:length(.)) 

# 6.3 执行点面拓扑关联 (Spatial Join Count)
# 计算每个六边形单元内包含的 Airbnb 房源总量
hex_grid$listing_count <- lengths(st_intersects(hex_grid, airbnb_sf_bng))

# 6.4 过滤低频统计噪声
# 采用均值过滤 (Mean filtering)，剔除零散区域，聚焦于高于平均分布密度的核心单元
hex_filtered <- hex_grid %>% 
  filter(listing_count >= mean(hex_grid$listing_count))

# 6.5 执行 Fisher-Jenks 自然断点分类
# 将连续的数值划分为 7 个离散层级，旨在最大化组间空间异质性
jenks_breaks <- classIntervals(hex_filtered$listing_count, n = 7, style = "fisher")$brks

# 6.6 坐标系逆向转换 (Visual Alignment)
# 重要：将生成的网格与底图转回 WGS 84，以纠正 BNG 投影在大尺度下的视觉偏转，对齐标准水平经纬网
hex_viz <- hex_filtered %>%
  mutate(count_class = cut(listing_count, breaks = jenks_breaks, include.lowest = TRUE, dig.lab = 10)) %>%
  st_transform(crs = 4326)

lnd_viz <- st_transform(lnd_bng, crs = 4326)

# ============================================================
# 7. [出图 2] 空间分布密度可视化 (Scientific Map: Hexagonal Binning)
# ============================================================
ggplot() +
  
  # --- 7.1 行政基底渲染 ---
  # 使用中性色调 (Neutral tones) 构建背景，减少视觉冗余
  geom_sf(data = lnd_viz, fill = "#F5F5F5", color = "#CCCCCC", linewidth = 0.5) +
  
  # --- 7.2 核心密度网格渲染 ---
  # 采用带白色描边的多边形呈现清晰的蜂窝拓扑结构
  geom_sf(data = hex_viz, 
          aes(fill = count_class), 
          color = "white",        
          linewidth = 0.1,        
          alpha = 0.9) +
  
  # --- 7.3 视觉通道配置 (Visual Channels) ---
  scale_fill_viridis_d(option = "magma", 
                       direction = -1,
                       name = "Number of Airbnb listing buildings \n(Natural Breaks)") + 
  
  # 维持底部横向布局，适配多分类图例的展示张力
  guides(fill = guide_legend(nrow = 1, 
                             title.position = "top", 
                             title.hjust = 0.5,
                             keywidth = unit(1.2, "cm"), 
                             keyheight = unit(0.3, "cm"))) +
  
  # --- 7.4 空间制图规范组件 ---
  annotation_north_arrow(location = "tl", which_north = "true", 
                         pad_x = unit(.2, "in"), pad_y = unit(.2, "in"),
                         style = north_arrow_nautical(text_col = "black", 
                                                      line_col = "black", 
                                                      fill = c("black", "white"))) +
  annotation_scale(location = "bl", width_hint = 0.25, 
                   style      = "ticks",
                   text_col   = "black", line_col = "black", 
                   pad_x      = unit(.2, "in"), pad_y = unit(.2, "in")) +
  
  # --- 7.5 视角约束与投影声明 ---
  # 使用 WGS 84 范围裁切，确保经纬网格线横平竖直
  coord_sf(xlim = st_bbox(lnd_viz)[c("xmin", "xmax")], 
           ylim = st_bbox(lnd_viz)[c("ymin", "ymax")],
           expand = TRUE) +
  
  # --- 7.6 学术主题定制 ---
  theme_void() +
  theme(
    plot.background  = element_rect(fill = "white", color = NA),
    panel.background = element_rect(fill = "white", color = NA),
    panel.grid.major = element_line(color = "grey90", linetype = "dotted", linewidth = 0.4),
    
    axis.text  = element_text(color = "grey50", size = 7), 
    legend.position   = "bottom",
    legend.background = element_rect(fill = "white", color = NA),
    legend.text       = element_text(color = "black", size = 8),
    legend.title      = element_text(color = "black", face = "bold", size = 9),
    
    plot.title    = element_text(color = "black", face = "bold", size = 16, hjust = 0.5, margin = margin(t = 15, b = 5)),
    plot.subtitle = element_text(color = "grey40", hjust = 0.5, margin = margin(b = 15)),
    plot.caption  = element_text(color = "grey60", hjust = 1, size = 8, margin = margin(t = 10)),
    plot.margin   = margin(10, 15, 10, 15) 
  ) +
  
  labs(
    title = "Airbnb Spatial Density in London",
    subtitle = paste0("Tessellation Scale: ~0.5 km² per Hexagon | Classed by Fisher-Jenks"),
    caption = "Data Source: Inside Airbnb & spData | Projection: WGS 84"
  )

以上代码部分由于需要在线下载数据，暂时使用 eval = F 暂停运行，最终出图效果如下两张地图

上述代码成果图1

上述代码成果图2

6.3 进阶代码挑战

进阶练习

本节不提供参考代码。任务目标是：获取城市空间数据（可沿用伦敦数据，但鼓励探索其他国际大都市或中国代表性城市，如北京），完成一组可用于研究报告或学术论文的“微观空间聚合与分面可视化”。重点训练 连续型变量映射 × 空间连接统计（Spatial Join） × 多维分面（Faceting） 的综合空间表达能力。

6.3.1 数据与研究问题

数据集：目标城市的开源空间数据（建议访问 Inside Airbnb 获取）

核心数据源与变量：

矢量面数据：目标城市的行政边界（可通过 neighbourhoods.geojson 读取）。
空间点数据：目标城市完整房源切片（listings.csv.gz）。
核心分析变量：
- 空间属性：longitude, latitude
- 评价属性：review_scores_value（性价比评分）或 review_scores_rating（综合评分）
- 分组属性：room_type（房源类型）或 host_is_superhost（是否为超赞房东）

研究问题：
- 在选定的大都市中，Airbnb 房客的综合满意度（评分）呈现何种微观空间分布特征？
- 当采用正方形网格（Square Tessellation）平铺城市空间时，不同属性组别的区域中位数评分分布格局是否存在显著的空间异质性？

6.3.2 可视化任务

任务 A｜微观空间点分布与离散化映射（原始点分布图）
- 数据获取与投影：读取目标城市的房源点集与行政边界，剔除评分字段存在缺失值（NA）的无效观测。需客观查明并指定该城市适用的局部投影坐标系（Projected CRS），完成精确投影对齐。
- 底图绘制：在基础行政边界上叠加房源点要素，放弃简单的连续映射，要求采用自然断点法（Natural Breaks / Fisher-Jenks）或分位数法（Quantile）将评分变量离散化后再进行分层设色。同时需妥善配置点径与透明度，以缓解高密度区域的过度绘制问题。

任务 B｜几何变换与空间聚合统计（正方形密铺与中位数运算）
- 网格划分：基于目标城市边界，生成物理边长为 500m（或近似面积 0.25 km²）的正方形空间网格。
- 空间频数计算与阈值清洗：为消除极值偏差（如某网格仅有 1 个房源且评分为满分），必须首先执行空间相交运算（st_intersects）计算每个网格内的房源总数（Count），并严格剔除房源数量少于特定阈值（如 5 个或 10 个）的微弱信号地块。
- 空间连接运算（Spatial Join）：对保留下来的有效网格，利用点面拓扑关联将网格属性挂载至房源点上，进而分组计算每个正方形单元内所包含房源的中位数评分（Median Score）。

任务 C｜高阶制图：分面网格特征图（Faceted Choropleth）
- 分组运算：在空间连接与聚合阶段，引入分类属性变量（如房源类型），分别计算不同组别在各有效网格内的评分中位数。
- 分面渲染：使用 facet_wrap() 将不同组别的正方形网格地图并排输出。
- 制图规范：
- 过滤无数据的网格（保持底图通透）。
- 采用适宜的连续型或离散型渐变色盘（如 Viridis 序列）。
- 图面需包含：严谨的标题、共享色带图例、比例尺、指北针及数据与投影说明。

6.3.3 技术要求

使用 R Markdown 完成，并 knit 输出为 HTML。
文档需包含：
1. 数据获取逻辑与清洗准则（需说明缺失值与离群值的处理依据）。
2. 作图代码（结构清晰，核心拓扑运算步骤需分段注释）。
3. 简短文字解读（2–4 句）：提炼该城市高分房源的核心空间聚集区，并对比不同分面组别间的空间分布差异。
代码注释要求：
- 必须说明选定局部投影坐标系（EPSG）的客观依据。
- 说明网格样本量过滤阈值的设定理由。

提示

正方形网格构建：在网格生成函数中指定相关参数（如 square = TRUE）即可生成正方形单元。
属性统计与空间关联范式：面对“求特定多边形内点要素中位数”的需求，单一的空间相交计数运算无法实现。建议遵循“空间关系传递 $\rightarrow$ 属性表分组聚合 $\rightarrow$ 几何对象挂载”的分析范式：首先利用空间连接操作将正方形网格的唯一标识符（ID）映射至每个微观房源点位上；随后在点数据集中，按网格 ID 与目标分组变量执行常规的分组汇总运算（Group-by Summary）；最终将计算得出的统计量以属性表连接（Attribute Join）的方式匹配回原始的正方形网格几何对象中。
投影陷阱预警：生成基于物理距离（米）的网格前，必须避开地理坐标系（如 WGS 84）。请查阅目标城市的平面投影（例如北京常采用 CGCS2000 投影系，具体视经度带而定）。

第四章：数据可视化 2: 空间制图

章节作者：刘贇喆 | 华东师范大学 社会发展学院

更新时间：2026-03-26

本章学习内容

学习导航

1.0 空间数据

1.0 空间数据导论

1.1 空间数据主要类型

1.2 空间数据的核心结构

用示例数据理解 geometry 列

关于示例数据 nc

逐行读 head(nc) 的输出

常见误解

1.3 矢量空间数据的几何类型

点（Point）：位置型对象

线（Line）：连接与路径型对象

面（Polygon）：边界与区域型对象

1.4 坐标参考系（Coordinate Reference System, CRS）

CRS 的两种常见类型：地理坐标系 vs 投影坐标系

为什么要统一 CRS？以及如何通过 sf 统一

常见 CRS 的例子

本节小结：课程定位与进一步学习

2.0 空间数据读写

2.0 空间数据读写

2.1 常见矢量数据格式

2.2 使用 sf 读写矢量空间数据

读写 .shp 文件

读写 .gpkg 文件

读写 .geojson 文件

从 csv 到 sf：read_csv() + st_as_sf() + 写出

格式说明：读入格式 ≠ 写出格式

2.3 按照地理键链接空间与非空间数据

地理键的常见形式

常用连接方式：left_join() 将指标并入空间边界

为何常用 left_join()？

3.0 空间数据处理

3.0 空间数据处理

3.1 常规属性数据处理

基于 geometry 的派生变量

3.2 基础空间处理入门

3.2.1 空间关系

3.2.2 空间连接

3.2.3 简单空间处理

裁剪（Clip）与相交（Intersection）

并集（Union）与溶解（Dissolve）

差集（Difference）与对称差（Symmetric Difference）

3.2.4 常用空间数据生成

缓冲区（Buffer）

重心（Centroid）

4.0 空间数据可视化

4.0 空间数据可视化

4.0 空间可视化指南

制图导航

制图基础

4.1 静态制图基础：sf × ggplot2::geom_sf()

基础制图框架

代码示例

综合：基于 OpenStreetMap 的多图层空间可视化

知识拓展：OpenStreetMap (OSM) 与 osmdata 包

分层设色

4.2 分层设色图（Choropleth map）

代码示例 1：spData::world 绘制人口密度分层设色图

代码示例 2：基于 tidycensus 的纽约市贫困率空间可视化

制图数据分级策略：为何选择 Fisher-Jenks 算法？

数据分类分色

4.3 数据分类分色

思考：为什么要把连续数据变成离散等级？

1. Fisher-Jenks 自然断点法 (style = "fisher")

2. 等量分位法 (style = "quantile")

3. 等间距法 (style = "equal")

延伸：更多数据分级与自定义策略”

底图绘制

4.4 底图：空间语境

交互式地图

4.5 交互式地图（概览）

5.0 制图规范与美化

5.0 制图规范与美化

5.0 规范与美化指南

规范导航

制图规范

章节作者：刘贇喆 | 华东师范大学社会发展学院

用示例数据理解 `geometry` 列

关于示例数据 `nc`

逐行读 `head(nc)` 的输出

为什么要统一 CRS？以及如何通过 `sf` 统一

2.2 使用 `sf` 读写矢量空间数据

读写 `.shp` 文件

读写 `.gpkg` 文件

读写 `.geojson` 文件

从 `csv` 到 `sf`：`read_csv()` + `st_as_sf()` + 写出

常用连接方式：`left_join()` 将指标并入空间边界

为何常用 `left_join()`？

基于 `geometry` 的派生变量

4.1 静态制图基础：`sf` × `ggplot2::geom_sf()`

代码示例 1：`spData::world` 绘制人口密度分层设色图

代码示例 2：基于 `tidycensus` 的纽约市贫困率空间可视化

1. Fisher-Jenks 自然断点法 (`style = "fisher"`)

2. 等量分位法 (`style = "quantile"`)

3. 等间距法 (`style = "equal"`)

数据说明：Inside Airbnb 与 `lnd` 面数据

几何拓展：如何推算六边形的 `cellsize`？