第四章：数据可视化 2: 空间制图
章节作者：刘贇喆 | 华东师范大学 社会发展学院
更新时间：2026-10-13
章节概述
本章基于 sf 与 ggplot2
构建空间可视化基础，串联矢量数据结构、坐标参考系（CRS）与核心几何运算（如空间连接与裁剪），打通从数据清洗到制图输出的完整工作流。
实践层面，本章聚焦计算社会科学中高频的分层设色图（Choropleth Map）与多图层综合制图。重点解析连续数据离散化的算法策略（如自然断点法），并引入网格聚合等进阶技巧，以科学化解空间要素的视觉拥挤。
规范层面，本章归纳了地图辅助要素的排版美学，并着重划定中国公开地图出图的合规红线与标准画法，旨在输出兼具科学严谨、政治合规与发表友好的专业出版级图件。
本章学习内容
1.0 空间数据
1.0 空间数据导论
基本概念
空间数据（spatial
data）记录对象的位置、形状，或某个位置上的观测值。空间可视化（geo-visualisation）将这些信息组织成可阅读的地图。本节先概述空间数据的常见形式；下图以矢量数据为例，示意几何与属性如何关联。本章以
sf 矢量对象为主线，学习坐标参考系、空间处理与 R
中的基础制图；栅格数据仅作简要介绍。
sf 矢量数据的几何与属性（示意）
读图范围：上图只展示
sf矢量要素表；图中的地名、坐标与比例尺用于说明结构，不用于实际量测。栅格数据则按像元网格组织，本节仅作概念介绍。
1.1 空间数据主要类型
空间数据常见的组织方式有以下两种：
- 矢量数据（vector）：以点、线、面等几何对象来表达空间实体，适合描述离散对象（discrete
features）的形状、边界与相互位置关系；
- 栅格数据（raster）：以规则网格单元（像元）组织空间信息，可记录连续数值，也可记录土地利用类型等类别。
矢量用点、线、面表示对象；栅格以规则像元网格组织信息。
本章重点是矢量数据。 行政区、道路、学校和兴趣点（POI）等研究对象常以点、线、面表示，也便于连接人口普查或调查数据。下文主要用
sf 处理矢量对象，并用 ggplot2::geom_sf()
绘图。
选学 · 栅格的像元网格如何存储信息？
栅格按行、列排列像元，每个像元覆盖一定的地面范围；在同一图层或波段中，每个像元记录一个数值或类别代码。像元大小称为空间分辨率。网格范围与坐标参考系共同确定像元在地图上的位置。
- 连续值：气温、降水量、夜间灯光强度等。
- 类别值：土地覆盖中的林地、农田、建设用地等，通常用代码及对应说明表示。
- 多波段：遥感影像可在同一位置记录多个波段的值。
栅格适合研究连续变化的现象，也能用于类别地图。它按行列定位像元并记录对应的值；sf
矢量表则是每行一个要素，在同一 geometry
列中保存各行的几何。
1.2 sf
矢量要素表的结构
基础 · 概念
下面讨论的是本章使用的 sf
矢量要素表。它保留普通表格的行列结构，并用一个几何列（geometry
column）记录对象的位置与形状。
- 属性（attributes）：
类似csv的普通列，用于描述对象“是什么”
（例如区域名称、人口、收入、病例数等）
- 几何（geometry）：
用于描述对象“在哪里/长什么样”
（例如点的位置、道路的线形、行政区的边界面）
在这里，一行对应一个要素，其属性与几何位于同一行。几何列常命名为
geometry，但名称也可以不同。sf 对象仍可进行
filter()、mutate()、select()
等表格操作。
用示例数据理解 geometry
列
基础 · 示例
下面使用 sf 包自带的示例数据 nc（North Carolina
counties）来直观看空间数据的表格结构。
library(sf)
# 1) 读取 sf 自带示例数据：North Carolina counties
# st_read() 读取矢量文件；安装 sf 只需在 Console 中进行一次
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（北卡罗来纳州）县界示例。每行对应一个县：名称、代码和统计指标等是属性，geometry
列保存县界的形状与坐标。下面通过打印结果辨认这些信息。
逐行读 head(nc)
的输出
head(nc)
的输出包含对象概览与前六行要素表。不同
sf 版本的打印文字可能略有差异。
1. 对象概览信息
先看输出顶部的四行信息：
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：
表示该数据采用 NAD27（North American Datum 1927） 地理坐标参考系；可用st_crs(nc)查看完整信息。CRS 的作用见 1.4。
2. 属性表
接下来的表格中，每行对应一个县，普通列给出属性，geometry
列保存边界。这里看到的是 sf 矢量要素表的结构。
geometry：
这是几何列，打印时通常会被简写为MULTIPOLYGON (((...)))。
此处无需逐个解读坐标值；geometry保存边界坐标及其组织方式，可用于绘图与空间计算。
区分属性值与几何计算
AREA、PERIMETER
是数据原有的属性值，未必是当前几何与 CRS
下重新计算的面积、周长。地图轮廓由 geometry
决定；需要新的面积或长度时，应核对 CRS 与计算方法后再计算。
1.3 矢量空间数据的几何类型
基础 · 概念
在矢量空间数据中，几何信息通常以 三类
基本对象表示：点、线、面。
它们也可以从 “几何维度” 来理解：
点是 0 维（仅表示位置），线是 1 维（沿某一方向延展的路径，具有长度），面是 2 维（具有封闭边界的区域，具有边界与面积）。
点、线、面在矢量数据中分别由以下坐标结构表示：
- 点：通常由一对坐标表示，如
(x, y)（如经纬度lon, lat）；
- 线：由一串按顺序连接的点构成，即
P1 → P2 → ... → Pn（最简单情形是两点连线）；
- 面：由闭合的边界点序列构成，即起点与终点相同，形成封闭边界
注：在复杂情形下还可能包含“空洞”或多个分离部分。
地图通常分别用点符号、线条和填充区域显示这三种几何。以下示例先用固定样式辨认几何类型，尚未将属性值映射为点的大小、线的宽度或面的颜色。
表示方式取决于研究尺度。 同一座学校在城市地图中可表示为点；研究校园用地时，则可表示为面。几何类型服务于研究问题，并非对象永远不变的属性。
以下三张图用浅色经纬网辅助定位：横轴表示经度，纵轴表示纬度，图框统一按 WGS 84（EPSG:4326）显示。塞纳河的源数据采用另一种 CRS，绘图时仅转换显示；CRS 的含义见 1.4。
点（Point）：位置型对象
点（point）表示在当前地图尺度上可视作一个位置的对象，例如 GPS 定位点、医院、地铁站和调查样本点。
下面用 spData 中的 cycle_hire_osm（伦敦共享单车站点）展示点要素及其基础地图。
library(sf)
library(spData)
library(ggplot2)
# 点数据：伦敦共享单车站点
data("cycle_hire_osm", package = "spData")
ggplot(cycle_hire_osm) +
geom_sf(colour = "#356F8A", size = 0.7, alpha = 0.8) +
labs(x = "经度", y = "纬度") +
theme_bw(base_size = 10)
伦敦共享单车站点的点要素；一个点表示一个站点，图框显示 WGS 84 经纬网（EPSG:4326）。
线（Line）：连接与路径型对象
线（line / linestring）表示连接关系或路径，常见于道路、河流、公交线路和迁移路径。线的形状由按顺序连接的坐标点决定。
下面用 spData 中的 seine
展示法国塞纳河及其支流（Marne、Yonne）的线要素与基础地图。
library(sf)
library(spData)
library(ggplot2)
# 线数据：塞纳河及其支流
data("seine", package = "spData")
ggplot(seine) +
geom_sf(colour = "#356F8A", linewidth = 0.7) +
# 仅改变地图上的显示坐标系，不修改 seine 数据。
coord_sf(crs = st_crs(4326), datum = st_crs(4326)) +
labs(x = "经度", y = "纬度") +
theme_bw(base_size = 10)
塞纳河及支流的线要素。源数据采用 EPSG:2154；图框按 WGS 84（EPSG:4326）显示经纬网。
面（Polygon）：边界与区域型对象
面（polygon / multipolygon）表示具有边界的区域，例如行政区、社区、土地利用地块和建筑轮廓。一个面要素也可能包含孔洞或多个分离部分。
下面用 spData 中的 world
绘制国家与地区边界，观察面要素如何呈现在地图上。
library(sf)
library(spData)
library(ggplot2)
# 面数据：国家与地区边界
data("world", package = "spData")
ggplot(world) +
geom_sf(fill = "#DDE7EC", colour = "#59768A", linewidth = 0.18) +
labs(x = "经度", y = "纬度") +
theme_bw(base_size = 10)
世界国家与地区边界的面要素及 WGS 84 经纬网（EPSG:4326）；经纬度直绘不保持面积。
1.4 坐标参考系（CRS）
坐标参考系（Coordinate Reference System, CRS）说明坐标如何对应到地球上的位置。经纬度或平面坐标若缺少 CRS，数值本身不足以确认地点和单位。
地图可叠放多个图层（layer），例如分别管理的站点、道路与行政区数据。当图层来自不同文件时，先用
st_crs() 检查各自的
CRS。绘图工具可以在显示时转换图层；开展空间连接、缓冲区、距离或面积计算前，则应明确输入
CRS、计算方法及所需单位，必要时将数据转换到适合研究区域的 CRS。
CRS 的两种常见类型：地理坐标系 vs 投影坐标系
入门阶段先区分两类常见 CRS：
地理坐标系（Geographic CRS / GCS）：
以经度与纬度来表达位置，数值单位通常是度（可用十进制度或度分秒表示）。这种坐标便于全球定位，但在距离、面积等度量计算上需要谨慎，因为“度”并不是线性距离单位。投影坐标系（Projected CRS / PCS）：
通过地图投影将地球表面转换为平面坐标，单位可为米或英尺等。投影都会产生变形；距离或面积是否适合计算，仍取决于投影、研究区域和测量目的。
地理坐标经投影转换为平面坐标；投影会产生变形（示意）。
图左的经度 λ、纬度 φ 标示位置 P；经过投影转换，同一位置在图右以平面坐标 X、Y 表示。中间的投影过程提醒：地图上的形状、面积、距离或方向可能发生变化。右侧的米制坐标与三类平面运算只作示意，实际量测须根据研究区域和目的选择合适的 CRS。
不同地图投影会改变经纬网（graticule）在平面上的形状。下面以同一份世界底图、相同的经纬网间距进行对照：
- WGS 84
经纬度直绘：经度作为横坐标、纬度作为纵坐标，经纬网呈水平和竖直线；
- Robinson 投影：经纬网弯曲、向边缘收敛，适合展示世界地图的整体形状，但不保持面积或距离。
两幅图均使用经线每 60°、纬线每 30° 的间距；形状差异来自坐标转换。
同一世界底图与相同经纬网间距下，WGS 84 经纬度直绘（上）与 Robinson 投影（下）的形态对比。数据：ggmapcn 世界底图。
查看并复制投影对比图的 R 代码
library(sf)
library(ggplot2)
library(ggmapcn)
library(patchwork)
# 1. 固定底图与网格间距，仅改变输出 CRS
crs_robin <- "+proj=robin +lon_0=0 +datum=WGS84 +units=m +no_defs"
p_wgs84 <- ggplot() +
geom_world(ocean_fill = "#EAF4FA", country_fill = "#E5EAF0",
coastline_color = "#5D86A1", international_boundary_color = "#526376") +
annotation_graticule(lon_step = 60, lat_step = 30,
line_color = "#A9BDCB", label_color = "#526779",
line_width = 0.25, label_size = 3.2, label_offset = 5) +
coord_sf(crs = st_crs(4326)) +
labs(title = "WGS 84（EPSG:4326）") +
theme_void() +
theme(plot.title = element_text(face = "bold", colour = "#234058", size = 13,
margin = margin(b = 8)),
plot.margin = margin(8, 16, 8, 16), panel.ontop = TRUE)
# 2. 以同样的视觉设置绘制 Robinson 投影
p_robin <- ggplot() +
geom_world(crs = crs_robin, ocean_fill = "#EAF4FA", country_fill = "#E5EAF0",
coastline_color = "#5D86A1", international_boundary_color = "#526376") +
annotation_graticule(crs = crs_robin, lon_step = 60, lat_step = 30,
line_color = "#A9BDCB", label_color = "#526779",
line_width = 0.25, label_size = 3.2, label_offset = 3e5) +
coord_sf(crs = crs_robin) +
labs(title = "Robinson 投影") +
theme_void() +
theme(plot.title = element_text(face = "bold", colour = "#234058", size = 13,
margin = margin(b = 8)),
plot.margin = margin(8, 16, 8, 16), panel.ontop = TRUE)
p_wgs84 / p_robin额外阅读 · 世界地图为何讨论等面积投影？
2026 年 9 月 4 日，联合国大会通过第 A/RES/80/307 号决议。决议鼓励在比较各洲相对面积的地图中使用等面积投影，包括 Equal Earth；同时承认没有一种投影能满足所有用途。这是关于地图如何呈现面积的倡议，并非要求所有地图采用同一投影，更不是更改空间数据的通用 CRS。上图的 Robinson 是用于世界概览的折中投影，不保持面积。联合国会议新闻稿提供了通过情况。
检查、声明与转换 CRS
st_crs() 用于查看
CRS；st_set_crs() 为已经知道原
CRS、但文件缺少记录的坐标补上声明，不会改变坐标数值；st_transform()
才会把坐标计算到目标 CRS。原 CRS
不明时应先查数据说明，不能猜一个 EPSG 编号直接赋给数据。
下面的 world 与 seine 分别采用 EPSG:4326 和
EPSG:2154。将河流数据转换到世界底图的 CRS 后，再核对编号：
library(sf)
library(spData)
data("world", package = "spData")
data("seine", package = "spData")
# 1. 查看原始 CRS
c(world = st_crs(world)$epsg, seine = st_crs(seine)$epsg)world seine
4326 2154 # 2. 真正转换坐标，并再次检查
seine_to_world <- st_transform(seine, st_crs(world))
c(world = st_crs(world)$epsg, seine_to_world = st_crs(seine_to_world)$epsg) world seine_to_world
4326 4326
ggplot2::geom_sf()可在绘图时自动转换图层以便共同显示。空间运算仍需检查 CRS 和方法；经纬度数据也可借助sf的测地计算处理距离、面积，不能仅凭“单位是米”判断结果准确。
常见 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 代码）。
查询其他地区的 CRS 时，优先核对数据提供方的元数据与 EPSG Registry；相同地区也可能因研究任务而选择不同投影。
本节回顾
先辨认数据结构，再检查 CRS。 矢量 sf
用几何列连接要素与属性；栅格用带空间参照的像元网格记录值。点、线、面的选择与研究尺度有关。地图叠加和空间量测前，确认源
CRS、目标 CRS 与所用方法是否适合研究区域。
2.0 空间数据读写
2.0 空间数据读写
上一节已经明确：矢量空间数据本质上是“属性表 + 几何列（
geometry）”。
本节聚焦于建立一个基于 sf 与 tidyverse
的空间数据读写工作流，目标是将不同来源的数据整理为可直接用于制图的
sf
输入，并支持后续复用与共享。该工作流主要包含以下四类操作：
读入：用
sf::st_read()读取常见矢量文件（如 Shapefile、GeoJSON、GeoPackage），得到sf对象；转换：将表格型数据（如
.csv中的x/y或lon/lat坐标列）转换为空间点数据（sf::st_as_sf()），并显式设置 CRS；连接：将外部表格数据与空间数据按键连接（attribute join），形成带有制图指标的“空间属性表”；
写出：将处理结果导出为常用空间格式（如
gpkg/geojson），以便复用、共享与后续制图。
2.1 常见空间数据存储方式
基础 · 概念
空间信息可以保存在矢量文件、带坐标的普通表格或空间数据库中。下图并列展示六种常见的存储方式；其中
CSV 坐标表需要补充已知 CRS 才能转为 sf 点对象，PostGIS
则是数据库扩展，并非文件格式。
- Shapefile（
.shp）：由同名的.shp、.shx、.dbf等文件共同组成；.prj通常记录 CRS。复制或分享时要保留整组文件。字段名长度和字段类型也有限制。
- GeoPackage（
.gpkg）：一个文件可容纳多个空间图层和属性表，适合组织项目数据。
- GeoJSON（
.geojson）：基于 JSON 的矢量数据格式，常用于网页展示和跨平台交换；普通.json文件不一定符合 GeoJSON 规范。
- CSV + 坐标列：普通表格可记录点的
lon/lat或x/y，但 CSV 本身不记录 CRS，也不直接保存sf的几何列。
- KML / GML：两种基于 XML 的空间数据格式；KML
常用于地图展示，GML 常见于标准化数据交换。KMZ 通常是压缩的 KML
文件。
- PostGIS：PostgreSQL 的空间数据库扩展，可存储空间表并支持空间查询；使用时需要数据库连接，而非单个文件路径。
六种常见空间数据存储方式
图内标题将不同载体一并概括：Shapefile、GeoPackage、GeoJSON、KML/GML 是矢量文件格式，CSV 是坐标表，PostGIS 是数据库扩展。PostGIS 中的数据版本仍需配合额外的记录与管理流程。窄屏阅读时可点击图片放大查看。
2.2 使用 sf
读写矢量空间数据
基础 · 操作
st_read() 读取矢量文件，st_write() 将
sf 对象写出。两者借助 GDAL 驱动支持
Shapefile、GeoPackage、GeoJSON 等格式。本节先读入 R
包附带的示例文件，再写到临时目录；重新读入后，用简短的结果表核对要素数和
CRS。
代码中的 system.file(..., package = ...) 用来定位已安装
R 包里的示例文件；tempdir() 返回本次 R
会话的临时目录。临时目录中的输出会在会话结束后被清理，实际项目应改用明确的保存路径。
读写 .shp 文件
常用
sf 包附带的 nc
是北卡罗来纳州县级边界。system.file() 找到其中的
.shp 文件，st_read()
会一并读取同名配套文件中的属性与
CRS。此处只覆盖专门创建的临时输出，原始文件应保留。
library(sf)
path_nc <- system.file("shape/nc.shp", package = "sf")
nc <- st_read(path_nc, quiet = TRUE)
# 写出后重读；Shapefile 会生成一组同名文件
out_shp <- file.path(tempdir(), "nc_out.shp")
st_write(nc, out_shp, delete_layer = TRUE, quiet = TRUE)
nc_out <- st_read(out_shp, quiet = TRUE)
data.frame(input_rows = nrow(nc), output_rows = nrow(nc_out),
same_crs = st_crs(nc_out) == st_crs(nc),
same_fields = identical(names(nc_out), names(nc))) input_rows output_rows same_crs same_fields
1 100 100 TRUE TRUE核对表中的 input_rows 与 output_rows
应相等；same_crs 和 same_fields 应为
TRUE。后两个示例也按相同方式查看写出结果。
读写 .gpkg 文件
GeoPackage
可以在同一文件中保存多个图层。一个矢量图层可理解为文件中的一张空间要素表：先用
st_layers() 查看图层名称，再以 layer =
读入所需图层。重复运行写出代码时，delete_layer = TRUE
只替换指定的输出图层，保留文件中的其他图层。
查看 GeoPackage 图层读写代码
library(sf)
library(spData)
path_gpkg <- system.file("shapes/world.gpkg", package = "spData")
st_layers(path_gpkg) # 先查看图层名称Driver: GPKG
Available layers:
layer_name geometry_type features fields crs_name
1 world Multi Polygon 177 10 WGS 84world_gpkg <- st_read(path_gpkg, layer = "world", quiet = TRUE)
out_gpkg <- file.path(tempdir(), "world_out.gpkg")
st_write(world_gpkg, out_gpkg, layer = "world_out",
delete_layer = TRUE, quiet = TRUE)
st_layers(out_gpkg)Driver: GPKG
Available layers:
layer_name geometry_type features fields crs_name
1 world_out Multi Polygon 177 10 WGS 84world_out <- st_read(out_gpkg, layer = "world_out", quiet = TRUE)
data.frame(input_rows = nrow(world_gpkg), output_rows = nrow(world_out),
same_crs = st_crs(world_out) == st_crs(world_gpkg),
same_fields = identical(names(world_out), names(world_gpkg))) input_rows output_rows same_crs same_fields
1 177 177 TRUE TRUE读写 .geojson 文件
重点 · 常用
这里读取 spData 附带的伦敦共享单车站点
GeoJSON，再写出为另一个 GeoJSON
文件。覆盖操作只针对专门创建的临时文件。
查看 GeoJSON 文件读写代码
library(sf)
library(spData)
path_geojson <- system.file("shapes/cycle_hire_osm.geojson", package = "spData")
cycle_hire_osm_gj <- st_read(path_geojson, quiet = TRUE)
out_geojson <- file.path(tempdir(), "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)
data.frame(input_rows = nrow(cycle_hire_osm_gj),
output_rows = nrow(cycle_out),
same_crs = st_crs(cycle_out) == st_crs(cycle_hire_osm_gj)) input_rows output_rows same_crs
1 532 532 TRUE从 csv 到
sf：read_csv() + st_as_sf() +
写出
重点 · 常用
点位 CSV
通常以每行一条记录、两列保存坐标。转换前要确认坐标列是数值、没有缺失，并从数据说明或原始对象核对源
CRS；列名 lon/lat 和数值范围只能辅助检查，不能证明
CRS。
为使示例无需额外下载，下面先核对 spData::cycle_hire_osm
的源 CRS 为 WGS 84（EPSG:4326），再从中生成一份真实的临时 CSV
文件并用 read_csv() 读回。CSV 本身不保存
CRS；st_as_sf() 使用前述已核对的 EPSG:4326
声明坐标，coords = c("lon", "lat") 表示经度（横坐标）在前、纬度（纵坐标）在后。st_transform()
随后才将几何坐标转换到英国国家格网（EPSG:27700）。最后写出 GeoPackage
并重新读入检查。
osm_id 是标识符，读入 CSV
时指定为字符型。remove = FALSE 会保留原始
lon/lat 属性列；投影转换后，它们仍是 WGS 84
经纬度，不能当作 EPSG:27700 的平面坐标。本例用
stopifnot() 检查源 CRS
和坐标列：条件不成立时停止运行，避免把错误坐标生成点位。实际数据还应核对经纬度范围，并抽查点位是否落在预期地区。
前例读取的是包内 GeoJSON 文件，这里使用的是包内同名 R 对象。两段代码分别用各自的输入核对写出结果，不以两者的要素数互相校验。
library(sf)
library(spData)
library(tidyverse)
# 1. 源 CRS 必须是 EPSG:4326；否则不能按 WGS 84 声明坐标
data("cycle_hire_osm", package = "spData")
stopifnot(st_crs(cycle_hire_osm)$epsg == 4326L)
coords <- st_coordinates(cycle_hire_osm)
cycle_csv_source <- tibble(
osm_id = as.character(cycle_hire_osm$osm_id),
name = as.character(cycle_hire_osm$name),
lon = coords[, 1], lat = coords[, 2]
)
path_csv <- file.path(tempdir(), "cycle_hire_osm.csv")
write_csv(cycle_csv_source, path_csv)
cycle_tbl <- read_csv(
path_csv,
col_types = cols(osm_id = col_character(), .default = col_guess()),
show_col_types = FALSE
)
# 2. 坐标须为有效数值；再声明源 CRS 并转换几何坐标
stopifnot(is.numeric(cycle_tbl$lon), is.numeric(cycle_tbl$lat),
all(is.finite(cycle_tbl$lon)), all(is.finite(cycle_tbl$lat)))
cycle_sf_wgs84 <- cycle_tbl %>%
st_as_sf(coords = c("lon", "lat"), crs = 4326, remove = FALSE)
cycle_sf_bng <- st_transform(cycle_sf_wgs84, 27700)
# 3. 写出 GeoPackage，重读并查看行数与 CRS
out_cycle_gpkg <- file.path(tempdir(), "cycle_hire_bng.gpkg")
st_write(cycle_sf_bng, out_cycle_gpkg, layer = "cycle_hire_bng",
delete_layer = TRUE, quiet = TRUE)
cycle_bng_out <- st_read(out_cycle_gpkg, layer = "cycle_hire_bng", quiet = TRUE)
tibble(csv_rows = nrow(cycle_tbl), source_epsg = st_crs(cycle_sf_wgs84)$epsg,
gpkg_rows = nrow(cycle_bng_out), gpkg_epsg = st_crs(cycle_bng_out)$epsg)# A tibble: 1 × 4
csv_rows source_epsg gpkg_rows gpkg_epsg
<int> <int> <int> <int>
1 540 4326 540 27700格式说明：读入格式 ≠ 写出格式
前三个示例展示同格式读写，CSV 示例则展示从普通表格生成 GeoPackage。更换文件格式与变换 CRS 是两项操作；需要目标 CRS 时应显式转换，并在读回后核对。这里的要素数、CRS 和字段检查只是基本检查，重要数据还应核对属性值和几何；目标格式对字段名及类型的支持可能不同。
- 如果只需将点数据保存为 CSV，可先用
st_drop_geometry()去掉几何列，并保留lon/lat等坐标列；再次转回sf时还需要另行记录的源 CRS。
- 对线、面等复杂几何，也可将几何保存为 WKT 文本列，并另行记录 CRS。普通 CSV 不会自动保留完整的空间对象结构；需要可靠保留几何和 CRS 时优先使用 GeoPackage。
实践建议 GeoPackage 适合项目内复用，也可由 QGIS、ArcGIS 等 GIS 软件读取；GeoJSON 常用于网页与跨平台交换。只有接收方或工作规范要求 Shapefile 时，才需额外导出整组
.shp配套文件。
2.3 按照地理键链接空间与非空间数据
重点 · 常用
空间边界数据（boundary
data）通常包含几何列和一些基础属性；研究需要的人口、经济等指标可能保存在另一张统计表中。依据两表共有的地理键（geographic
key），可以把指标并入边界对象，供后续制图和分析。该操作与第二章的表格连接相同，通常使用
dplyr::left_join()。
下图用美国州代码与人口密度说明连接过程。后面的代码沿用同一空间层级，改用美国人口调查局的 2025 年州级边界和 2025 年 7 月 1 日人口估计完成实际连接；图中的数字仅用于流程示意。
按州代码连接边界与人口密度（流程示意）
图中 geo_id
对应州代码，右侧示意连接后的分层设色图；图内数值并非下方代码的实际结果。窄屏阅读时可点击图片放大查看。
地理键的常见形式
基础 · 概念
地理键（geographic key / identifier）用于标识空间单元。若边界表每行代表一个地理单元，键应在边界表中唯一；要维持连接后每个边界一行，指标表中的同一键也应至多出现一次。
连接前还应确认两表对应相同的地理层级，并核对统计时期、行政区划或边界版本及键的字段类型。有前导零的编码应按字符读取，避免读成数值后丢失前导零。
常见地理键形式主要包括两类：
地理编号/代码（code）：
以标准化编码唯一标识空间单元，通常更稳定、更适合作为连接键。
例如：英国统计地理体系中的OA / LSOA / MSOA代码、美国州/县/普查区的GEOID，以及各类行政区划代码。邮编（postcode）属于邮政地理体系；与统计区边界连接时，应先核对空间单元是否一致，必要时使用对应关系表。地名字段（name）：
以名称描述空间单元，例如“上海市·浦东新区”“某某区/县”“某某州”等。
名称字段更直观，但在连接时更容易受到同名、别名、拼写差异、大小写、空格与符号差异等因素影响。
建议 若同时提供
code与name，通常优先使用经过核对的code作为连接键；name更适合作为地图标注或图例展示字段。
常用连接方式：left_join()
将指标并入空间边界
常用
sf 对象的属性部分遵循数据框结构，因此可直接使用
dplyr::left_join()
把外部指标并入边界对象。该函数比较的是两表的键值，不会检查几何位置是否重合。右表键唯一时，每个左表边界仍只占一行，几何列随该行保留；若右表同一键出现多次，连接后可能复制边界行，制图前必须检查。
基本写法
地理键同名（如，地理键所在的那一列都叫做
'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()
会保留全部边界；未匹配单元的指标为
NA，地图可用缺失色标识。应先查看未匹配项的原因，不能直接把
NA 当作零。
相对地，inner_join()
仅保留两表都能匹配到的空间单元：未匹配的边界将被剔除，地图上对应区域也会“消失”。该方式适用于明确需要排除缺失单元的场景，应结合研究目的选择。
示范案例
这里使用两份独立的近期公开数据：2025
年州级制图边界和 Vintage
2025 州人口估计，均来自美国人口调查局。POPESTIMATE2025
是 2025 年 7 月 1
日的人口估计，并非人口普查计数。课程项目附带两个原始文件；单独复制代码而没有本地文件时，代码会下载相同版本的官方文件。
边界的 GEOID 与人口表的 STATE
都是两位州级代码，须按字符读取以保留前导零。人口 CSV
还包含全国、区域和分区汇总行，因此先筛选州级的
SUMLEV == "040"。anti_join()
分别找出两侧未匹配的代码；relationship = "one-to-one"
在同一代码对应多行时会报错，避免连接后重复边界。
连接后用人口估计除以边界文件的 ALAND（陆地面积，平方米），得到每平方公里陆地面积的估计人口。此处边界经过制图概化，计算密度应使用文件提供的统计面积字段，不从概化后的地图形状反算面积。
先运行数据读取代码，再运行下方的键值检查与连接代码。读取步骤默认折叠，便于聚焦本节的连接逻辑。
查看 2025 年官方数据读取代码
library(sf)
library(tidyverse)
# 1. 读取教材快照；单独运行代码时下载同版本的官方原始文件
local_dir <- file.path("data", "census-2025")
state_zip <- file.path(local_dir, "cb_2025_us_state_20m.zip")
pop_csv <- file.path(local_dir, "NST-EST2025-ALLDATA.csv")
if (!file.exists(state_zip)) {
state_zip <- file.path(tempdir(), basename(state_zip))
if (!file.exists(state_zip)) download.file(
"https://www2.census.gov/geo/tiger/GENZ2025/shp/cb_2025_us_state_20m.zip",
state_zip, mode = "wb", quiet = TRUE
)
}
if (!file.exists(pop_csv)) {
pop_csv <- file.path(tempdir(), basename(pop_csv))
if (!file.exists(pop_csv)) download.file(
"https://www2.census.gov/programs-surveys/popest/datasets/2020-2025/state/totals/NST-EST2025-ALLDATA.csv",
pop_csv, mode = "wb", quiet = TRUE
)
}
# 2. Shapefile 是一组文件：解压后读入；人口表只保留州级记录
shape_dir <- file.path(tempdir(), "census_states_2025")
dir.create(shape_dir, showWarnings = FALSE)
unzip(state_zip, exdir = shape_dir)
states_2025 <- st_read(
file.path(shape_dir, "cb_2025_us_state_20m.shp"), quiet = TRUE
)
boundary_sf <- states_2025 %>%
select(geo_id = GEOID, state_name = NAME, land_m2 = ALAND)
# 州代码按字符读取，保留 "01" 等前导零
attr_tbl <- read_csv(
pop_csv,
col_types = cols(SUMLEV = col_character(), STATE = col_character(),
POPESTIMATE2025 = col_double(), .default = col_skip()),
show_col_types = FALSE
) %>%
filter(SUMLEV == "040") %>%
transmute(geo_id = STATE, population_2025 = POPESTIMATE2025)# 3. 分别找出未匹配的代码，再做一对一属性连接
boundary_keys <- boundary_sf %>% st_drop_geometry() %>% select(geo_id)
unmatched_boundary <- anti_join(boundary_keys, attr_tbl, by = "geo_id")
unused_indicator <- anti_join(attr_tbl, boundary_keys, by = "geo_id")
states_joined <- boundary_sf %>%
left_join(attr_tbl, by = "geo_id", relationship = "one-to-one") %>%
mutate(density_2025 = population_2025 / (land_m2 / 1e6))
# 结果表中未匹配、未使用和缺失数量均应为 0
tibble(
check = c("边界", "人口表", "连接后", "边界未匹配", "人口表未使用", "人口值缺失"),
count = c(nrow(boundary_sf), nrow(attr_tbl), nrow(states_joined),
nrow(unmatched_boundary), nrow(unused_indicator),
sum(is.na(states_joined$population_2025)))
)# A tibble: 6 × 2
check count
<chr> <int>
1 边界 52
2 人口表 52
3 连接后 52
4 边界未匹配 0
5 人口表未使用 0
6 人口值缺失 0这份固定快照应得到 52 个州级单元（50 州、华盛顿特区和波多黎各），连接后仍为 52
行，未匹配与人口值缺失均为
0。若使用其他数据时检查表出现非零值，应先核对地理键和统计范围。新变量
density_2025
是基于人口估计与陆地面积计算的描述性指标；解读地图时仍需注明人口估计日期与面积口径。
3.0 空间数据处理
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() 适用于多条件分支（“分组规则表”写法）。常用
示例代码
沿用 2.3 节连接得到的 states_joined，以美国 2025
年州级边界和 2025 年 7 月 1
日人口估计演示属性处理。density_2025
已按人口估计除以 ALAND 陆地面积计算；这里先用
case_when() 区分州、华盛顿特区和波多黎各，再筛选 50
个州，分别查看人口密度最高与最低的州。单独运行本例时，须先运行 2.3
节的数据读取与连接代码。
library(sf)
library(tidyverse)
# sf 的 select() 会保留几何列，后续仍可制图
state_profile <- states_joined %>%
select(geo_id, state_name, population_2025, density_2025) %>%
mutate(
unit_type = case_when(
geo_id == "11" ~ "华盛顿特区",
geo_id == "72" ~ "波多黎各",
TRUE ~ "州"
),
population_million = population_2025 / 1e6
)
# 排序前只取 50 个州；去掉几何列使结果表更易阅读
density_table <- state_profile %>%
filter(unit_type == "州") %>%
st_drop_geometry() %>%
select(state_name, population_million, density_2025)
density_table %>% arrange(desc(density_2025)) %>% slice_head(n = 6) state_name population_million density_2025
1 New Jersey 9.548215 501.2333
2 Rhode Island 1.114521 416.2123
3 Massachusetts 7.154084 354.0841
4 Connecticut 3.688496 294.0892
5 Maryland 6.265347 249.1063
6 Delaware 1.059952 210.0284 state_name population_million density_2025
1 Alaska 0.737270 0.4981913
2 Wyoming 0.588753 2.3413554
3 Montana 1.144694 3.0365369
4 North Dakota 0.799358 4.4733242
5 South Dakota 0.935094 4.7625779
6 New Mexico 2.125498 6.7648248两张结果表分别显示人口密度最高和最低的六个州。筛选前的
state_profile 仍是包含 52 个州级单元的 sf
对象；st_drop_geometry()
只用于打印属性结果，不改变原数据。密度的分母是官方陆地面积字段，并非概化地图边界的几何面积。
基于 geometry
的派生变量
重点
sf 的关键特征在于：几何列 geometry
可直接参与计算；因此 mutate()
不仅能生成普通属性变量，也能生成基于几何度量的空间指标。常见操作包括：
- 面要素：
st_area()计算面积；
- 边界/线要素：
st_length()计算长度（对面数据通常先取边界再计算周长）；
- 派生指标：将属性量与几何度量结合形成密度类指标（“属性 ÷ 面积/长度”）。
注意
st_area()/st_length()的单位与计算方式受 CRS 影响。量测时须选择适用于研究地区和量测任务的投影；仅看到坐标单位为“米”，不能保证面积或长度准确。
示例代码
下面使用 sf 自带的北卡罗来纳州县界数据
nc，转换到适用于当地的 NAD83 / North
Carolina（EPSG:32119）
后计算面积和周长。该投影的坐标单位为米，计算结果再换算为平方千米和千米。
library(sf)
library(tidyverse)
path_nc <- system.file("shape/nc.shp", package = "sf")
nc <- st_read(path_nc, quiet = TRUE)
nc_geom <- nc %>%
select(NAME) %>%
st_transform(32119) %>%
mutate(
area_km2 = as.numeric(st_area(geometry)) / 1e6,
# 面的周长先取边界，再计算边界长度
perimeter_km = as.numeric(st_length(st_boundary(geometry))) / 1e3
)
nc_geom %>%
st_drop_geometry() %>%
select(NAME, area_km2, perimeter_km) %>%
head(6) NAME area_km2 perimeter_km
1 Ashe 1137.5901 141.6782
2 Alleghany 611.1970 119.9423
3 Surry 1423.7288 160.5110
4 Currituck 694.6611 301.5390
5 Northampton 1520.9918 211.9729
6 Hertford 967.8553 160.90393.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 是否“只接触边界”- 典型理解：两个面共享边界/顶点，但内部不重叠（没有共享面积）。
- 典型理解：两个面共享边界/顶点，但内部不重叠（没有共享面积）。
这些函数返回满足条件的对象配对，不会裁切几何；需要新的交叠范围时应使用后文的几何运算。
这类函数的输出本质上是在回答“哪些对象彼此满足某种空间关系”。它们通常被用来：
- 筛选：保留满足空间关系的要素；
- 匹配：为每个对象找到与之对应的对象；
- 连接：把匹配到的属性并入数据中（空间连接）。
包含、相交与边界接触的空间关系（示意）。
3.2.2 空间连接
重点 · 常用
空间连接（spatial
join）根据几何之间的关系，将 y 图层的属性并入
x 图层，输出保留 x 的几何。与 2.3
节按地理键连接不同，它以“位于内部”“相交”等位置关系匹配对象；连接本身不会裁切线段或计算交叠面积。
常见写法是 st_join(x, y, join = st_within) 或
st_join(x, y, join = st_intersects)。连接前须确认两个图层的
CRS 一致。同一个 x 要素若匹配多个 y
要素，结果中可能出现多行，因此汇总前应检查重复匹配。
按位置关系把区域属性连接到点要素（示意）。
典型案例
点 → 面：
将点位归属到所在的面单元（经典point-in-polygon问题）。
常用于将 POI/事件点“分配”到行政区或规则网格，再进一步计算计数、密度或比率并制图。 常用线 → 面：
可用空间连接识别一条线与哪些面相交；若需计算面内线段长度，应先做几何相交或裁剪，再在合适的投影下量测。面 ↔︎ 面：
可用空间连接识别边界间的对应关系；若需交叠面积或面积权重，应另做叠置运算。简单的“面积权重”并不自动等于人口等属性的真实分布。
示例：point in polygon问题
示例
下面构造一个矩形面作为“区域边界”，再随机生成一批点位。
目标是：用 st_join()
将点位归属到该面中，并进一步统计“面内点的数量”。
该流程与实践中“事件点/POI → 行政区/缓冲区”一致。
library(sf)
library(tidyverse)
# 1. 构造同一 CRS 下的矩形区域和 30 个示例点
region <- st_sf(
region_id = "A",
geometry = st_sfc(st_polygon(list(rbind(
c(0, 0), c(10, 0), c(10, 10), c(0, 10), c(0, 0)
))), crs = 4326)
)
set.seed(123)
points_sf <- tibble(
point_id = 1:30,
lon = runif(30, -2, 12),
lat = runif(30, -2, 12)
) %>%
st_as_sf(coords = c("lon", "lat"), crs = 4326)
# 2. 把区域编号连接到点；未匹配的点得到 NA
points_joined <- st_join(points_sf, region, join = st_within)
points_joined %>% st_drop_geometry() %>% count(region_id)# A tibble: 2 × 2
region_id n
<chr> <int>
1 A 15
2 <NA> 15# 3. 先按区域编号计数，再接回原来的面，保留其几何
point_counts <- points_joined %>%
st_drop_geometry() %>%
filter(!is.na(region_id)) %>%
count(region_id, name = "n_points")
region_count <- region %>%
left_join(point_counts, by = "region_id") %>%
mutate(n_points = replace_na(n_points, 0L))
region_count %>% st_drop_geometry() region_id n_points
1 A 15ggplot() +
geom_sf(data = region_count, fill = "#DCECEF", colour = "#537687") +
geom_sf(data = points_joined,
aes(colour = if_else(is.na(region_id), "面外", "面内")),
size = 1.9, alpha = 0.85) +
scale_colour_manual(values = c("面内" = "#197A85", "面外" = "#C37350")) +
labs(x = "经度（示意）", y = "纬度（示意）", colour = NULL) +
theme_bw(base_size = 10) +
theme(panel.grid.minor = element_blank(), legend.position = "right")
点与区域的空间连接：面内点获得区域编号，面外点未匹配。
这里的 region_count 仍是可制图的 sf
面对象。示例只有一个面，点不会重复匹配；若面彼此重叠，同一点可能出现在多个匹配结果中，不能直接把连接后的行数当作独立点数。
3.2.3 简单空间处理
几何运算会根据输入对象生成新的空间范围。下图保留原有的操作概览；随后用同一对示例面
A、B 运行 sf 代码，比较运算前后的结果。关系判断函数（如
st_intersects()）只返回匹配关系，st_intersection()
才生成交叠几何。
相交、并集与裁剪等空间处理（示意）。
相交与裁剪
st_intersection(x, y) 保留 x 与 y 的共同部分。若把 y
看作研究区边界，这一步就是“以 y 裁剪
x”；例如裁剪道路时，结果是边界内的线段，并非仅给整条道路贴上区域编号。两个面若仅共享边或顶点，相交结果也可能是线或点，不能一律解释为重叠面积。
并集与溶解
st_union() 合并几何覆盖范围。若 x、y
各只有一个几何，st_union(x, y)
可得到两者的整体范围；若各图层含多个要素，先分别用
st_union(st_geometry(x)) 与
st_union(st_geometry(y))
汇合，再对两者求并集，避免把逐要素运算误读为“整个图层的并集”。
溶解（dissolve）是在同一图层中按组消除内部边界。对于
sf 对象，group_by(group) %>% summarise(...)
默认已联合每组的几何；应同时确定数值属性的汇总方式，不必再对结果整体调用
st_union()，否则会丢失分组。
差集与对称差
st_difference(x, y) 保留 x 中未被 y
覆盖的部分，方向不能颠倒；st_sym_difference(x, y)
则保留双方各自独有的部分。若 y
含多个要素且目标是扣除其整体覆盖范围，可先写
st_union(st_geometry(y))，再计算差集。
下面的五个面板使用同一比例尺。虚线保留原始 A、B 边界，帮助辨认每种运算改变了哪一部分；该图由折叠区中的 R 代码实际生成。
两个示例面在相交、并集、差集与对称差运算前后的范围。
查看并复制空间叠置图的 R 代码
library(sf)
library(ggplot2)
# 1. 在同一米制投影下构造相互重叠的示例面
area_a <- st_sf(
id = "A",
geometry = st_as_sfc(st_bbox(
c(xmin = 300000, ymin = 3450000, xmax = 304000, ymax = 3454000),
crs = st_crs(32651)
))
)
area_b <- st_sf(
id = "B",
geometry = st_as_sfc(st_bbox(
c(xmin = 302000, ymin = 3452000, xmax = 306000, ymax = 3456000),
crs = st_crs(32651)
))
)
a_geom <- st_union(st_geometry(area_a))
b_geom <- st_union(st_geometry(area_b))
# 2. 计算结果；两个输入各只有一个面
stages <- c("处理前", "相交／裁剪", "并集", "差集 A－B", "对称差")
overlay_sf <- st_sf(
stage = factor(c(stages[1], stages[1], stages[-1]), levels = stages),
kind = factor(c("区域 A", "区域 B", rep("运算结果", 4)),
levels = c("区域 A", "区域 B", "运算结果")),
geometry = st_sfc(
st_geometry(area_a)[[1]], st_geometry(area_b)[[1]],
st_intersection(a_geom, b_geom)[[1]],
st_union(a_geom, b_geom)[[1]],
st_difference(a_geom, b_geom)[[1]],
st_sym_difference(a_geom, b_geom)[[1]],
crs = 32651
)
)
plot_family <- "sans"
if (.Platform$OS.type == "windows") {
grDevices::windowsFonts(c4_zh = grDevices::windowsFont("Microsoft YaHei"))
plot_family <- "c4_zh"
}
ggplot(overlay_sf) +
geom_sf(aes(fill = kind), colour = "#384C59", alpha = 0.80,
linewidth = 0.35) +
geom_sf(data = area_a, fill = NA, colour = "#287F91",
linetype = "dashed", linewidth = 0.45) +
geom_sf(data = area_b, fill = NA, colour = "#B56D42",
linetype = "dashed", linewidth = 0.45) +
facet_wrap(~stage, ncol = 3) +
scale_fill_manual(values = c("区域 A" = "#5E9FAE",
"区域 B" = "#DEA46C",
"运算结果" = "#827CB2")) +
coord_sf(datum = NA) +
labs(fill = NULL) +
theme_bw(base_size = 10, base_family = plot_family) +
theme(panel.grid = element_blank(),
axis.text = element_blank(), axis.ticks = element_blank(),
axis.title = element_blank(),
strip.background = element_rect(fill = "#EFF3F5", colour = "#CAD5DA"),
strip.text = element_text(face = "bold"),
legend.position = "bottom")3.2.4 常用空间数据生成
缓冲区（Buffer）
重点 · 常用
st_buffer()
根据给定距离，在点、线或面的周围生成新的面几何。点缓冲形成近似圆形范围，线缓冲形成沿线走廊；面缓冲可向外扩张。在合适的平面投影下，对面使用负距离还可向内收缩。本例的缓冲区表示平面直线距离内的邻近范围，不能直接当作路网通行时间或服务可达性。
距离与 CRS 缓冲距离的单位不能只凭坐标轴上的“度”或“米”判断：
sf在默认 s2 球面计算下，EPSG:4326 的数值dist按米解释；投影数据则按其坐标单位解释。本例选用适合北卡罗来纳州的米制投影 EPSG:32119，并把 5、10、20 km 写为 5000、10000、20000 米。
点、线、面的缓冲区与面内缩（示意）。
下面以 sf 包自带的 nc 县界数据中的 Wake
县为位置参考，构造一个示例点、一条示例线和一个示例面。三个面板采用相同尺度；每个面板同时显示
5、10、20 km
缓冲距离。为避免半透明面相叠后混色，图中把相邻缓冲区相减，显示为
0–5、5–10、10–20 km
距离带。这些对象仅用于演示，缓冲范围没有按县界裁剪。
左右滑动可查看完整的三个面板。
点、线、面的 5、10、20 km 缓冲距离带；玫红色表示原始几何。
查看并复制缓冲区对照图的 R 代码
library(sf)
library(ggplot2)
# 1. 选用米制投影，并在 Wake 县内构造三种教学几何
buffer_nc <- st_read(system.file("shape/nc.shp", package = "sf"), quiet = TRUE)
buffer_wake <- st_transform(buffer_nc[buffer_nc$NAME == "Wake", ], 32119)
buffer_center <- as.numeric(st_coordinates(
st_point_on_surface(st_geometry(buffer_wake))
)[1, 1:2])
buffer_point <- st_point(buffer_center)
buffer_line <- st_linestring(rbind(
buffer_center + c(-9000, -6000),
buffer_center + c(-2000, 4000),
buffer_center + c(9500, 6500)
))
buffer_polygon <- st_polygon(list(rbind(
buffer_center + c(-6000, -4000),
buffer_center + c(5000, -4000),
buffer_center + c(7000, 3000),
buffer_center + c(-3000, 6500),
buffer_center + c(-6000, -4000)
)))
buffer_source <- st_sf(
feature_type = factor(c("点", "线", "面"), levels = c("点", "线", "面")),
geometry = st_sfc(buffer_point, buffer_line, buffer_polygon, crs = 32119)
)
# 2. 先计算三种缓冲半径，再相减形成不重叠的距离带
buffer_band_levels <- c("0–5 km", "5–10 km", "10–20 km")
buffer_bands <- do.call(rbind, lapply(seq_len(nrow(buffer_source)), function(i) {
source_geometry <- st_geometry(buffer_source)[i]
buffer_5 <- st_buffer(source_geometry, dist = 5000)
buffer_10 <- st_buffer(source_geometry, dist = 10000)
buffer_20 <- st_buffer(source_geometry, dist = 20000)
st_sf(
feature_type = buffer_source$feature_type[i],
distance_band = factor(buffer_band_levels, levels = buffer_band_levels),
geometry = st_sfc(
buffer_5[[1]],
st_difference(buffer_10, buffer_5)[[1]],
st_difference(buffer_20, buffer_10)[[1]],
crs = 32119
)
)
}))
# 3. 相同色序与比例尺用于全部面板；原始几何放在最上层
buffer_band_cols <- c(
"0–5 km" = "#276F98",
"5–10 km" = "#55A8A8",
"10–20 km" = "#A5CEB9"
)
buffer_plot_family <- "sans"
if (.Platform$OS.type == "windows") {
grDevices::windowsFonts(c4_zh = grDevices::windowsFont("Microsoft YaHei"))
buffer_plot_family <- "c4_zh"
}
ggplot() +
geom_sf(data = buffer_wake, fill = "#FAFBFC", colour = "#C6D2D8",
linewidth = 0.30) +
geom_sf(data = buffer_bands, aes(fill = distance_band),
colour = "#547987", linewidth = 0.32, alpha = 0.74) +
geom_sf(data = buffer_source[buffer_source$feature_type == "点", ],
colour = "#A53F5D", size = 2.6) +
geom_sf(data = buffer_source[buffer_source$feature_type == "线", ],
colour = "#A53F5D", linewidth = 1.05) +
geom_sf(data = buffer_source[buffer_source$feature_type == "面", ],
fill = "#A53F5D", colour = "#78304A", linewidth = 0.50) +
facet_wrap(~feature_type, nrow = 1) +
scale_fill_manual(name = "距原几何", values = buffer_band_cols,
breaks = buffer_band_levels) +
scale_x_continuous(
breaks = buffer_center[1] + c(-20000, 0, 20000),
labels = function(x) round((x - buffer_center[1]) / 1000)
) +
scale_y_continuous(
breaks = buffer_center[2] + c(-20000, 0, 20000),
labels = function(y) round((y - buffer_center[2]) / 1000)
) +
coord_sf(
xlim = buffer_center[1] + c(-33000, 33000),
ylim = buffer_center[2] + c(-33000, 33000),
datum = NA, expand = FALSE
) +
labs(x = "相对中心的东西距离（km）", y = "相对中心的南北距离（km）") +
theme_bw(base_size = 12, base_family = buffer_plot_family) +
theme(
panel.grid.major = element_line(colour = "#E9EFF2", linewidth = 0.25),
panel.grid.minor = element_blank(),
axis.title = element_text(size = 10, colour = "#3D5360"),
axis.text = element_text(size = 9.5, colour = "#4C616C"),
strip.background = element_rect(fill = "#EEF3F5", colour = "#C6D2D8"),
strip.text = element_text(face = "bold", size = 11, colour = "#2E4958"),
legend.position = "bottom",
legend.title = element_text(face = "bold", size = 10),
legend.text = element_text(size = 10),
legend.key.width = grid::unit(1.6, "lines"),
panel.spacing = grid::unit(0.55, "lines")
)重心（Centroid）
常用
st_centroid() 可从面或线几何计算几何质心（centroid），常用于标注或表示位置。凹形、多部件或带孔多边形的质心可能位于面外；需要保证点位落在面内时，可用
st_point_on_surface()。
面质心在面内与面外的两种情况（示意）。
下面先将北卡罗来纳州各县边界转换到适合当地的平面投影，再把面几何替换为质心点，同时保留县名等属性。对几何列单独计算，可避免把整张属性表传入
st_centroid() 时出现“属性在几何范围内保持恒定”的提示。
library(sf)
library(ggplot2)
nc <- st_read(system.file("shape/nc.shp", package = "sf"), quiet = TRUE)
nc <- st_transform(nc, 32119)
nc_cent <- nc
st_geometry(nc_cent) <- st_centroid(st_geometry(nc))
# 几何从面变为点，县名等属性仍保留
data.frame(
county = head(nc_cent$NAME),
geometry_type = as.character(head(st_geometry_type(nc_cent)))
) county geometry_type
1 Ashe POINT
2 Alleghany POINT
3 Surry POINT
4 Currituck POINT
5 Northampton POINT
6 Hertford POINTggplot() +
geom_sf(data = nc, fill = "#EFF3F5", colour = "#8296A0",
linewidth = 0.25) +
geom_sf(data = nc_cent, colour = "#287F91", size = 1.0) +
theme_bw(base_size = 10) +
theme(panel.grid.minor = element_blank(),
axis.title = element_blank())
北卡罗来纳州各县边界及其几何质心。数据：sf 内置 nc。
4.0 空间数据可视化
4.0 空间数据可视化
4.0 空间可视化指南
制图基础
4.1 静态制图基础：sf ×
ggplot2::geom_sf()
下图以英国区域分层设色图为例，串联原始文件、sf
对象、属性整理、geom_sf()
设色与地图整饰。图中的数值用于说明步骤，后面的代码示例使用各自注明的数据。
sf
负责记录矢量几何（点、线、面）、属性与坐标参考系（CRS）；ggplot2
负责图层组织和图形呈现。geom_sf() 绘制 sf
图层，coord_sf()
控制地图的显示范围和显示坐标系。coord_sf(crs = ...)
不会改写原始数据的坐标；需要转换数据本身时使用
st_transform()。
因此，静态地图的基本写法可沿用第三章已掌握的 ggplot2 结构：以
ggplot() 建立画布，通过 + 逐层叠加
geom_sf()、coord_sf()、labs() 与
theme_*() 等组件。整体流程可概括为：数据准备 →
图层叠加 → 坐标控制 → 主题布局 → 导出输出。
以英国区域为例的 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)
按范围和标签查询水域、绿地、道路与建筑等矢量要素。以下案例用于演示图层叠加；在线查询结果会随
OSM 编辑及服务状态变化。
知识拓展：OpenStreetMap (OSM) 与 osmdata 包
OpenStreetMap (OSM) 是由志愿者协作维护的开放地理数据库，包含建筑、道路、兴趣点等要素。不同地区的覆盖率、位置精度和属性完整度并不相同，使用前应检查所需标签与研究区的数据质量。
osmdata 通过 Overpass API 查询 OSM
要素，并可将结果转为 sf
对象。按地名查找外接矩形时还依赖在线地名解析服务；外接矩形并不等于行政边界。
📚 推荐参考资料与进阶指南：
- OSM Map Features ：https://wiki.openstreetmap.org/wiki/Map_features
- osmdata 官方教程：https://docs.ropensci.org/osmdata/
OSM 以键值标签（tag）描述要素。查询前先查阅上面的标签文档，并抽查当地实际记录；同一类设施可能使用不同标签，标签缺失也不等于设施不存在。
在完成要素提取后，依托 ggplot2
的图层叠加逻辑即可构建综合地图。
本案例在制图过程中结合了以下几个关键点：
1.
多图层叠加：分别显示水域、绿地与建筑，并用颜色区分类别；
2. 地图辅助要素：用
ggspatial 添加比例尺和指北针；
3.
经纬网与范围：用 coord_sf()
限定显示范围，以细线提供坐标参照。
运行提示：以下代码需要连接地名解析服务和 Overpass API；网络慢、请求限流或服务暂不可用时可能失败。课堂页面展示的是预先保存的示意图，代码设置为不在 Knit 时联网运行。复制代码重新运行所得 OSM 要素和画面可能与示意图不同；研究使用时应记录查询日期，并注明 © OpenStreetMap contributors。
代码 · 示例3
# 如果未安装，请先运行 install.packages("osmdata")
library(tidyverse)
library(sf)
library(osmdata) # 用于通过 API 实时获取 OpenStreetMap 矢量数据
library(ggspatial)
# ============================================================
# 1. 定义查询范围
# ============================================================
# 根据地名获取外接矩形；它仅限定查询范围，并非行政边界
bbox_manhattan <- getbb("Lower Manhattan, New York")
# ============================================================
# 2. 构建查询并分别获取各个图层的数据
# ============================================================
# 以下会向 Overpass API 发送 3 次请求；耗时取决于网络和服务状态
# 水系 (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
# 某类要素可能未返回面几何；空图层不添加到地图
polygon_layer <- function(data, label) {
if (is.null(data) || nrow(data) == 0L) return(NULL)
geom_sf(data = data, aes(fill = label), colour = NA)
}
# ============================================================
# 4. 图层叠加与城市综合地图绘制
# ============================================================
ggplot() +
polygon_layer(water_poly, "Water") +
polygon_layer(water_multi, "Water") +
polygon_layer(parks_poly, "Parks") +
polygon_layer(buildings_poly, "Buildings") +
# 指北针与比例尺依赖 ggspatial
annotation_north_arrow(location = "tl", which_north = "true",
pad_x = grid::unit(0.1, "in"), pad_y = grid::unit(0.1, "in"),
style = north_arrow_fancy_orienteering()) +
# br 为右下角
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(
panel.background = element_rect(fill = "#f2efe9", colour = NA),
panel.grid.major = element_line(color = "#c8c8c4", linetype = "dotted", linewidth = 0.25),
panel.grid.minor = element_blank(),
panel.border = element_rect(color = "grey35", fill = NA, linewidth = 0.6),
axis.text = element_text(color = "grey40", size = 8),
axis.ticks = element_line(color = "grey40"),
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),
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 = "Water, parks and buildings in an OSM query area",
caption = "Data: © OpenStreetMap contributors; queried via osmdata",
x = NULL, y = NULL
)
预存的 OSM 多图层示意图；重新查询所得要素可能变化。
分层设色
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）往往会被面积或人口规模主导，解释需谨慎。
代码示例 1：spData::world
绘制人口密度分层设色图
示例代码
下述示例以 spData::world
的国家边界为制图单元，利用其中的 2014 年人口
pop 与面积 area_km2（km²）构造人口密度
pop_density = pop / area_km2（人/平方千米），并将该指标映射到
fill。这是一份包内固定年份的示例数据。
世界人口密度：spData::world 的 2014 年人口除以面积；色标使用对数刻度，灰色表示缺失。不显示南极洲。
查看并复制世界人口密度图代码
library(spData)
library(dplyr)
library(ggplot2)
# 1. 计算人/平方千米；南极洲不参与人口密度比较。
data("world", package = "spData")
world_density <- world %>%
filter(name_long != "Antarctica") %>%
mutate(pop_density = if_else(pop > 0 & area_km2 > 0,
pop / area_km2, NA_real_))
# 2. 对数色标拉开偏态数值；缺失区域直接由 na.value 着灰色。
ggplot(world_density) +
geom_sf(aes(fill = pop_density), colour = "#F7F9FB", linewidth = 0.12) +
scale_fill_viridis_c(
option = "magma", direction = -1, trans = "log10",
breaks = c(0.1, 1, 10, 100, 1000),
labels = c("0.1", "1", "10", "100", "1,000"),
na.value = "#D9DEE4", name = "人口密度\n(人/km²)"
) +
coord_sf(ylim = c(-60, 85), expand = FALSE) +
labs(x = NULL, y = NULL) +
theme_minimal(base_size = 10) +
theme(
panel.background = element_rect(fill = "#FAFBFC", colour = NA),
panel.grid.major = element_line(colour = "#E6EBEF", linewidth = 0.22),
panel.grid.minor = element_blank(),
axis.text = element_text(colour = "#566270", size = 8.5),
legend.title = element_text(face = "bold", size = 9),
legend.text = element_text(size = 8.5),
plot.margin = margin(5, 8, 5, 5)
)代码示例 2：基于 ACS 固定快照的纽约市贫困率空间可视化
拓展 · 专业 · 常用
本示例使用美国人口普查局的 ACS 2020—2024
年五年估计数据，以纽约市的 Census Tract
为单元绘制贫困率分层设色图。ACS
是抽样调查，图中的比率由人口估计值计算，不能视作每个普查区的精确计数。
- ACS 是美国最权威的官方社会经济统计体系之一，数据口径统一且公开，广泛应用于社区尺度的空间差异分析与制图。
- Census Tract 是美国人口普查局设置的统计分区，设计人口通常约为 4,000 人；其边界并不等同于居民实际认知的邻里边界。
核心处理流程：
- 数据准备：
读取课程提供的 ACS 普查区固定快照；在线获取与更新方法见下方折叠代码。 - 指标计算：
以贫困人口估计值除以可判定贫困状态的人口估计值，得到poverty_rate。 - 离散分级：
使用classInt计算 Fisher–Jenks 断点，将贫困率划为若干等级。 - 制图输出：
以geom_sf()映射填充色；没有估计值及分母过小的区域单独编码，图例只显示本期实际出现的类别。
尺度提示
tidycensus可按州、县、Census Tract、Block Group 等尺度查询；Census Block 比 Block Group 更细。不同尺度的估计值和地图形状不能直接当作相同对象比较。
纽约市普查区贫困率：ACS 2020—2024 年五年估计，按 Fisher–Jenks 分为六级；灰色为分母≤100 的区域。
查看并复制纽约贫困率地图代码
数据准备与分级
library(sf)
library(dplyr)
library(ggplot2)
library(classInt)
library(ggspatial)
# 1. 读取固定快照；分母≤100 的普查区不参与率值分级。
nyc_tracts <- st_read(
"data/acs-nyc/nyc_acs_tracts_2020_2024.gpkg",
layer = "nyc_acs_tracts", quiet = TRUE
)
acs_period <- unique(nyc_tracts$acs_period)
nyc_tracts_map <- nyc_tracts %>%
mutate(
data_status = case_when(
!is.finite(totE) | !is.finite(povE) ~ "无可用估计",
totE <= 100 ~ "分母≤100",
TRUE ~ "可计算"
),
poverty_rate = if_else(data_status == "可计算", povE / totE, NA_real_)
) %>%
st_transform(2263) # 纽约本地投影，供地图与比例尺使用
# 2. 只用可计算的率值确定六级 Fisher–Jenks 断点。
poverty_values <- nyc_tracts_map$poverty_rate
poverty_values <- poverty_values[is.finite(poverty_values)]
poverty_breaks <- classInt::classIntervals(
poverty_values, n = 6, style = "fisher", warnSmall = FALSE
)$brks
poverty_labels <- paste0(
scales::label_percent(accuracy = 0.1)(head(poverty_breaks, -1)),
" – ",
scales::label_percent(accuracy = 0.1)(tail(poverty_breaks, -1))
)
nyc_tracts_map <- nyc_tracts_map %>%
mutate(
poverty_class = cut(poverty_rate, breaks = poverty_breaks,
include.lowest = TRUE, labels = poverty_labels),
poverty_class = if_else(data_status == "可计算",
as.character(poverty_class), data_status),
poverty_class = factor(
poverty_class,
levels = c(poverty_labels, "分母≤100", "无可用估计")
)
)
# 3. 固定颜色、图例顺序与纽约市外轮廓。
poverty_colors <- stats::setNames(
viridis::viridis(6, option = "mako", direction = -1), poverty_labels
)
poverty_colors <- c(poverty_colors,
"分母≤100" = "#CED4D8", "无可用估计" = "#858F97")
poverty_legend <- c(rev(poverty_labels), "分母≤100", "无可用估计")
nyc_outline <- st_sf(geometry = st_union(st_geometry(nyc_tracts_map)))地图绘制
ggplot(nyc_tracts_map) +
geom_sf(aes(fill = poverty_class), colour = NA) +
geom_sf(data = nyc_outline, fill = NA,
colour = "#24313A", linewidth = 0.35) +
scale_fill_manual(
values = poverty_colors, breaks = poverty_legend,
drop = TRUE, name = "贫困率\n(Fisher–Jenks)"
) +
ggspatial::annotation_scale(
location = "br", width_hint = 0.22,
height = grid::unit(0.13, "cm"),
pad_x = grid::unit(0.25, "cm"),
pad_y = grid::unit(0.20, "cm")
) +
# 左上留白放指北针，避免遮挡 Staten Island 的普查区。
ggspatial::annotation_north_arrow(
location = "tl", which_north = "true",
style = ggspatial::north_arrow_minimal(
line_col = "#24313A", text_col = "#24313A", text_size = 9
),
height = grid::unit(0.85, "cm"),
width = grid::unit(0.85, "cm"),
pad_x = grid::unit(0.35, "cm"),
pad_y = grid::unit(0.35, "cm")
) +
coord_sf(expand = TRUE) +
labs(x = NULL, y = NULL) +
theme_bw(base_size = 10) +
theme(
panel.background = element_rect(fill = "#FCFDFE", colour = NA),
panel.grid.major = element_line(colour = "#E6EBEF", linewidth = 0.22),
panel.grid.minor = element_blank(),
legend.key.height = grid::unit(0.52, "cm"),
legend.title = element_text(face = "bold", size = 9),
legend.text = element_text(size = 8.5),
plot.margin = margin(5, 8, 5, 5)
)读图提示 “分母≤100” 是本例为避免极小分母造成不稳定比率而设的教学阈值，并非 Census 官方缺失编码。区域面积也不代表人口规模；相邻普查区颜色接近时，还应结合 ACS 的误差范围判断差异。
首次使用 Census API：申请并保存密钥
课程地图直接读取本地数据，因此离线 Knit 不需要密钥。若要自行下载新年份数据，按以下步骤申请；密钥只保存在本机，不写进 Rmd，也不要上传到公开仓库。
02 填写申请表
填写组织名称与邮箱地址，同意服务条款，然后提交申请。
Census API key
申请表的组织、邮箱和提交按钮
03 查收邮件并在本机保存：打开申请邮箱中的 Census Bureau 邮件，按邮件提示完成密钥激活。在 RStudio Console 将下方占位文字换成自己的密钥，运行一次后重启 R 会话。教材代码使用环境变量读取，不展示密钥内容。
# 仅在 Console 运行；不要把真实密钥保存进课程 Rmd
tidycensus::census_api_key("替换为邮件中的密钥", install = TRUE)
# 重启 R 会话后检查：只显示 TRUE/FALSE，不显示密钥
nzchar(Sys.getenv("CENSUS_API_KEY"))截图来源：美国人口普查局公开申请页面；邮件画面因包含私人邮箱与密钥，不在教材中展示。
选学：用个人 API key 在线获取新年份 ACS 数据
下方代码只在主动更新数据时运行。个人密钥从本机环境变量读取；请求得到的新数据不会覆盖本章固定快照。变量编码与年份应以 Census 当期目录核对。
api_key <- Sys.getenv("CENSUS_API_KEY")
if (!nzchar(api_key)) stop("请先在 Console 保存并激活个人 Census API key")
nyc_acs_new <- tidycensus::get_acs(
geography = "tract", state = "NY",
county = c("Bronx", "Kings", "New York", "Queens", "Richmond"),
variables = c(tot = "B17001_001", pov = "B17001_002",
med = "B19013_001"),
survey = "acs5", year = 2024, geometry = TRUE,
output = "wide", key = api_key
)
head(sf::st_drop_geometry(nyc_acs_new))制图数据分级策略：为何选择 Fisher-Jenks 算法？
下一节详讲
分层设色图可以使用连续色带，也可以使用离散等级。对偏态数据，线性连续色带有时会使多数地区的颜色过于接近；此时可以比较转换色标与分级设色的效果。
纽约示例采用 Fisher–Jenks
自然断点法，沿排序后的数值寻找使各组内部差异较小的断点。断点取决于当前数据；不同年份的地图若需直接比较，应使用共同的断点。
- 原理机制：计算断点
breaks（阈值序列），将连续数值划分为离散区间，同一区间映射为同一颜色； - 实现路径：调用
classInt::classIntervals()计算断点，再用cut()将连续变量转为等级变量。
建议 分级数量不宜过多，通常控制在 5–7 级，以在视觉区分度与认知负担之间取得平衡。下一小节将进一步介绍其他常用分类方法。
代码示例 3：日本大阪府主要城市的人口密度
本例把日本国土交通省 N03（2024）的市町村边界与日本统计年鉴 2026 表 2-4中的 2024 年 1 月 1 日城市人口按市名连接。大阪府共有 43 个市町村；年鉴覆盖其中 33 个市，其余 9 町和 1 村显示为灰色。人口密度以快照中由边界计算的面积（km²）为分母，是教学计算值，并非官方发布的密度指标。数据来源与快照说明。
这里用连续色带显示密度，并描出人口至少 30 万的 7 个城市。这个阈值仅用于本图突出主要城市，不代表日本官方的城市等级。与纽约示例相比，制图单元从普查区变为市界，色标也从六级离散色块变为连续梯度。
大阪府城市人口密度（2024）；灰色为年鉴未收录的町村。
查看并复制大阪府分层设色图代码
本地数据连接
library(sf)
library(dplyr)
library(readr)
library(ggplot2)
library(ggspatial)
library(ggrepel)
# 1. 按市町村名连接边界与人口、面积；无人口记录的町村保留为 NA。
osaka_boundary <- st_read(
"data/japan-osaka/osaka_municipalities_2024.gpkg",
layer = "municipalities", quiet = TRUE
)
osaka_population <- read_csv(
"data/japan-osaka/osaka_population_2024.csv",
show_col_types = FALSE
)
osaka_map <- osaka_boundary %>%
left_join(
osaka_population %>%
select(municipality_ja, municipality_en, population_2024, area_km2),
by = "municipality_ja"
) %>%
st_transform(32653) # UTM 53N：供地图与比例尺使用
# 2. 人口 ÷ 边界面积（km²）得到密度；30 万人是本图的标注阈值。
osaka_map <- osaka_map %>%
mutate(
density_per_km2 = round(population_2024 / area_km2, 1),
major_city = !is.na(population_2024) & population_2024 >= 300000
)
# 3. 只标注主要城市；标签使用年鉴中的罗马字城市名。
osaka_major <- osaka_map %>% filter(major_city)
osaka_points <- st_point_on_surface(osaka_major)
osaka_labels <- bind_cols(
st_drop_geometry(osaka_points),
as.data.frame(st_coordinates(osaka_points))
) %>%
mutate(city_label = sub("-shi$", "", municipality_en))
osaka_outline <- st_sf(geometry = st_union(st_geometry(osaka_map)))地图绘制
ggplot(osaka_map) +
geom_sf(aes(fill = density_per_km2),
colour = "#FFFFFF", linewidth = 0.18) +
geom_sf(data = osaka_major, fill = NA,
colour = "#25384C", linewidth = 0.46) +
geom_sf(data = osaka_outline, fill = NA,
colour = "#485B68", linewidth = 0.36) +
ggrepel::geom_label_repel(
data = osaka_labels,
aes(x = X, y = Y, label = city_label),
inherit.aes = FALSE, seed = 42,
size = 2.6, fontface = "bold", colour = "#25384C",
fill = "#FFFFFF", alpha = 0.94, label.size = 0.16,
box.padding = 0.22, point.padding = 0,
min.segment.length = 0, segment.colour = "#667783",
max.overlaps = Inf
) +
scale_fill_viridis_c(
option = "cividis", direction = -1,
labels = scales::label_comma(), na.value = "#E5EAEE",
name = "人口密度\n(人/km²)"
) +
ggspatial::annotation_scale(
location = "bl", width_hint = 0.20,
pad_x = grid::unit(0.24, "cm"),
pad_y = grid::unit(0.22, "cm")
) +
coord_sf(expand = TRUE) +
labs(x = NULL, y = NULL) +
theme_minimal(base_size = 10) +
theme(
panel.background = element_rect(fill = "#FBFCFD", colour = NA),
panel.grid.major = element_line(colour = "#E5EBF0", linewidth = 0.22),
panel.grid.minor = element_blank(),
axis.text = element_text(colour = "#586774", size = 8.5),
legend.title = element_text(face = "bold", size = 9),
legend.text = element_text(size = 8.5),
plot.margin = margin(6, 10, 6, 6)
)读图提示 灰色表示年鉴未收录的町村，不表示人口为零。人口密度能比较不同面积城市的居住人口集中程度；市界内部的分布仍可能很不均匀。
数据分类分色
4.3 数据分类分色：同一县级地图，不同断点
重要 · 常用
连续色带将数值逐步映射为颜色。若读图任务需要按明确区间识别地区，可以先把数值划为少数等级，再让每级对应一种颜色。本节用同一批县级单元比较三种断点的效果。
分级的用途与边界
分级会合并同一级内的差异；断点两侧十分接近的数值，也可能显示为两种颜色。若有政策或研究上明确的阈值，应优先使用有依据的固定断点。比较不同年份的绝对变化时，各年份还须使用同一组断点。
本例使用美国人口调查局的2025
年县级制图边界和2025
年县级人口估计。人口密度由 2025 年 7 月 1 日人口估计 ÷
官方陆地面积 ALAND
得到，单位为人/km²。地图仅保留美国本土连续 48 州，共 3,108
个县及县级等价单元；不纳入阿拉斯加、夏威夷、华盛顿特区、波多黎各及其他美国属地。这里的面积采用边界文件附带的官方字段，不从概化后的地图形状重新计算。
县级密度差异很大：大多数县处于较低区间，少数高度城市化县位于长尾。若直接按原始密度分级，地图容易由最低一级占据大部分面积。为让三种方法的差异可读，先对每个县计算
log10(1 + 人口密度)，再在这组变换后的数值上分别求五级断点。变换参与分级计算；图例则把断点逆变换回人/km²，便于按原始单位读图。对数变换保留密度的高低顺序，并压缩高值间距，但不会消除区内差异或改变县界。
连续参照图和三种分级图始终使用同一批县级单元、同一投影与范围，以及同方向的
mako
色带：浅薄荷绿表示较低密度，深紫色表示较高密度。三幅分级图只更换断点；跨方法的同一种颜色可能代表不同的人/km²区间，需分别阅读图下方的区间说明。左侧直方图使用固定的对数宽度分箱，并按右侧地图所用等级堆叠着色；虚线是该方法的四个内部断点。横轴虽标出原始人/km²数值，柱宽对应的是相同的
log10(1 + 密度)
宽度，不能按相同的人/km²跨度理解。
连续色带参照
先用同一批县级数据画一张连续色带地图。颜色位置按照
log10(1 + 密度)
计算，图例仍标原始人/km²。对数色标让低密度县之间的差异比线性色标更容易辨认；颜色连续变化，不预先合并成五级。
连续色带参照：美国本土 48 州的县级人口密度（2025 年 7 月 1 日人口估计）。
读图提示：以下三组图用于比较分级方法，密度区间与县级单元数已列在各组图下。方向和距离可参照上方连续地图，因此不在每幅分级地图内重复放置指北针、比例尺和图内图例。比例尺仅供其所在位置附近参考。
1. Fisher–Jenks
自然断点法（style = "fisher"）
Fisher–Jenks
在排序后的对数密度中寻找断点，使同级数值尽量接近。它适合观察这批县级单元的聚集与间隙，断点由本次数据决定。左图的彩色柱段与右图的县级颜色来自同一份分级结果；竖向虚线标出四个内部断点。全国县级数据超过
classInt 的默认抽样阈值，代码设置
largeN = Inf，让自然断点使用全量数据且结果可复现。
A · 对数密度分布与 Fisher–Jenks 断点
B · 同一分级在县级地图上的颜色
- 0.03–3.9 人/km²557 个县级单元
- 3.9–15.7 人/km²894 个县级单元
- 15.7–54.4 人/km²931 个县级单元
- 54.4–269.6 人/km²514 个县级单元
- 269.6–28,370 人/km²212 个县级单元
2.
等量分位法（style = "quantile"）
分位数法按排序位置划分，使每级县级单元数尽量接近，适合比较地区在当前样本中的相对位置。由于对数变换保留原始密度的高低顺序，这里的分位数分组与直接对原始密度求分位数相同；图例中的数值仍须读取本组断点。各级单元数接近，不意味着地图上各色覆盖的面积接近。
A · 对数密度分布与分位数断点
B · 相同县级地图按分位数重新设色
- 0.03–4.7 人/km²622 个县级单元
- 4.7–12.3 人/km²621 个县级单元
- 12.3–25.0 人/km²622 个县级单元
- 25.0–66.4 人/km²621 个县级单元
- 66.4–28,370 人/km²622 个县级单元
3.
等间距法（style = "equal"）
等间距法把对数密度的最小值到最大值划为五段相同跨度。对数尺度上的相等间距，逆变换后对应的是逐级扩大的原始人/km²区间，因此下方图例不会显示五个相同的数值跨度。这种做法让高密度长尾不至于把绝大多数县都压进最低一级，也说明分级前的数据变换会改变等间距法的含义。若研究要求原始密度每级固定增加相同的人/km²，应直接在未变换的密度上设置断点，并接受地图可能较为单色的结果。
A · 对数密度分布与等间距断点
B · 相同县级地图按对数等间距重新设色
- 0.03–7.0 人/km²808 个县级单元
- 7.0–60.5 人/km²1626 个县级单元
- 60.5–474.2 人/km²548 个县级单元
- 474.2–3,670.7 人/km²117 个县级单元
- 3,670.7–28,370 人/km²9 个县级单元
三组图中的县级单元、直方图分箱、地图范围和颜色方向都保持不变。图例端点经过四舍五入，实际分级使用未取整的断点。由浅到深的五级单元数如下。颜色覆盖的地图面积不能代替该级的单元数，同色县的原始密度也不一定相近。
| 分级方法 | 由浅到深五级的县级单元数 |
|---|---|
| Fisher–Jenks | 557 / 894 / 931 / 514 / 212 |
| 分位数 | 622 / 621 / 622 / 621 / 622 |
| 对数等间距 | 808 / 1626 / 548 / 117 / 9 |
查看并复制：读取县级数据、三种分级与绘图代码
stopifnot()
用于核对县界、人口和密度计算结果；若数据不匹配或出现无效数值，代码会停止，避免把缺失值带入地图。
一次读取与密度构建
library(sf)
library(dplyr)
library(readr)
library(ggplot2)
library(classInt)
# 1. 读取同一版官方快照。课程附带本地文件；独立复制代码时可按固定网址下载。
get_class_43_snapshot <- function(local_path, source_url) {
if (file.exists(local_path)) return(local_path)
target <- file.path(tempdir(), basename(local_path))
if (!file.exists(target)) {
download.file(source_url, target, mode = "wb", quiet = TRUE)
}
target
}
county_zip <- get_class_43_snapshot(
"data/census-county-2025/cb_2025_us_county_5m.zip",
"https://www2.census.gov/geo/tiger/GENZ2025/shp/cb_2025_us_county_5m.zip"
)
state_zip <- get_class_43_snapshot(
"data/census-2025/cb_2025_us_state_20m.zip",
"https://www2.census.gov/geo/tiger/GENZ2025/shp/cb_2025_us_state_20m.zip"
)
population_csv <- get_class_43_snapshot(
"data/census-county-2025/co-est2025-alldata.csv",
"https://www2.census.gov/programs-surveys/popest/datasets/2020-2025/counties/totals/co-est2025-alldata.csv"
)
county_dir <- file.path(tempdir(), "class_43_county_shape")
state_dir <- file.path(tempdir(), "class_43_state_shape")
dir.create(county_dir, showWarnings = FALSE)
dir.create(state_dir, showWarnings = FALSE)
utils::unzip(county_zip, exdir = county_dir)
utils::unzip(state_zip, exdir = state_dir)
county_boundary <- st_read(
file.path(county_dir, "cb_2025_us_county_5m.shp"), quiet = TRUE
)
state_boundary <- st_read(
file.path(state_dir, "cb_2025_us_state_20m.shp"), quiet = TRUE
)
# 州界快照确定本土 48 州的代码；县界另含其他属地，不能只排除四个代码。
class_43_state_outline <- state_boundary %>%
filter(!GEOID %in% c("02", "15", "11", "72"))
statefp_48 <- class_43_state_outline$GEOID
stopifnot(length(statefp_48) == 48L, !anyDuplicated(statefp_48))
county_population <- read_csv(
population_csv,
col_types = cols(
SUMLEV = col_character(), STATE = col_character(),
COUNTY = col_character(), POPESTIMATE2025 = col_double(),
.default = col_skip()
),
locale = locale(encoding = "Latin1"),
show_col_types = FALSE
) %>%
filter(SUMLEV == "050") %>%
transmute(
GEOID = paste0(STATE, COUNTY),
population_2025 = POPESTIMATE2025
)
county_boundary_48 <- county_boundary %>%
filter(STATEFP %in% statefp_48)
unmatched_counties <- anti_join(
st_drop_geometry(county_boundary_48) %>% select(GEOID),
county_population %>% select(GEOID), by = "GEOID"
)
stopifnot(nrow(county_boundary_48) == 3108L, nrow(unmatched_counties) == 0L)
class_43_counties <- county_boundary_48 %>%
left_join(county_population, by = "GEOID", relationship = "one-to-one") %>%
mutate(
density_2025 = population_2025 / (as.numeric(ALAND) / 1e6),
log_density = log10(1 + density_2025)
) %>%
st_transform(5070) # 美国本土等面积投影；三种方法使用同一幅边界。
class_43_state_outline <- st_transform(class_43_state_outline, 5070)
stopifnot(
nrow(class_43_counties) == 3108L,
all(is.finite(class_43_counties$log_density)),
all(class_43_counties$population_2025 > 0),
all(as.numeric(class_43_counties$ALAND) > 0)
)
# 精确配色：viridisLite::mako(n, direction = -1, begin = 0.12, end = 0.88)。
# 低值浅薄荷绿、高值深紫色；六位色值同时供 ggplot 图形和 HTML 图例使用。
class_43_palette <- function(n) {
substr(viridisLite::mako(
n, direction = -1, begin = 0.12, end = 0.88
), 1, 7)
}
class_43_colors <- stats::setNames(
class_43_palette(5), as.character(1:5)
)通用分级与绘图函数
# 2. 所有方法均在 log10(1 + 密度) 上分级，再把断点逆变换为原始单位。
classify_county_density <- function(style) {
breaks_log <- classInt::classIntervals(
class_43_counties$log_density,
n = 5, style = style, largeN = Inf
)$brks # Fisher 在约 3100 个不同数值上使用全量数据，避免默认随机抽样。
stopifnot(length(breaks_log) == 6L, all(diff(breaks_log) > 0))
map <- class_43_counties %>%
mutate(level = factor(
cut(log_density, breaks = breaks_log, right = TRUE,
include.lowest = TRUE, labels = FALSE),
levels = 1:5
))
stopifnot(!anyNA(map$level))
list(
map = map,
breaks_log = breaks_log,
breaks_density = 10^breaks_log - 1,
counts = tabulate(as.integer(map$level), nbins = 5)
)
}
# 所有方法共用对数轴上的 0.2 宽分箱：图中每根柱覆盖同样的对数距离。
class_43_hist_breaks <- seq(0, 4.6, by = 0.2)
class_43_hist_ticks <- c(0, 10, 100, 1000, 10000)
class_43_hist_max <- max(hist(
class_43_counties$log_density,
breaks = class_43_hist_breaks, plot = FALSE, right = TRUE
)$counts)
plot_class_hist <- function(result) {
ggplot(st_drop_geometry(result$map),
aes(x = log_density, fill = level)) +
geom_histogram(
breaks = class_43_hist_breaks, closed = "right",
position = "stack", colour = "white", linewidth = 0.15
) +
# 白色细底线使深色柱中的红色断点仍可辨认。
geom_vline(
xintercept = result$breaks_log[2:5],
linetype = "longdash", colour = "white", linewidth = 1.15
) +
geom_vline(
xintercept = result$breaks_log[2:5],
linetype = "longdash", colour = "#9F3D59", linewidth = 0.55
) +
scale_fill_manual(values = class_43_colors, drop = FALSE, guide = "none") +
scale_x_continuous(
limits = range(class_43_hist_breaks),
breaks = log10(1 + class_43_hist_ticks),
labels = c("0", "10", "100", "1,000", "10,000"),
expand = c(0, 0)
) +
scale_y_continuous(
limits = c(0, ceiling(class_43_hist_max * 1.1)),
breaks = scales::breaks_pretty(n = 4),
expand = c(0, 0)
) +
labs(x = "人口密度（人/km²；对数展示）", y = "县级单元数") +
theme_classic(base_size = 11.5, base_family = "sans") +
theme(
axis.line = element_line(colour = "#52636D", linewidth = 0.32),
axis.text = element_text(colour = "#4E5D66", size = 11.5),
axis.title = element_text(colour = "#263943", size = 12),
panel.grid.major.y = element_line(colour = "#EDF1F3", linewidth = 0.24),
plot.margin = margin(10, 9, 9, 9)
)
}
plot_class_map <- function(result) {
ggplot() +
geom_sf(
data = result$map, aes(fill = level),
colour = "#E2E8E8", linewidth = 0.06
) +
geom_sf(
data = class_43_state_outline,
fill = NA, colour = "#64747C", linewidth = 0.28
) +
scale_fill_manual(
values = class_43_colors, limits = as.character(1:5),
drop = FALSE, guide = "none"
) +
coord_sf(crs = st_crs(5070), datum = NA, expand = FALSE) +
theme_void(base_size = 11.5, base_family = "sans") +
theme(
panel.background = element_rect(fill = "#FBFCFD", colour = NA),
plot.background = element_rect(fill = "white", colour = NA),
plot.margin = margin(10, 9, 9, 9)
)
}
# 图例展示逆变换后的原始密度；显示值保留一位，实际分类保留完整精度。
print_class_key <- function(result, method_name) {
format_density <- scales::label_number(accuracy = 0.1, big.mark = ",")
lower <- format_density(head(result$breaks_density, -1))
upper <- format_density(tail(result$breaks_density, -1))
lower[1] <- scales::label_number(accuracy = 0.01)(
result$breaks_density[1]
)
upper[5] <- scales::label_number(accuracy = 1, big.mark = ",")(
result$breaks_density[6]
)
items <- vapply(seq_len(5), function(i) {
sprintf(paste0(
'<li><span class="c4-class-swatch" style="background:%s" aria-hidden="true"></span>',
'<span class="c4-class-interval">%s–%s 人/km²</span>',
'<small>%d 个县级单元</small></li>'
), class_43_colors[i], lower[i], upper[i], result$counts[i])
}, character(1))
cat(sprintf(
'<ul class="c4-class-key" aria-label="%s：五级人口密度区间和各级县级单元数">%s</ul>',
method_name, paste(items, collapse = "")
))
}分别生成三种分级结果
class_43_fisher <- classify_county_density("fisher")
class_43_quantile <- classify_county_density("quantile")
class_43_equal <- classify_county_density("equal")绘制连续色带参照图
ggplot() +
geom_sf(
data = class_43_counties, aes(fill = log_density),
colour = "#E2E8E8", linewidth = 0.06
) +
geom_sf(
data = class_43_state_outline,
fill = NA, colour = "#64747C", linewidth = 0.28
) +
ggspatial::annotation_scale(
location = "bl", width_hint = 0.16, style = "ticks",
unit_category = "metric", text_cex = 0.7,
line_width = 0.55,
pad_x = grid::unit(0.3, "cm"),
pad_y = grid::unit(0.35, "cm")
) +
ggspatial::annotation_north_arrow(
location = "tl", which_north = "true",
height = grid::unit(0.75, "cm"),
width = grid::unit(0.75, "cm"),
pad_x = grid::unit(0.3, "cm"),
pad_y = grid::unit(0.3, "cm"),
style = ggspatial::north_arrow_minimal(
text_size = 8, line_width = 0.55,
line_col = "#263943", fill = "#263943", text_col = "#263943"
)
) +
scale_fill_gradientn(
colours = class_43_palette(256), name = "人/km²",
limits = c(0, max(class_43_counties$log_density)),
breaks = log10(1 + c(0, 10, 100, 1000, 10000)),
labels = c("0", "10", "100", "1,000", "10,000")
) +
coord_sf(crs = st_crs(5070), datum = NA, expand = FALSE) +
theme_void(base_size = 11.5, base_family = "sans") +
theme(
panel.background = element_rect(fill = "#FBFCFD", colour = NA),
plot.background = element_rect(fill = "white", colour = NA),
legend.title = element_text(size = 10, face = "bold"),
legend.text = element_text(size = 9),
plot.margin = margin(9, 12, 9, 9)
)运行 plot_class_hist(class_43_fisher) 与
plot_class_map(class_43_fisher) 可查看第一组；将对象名改为
class_43_quantile 或
class_43_equal，即可查看另外两组。代码应在第四章目录运行，使相对路径能找到本地快照。
延伸：固定阈值与其他分级方法
标准差法可表达相对均值的偏离；自定义分级适用于有理论、政策或专业依据的固定阈值。比较多个年份的绝对密度时，可用
style = "fixed" 与 fixedBreaks = c(...)
设置一套共同断点。若在对数密度上分级，应先把原始单位的固定阈值转换成
log10(1 + 阈值)，并在图例中再换算回来。
参考资料：《Geocomputation with R》 的地图制作章节。
底图绘制
4.4 底图：空间语境
常用
地图可以用行政边界、道路等矢量图层提供位置参照，也可以叠加在线瓦片底图（basemap tiles）显示街道和地名。先用本地数据演示“边界＋点位”的离线参照图，再看在线瓦片提供的街区细节。
底图应服务于主题信息。瓦片的标签、颜色与更新时间由提供方决定；其内容也可能与叠加的 OSM 点位并非同一时刻。
制图思路
叠加在线瓦片时，先确保主题点位清晰，再选择合适的底图样式：
- 色彩退隐（Color
Muting）：
若主题图层已用颜色编码类别，宜选择低饱和度瓦片，避免道路与建筑色彩干扰点位识别。 - 标签克制（Label
Control）：
街道与地名标签过密时，可选标签较少的瓦片样式，并检查缩放后点位是否仍然可见。 - 暗色底图（Dark
Basemap）：
暗色瓦片（如 Carto Dark Matter）适合与浅色点位搭配；需要同时检查比例尺、坐标文字和图例的对比度。这里的底图样式与网页的浅色／深色阅读模式是两个独立设置。
离线示例：行政区边界作为底图
伦敦行政区边界和共享单车站点都由 spData
随包提供。图中的浅灰色行政区就是离线矢量底图，红色站点是叠加其上的主题图层；运行时无需下载在线瓦片。边界提供区域位置参照，街道与地名则留给后面的瓦片示例展示。
伦敦行政区底图与共享单车站点（spData；站点数据源：OpenStreetMap）。
查看并复制离线底图的 R 代码
library(sf)
library(spData)
library(ggplot2)
data("lnd", package = "spData")
data("cycle_hire_osm", package = "spData")
# 两层统一 CRS；外轮廓用于区分整个伦敦范围与区内边界
cycle_hire_lnd <- st_transform(cycle_hire_osm, st_crs(lnd))
lnd_outline <- st_sf(geometry = st_union(st_geometry(lnd)))
ggplot() +
geom_sf(data = lnd, fill = "#F1F5F6", colour = "#A4B6BE", linewidth = 0.25) +
geom_sf(data = lnd_outline, fill = NA,
colour = "#667F8B", linewidth = 0.48) +
geom_sf(data = cycle_hire_lnd, aes(fill = "共享单车站点"),
shape = 21, colour = "white", stroke = 0.19,
size = 1.38, alpha = 0.88) +
scale_fill_manual(NULL, values = c("共享单车站点" = "#B84B63"),
guide = guide_legend(override.aes = list(size = 3.2, alpha = 1))) +
ggspatial::annotation_scale(
location = "bl", width_hint = 0.19, style = "ticks",
line_width = 0.35,
pad_x = grid::unit(0.26, "cm"),
pad_y = grid::unit(0.28, "cm")
) +
ggspatial::annotation_north_arrow(
location = "tr", which_north = "true",
height = grid::unit(0.7, "cm"), width = grid::unit(0.7, "cm"),
pad_x = grid::unit(0.2, "cm"), pad_y = grid::unit(0.2, "cm"),
style = ggspatial::north_arrow_minimal
) +
coord_sf(datum = NA, expand = TRUE) +
labs(x = NULL, y = NULL) +
theme_bw(base_size = 10, base_family = "sans") +
theme(
panel.background = element_rect(fill = "#FCFDFD", colour = NA),
panel.grid = element_blank(),
axis.text = element_blank(), axis.ticks = element_blank(),
legend.position = "right",
legend.text = element_text(colour = "#344B58", size = 9),
legend.key = element_rect(fill = "white", colour = NA),
panel.border = element_rect(colour = "#A7B7BE", fill = NA, linewidth = 0.4),
plot.margin = margin(6, 10, 6, 6)
)选择瓦片来源
本节用 ggspatial::annotation_map_tile()
添加在线底图：先放瓦片层，再叠加 geom_sf()
绘制的主题要素。type 选择瓦片服务与样式，zoom
控制街道细节，alpha 调整整张瓦片的透明度；地图范围由数据或
coord_sf() 决定。alpha
不会分别改变水体、公园等单个对象的颜色。
| 瓦片提供方 | 本节选择 | 使用前确认 |
|---|---|---|
| CARTO | 东京使用浅色 cartolight；上海使用深色
cartodark |
需 CARTO 底图密钥；标明 CARTO 与 OpenStreetMap 数据来源 |
| OpenStreetMap Deutschland（FOSSGIS） | 伦敦使用浅色街道瓦片 | 无需密钥；有请求限制，不能保证各地持续可达；标明瓦片与数据来源 |
伦敦示例先用 rosm::register_tile_source() 注册 HTTPS
瓦片地址，再在 type 中选择它。三个示例都用
cachedir
将瓦片暂存于系统临时目录，避免中文项目路径干扰缓存读写。
其他选择包括 tmap::tm_basemap()（适合在静态与交互地图之间切换）和
maptiles（获取并管理瓦片）。Esri、MapTiler
等也提供底图，但需分别核对接口、密钥及授权条款。保存、分发与大量请求瓦片应遵守所选服务的条款；程序支持缓存，不等于获得离线分发许可。
运行 CARTO 示例前：先在 CARTO 底图页面申请个人密钥，并更新
ggspatial至支持密钥的 1.1.11 或更新版本。在 R Console 运行file.edit("~/.Renviron")，添加一行CARTO_API_KEY=个人密钥，保存并重启 R。密钥只保存在个人电脑上，不写进.Rmd或共享脚本。伦敦示例使用的 OpenStreetMap Deutschland 瓦片无需此密钥。
国内课堂提示：境外瓦片与 OSM 查询的可达性会随网络变化；预存图片可供离线阅读，在线代码需先测试。若改用高德等国内服务，应核对密钥、署名及 GCJ-02 坐标系，不能将 WGS 84 的 OSM 点位直接叠加。
示例 1 · 东京咖啡馆：浅色瓦片
先用地名取得东京千代田区周边的外接矩形，再查询其中以点要素记录的咖啡馆。该矩形并非行政区边界，面要素记录的店铺也未纳入。浅色瓦片放在点位之下，保留街道参照。地名解析、OSM Overpass 与 CARTO 瓦片均需联网。
东京咖啡馆与浅色街道瓦片（历史快照）。
读图与署名：图片是历史快照；重新查询 OSM 后，点位与类别数量可能变化。图中的 “Independent Cafe” 只能解读为品牌未记录，“Other Chains” 表示其他已记录品牌。瓦片：© CARTO、© OpenStreetMap contributors；点位：© OpenStreetMap contributors。
查看并复制东京瓦片地图的 R 代码
library(tidyverse)
library(sf)
library(osmdata)
library(ggspatial)
# 1. 获取点位；地名解析返回矩形范围，不是行政区边界
chiyoda_bbox <- getbb("Chiyoda, Tokyo, Japan")
chiyoda_query <- opq(bbox = chiyoda_bbox) %>%
add_osm_feature(key = "amenity", value = "cafe")
chiyoda_cafes <- osmdata_sf(chiyoda_query)$osm_points
if (is.null(chiyoda_cafes) || nrow(chiyoda_cafes) == 0L) {
stop("未取得咖啡馆点位，请检查网络或 OSM 查询。")
}
# 2. 整理 OSM 品牌标签；缺失标签不能解释为独立经营
for (tag in c("brand", "brand:en")) {
if (!tag %in% names(chiyoda_cafes)) chiyoda_cafes[[tag]] <- NA_character_
}
chiyoda_cafes <- chiyoda_cafes %>%
mutate(brand_name = coalesce(
na_if(str_squish(brand), ""),
na_if(str_squish(`brand:en`), "")
))
top_brands <- chiyoda_cafes %>%
st_drop_geometry() %>%
filter(!is.na(brand_name)) %>%
count(brand_name, sort = TRUE) %>%
slice_head(n = 6) %>%
pull(brand_name)
chiyoda_cafes <- chiyoda_cafes %>%
mutate(
cafe_type = case_when(
is.na(brand_name) ~ "Brand not recorded",
brand_name %in% top_brands ~ brand_name,
TRUE ~ "Other recorded brands"
),
cafe_type = factor(cafe_type, levels = c(
top_brands, "Other recorded brands", "Brand not recorded"
))
)
# 3. 先绘瓦片，再叠加点；CARTO 密钥从个人环境变量读取
if (!nzchar(Sys.getenv("CARTO_API_KEY"))) {
stop("请先在个人 .Renviron 中设置 CARTO_API_KEY。")
}
tile_cache <- file.path(tempdir(), "chapter4_tiles")
dir.create(tile_cache, showWarnings = FALSE)
ggplot() +
annotation_map_tile(
type = "cartolight", zoom = 14, alpha = 0.90,
cachedir = tile_cache, quiet = TRUE
) +
geom_sf(data = chiyoda_cafes, aes(colour = cafe_type),
size = 1.1, alpha = 0.90) +
scale_colour_brewer(palette = "Set1", name = "OSM brand tag") +
guides(colour = guide_legend(override.aes = list(size = 2.8, alpha = 1))) +
annotation_scale(location = "bl", width_hint = 0.25,
style = "ticks", text_cex = 0.8) +
annotation_north_arrow(
location = "br", which_north = "true",
style = north_arrow_fancy_orienteering(
fill = c("gray40", "white"), line_col = "gray20"
)
) +
coord_sf(expand = FALSE) +
labs(
title = "Cafes in central Tokyo",
subtitle = "Chiyoda query area | OSM point features",
x = "Longitude", y = "Latitude",
caption = "Basemap: CARTO / OpenStreetMap contributors | POI: OpenStreetMap contributors"
) +
theme_bw(base_size = 10) +
theme(
legend.position = "right",
legend.title = element_text(face = "bold"),
legend.key.size = grid::unit(0.5, "cm"),
panel.grid.major = element_line(colour = "grey88", linetype = "dashed"),
panel.grid.minor = element_blank(),
plot.title = element_text(face = "bold", size = 13),
plot.subtitle = element_text(colour = "grey35")
)示例 2 · 上海便利店：深色瓦片
上海示例直接给出陆家嘴与外滩附近的经纬度边界框，查询其中以点要素记录的便利店。该范围约 4.8 × 4.4 公里，并非行政区界。图中的深色瓦片衬托浅色点位；图例数字表示各类点位数，不能解释为区域内全部便利店的数量。
上海便利店与深色街道瓦片（历史快照）。
读图与署名：图片是历史快照；图中的 “Independent Store” 表示品牌未记录，“Other Chains” 表示其他已记录品牌。重新查询会改变点位与类别数量。瓦片：© CARTO、© OpenStreetMap contributors；点位：© OpenStreetMap contributors。
查看并复制上海瓦片地图的 R 代码
library(tidyverse)
library(sf)
library(osmdata)
library(ggspatial)
# 1. 查询固定矩形内的便利店点；四个数依次为西、南、东、北
lujiazui_bbox <- c(121.47, 31.22, 121.52, 31.26)
shanghai_query <- opq(bbox = lujiazui_bbox) %>%
add_osm_feature(key = "shop", value = "convenience")
shanghai_stores <- osmdata_sf(shanghai_query)$osm_points
if (is.null(shanghai_stores) || nrow(shanghai_stores) == 0L) {
stop("未取得便利店点位，请检查网络或 OSM 查询。")
}
# 2. 品牌标签可能缺失；统计类别点数用于图例
for (tag in c("brand", "brand:en")) {
if (!tag %in% names(shanghai_stores)) shanghai_stores[[tag]] <- NA_character_
}
shanghai_stores <- shanghai_stores %>%
mutate(brand_name = coalesce(
na_if(str_squish(brand), ""),
na_if(str_squish(`brand:en`), "")
))
top_brands <- shanghai_stores %>%
st_drop_geometry() %>%
filter(!is.na(brand_name)) %>%
count(brand_name, sort = TRUE) %>%
slice_head(n = 3) %>%
pull(brand_name)
shanghai_stores <- shanghai_stores %>%
mutate(
store_type = case_when(
is.na(brand_name) ~ "Brand not recorded",
brand_name %in% top_brands ~ brand_name,
TRUE ~ "Other recorded brands"
),
store_type = factor(store_type, levels = c(
top_brands, "Other recorded brands", "Brand not recorded"
))
) %>%
add_count(store_type, name = "store_n") %>%
mutate(store_label = paste0(store_type, " (", store_n, ")"))
store_levels <- shanghai_stores %>%
st_drop_geometry() %>%
distinct(store_type, store_label) %>%
arrange(store_type) %>%
pull(store_label)
shanghai_stores <- shanghai_stores %>%
mutate(store_label = factor(store_label, levels = store_levels))
# 3. 深色瓦片位于点图层之下，比例尺和指北针用浅色
if (!nzchar(Sys.getenv("CARTO_API_KEY"))) {
stop("请先在个人 .Renviron 中设置 CARTO_API_KEY。")
}
tile_cache <- file.path(tempdir(), "chapter4_tiles")
dir.create(tile_cache, showWarnings = FALSE)
ggplot() +
annotation_map_tile(
type = "cartodark", zoom = 15, alpha = 0.90,
cachedir = tile_cache, quiet = TRUE
) +
geom_sf(data = shanghai_stores, aes(colour = store_label),
size = 2, alpha = 0.90) +
scale_colour_brewer(palette = "Set2", name = "OSM brand tag (points)") +
guides(colour = guide_legend(override.aes = list(size = 3, alpha = 1))) +
annotation_scale(location = "bl", width_hint = 0.25,
style = "ticks", text_cex = 0.8,
text_col = "grey90", line_col = "grey90") +
annotation_north_arrow(
location = "br", which_north = "true",
style = north_arrow_fancy_orienteering(
fill = c("grey50", "white"),
line_col = "grey90", text_col = "grey90"
)
) +
coord_sf(expand = FALSE) +
labs(
title = "Convenience stores in central Shanghai",
subtitle = "Lujiazui and the Bund query area | OSM point features",
x = "Longitude", y = "Latitude",
caption = "Basemap: CARTO / OpenStreetMap contributors | POI: OpenStreetMap contributors"
) +
theme_bw(base_size = 10) +
theme(
panel.background = element_rect(fill = "#25282D", colour = NA),
panel.grid.major = element_line(colour = "grey35", linetype = "dashed"),
panel.grid.minor = element_blank(),
legend.position = "right",
legend.title = element_text(face = "bold"),
legend.key.size = grid::unit(0.5, "cm"),
plot.title = element_text(face = "bold", size = 13),
plot.subtitle = element_text(colour = "grey35")
)示例 3 · 伦敦南岸：500 米缓冲区与浅色瓦片
以 Waterloo Bridge, South Bank
共享单车站为中心，在适合伦敦的米制投影 EPSG:27700
中建立 500
米缓冲区，再把邻近站点叠加到街道瓦片上。缓冲区表示地图上的直线距离，不能代表沿道路步行
500 米可到达的范围。spData 随包提供固定的历史 OSM
站点快照，点位数据无需在线查询；显示瓦片时仍需联网。
伦敦 Waterloo Bridge 共享单车站周边 500 米缓冲区与站点（spData 历史快照）。
读图与署名：此版
spData快照中，缓冲区内有 11 个其他站点。瓦片透明度统一设为0.56，以突出缓冲边界与点位；水体、公园没有单独调色。站点数据：© OpenStreetMap contributors（经spData整理）；瓦片：© OpenStreetMap contributors、© OpenStreetMap Deutschland / FOSSGIS。
查看并复制伦敦缓冲区瓦片地图的 R 代码
library(sf)
library(spData)
library(ggplot2)
library(ggspatial)
# 1. 注册可通过 HTTPS 访问的浅色瓦片源
osm_de <- rosm::as.tile_source(
"https://tile.openstreetmap.de/tiles/osmde/${z}/${x}/${y}.png",
attribution = "Map data © OpenStreetMap contributors; tiles © OpenStreetMap Deutschland / FOSSGIS"
)
rosm::register_tile_source(osmde = osm_de)
tile_cache <- file.path(tempdir(), "chapter4_tiles")
dir.create(tile_cache, showWarnings = FALSE)
# 2. 先转到米制投影，再建立 500 米缓冲区并筛选邻近站点
data("cycle_hire_osm", package = "spData")
stations_m <- st_transform(cycle_hire_osm, 27700)
centre_m <- stations_m[which(stations_m$name == "Waterloo Bridge, South Bank"), ]
if (nrow(centre_m) != 1L) stop("未找到唯一的 Waterloo Bridge 站点。")
buffer_500m <- st_buffer(centre_m, dist = 500)
nearby_m <- stations_m[lengths(st_intersects(stations_m, buffer_500m)) > 0L, ]
nearby_m$station_type <- ifelse(
nearby_m$osm_id == centre_m$osm_id,
"Centre station", "Other stations"
)
nearby_m$station_type <- factor(
nearby_m$station_type,
levels = c("Centre station", "Other stations")
)
map_bbox <- st_bbox(buffer_500m)
padding_m <- 300
# 3. 瓦片在最底层；半透明缓冲区与站点依次叠加
ggplot() +
annotation_map_tile(
type = "osmde", zoom = 15, alpha = 0.56,
cachedir = tile_cache, quiet = TRUE
) +
geom_sf(data = buffer_500m, fill = "#F5B975",
colour = "#C76143", alpha = 0.13,
linewidth = 0.8, linetype = "22") +
geom_sf(
data = nearby_m,
aes(fill = station_type, shape = station_type, size = station_type),
colour = "white", stroke = 0.65
) +
scale_fill_manual(NULL, values = c(
"Centre station" = "#BB495B", "Other stations" = "#207B89"
)) +
scale_shape_manual(NULL, values = c(
"Centre station" = 23, "Other stations" = 21
)) +
scale_size_manual(values = c(
"Centre station" = 4.2, "Other stations" = 2.7
), guide = "none") +
guides(fill = guide_legend(override.aes = list(size = 4.4))) +
annotation_scale(location = "bl", width_hint = 0.20, style = "ticks") +
annotation_north_arrow(
location = "tr", which_north = "true",
style = north_arrow_minimal,
height = grid::unit(0.7, "cm"), width = grid::unit(0.7, "cm")
) +
coord_sf(
crs = st_crs(27700),
xlim = c(map_bbox["xmin"] - padding_m, map_bbox["xmax"] + padding_m),
ylim = c(map_bbox["ymin"] - padding_m, map_bbox["ymax"] + padding_m),
expand = FALSE, datum = NA
) +
labs(
title = "Cycle hire within 500 m of Waterloo Bridge",
subtitle = "London South Bank | Straight-line buffer",
caption = "Map data: OpenStreetMap contributors | Tiles: OSM.de / FOSSGIS"
) +
theme_minimal(base_size = 10) +
theme(
panel.grid = element_blank(),
axis.title = element_blank(), axis.text = element_blank(),
legend.position = "right",
legend.text = element_text(size = 9, colour = "#344B58"),
legend.key = element_rect(fill = "white", colour = NA),
plot.title = element_text(face = "bold", size = 13),
plot.subtitle = element_text(colour = "#52636E"),
plot.caption = element_text(colour = "#52636E", size = 7.5)
)交互式地图
4.5 交互式地图（概览）
拓展 · 进阶
交互式地图适合网页报告与数据探索。读者可以缩放、漫游，并通过悬停或点击查看具体对象。本节用两张地图分别展示区域属性与地点连接：曼哈顿普查区的家庭收入，以及纽约肯尼迪机场（JFK）至美国本土主要目的机场的航班联系。两张地图都先呈现结果，示例代码可在图下展开。
第一张地图沿用分层设色的“指标 → 分级区间 → 颜色”映射；第二张地图把汇总后的航班数映射为连线宽度。交互功能补充具体数值，不改变地图变量与符号的含义。
本节以 leaflet 作为示例框架：
leaflet提供基于 HTML 的交互地图容器，可直接在 R Markdown knit 后输出为可交互页面；
sf管理边界与连线几何，供 Leaflet 使用的坐标统一为EPSG:4326（经纬度）；addPolygons()绘制面，addPolylines()绘制连接；label、popup分别提供悬停标签与点击详情；addTiles()在数据图层下方加载 OpenStreetMap Deutschland 瓦片。服务无需 API 密钥，但打开网页时仍需联网；两张地图均保留瓦片署名。
示例 1 · 曼哈顿家庭收入
地图沿用 4.2 的 ACS 本地快照，只显示曼哈顿，以控制网页负载。桌面端悬停可查看普查区标签；点击或轻触可查看估计值及误差范围。灰色表示收入估计缺失。
指标口径 家庭收入中位数来自 ACS 2020—2024 年五年估计，以 2024 年通胀调整后的美元表示。弹窗中的 MOE 是 90% 置信水平下的误差范围；收入缺失的普查区单独用灰色表示。
查看并复制曼哈顿交互地图的 R 代码
# 交互式分层设色图：曼哈顿家庭收入中位数
library(sf)
library(tidyverse)
library(leaflet)
library(classInt)
library(RColorBrewer)
library(htmltools)
# 1. 与 4.2 共用固定快照，仅取曼哈顿（GEOID 中的县码 061）
nyc_acs <- st_read(
"data/acs-nyc/nyc_acs_tracts_2020_2024.gpkg",
layer = "nyc_acs_tracts", quiet = TRUE
)
acs_period <- unique(nyc_acs$acs_period)
inc_tract <- nyc_acs %>%
filter(substr(GEOID, 3, 5) == "061") %>%
transmute(
GEOID,
NAME,
med_hh_income,
moe = income_moe
) %>%
st_transform(4326)
n_class <- 7
# ------------------------------------------------------------
# 3) Fisher 分级（classInt）
# - 仅用有效值生成断点，避免 NA/Inf 干扰
# ------------------------------------------------------------
x_ok <- inc_tract$med_hh_income
x_ok <- x_ok[is.finite(x_ok)]
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 = "#8C8C8C"
)
# ------------------------------------------------------------
# 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 (90%): ",
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)] / 1000, format = "f", digits = 1),
"–",
formatC(brks[-1] / 1000, format = "f", digits = 1)
)
legend_labels_rev <- c(rev(legend_labels), "无可用估计")
legend_colors_rev <- c(rev(pal_vec), "#8C8C8C")
# ------------------------------------------------------------
# 7) leaflet 交互地图
# ------------------------------------------------------------
m <- leaflet(inc_tract2, width = "100%", height = 520) %>%
# 先加载街道底图，再叠加半透明普查区；瓦片需在浏览器中联网获取
addTiles(
urlTemplate = "https://tile.openstreetmap.de/{z}/{x}/{y}.png",
attribution = paste0(
'© <a href="https://www.openstreetmap.org/copyright">',
'OpenStreetMap contributors</a> · ',
'<a href="https://openstreetmap.de/">OpenStreetMap Deutschland</a>'
)
) %>%
addPolygons(
fillColor = pal(inc_tract2$med_hh_income),
fillOpacity = 0.70,
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,
popupOptions = popupOptions(maxWidth = 220, autoPan = TRUE),
# 高亮（hover 时描边变亮）
highlightOptions = highlightOptions(
weight = 1.0,
color = "#ffd54f",
bringToFront = TRUE
)
) %>%
addLegend(
colors = legend_colors_rev,
labels = legend_labels_rev,
title = paste0("Median HH income<br>(USD, thousands)<br>ACS ", acs_period),
opacity = 0.9,
position = "bottomright"
)
# 按所选普查区设定初始视窗
map_bbox <- st_bbox(inc_tract2)
m <- m %>% fitBounds(map_bbox[["xmin"]], map_bbox[["ymin"]],
map_bbox[["xmax"]], map_bbox[["ymax"]])
# 在 R Script 里务必显式打印对象，确保 Viewer 渲染
m示例 2 · JFK 出发航线
nycflights23 收录 2023
年纽约地区机场的航班记录。本例筛选 JFK
有实际起飞时间的航班，汇总至美国本土班次最多的 20
个目的机场。连线宽度与终点圆点大小表示航班记录次数；悬停或点击可查看机场代码与具体班次。
运行前：若尚未安装数据包，先在 R Console 执行
install.packages("nycflights23")。航班数据已随包提供，无需在 Knit 时在线下载。
读图口径：图中 20 条连接共对应 78,837 班实际起飞记录，表示 JFK 与目的机场之间的联系，不代表乘客人数。两机场之间的直线仅为示意连接，并非真实飞行轨迹。航班记录来自美国交通统计局（经
nycflights23整理），机场坐标源于 OpenFlights。
查看并复制 JFK 航线地图的 R 代码
# 2023 年 JFK 出发至美国本土主要目的机场的交互式航线图
library(dplyr)
library(sf)
library(leaflet)
# 1. 仅统计有实际起飞时间的航班；每个目的机场汇总为一条连接。
airports <- nycflights23::airports %>%
select(faa, lat, lon)
jfk <- airports %>% filter(faa == "JFK")
routes <- nycflights23::flights %>%
filter(origin == "JFK", !is.na(dep_time)) %>%
count(dest, name = "flight_count") %>%
inner_join(airports, by = c("dest" = "faa")) %>%
filter(between(lon, -125, -66), between(lat, 24, 50)) %>%
slice_max(flight_count, n = 20, with_ties = FALSE) %>%
mutate(
line_weight = scales::rescale(flight_count, to = c(2.2, 5.8)),
point_radius = scales::rescale(flight_count, to = c(4.5, 9)),
tooltip = paste0("JFK → ", dest, "：", format(flight_count, big.mark = ","), " 班"),
popup = paste0(
"<strong>JFK → ", dest, "</strong><br>",
"2023 年实际起飞：", format(flight_count, big.mark = ","), " 班"
)
)
# 可用于核对图面范围：20 个目的机场、78,837 班实际起飞航班。
route_summary <- routes %>%
summarise(destinations = n(), mapped_flights = sum(flight_count))
# 2. 两机场坐标组成一条示意直线；线宽对应实际起飞班次。
route_geometry <- st_sfc(
lapply(seq_len(nrow(routes)), function(i) {
st_linestring(rbind(
c(jfk$lon, jfk$lat),
c(routes$lon[i], routes$lat[i])
))
}),
crs = 4326
)
route_lines <- st_sf(routes, geometry = route_geometry)
legend_html <- paste0(
'<div style="background:#ffffffee;color:#243642;padding:9px 12px;',
'border-radius:7px;box-shadow:0 1px 8px #0002;font-size:12px;line-height:1.6">',
'<strong>JFK 出发 · 2023</strong><br>',
'<span style="color:#C65B3C">●</span> JFK 起点 ',
'<span style="color:#155A78">●</span> 目的机场<br>',
'<span style="display:inline-block;width:20px;border-top:3px solid #276782;',
'vertical-align:middle"></span> 线越粗，航班越多<br>',
'前 20 个目的机场 · ', format(route_summary$mapped_flights, big.mark = ","), ' 班',
'</div>'
)
# 3. 底图、航线和机场依次叠加；在线打开页面时加载瓦片。
route_map <- leaflet(
width = "100%", height = 500,
options = leafletOptions(minZoom = 3, maxZoom = 9)
) %>%
addTiles(
urlTemplate = "https://tile.openstreetmap.de/{z}/{x}/{y}.png",
attribution = paste0(
'© <a href="https://www.openstreetmap.org/copyright">',
'OpenStreetMap contributors</a> · ',
'<a href="https://openstreetmap.de/">OpenStreetMap Deutschland</a>'
),
options = tileOptions(opacity = 0.62)
) %>%
addPolylines(
data = route_lines,
color = "#276782", weight = ~line_weight, opacity = 0.72,
label = ~tooltip, popup = ~popup,
highlightOptions = highlightOptions(color = "#C65B3C", weight = 6)
) %>%
addCircleMarkers(
data = routes, lng = ~lon, lat = ~lat,
radius = ~point_radius,
color = "#FFFFFF", weight = 1.2,
fillColor = "#155A78", fillOpacity = 0.92,
label = ~tooltip, popup = ~popup
) %>%
addCircleMarkers(
data = jfk, lng = ~lon, lat = ~lat,
radius = 10, color = "#FFFFFF", weight = 2,
fillColor = "#C65B3C", fillOpacity = 1,
label = "JFK · 纽约起点",
popup = "<strong>JFK · 纽约起点</strong><br>图中连线仅为前 20 个本土目的机场。"
) %>%
addControl(legend_html, position = "bottomleft") %>%
# 按所有机场自动确定初始视窗，窄屏也能看到起点和目的地。
fitBounds(
lng1 = min(c(jfk$lon, routes$lon)) - 1.5,
lat1 = min(c(jfk$lat, routes$lat)) - 1.5,
lng2 = max(c(jfk$lon, routes$lon)) + 1.5,
lat2 = max(c(jfk$lat, routes$lat)) + 1.5
)
route_map进阶提示
leaflet的地图由浏览器渲染，交互依赖 HTML、CSS 与 JavaScript。两例分别使用addPolygons()和addPolylines()；需要发布为静态图片时，应另行检查图例、标注与数据来源是否能脱离悬停信息独立阅读。
5.0 制图规范与美化
5.0 制图规范与美化
5.0 规范与美化指南
制图规范
5.1.1 通用地图出图规范
地图是对空间证据的表达。出图前可依次检查范围与投影、变量与图例、方向与距离、来源与注记。辅助元素按读图需要选择，不必每张地图都堆齐。
读图顺序 先看“在哪里、是什么数据”，再看“颜色和符号表示什么”，最后看“这些差异能否比较”。
- 范围与投影：研究区边界应与问题及数据覆盖范围一致。需要交代位置时，可加入定位小图（locator map）；涉及距离或面积计算时，应选用适合研究区的投影。闵行区研究可用“闵行区—上海市—中国”的位置关系帮助读者定位。
- 变量与图例：写清变量、单位和符号含义。分级设色图宜在图例或题注交代分级方法与阈值；同一指标跨图比较时，尽量使用相同阈值与色标范围。竖式顺序图例可按从下到上由低到高排列。
- 缺失值：缺失不是零值。用有别于数值色带的颜色或纹理表示缺失区域，并在图例中说明，避免误读为低值或背景。
- 方向与距离：局部地图若需要估计距离，宜加入与投影及图幅相适应的比例尺。全球图或变形随位置变化较大的地图，不宜用一根比例尺暗示全图距离处处相同。指北针、经纬网按需要选用，地图旋转时尤其应让方向可判断。
- 来源与注记：在图面或题注说明专题数据的来源、年份和必要的统计口径；使用外部瓦片时，按服务条款保留底图署名。优先标出解释地图所必需的地名和数值，避免遮挡主题图层。
- 跨图一致性：同一组地图尽量统一字体、色标逻辑、边界线和缺失值表示。指标不同、统计口径不同的地图，即使使用同一色系，也不能仅凭颜色直接比较数值。
题注与图面分工：题注已交代的主题、数据来源和年份，通常无需在图内重复；图内保留读图必需的图例、标签与注记。
5.1.2 中国地图的规范表达
中国地图的制图检查可沿着底图来源、图面内容、使用方式三个问题展开。底图的审图号、专题数据的来源和新制作地图的审核状态分别记录，不能相互替代。
检查顺序 先确认底图与审图号 → 再核对版图、附图和资料缺失 → 最后按课堂练习、直接使用或公开发表的用途检查要求。
01 · 底图有来源
优先从自然资源部标准地图服务系统查找适合图幅的标准地图，保存图名、审图号和下载日期。第三方 GeoJSON 可用于练习空间连接与排版；只有文件名或网页说明，尚不足以证明其地图内容已经审核。
02 · 内容按图幅核对
中国全图应完整表达大陆、海南岛、台湾岛、南海诸岛，以及按比例尺和图幅应表示的重要岛屿。南海诸岛可以放在主图，也可以用附图表示；将附图安排在右下角通常便于阅读，但仍应以不遮挡主图、内容清楚为准。公开图还应避免标出依法不应公开的敏感信息。
03 · 缺失值明示
按行政区设色时，缺失值使用独立于数值色带的颜色。台湾省没有同口径资料时，应在地图适当位置注明“台湾省资料暂缺”。香港、澳门等地区若未纳入统计来源，也须说明，不能当作低值填色。
04 · 坐标与读图要素
EPSG:4490 表示 CGCS2000
的地理坐标系；应先核对源图层实际使用的 CRS，再决定是否调用
st_transform()。坐标转换不能证明边界来源或地图已通过审核。比较面积时可选适合研究区的等面积投影。比例尺、经纬网和指北针按阅读需要配置；大范围地图上的单一比例尺尤其要说明适用位置。
附图排版建议：横向图幅可在主图右下侧留出南海诸岛附图的位置。正式地图还须对照《公开地图内容表示规范》核查附图范围、重要岛屿和界线表示，确保附图清晰且不遮挡主图。
示例：2020 年县级常住人口
下面用第七次全国人口普查年鉴的整理数据演示县级人口与行政区边界的连接。地图叠加县、市、省三级边界，以五级
mako
色系表示人口从少到多；主图与南海诸岛附图共用分级和颜色，岛礁直接按县级面着色。灰色在图例中单列，表示没有可直接对应的县级数值，并不代表人口为零。
浅蓝灰背景表示海域。上海附近的长江口和杭州湾是水域，沿海零散的小岛来自原始边界，并非人口连接后漏绘的地块；县级细线在全国比例尺下适当减淡，以便辨认海岸。
分级断点用 Fisher–Jenks 方法在
log10(人口)
上计算，再以原始人数（万人）标注图例。图例阈值保留两位小数，实际分组使用未取整的断点。这样，人口极高的少数县级单位不会把其余大部分县挤在同一档；区间左开右闭，最低档包含最小值。
教学预览：数据年份与审图号
人口为 2020 年数据，边界来自本地 2025 年行政区划数据包；图上的 GS（2024）0650号标识的是边界底图，并非叠加人口后这张专题图自身的审图号。源 SHP 以 WGS84（EPSG:4326）存储，绘图时在 WGS84 基准的 Albers 等面积投影下显示。
核对更名、改码与撤并后，2,844/2,891 个县级边界接上人口数；余下 47 个保留灰色，主要涉及港澳台资料未覆盖及跨年区划无法稳妥对齐。该图供课堂学习与预览，公开发表时应按实际用途核查地图审核要求。
复现用的县级人口表与三沙人口分项表可在 DAWN 课程数据页下载。两个
CSV 均放入本章 data/
文件夹；边界数据需另行按来源授权取得，并将包含省、市、县、国界线和十段线
SHP 的文件夹路径设为
CHINA_DIVISION_DIR。数据页目前提供人口表，不提供这套边界文件。
2020 年县级常住人口五级分层设色（教学预览）。
查看 R 代码：数据连接、设色与南海附图
代码读取省、市、县边界与两张人口表，核对县级代码及跨年区划变化，再分别生成主图和南海附图。运行前需安装脚本加载的
R 包，并将本地边界数据包路径设为
CHINA_DIVISION_DIR。完整脚本见 data/china_county_population_2020_local.R。
# 第 4 章 5.1.2 本地制图预览：2020 年县级常住人口。
# 2020 年人口来自 Lei Dong 团队整理的第七次人口普查年鉴数据；
# 省、市、县边界来自本地 SHP；边界底图审图号：GS（2024）0650号。
# 本脚本生成教学预览；新增人口专题内容后的图面需按实际用途核对审核要求。
Sys.setlocale("LC_CTYPE", ".UTF-8")
library(cowplot)
library(classInt)
library(dplyr)
library(ggplot2)
library(readr)
library(sf)
library(viridisLite)
shape_dir <- Sys.getenv("CHINA_DIVISION_DIR")
if (!dir.exists(shape_dir)) {
stop("请先将边界 SHP 文件夹路径设为 CHINA_DIVISION_DIR。")
}
read_shape <- function(stem) {
st_read(file.path(shape_dir, paste0(stem, ".shp")), quiet = TRUE)
}
# 1. 几何、代码与普查表只读取一次。
county <- read_shape("中国_县")
city <- read_shape("中国_市")
province <- read_shape("中国_省")
border_raw <- read_shape("国界线")
border <- border_raw %>% filter(name == "国界线-已定国界")
border_undefined <- border_raw %>% filter(name == "国界线-未定国界")
ten_dash <- read_shape("十段线")
layers <- list(county, city, province, border_raw, ten_dash)
if (any(vapply(layers, function(x) {
is.na(st_crs(x)) || is.na(st_crs(x)$epsg) || st_crs(x)$epsg != 4326
}, logical(1)))) {
stop("图层 CRS 与预期的 EPSG:4326 不一致，请检查源文件。")
}
census <- read_csv(
"data/china_county_population_2020.csv",
col_types = cols(
county_code = col_character(),
province = col_character(),
city = col_character(),
county_name = col_character(),
population_2020 = col_double()
),
show_col_types = FALSE
) %>%
select(county_code, census_name = county_name, population_2020)
county_key <- county %>%
st_drop_geometry() %>%
transmute(county_code = substr(as.character(gb), 4, 9), map_name = name)
if (anyDuplicated(county_key$county_code) || anyDuplicated(census$county_code)) {
stop("县级代码有重复，不能执行一对一连接。")
}
# 2. 默认要求六位县码与名称同时一致；区划更名和撤并只使用下列逐条核实的映射。
# 原始人口来自 data/census_county_2020_v1.xlsx 的“七普（年鉴）”工作表。
# 同码异名 8 项与一对一改码 12 项均显式列出旧名，防止“南沙区”等重名误连。
# 撤并依据：福建省政府闽政文〔2021〕49号、洛阳市政府历史沿革。
# 三沙依据：民政部 2020-04-18 公告；西沙区管辖西沙、代管中沙，南沙区管辖南沙。
# https://www.fujian.gov.cn/zwgk/zxwj/szfwj/202102/t20210203_5530849.htm
# https://www.ly.gov.cn/mlly/lsyg/
# https://mzj.panjin.gov.cn/2020_04/20_14/content-63560.html
exact_population <- county_key %>%
inner_join(census, by = "county_code") %>%
filter(map_name == census_name) %>%
transmute(county_code, population_2020, match_type = "代码及名称一致")
renamed_units <- tibble::tribble(
~county_code, ~map_name, ~source_code, ~source_name, ~match_type,
"632301", "同仁市", "632301", "同仁县", "同码异名",
"520204", "水城区", "520204", "水城县", "同码异名",
"421088", "监利市", "421088", "监利县", "同码异名",
"540422", "米林市", "540422", "米林县", "同码异名",
"320614", "海门区", "320614", "海门市", "同码异名",
"340212", "繁昌区", "340212", "繁昌县", "同码异名",
"340210", "湾沚区", "340210", "芜湖县", "同码异名",
"540530", "错那市", "540530", "错那县", "同码异名",
"520581", "黔西市", "520522", "黔西县", "一对一改码",
"450181", "横州市", "450127", "横县", "一对一改码",
"513402", "会理市", "513425", "会理县", "一对一改码",
"610981", "旬阳市", "610928", "旬阳县", "一对一改码",
"350604", "龙海区", "350681", "龙海市", "一对一改码",
"350605", "长泰区", "350625", "长泰县", "一对一改码",
"410307", "偃师区", "410381", "偃师市", "一对一改码",
"654203", "沙湾市", "654223", "沙湾县", "一对一改码",
"532302", "禄丰市", "532331", "禄丰县", "一对一改码",
"610305", "凤翔区", "610322", "凤翔县", "一对一改码",
"431181", "祁阳市", "431121", "祁阳县", "一对一改码",
"350405", "沙县区", "350427", "沙县", "一对一改码"
)
if (nrow(dplyr::anti_join(renamed_units, county_key,
by = c("county_code", "map_name"))) != 0) {
stop("更名映射与当前县界代码或名称不符。")
}
renamed_population <- renamed_units %>%
left_join(census %>%
transmute(source_code = county_code,
observed_name = census_name, population_2020),
by = "source_code")
if (any(is.na(renamed_population$population_2020)) ||
any(renamed_population$source_name != renamed_population$observed_name)) {
stop("更名映射与 2020 人口表的代码或名称不符。")
}
renamed_population <- renamed_population %>%
select(county_code, population_2020, match_type)
# 两个新设区是完整撤并：分别合计 2020 年前身单元的人口。
merged_parts <- tibble::tribble(
~county_code, ~map_name, ~source_code, ~source_name,
"410308", "孟津区", "410322", "孟津县",
"410308", "孟津区", "410306", "吉利区",
"350404", "三元区", "350402", "梅列区",
"350404", "三元区", "350403", "三元区"
)
if (nrow(dplyr::anti_join(merged_parts, county_key,
by = c("county_code", "map_name"))) != 0) {
stop("撤并映射与当前县界代码或名称不符。")
}
merged_parts <- merged_parts %>%
left_join(census %>%
transmute(source_code = county_code,
observed_name = census_name, population_2020),
by = "source_code")
if (any(is.na(merged_parts$population_2020)) ||
any(merged_parts$source_name != merged_parts$observed_name)) {
stop("撤并映射与 2020 人口表的代码或名称不符。")
}
merged_population <- merged_parts %>%
group_by(county_code) %>%
summarise(population_2020 = sum(population_2020),
match_type = "撤并单元合计", .groups = "drop")
# 三沙三个分项在原始年鉴表中没有六位区划码，因此另存为核对表。
# 西沙 1,549 + 中沙 145 = 西沙区 1,694；南沙 639 = 南沙区；合计须等于三沙市 2,333。
sansa_raw <- read_csv(
"data/china_sansha_population_2020.csv",
col_types = cols(
record_type = col_character(), source_name = col_character(),
source_code = col_character(), population_2020 = col_double()
),
show_col_types = FALSE
)
sansa_parts <- sansa_raw %>%
filter(record_type == "群岛分项") %>%
select(source_name, source_code, population_2020)
sansa_total <- sansa_raw %>%
filter(record_type == "市级合计", source_code == "460300",
source_name == "三沙市") %>%
pull(population_2020)
if (nrow(sansa_parts) != 3 || any(!is.na(sansa_parts$source_code)) ||
any(is.na(sansa_parts$population_2020)) ||
length(sansa_total) != 1 || sum(sansa_parts$population_2020) != sansa_total) {
stop("三沙人口分项与原始年鉴表的三沙市总数不一致。")
}
sansa_mapping <- tibble::tribble(
~county_code, ~map_name, ~source_name,
"460321", "西沙区", "西沙群岛",
"460321", "西沙区", "中沙群岛的岛礁及其海域",
"460322", "南沙区", "南沙群岛"
)
if (nrow(dplyr::anti_join(sansa_mapping, county_key,
by = c("county_code", "map_name"))) != 0) {
stop("三沙映射与当前县界代码或名称不符。")
}
sansa_population_parts <- sansa_mapping %>%
left_join(sansa_parts, by = "source_name")
if (any(is.na(sansa_population_parts$population_2020))) {
stop("三沙映射缺少原始年鉴表人口分项。")
}
sansa_population <- sansa_population_parts %>%
group_by(county_code) %>%
summarise(population_2020 = sum(population_2020),
match_type = "三沙分项合计", .groups = "drop")
population_lookup <- bind_rows(
exact_population, renamed_population, merged_population, sansa_population
)
if (anyDuplicated(population_lookup$county_code) ||
any(is.na(population_lookup$population_2020)) ||
any(population_lookup$population_2020 <= 0)) {
stop("人口映射重复、缺失或含非正数，不能用于对数分级。")
}
invalid_count <- sum(!st_is_valid(county))
county <- county %>%
mutate(county_code = substr(as.character(gb), 4, 9)) %>%
left_join(population_lookup, by = "county_code") %>%
mutate(geometry = st_make_valid(geometry))
mainland <- !substr(county$county_code, 1, 2) %in% c("71", "81", "82")
cat("县界总数：", nrow(county), "；非港澳台代码县界：", sum(mainland),
"；有值：", sum(mainland & !is.na(county$population_2020)),
"；修复前无效几何：", invalid_count, "\n")
print(county %>% st_drop_geometry() %>% filter(mainland) %>%
count(match_type, name = "县界数"))
cat("仍缺人口的非港澳台代码县界：\n")
print(county %>% st_drop_geometry() %>%
filter(mainland & is.na(population_2020)) %>%
select(county_code, name), row.names = FALSE)
# 3. 在 log10(人口) 上求 Fisher–Jenks 断点，再按原始人数分组。
# 少数人口极高的县会拉开原值分级；对数分级让地图主体更易辨认。
population_values <- county$population_2020[!is.na(county$population_2020)]
class_count <- 5L
log_breaks <- classInt::classIntervals(
log10(population_values), n = class_count, style = "fisher",
largeN = Inf # 固定使用全部有值县，不进行抽样
)$brks
population_breaks <- 10^log_breaks
population_breaks[c(1, class_count + 1L)] <- range(population_values)
if (any(diff(population_breaks) <= 0)) {
stop("人口分级断点没有严格递增，请检查数据。")
}
# 实际分组使用计算所得断点；图例仅把人数格式化为两位“万人”。
inner_wan <- population_breaks[2:class_count] / 10000
class_labels <- c(
sprintf("≤%.2f 万", inner_wan[1]),
vapply(seq_len(class_count - 2L), function(i) {
sprintf(">%.2f–%.2f 万", inner_wan[i], inner_wan[i + 1L])
}, character(1)),
sprintf(">%.2f 万", tail(inner_wan, 1))
)
missing_label <- "无同口径数据"
class_levels <- c(class_labels, missing_label)
population_class <- as.character(cut(
county$population_2020, breaks = population_breaks,
labels = class_labels, include.lowest = TRUE, right = TRUE
))
population_class[is.na(population_class)] <- missing_label
county$population_class <- factor(population_class, levels = class_levels)
print(table(county$population_class))
# 仍沿用 mako 色系，只取较亮的五级；缺失值单列为浅灰色。
missing_fill <- "#BBC3C7"
class_colours <- c(
setNames(viridisLite::mako(class_count, begin = 0.42, end = 0.95,
direction = -1), class_labels),
setNames(missing_fill, missing_label)
)
# 4. 县/市/省界由细到粗；主图和南海附图共用类别及颜色。
map_crs <- st_crs(paste(
"+proj=aea +lat_1=25 +lat_2=47 +lat_0=0 +lon_0=105",
"+datum=WGS84 +units=m +no_defs"
))
font_family <- "sans" # Windows 绘图设备会为中文选择系统回退字体
map_theme <- theme_minimal(base_size = 9, base_family = font_family) +
theme(
panel.grid.major = element_line(colour = "#D5E3E9", linewidth = 0.18),
panel.grid.minor = element_blank(),
panel.background = element_rect(fill = "#EAF3F6", colour = NA),
plot.background = element_rect(fill = "white", colour = NA),
axis.title = element_blank(),
axis.text = element_text(colour = "#53646B", size = 7.7),
legend.title = element_text(colour = "#23363E", face = "bold", size = 10),
legend.text = element_text(colour = "#485B64", size = 9),
legend.key.spacing.y = unit(1.3, "mm"),
plot.margin = margin(4, 4, 4, 4)
)
boundary_layers <- list(
geom_sf(data = city, fill = NA, colour = "#5F7077", linewidth = 0.10,
alpha = 0.50, show.legend = FALSE),
geom_sf(data = province, fill = NA, colour = "#2D4049", linewidth = 0.26,
show.legend = FALSE),
geom_sf(data = border, colour = "#21343D", linewidth = 0.48,
show.legend = FALSE),
geom_sf(data = border_undefined, colour = "#21343D", linewidth = 0.48,
linetype = "dashed", show.legend = FALSE),
geom_sf(data = ten_dash, colour = "#21343D", linewidth = 0.52,
show.legend = FALSE)
)
main_map <- ggplot() +
geom_sf(data = county, aes(fill = population_class),
colour = "#8A989D80", linewidth = 0.035, show.legend = TRUE) +
boundary_layers +
annotate("label", x = 128.2, y = 24.0,
label = "台湾省\n资料暂缺", family = font_family,
colour = "#4B5A60", fill = "white", size = 2.7,
linewidth = 0.18, label.padding = unit(0.10, "lines")) +
scale_fill_manual(
values = class_colours, breaks = class_levels, drop = FALSE,
name = "县级常住人口\n2020 年",
guide = guide_legend(
keyheight = unit(5.1, "mm"), keywidth = unit(5.1, "mm"),
override.aes = list(colour = "#71818A", linewidth = 0.22)
)
) +
coord_sf(crs = map_crs, default_crs = st_crs(4326),
xlim = c(73, 136), ylim = c(16, 57), expand = FALSE) +
map_theme
map_legend <- cowplot::get_legend(main_map)
south_map <- ggplot() +
geom_sf(data = county, aes(fill = population_class),
colour = "#8A989D80", linewidth = 0.055) +
geom_sf(data = province, fill = NA, colour = "#60747D", linewidth = 0.20) +
geom_sf(data = ten_dash, colour = "#25414E", linewidth = 0.52) +
scale_fill_manual(values = class_colours, breaks = class_levels,
drop = FALSE, guide = "none") +
coord_sf(crs = map_crs, default_crs = st_crs(4326),
xlim = c(105, 125), ylim = c(2, 26), expand = FALSE) +
labs(title = "南海诸岛", subtitle = "岛礁与主图共用人口分级") +
map_theme +
theme(
plot.title = element_text(face = "bold", size = 8.5, colour = "#23363E"),
plot.subtitle = element_text(size = 7.2, colour = "#50616A"),
panel.border = element_rect(fill = NA, colour = "#84969F", linewidth = 0.38),
axis.text = element_blank(), axis.ticks = element_blank(),
panel.grid = element_blank(),
plot.margin = margin(3, 3, 3, 3)
)
preview <- ggdraw() +
draw_plot(main_map + theme(legend.position = "none"),
x = 0.005, y = 0.075, width = 0.78, height = 0.91) +
draw_grob(map_legend, x = 0.79, y = 0.55, width = 0.20, height = 0.40) +
draw_plot(south_map, x = 0.79, y = 0.13, width = 0.19, height = 0.38) +
draw_label(
"人口：第七次人口普查（2020） · 边界底图审图号：GS（2024）0650号",
x = 0.022, y = 0.030, hjust = 0, size = 8.1,
fontfamily = font_family, colour = "#52646C"
)
ggsave("images/china_county_population_2020_local.png", preview,
width = 10.2, height = 6.5, dpi = 300, bg = "white",
device = ragg::agg_png)另一种底图：官方在线瓦片
天地图提供可缩放的在线底图。它适合查看道路、地名和周边地理背景；瓦片署名与底图审图号应按所用服务标注。底图本身的审图号不等于后来叠加的专题图审图号。下例只加载矢量瓦片及注记，不叠加未经核验的省级界线。
运行前准备：从天地图官网进入开发服务，申请浏览器端 Key，并在本机 R Console 设置
TIANDITU_KEY环境变量。个人 Key 不写入 Rmd 或公开仓库；发布交互网页时，Key 仍会随瓦片请求在浏览器中可见，应使用专门的、受服务方允许且限制来源的 Key。
查看 R 代码：天地图浅色瓦片底图
library(leaflet)
library(magrittr)
tile_key <- Sys.getenv("TIANDITU_KEY")
if (!nzchar(tile_key)) {
stop("请先在本机设置 TIANDITU_KEY，再运行瓦片示例。")
}
# vec_w 为浅色矢量底图；cva_w 为地名注记。
tile_url <- function(layer) {
paste0(
"https://t{s}.tianditu.gov.cn/", layer,
"/wmts?SERVICE=WMTS&REQUEST=GetTile&VERSION=1.0.0",
"&LAYER=", substr(layer, 1, 3),
"&STYLE=default&TILEMATRIXSET=w&FORMAT=tiles",
"&TILEMATRIX={z}&TILEROW={y}&TILECOL={x}&tk=", tile_key
)
}
leaflet(options = leafletOptions(zoomControl = TRUE)) %>%
addTiles(
tile_url("vec_w"),
options = tileOptions(subdomains = as.character(0:7), opacity = 0.68),
attribution = "底图 © 天地图"
) %>%
addTiles(
tile_url("cva_w"),
options = tileOptions(subdomains = as.character(0:7)),
attribution = "注记 © 天地图"
) %>%
setView(lng = 104, lat = 35, zoom = 4) %>%
addScaleBar(position = "bottomleft", options = scaleBarOptions(imperial = FALSE))延伸阅读：规范原文与参考地图
- 《公开地图内容表示规范》（2023）列出中国全图、南海诸岛、台湾省等内容的表示要求。
- 《地图审核管理规定》说明地图审核范围与免审情形。
- 自然资源部标准地图服务系统提供带审图号的标准地图；下载时应核对所选图幅的实际审图号。
- 《地理科学进展》的一幅中国专题图在图注中列出所用标准地图的审图号和底图处理说明，可作为来源写法的参照；请在期刊网页查看原图。
地图美化
5.2 地图美化
地图的视觉设计应帮助读者先看到研究变量，再看到必要的地理背景。颜色、线宽和版式可调整；范围、数值、缺失值与边界含义不能因美化而变得含糊。
先后顺序 确认数据与地图内容正确 → 安排视觉层级 → 在最终网页尺寸检查图例、注记和局部细节。
视觉层级。 专题数据通常比路网、行政边界和地形底色更醒目。底图可选低饱和度色系，边界线保持可辨但不过重；当边界本身是研究对象时，应提高它的视觉权重。图例与比例尺的位置由图幅留白和读图路径决定，无须固定在某一角落。
色盘选择。 表示数量或强度时可用由浅到深的顺序色盘；围绕零值或其他有意义的参照值比较增减时，可用发散色盘；区分无序类别时，应选定性色盘。缺失值另设颜色，并检查图例与底图的对比度。色盘类型、R 配色包与预览方法可回看第三章 5.0：可视化美化与配色选择。
一致性。 同一研究中的地图沿用一致的字体、字号层级、数值色标和缺失值符号。浅色与深色网页中都要检查图内文字和边界对比度；瓦片底图的署名不能被图层或装饰遮盖。
案例：两张专题地图如何合成一张双变量地图
以美国本土 48 州的 3,108 个县及县级等价单元为例，先分别绘制 2025 年人口密度（A）与 2020—2025 年人口变化率（B）。人口密度用估计人口除以 Census 提供的土地面积；变化率使用同一版人口估计中的两个年份。两图保持相同的县界、范围和灰色底图，各自按县级数值的四分位分成四档，便于先读懂单项指标。
A · 2025 年县级人口密度
B · 2020—2025 年县级人口变化率
两图并列时，浅灰至蓝色只表示密度等级，浅灰至珊瑚色只表示变化率等级。点击图片可放大查看县界、图例和局部地理细节。
查看并复制：共同数据准备及 A、B 两张地图的代码
代码读取本章 data/ 下的 Census
固定快照，计算两个指标，再用相同范围和制图元素绘制
A、B。可从课程根目录、第四章目录或脚本所在目录运行；在其他位置运行时，先设置
CHAPTER4_DIR。在线 Esri 底图需要联网。
invisible(Sys.setlocale("LC_CTYPE", ".UTF-8"))
library(sf)
library(dplyr)
library(readr)
library(ggplot2)
library(patchwork)
library(ggspatial)
# 1. 读取第四章已经保存的 Census 2025 固定快照 ------------------------------
chapter_dir <- Sys.getenv("CHAPTER4_DIR", "")
if (!nzchar(chapter_dir)) {
candidates <- c(".", "第4章", "..")
has_county_data <- vapply(candidates, function(path) {
file.exists(file.path(path, "data/census-county-2025/cb_2025_us_county_5m.zip"))
}, logical(1))
if (!any(has_county_data)) {
stop("找不到第四章 data 目录；请先设置 CHAPTER4_DIR 为第四章目录。")
}
chapter_dir <- candidates[which(has_county_data)[1]]
}
chapter_dir <- normalizePath(chapter_dir, winslash = "/", mustWork = TRUE)
output_dir <- file.path(chapter_dir, "bivariate_county_preview")
dir.create(output_dir, showWarnings = FALSE, recursive = TRUE)
county_zip <- file.path(chapter_dir, "data/census-county-2025/cb_2025_us_county_5m.zip")
state_zip <- file.path(chapter_dir, "data/census-2025/cb_2025_us_state_20m.zip")
population_csv <- file.path(chapter_dir, "data/census-county-2025/co-est2025-alldata.csv")
input_files <- c(county_zip, state_zip, population_csv)
if (!all(file.exists(input_files))) {
stop("缺少 Census 固定快照：", paste(input_files[!file.exists(input_files)], collapse = "; "))
}
# SHP 是一组配套文件；解压到本次 R 会话的临时目录，不向章节目录留过程文件。
shape_dir <- file.path(tempdir(), "bivariate_county_shapes")
county_dir <- file.path(shape_dir, "county")
state_dir <- file.path(shape_dir, "state")
dir.create(county_dir, recursive = TRUE, showWarnings = FALSE)
dir.create(state_dir, recursive = TRUE, showWarnings = FALSE)
unzip(county_zip, exdir = county_dir)
unzip(state_zip, exdir = state_dir)
county_boundary <- st_read(
file.path(county_dir, "cb_2025_us_county_5m.shp"), quiet = TRUE
)
state_boundary <- st_read(
file.path(state_dir, "cb_2025_us_state_20m.shp"), quiet = TRUE
)
# 主图只画本土连续 48 州；放大图另纳入华盛顿特区。
states_48 <- state_boundary %>%
filter(!GEOID %in% c("02", "15", "11", "72")) %>%
st_transform(4326)
state_codes <- states_48$GEOID
states_with_dc <- state_boundary %>%
filter(GEOID %in% c(state_codes, "11")) %>%
st_transform(4326)
population <- read_csv(
population_csv,
col_types = cols(
SUMLEV = col_character(), STATE = col_character(), COUNTY = col_character(),
POPESTIMATE2020 = col_double(), POPESTIMATE2025 = col_double(),
.default = col_skip()
),
locale = locale(encoding = "Latin1"),
show_col_types = FALSE
) %>%
filter(SUMLEV == "050") %>%
transmute(
GEOID = paste0(STATE, COUNTY),
population_2020 = POPESTIMATE2020,
population_2025 = POPESTIMATE2025
)
# 2. 构建两个比率指标，并核查县界与人口表 ------------------------------
counties_with_dc <- county_boundary %>%
filter(STATEFP %in% c(state_codes, "11")) %>%
left_join(population, by = "GEOID", relationship = "one-to-one") %>%
mutate(
# ALAND 是 Census 提供的土地面积（平方米）；不用概化后的地图面积。
density_2025 = population_2025 / (as.numeric(ALAND) / 1e6),
# 两年人口均是 Vintage 2025 中 7 月 1 日的估计值。
growth_pct = 100 * (population_2025 / population_2020 - 1)
) %>%
st_transform(4326)
counties <- counties_with_dc %>% filter(STATEFP %in% state_codes)
if (nrow(states_48) != 48L || nrow(counties) != 3108L ||
nrow(counties_with_dc) != 3109L ||
sum(counties_with_dc$GEOID == "11001") != 1L ||
anyDuplicated(counties_with_dc$GEOID) ||
any(!is.finite(counties_with_dc$density_2025)) ||
any(!is.finite(counties_with_dc$growth_pct)) ||
any(counties_with_dc$population_2020 <= 0) ||
any(as.numeric(counties_with_dc$ALAND) <= 0)) {
stop("县级行数、DC 单元、地理键、人口或面积未通过核查，请检查输入快照。")
}
# 两个变量分别按主图 3,108 个县级单元的四分位分档；DC 沿用断点。
quartile_probs <- c(0.25, 0.50, 0.75)
quartile_levels <- paste0("Q", 1:4)
density_breaks <- unname(quantile(counties$density_2025, quartile_probs))
growth_breaks <- unname(quantile(counties$growth_pct, quartile_probs))
classify_quartile <- function(values, breaks) {
cut(values, breaks = c(-Inf, breaks, Inf),
labels = quartile_levels, include.lowest = TRUE)
}
counties_with_dc <- counties_with_dc %>%
mutate(
density_class = classify_quartile(density_2025, density_breaks),
growth_class = classify_quartile(growth_pct, growth_breaks),
density_rank = as.integer(density_class),
growth_rank = as.integer(growth_class),
bivariate_class = paste0(growth_rank, "-", density_rank)
)
counties <- counties_with_dc %>% filter(STATEFP %in% state_codes)
# 沿用 GrPink 的九个锚点，在 Lab 色彩空间插值为 4 × 4；
# 横向为变化率（灰 → 珊瑚），纵向为人口密度（灰 → 蓝）。
palette_anchors <- matrix(c(
"#e8e8e8", "#e4acac", "#c85a5a",
"#b0d5df", "#ad9ea5", "#985356",
"#64acbe", "#627f8c", "#574249"
), nrow = 3, byrow = TRUE)
palette_rows <- t(apply(palette_anchors, 1, function(colours) {
grDevices::colorRampPalette(colours, space = "Lab")(4)
}))
palette_grid <- apply(palette_rows, 2, function(colours) {
grDevices::colorRampPalette(colours, space = "Lab")(4)
})
bivariate_palette <- setNames(
as.vector(palette_grid),
unlist(lapply(1:4, function(growth) paste0(growth, "-", 1:4)))
)
density_palette <- setNames(palette_grid[, 1], quartile_levels)
growth_palette <- setNames(palette_grid[1, ], quartile_levels)
# 3. Esri 灰色底图与三张主图共用的制图元素 -------------------------
# World Light Gray Base 保留浅灰地理背景；不叠加会增加文字的 Reference 层。
# 缓存只保存在本次 R 会话临时目录，不随本项目分发瓦片文件。
tile_cache <- file.path(tempdir(), "bivariate_esri_gray_tiles")
dir.create(tile_cache, showWarnings = FALSE, recursive = TRUE)
rosm::register_tile_source(
esri_gray = rosm::as.tile_source(
"https://server.arcgisonline.com/ArcGIS/rest/services/Canvas/World_Light_Gray_Base/MapServer/tile/${z}/${y}/${x}",
attribution = "Esri, HERE, Garmin, OpenStreetMap contributors, and the GIS user community",
extension = "jpg", max_zoom = 16
)
)
tile_credit <- "底图：Esri Light Gray Canvas © Esri、HERE、Garmin、OpenStreetMap 贡献者等；专题数据：U.S. Census Bureau"
county_alpha <- 0.82
font_family <- "Microsoft YaHei"
# 以 EPSG:4326 显示：经线竖直、纬线水平，真北朝向页面正上方。
# 密度仍使用 Census ALAND 计算，不使用经纬度坐标下的多边形面积。
x_limits <- c(-125.5, -66.0)
y_limits <- c(23.7, 50.5)
map_graticule <- st_graticule(
st_as_sfc(st_bbox(c(xmin = x_limits[1], ymin = y_limits[1],
xmax = x_limits[2], ymax = y_limits[2]),
crs = st_crs(4326))),
crs = st_crs(4326), datum = st_crs(4326),
lon = seq(-120, -70, by = 10), lat = seq(25, 50, by = 5)
)
# 经纬度图的东西向比例随纬度改变。按下方比例尺所在的 25°N 附近，
# 用 WGS84 椭球上的实际距离求出 500 / 1000 km 对应的经度长度。
scale_lat <- 25.3
scale_start <- -123.5
distance_km <- function(delta_lon) {
points <- st_sfc(st_point(c(scale_start, scale_lat)),
st_point(c(scale_start + delta_lon, scale_lat)),
crs = 4326)
as.numeric(st_distance(points[1], points[2])) / 1000
}
lon_at_km <- function(km) {
scale_start + uniroot(function(delta) distance_km(delta) - km,
interval = c(0, 20))$root
}
scale_middle <- lon_at_km(500)
scale_end <- lon_at_km(1000)
scale_bar_layers <- list(
annotate("rect", xmin = scale_start, xmax = scale_middle,
ymin = scale_lat, ymax = scale_lat + 0.28,
fill = "#283944", colour = "#283944", linewidth = 0.35),
annotate("rect", xmin = scale_middle, xmax = scale_end,
ymin = scale_lat, ymax = scale_lat + 0.28,
fill = "white", colour = "#283944", linewidth = 0.35),
annotate("text", x = c(scale_start, scale_middle, scale_end),
y = scale_lat - 0.38, label = c("0", "500", "1000 km"),
size = 2.6, family = font_family, colour = "#34454e"),
annotate("text", x = (scale_start + scale_end) / 2,
y = scale_lat + 0.72, label = "约 25°N 处距离",
size = 2.35, family = font_family, colour = "#51616b")
)
one_decimal <- function(value) formatC(value, format = "f", digits = 1)
signed_pct <- function(value) sprintf("%+.1f%%", value)
density_labels <- c(
paste0("Q1 ≤", one_decimal(density_breaks[1])),
paste0("Q2 ", one_decimal(density_breaks[1]), "–", one_decimal(density_breaks[2])),
paste0("Q3 ", one_decimal(density_breaks[2]), "–", one_decimal(density_breaks[3])),
paste0("Q4 >", one_decimal(density_breaks[3]))
)
growth_labels <- c(
paste0("Q1 ≤", signed_pct(growth_breaks[1])),
paste0("Q2 ", signed_pct(growth_breaks[1]), "–", signed_pct(growth_breaks[2])),
paste0("Q3 ", signed_pct(growth_breaks[2]), "–", signed_pct(growth_breaks[3])),
paste0("Q4 >", signed_pct(growth_breaks[3]))
)
map_theme <- theme_minimal(base_family = font_family, base_size = 11) +
theme(
panel.grid = element_blank(), axis.title = element_blank(),
axis.text = element_text(size = 8.8, colour = "#586c77"),
axis.ticks = element_line(colour = "#8b999f", linewidth = 0.3),
axis.ticks.length = grid::unit(2, "pt"),
panel.background = element_rect(fill = "#fbfbfa", colour = NA),
panel.border = element_rect(fill = NA, colour = "#d8dcde", linewidth = 0.35),
plot.background = element_rect(fill = "white", colour = NA),
plot.title.position = "panel",
plot.title = element_text(size = 14, face = "bold", colour = "#26343e",
margin = margin(b = 4)),
plot.subtitle = element_text(size = 9, colour = "#64717b",
margin = margin(b = 8)),
legend.position = "bottom", legend.direction = "horizontal",
legend.title = element_text(size = 9, face = "bold", colour = "#33424d"),
legend.text = element_text(size = 8.5, colour = "#45525d"),
legend.background = element_rect(fill = "white", colour = NA),
plot.caption = element_text(size = 7.4, colour = "#71818a",
hjust = 1, margin = margin(t = 6)),
plot.margin = margin(10, 12, 9, 12)
)
map_layers <- function(fill_column, title, subtitle,
north_location = "br") {
ggplot(counties) +
annotation_map_tile(
type = "esri_gray", zoom = 4, alpha = 0.55,
cachedir = tile_cache, quiet = TRUE
) +
geom_sf(aes(fill = .data[[fill_column]]),
alpha = county_alpha, colour = "#f6f7f7", linewidth = 0.035) +
geom_sf(data = map_graticule, colour = "#718993",
alpha = 0.42, linewidth = 0.16, inherit.aes = FALSE) +
geom_sf(data = states_48, fill = NA, colour = "#5f6f79",
linewidth = 0.22, inherit.aes = FALSE) +
scale_x_continuous(
breaks = seq(-120, -70, by = 10),
labels = function(value) paste0(abs(value), "°W")
) +
scale_y_continuous(
breaks = seq(25, 50, by = 5),
labels = function(value) paste0(value, "°N")
) +
coord_sf(crs = st_crs(4326), default_crs = st_crs(4326),
datum = st_crs(4326), xlim = x_limits, ylim = y_limits,
expand = FALSE) +
scale_bar_layers +
annotation_north_arrow(
location = north_location, which_north = "true",
height = grid::unit(0.52, "in"), width = grid::unit(0.44, "in"),
pad_x = grid::unit(0.20, "in"),
pad_y = grid::unit(if (north_location == "br") 0.43 else 0.18, "in"),
style = north_arrow_fancy_orienteering(
line_col = "#263b45", fill = c("white", "#263b45"),
text_col = "#263b45", text_family = font_family, text_size = 9
)
) +
labs(title = title, subtitle = subtitle, caption = tile_credit) +
map_theme
}
four_class_guide <- guide_legend(
title.position = "top", label.position = "bottom", nrow = 1,
keywidth = grid::unit(0.54, "in"), keyheight = grid::unit(0.16, "in"),
override.aes = list(colour = NA)
)
map_density <- map_layers(
"density_class", "A 2025 年人口密度",
"县级人口 ÷ 官方土地面积；颜色由浅到深表示相对等级"
) +
scale_fill_manual(
values = density_palette, breaks = quartile_levels,
labels = density_labels, name = "人 / km²",
guide = four_class_guide, na.value = "#c7cbce"
)
map_growth <- map_layers(
"growth_class", "B 2020—2025 年人口变化率",
"变化率 =（2025 年人口 ÷ 2020 年人口 − 1）× 100%"
) +
scale_fill_manual(
values = growth_palette, breaks = quartile_levels,
labels = growth_labels, name = "变化率",
guide = four_class_guide, na.value = "#c7cbce"
)
if (interactive()) {
print(map_density)
print(map_growth)
}将两个四档结果组合，每个县会落入 4 × 4 种组合之一。C 图以一格颜色同时表示两项指标：色标横向读取变化率分位，纵向读取密度分位；下方四个局部图帮助查看纽约、华盛顿、洛杉矶与休斯敦周边的县级细节。
C · 两项县级指标的双变量地图与四个局部视窗
查看并复制：C 图的双变量配色、局部图与导出代码
先运行上方 A、B 的共同准备代码，再运行以下代码。C 图沿用相同的县界、分位断点、色盘和底图；四个局部图还沿用主图断点，因此可以直接比较颜色含义。
# 放大范围是城市周边的县级地图；点仅用于定位城市中心。
city_windows <- tibble::tribble(
~code, ~city_name, ~lon, ~lat, ~xmin, ~xmax, ~ymin, ~ymax,
"1", "纽约周边 · New York", -74.006, 40.713, -75.1, -72.7, 40.0, 41.3,
"2", "华盛顿周边 · Washington, DC", -77.037, 38.895,
-78.2, -75.8, 38.1, 39.5,
"3", "洛杉矶周边 · Los Angeles", -118.244, 34.052,
-119.7, -116.5, 33.15, 34.9,
"4", "休斯敦周边 · Houston", -95.370, 29.760,
-96.4, -94.3, 29.1, 30.5
)
window_geometries <- lapply(seq_len(nrow(city_windows)), function(i) {
st_as_sfc(st_bbox(c(xmin = city_windows$xmin[i],
ymin = city_windows$ymin[i],
xmax = city_windows$xmax[i],
ymax = city_windows$ymax[i]),
crs = st_crs(4326)))[[1]]
})
city_boxes <- st_sf(
city_windows,
geometry = do.call(st_sfc, c(window_geometries, list(crs = 4326)))
)
map_bivariate <- map_layers(
"bivariate_class", "C 人口密度 × 人口变化率",
"16 种颜色标示两项县级指标的空间组合；方框对应下方四个局部视窗"
) +
geom_sf(data = city_boxes, fill = NA, colour = "#263b45",
linewidth = 0.55, inherit.aes = FALSE) +
geom_label(data = city_windows,
aes(x = xmin, y = ymax, label = code),
inherit.aes = FALSE, hjust = 0, vjust = 1,
size = 2.5, family = font_family, fontface = "bold",
colour = "#263b45", fill = "white", linewidth = 0.18,
label.padding = grid::unit(0.07, "lines")) +
scale_fill_manual(values = bivariate_palette, guide = "none",
na.value = "#c7cbce") +
labs(caption = NULL) +
theme(legend.position = "none",
axis.text = element_text(size = 9.8, colour = "#586c77"))
# 图例两轴显示县级分位区间，不把分位百分比误写成实际变化率。
legend_data <- expand.grid(growth_rank = 1:4, density_rank = 1:4) %>%
mutate(bivariate_class = paste0(growth_rank, "-", density_rank))
legend_note <- paste0(
"密度断点 ", paste(one_decimal(density_breaks), collapse = " / "), " 人/km²\n",
"变化率断点 ", paste(signed_pct(growth_breaks), collapse = " / ")
)
map_legend <- ggplot(legend_data, aes(growth_rank, density_rank,
fill = bivariate_class)) +
geom_tile(width = 0.98, height = 0.98, colour = "white", linewidth = 0.55) +
scale_fill_manual(values = bivariate_palette, guide = "none") +
scale_x_continuous(breaks = 1:4,
labels = c("0–25", "25–50", "50–75", "75–100"),
limits = c(0.5, 4.5), expand = c(0, 0)) +
scale_y_continuous(breaks = 1:4,
labels = c("0–25", "25–50", "50–75", "75–100"),
limits = c(0.5, 4.5), expand = c(0, 0)) +
coord_equal() +
labs(title = "双变量色标 · 4 × 4", subtitle = "两轴均为县级分位区间（%）",
x = "变化率分位 →", y = "密度分位 ↑", caption = legend_note) +
theme_minimal(base_family = font_family, base_size = 10) +
theme(
panel.grid = element_blank(), plot.background = element_rect(fill = "white", colour = NA),
plot.title = element_text(face = "bold", size = 10.5, colour = "#26343e"),
plot.subtitle = element_text(size = 8, colour = "#64717b"),
axis.title = element_text(size = 8.5, colour = "#354550"),
axis.text = element_text(size = 7.5, colour = "#45525d"),
axis.ticks = element_blank(),
plot.caption = element_text(size = 7.3, colour = "#5f6b75", hjust = 0,
lineheight = 1.3, margin = margin(t = 8)),
plot.margin = margin(28, 8, 12, 3)
)
# 简洁的局部比例尺：按各图左下角纬度求实际 50 km 所对应的经度长度。
local_scale_bar <- function(city, km = 50) {
width <- city$xmax - city$xmin
height <- city$ymax - city$ymin
x0 <- city$xmin + 0.07 * width
y0 <- city$ymin + 0.105 * height
distance_at <- function(delta_lon) {
points <- st_sfc(st_point(c(x0, y0)),
st_point(c(x0 + delta_lon, y0)), crs = 4326)
as.numeric(st_distance(points[1], points[2])) / 1000
}
x1 <- x0 + uniroot(function(delta_lon) distance_at(delta_lon) - km,
interval = c(0, 2))$root
tick_height <- 0.022 * height
list(
annotate("rect", xmin = x0 - 0.03 * width,
xmax = x1 + 0.03 * width,
ymin = y0 - 0.050 * height, ymax = y0 + 0.095 * height,
fill = "white", alpha = 0.82, colour = NA),
annotate("segment", x = x0, xend = x1, y = y0, yend = y0,
linewidth = 0.55, colour = "#263b45"),
annotate("segment", x = c(x0, x1), xend = c(x0, x1),
y = y0 - tick_height, yend = y0 + tick_height,
linewidth = 0.55, colour = "#263b45"),
annotate("text", x = (x0 + x1) / 2, y = y0 + 0.055 * height,
label = paste0(km, " km"), size = 2.6,
family = font_family, colour = "#263b45")
)
}
make_city_inset <- function(i) {
city <- city_windows[i, ]
window <- city_boxes[i, ]
nearby_counties <- st_filter(counties_with_dc, window)
nearby_states <- st_filter(states_with_dc, window)
city_point <- st_as_sf(city, coords = c("lon", "lat"), crs = 4326)
inset_graticule <- st_graticule(
window, crs = st_crs(4326), datum = st_crs(4326),
lon = pretty(c(city$xmin, city$xmax), n = 3),
lat = pretty(c(city$ymin, city$ymax), n = 3)
)
ggplot(nearby_counties) +
annotation_map_tile(
type = "esri_gray", zoom = 8, alpha = 0.55,
cachedir = tile_cache, quiet = TRUE
) +
geom_sf(aes(fill = bivariate_class),
alpha = county_alpha, colour = "#f6f7f7", linewidth = 0.13) +
geom_sf(data = inset_graticule, colour = "#748b95",
linewidth = 0.16, alpha = 0.48, inherit.aes = FALSE) +
geom_sf(data = nearby_states, fill = NA, colour = "#536873",
linewidth = 0.35, inherit.aes = FALSE) +
geom_sf(data = city_point, shape = 21, fill = "#263b45",
colour = "white", stroke = 0.45, size = 2,
inherit.aes = FALSE) +
local_scale_bar(city) +
scale_fill_manual(values = bivariate_palette, guide = "none",
na.value = "#c7cbce") +
scale_x_continuous(
breaks = pretty(c(city$xmin, city$xmax), n = 3),
labels = function(value) paste0(abs(value), "°W")
) +
scale_y_continuous(
breaks = pretty(c(city$ymin, city$ymax), n = 3),
labels = function(value) paste0(value, "°N")
) +
coord_sf(crs = st_crs(4326), datum = st_crs(4326),
xlim = c(city$xmin, city$xmax),
ylim = c(city$ymin, city$ymax), expand = FALSE) +
labs(title = paste0(city$code, " ", city$city_name)) +
theme_minimal(base_family = font_family, base_size = 8) +
theme(
panel.grid = element_blank(),
panel.background = element_rect(fill = "#fbfbfa", colour = NA),
panel.border = element_rect(fill = NA, colour = "#cbd3d7",
linewidth = 0.4),
axis.title = element_blank(),
axis.text = element_text(size = 8.1, colour = "#526873"),
axis.ticks = element_line(colour = "#9aabb2", linewidth = 0.25),
axis.ticks.length = grid::unit(1.7, "pt"),
plot.title = element_text(face = "bold", size = 10,
colour = "#26343e", margin = margin(b = 4)),
plot.margin = margin(2, 8, 2, 8)
)
}
city_insets <- wrap_plots(
lapply(seq_len(nrow(city_windows)), make_city_inset), ncol = 4
)
bivariate_plate <-
(map_bivariate + map_legend + plot_layout(widths = c(4.8, 0.95))) /
city_insets +
plot_layout(heights = c(2.3, 1)) +
plot_annotation(
caption = paste0(
"城市点仅用于定位；放大图颜色表示县及县级等价单元，不代表城市本身的指标。 ",
tile_credit
),
theme = theme(
plot.caption = element_text(family = font_family, size = 8.5,
colour = "#5e6e77", hjust = 0,
margin = margin(t = 6)),
plot.background = element_rect(fill = "white", colour = NA)
)
)
if (interactive()) print(bivariate_plate)
# 4. 固定尺寸导出，保存分类数据和文字说明 ------------------------------
save_preview <- function(plot, filename, width, height) {
target <- file.path(output_dir, filename)
temp_file <- tempfile(pattern = "map_preview_", tmpdir = output_dir,
fileext = ".png")
on.exit(unlink(temp_file), add = TRUE)
ggsave(
filename = temp_file, plot = plot,
device = ragg::agg_png, width = width, height = height,
units = "in", dpi = 240, bg = "white"
)
if (file.exists(target) && unlink(target) != 0L) {
stop("无法替换旧图片：", target, "；请关闭占用该文件的预览窗口后重试。")
}
if (!file.rename(temp_file, target)) {
stop("无法保存新图片：", target)
}
}
save_preview(map_density, "01_density.png", 10.8, 6.7)
save_preview(map_growth, "02_growth.png", 10.8, 6.7)
save_preview(bivariate_plate, "03_bivariate.png", 13.6, 8.8)
# 正式网页引用固定文件名；一次绘图后同步更新三张网页图片。
formal_images <- file.path(chapter_dir, "images")
dir.create(formal_images, showWarnings = FALSE, recursive = TRUE)
formal_names <- c(
"01_density.png" = "c4_s5_bivariate_density.png",
"02_growth.png" = "c4_s5_bivariate_growth.png",
"03_bivariate.png" = "c4_s5_bivariate_joint.png"
)
for (source_name in names(formal_names)) {
copied <- file.copy(
from = file.path(output_dir, source_name),
to = file.path(formal_images, formal_names[[source_name]]),
overwrite = TRUE
)
if (!copied) stop("无法更新网页图片：", formal_names[[source_name]])
}读图与复现：色标上的“0—25%”等是县在当前样本中的分位区间，不是人口变化率的实际百分数；实际断点写在图例下方。16 色能呈现两项指标的空间组合，也增加了辨色负担，阅读时应配合二维色标和局部放大图。C 图显示空间共现，不能据此推断两项指标存在因果关系。地图使用的 Census 本地快照见本章
data/，重绘时需要联网获取 Esri 底图。
空间拥挤时的四种处理方式
大量点位或轨迹重叠时，直接加大颜色对比通常不能解决遮挡。下面四种方式分别保留或概括了不同的信息，选用时应说明它改变了什么。
01 · 透明度叠加
降低点或线的不透明度，让重叠区域逐渐显现。它保留单个要素位置，但亮度受符号大小和叠加顺序影响，适合定性查看聚集，不直接等于密度估计。
透明度叠加：重叠区域更醒目。
02 · 既有空间单元聚合
将点归入行政区或统计区，展示每个单元的数量或比例。若研究问题涉及人口风险或服务覆盖，应根据问题选择合适的人口、面积等分母；总量、面积密度和人均值回答不同问题。
按既有空间单元汇总点位。
03 · 规则网格聚合
按等尺寸方格或六边形汇总点位。网格便于比较局部差异，但单元大小和起点会影响图面；若各网格有效面积不同，还须说明是否进行了面积归一。
规则网格将点位概括为单元统计值。
04 · 核密度估计
将离散点位转为连续的平滑表面，便于观察聚集中心。结果受带宽、边界处理和点位分布影响；热点轮廓是估计结果，不是精确边界。
核密度估计的连续表面示意。
选法提示：要保留具体位置，先考虑透明度；要比较行政或统计单元，使用既有单元聚合；要比较规则空间位置，使用网格；要描述连续聚集趋势，使用核密度。四种方法可并用，但图注要说清各图的统计对象。
延伸阅读：地图设计与 R 制图资料
- Geocomputation with R：Making maps 与 ggplot2 地图章节提供 R 制图示例。
- ColorBrewer 2.0可筛选顺序、发散和定性色盘，并查看色觉与打印适配提示。
- DataShine和 London: The Information Capital展示城市数据地图的表达方式；#30DayMapChallenge可供观察不同题材的版式。
6.0 章节练习与测评
6.0 章节练习与测评
本节分为三个环节：先对照两张伦敦地图运行空间制图示例，再从城市房源、道路噪声与绿色基础设施三个案例中选择一种完成专题地图，最后用六道自测题回顾本章主要知识点。
任务 1：跟随示例运行
示例使用 Inside Airbnb 2025 年 9 月 14
日的伦敦房源快照与 spData::lnd
大伦敦行政区边界。选用约 14.6 MB 的精简房源表（visualisations/listings.csv）；它与约
50 MB 的完整压缩表含有相同的房源
ID、坐标和房型，足以复现本节两图。表中一行对应一条房源记录，经纬度确定点的位置，room_type
区分房源类型。首次运行需要联网下载；若数据源以后调整目录，应从 Inside Airbnb
数据页查找新的伦敦快照，并核对图中日期与房源数。
- 图 1：点分布图。 房源类型用颜色区分。重叠点会增强局部亮度，因此亮度只能辅助识别聚集区域，不能直接当作精确密度。
- 图 2：六边形网格图。 在英国国家格网坐标系（EPSG:27700）中生成约 0.5 km² 的六边形，仅保留与伦敦边界相交的网格，再统计各格房源数并按自然断点分级。图中只显示房源数不低于这些网格平均值的单元；市界内的空白区域不能读作“没有房源”。
下面先展示预存的参考图。它们来自上述日期的数据快照，不是本页 Knit 时重新下载并计算的结果。将完整代码复制到自己的 Rmd 后运行，再比较边界、图例、点群与网格分布。
图 2：六边形网格中的房源数
伦敦约 0.5 km² 六边形网格的房源数：仅显示达到与伦敦相交网格平均房源数的单元，颜色按 Fisher–Jenks 自然断点分为七级。
运行顺序：先加载依赖包并读取房源数据，再检查坐标缺失值与两份数据的
CRS；完成点图后，将数据转换到 EPSG:27700
生成网格并统计房源数。最后核对网格面积、筛选阈值及图例分级。install.packages()
只需在 Console 中按需运行一次；Rmd 内保留 library()
调用。
查看并复制两张地图的完整 R 代码
# ============================================================
# 系统环境配置与核心依赖库加载
# ============================================================
library(tidyverse)
library(sf)
library(ggspatial)
library(spData) # 提供大伦敦行政区划 (lnd) 边界数据集
library(classInt) # 执行 Fisher-Jenks 自然断点分类
# ============================================================
# 1. 读取精简房源表：保留本例所需的坐标与房型字段
# ============================================================
# 精简表仍是一行一房源，下载量比完整压缩表小。
airbnb_url <- "https://data.insideairbnb.com/united-kingdom/england/london/2025-09-14/visualisations/listings.csv"
london_airbnb <- read_csv(airbnb_url)
# ============================================================
# 2. 空间对象实例化与数据清洗 (Spatial Object Instantiation)
# ============================================================
airbnb_sf <- london_airbnb %>%
# 剔除坐标字段缺失的无效观测值，保障底层几何转换的稳健性
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 在米制投影下生成约 0.5 km² 的六边形，并限定研究范围
hex_grid <- st_make_grid(lnd_bng, cellsize = 759.83, square = FALSE) %>%
st_sf(hex_id = seq_along(.))
london_outline <- st_union(lnd_bng)
hex_grid <- hex_grid[lengths(st_intersects(hex_grid, london_outline)) > 0, ]
# 6.3 执行点面拓扑关联 (Spatial Join Count)
# 计算每个六边形单元内包含的 Airbnb 房源总量
hex_grid$listing_count <- lengths(st_intersects(hex_grid, airbnb_sf_bng))
# 6.4 以与伦敦相交的全部网格为参照计算平均值
# 边缘六边形保留完整形状；格内房源数是计数，不是按有效面积校正的密度。
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 使用 WGS 84 展示地图；米制网格已经在 EPSG:27700 中完成计算
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 listings \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 Listings per Hexagon in London",
subtitle = paste0("Tessellation Scale: ~0.5 km² per Hexagon | Classed by Fisher-Jenks"),
caption = "Data Source: Inside Airbnb & spData | Projection: WGS 84"
)运行后核对
- 核对点数据坐标是否为 WGS 84、
lnd边界是否与点位对齐；图 1 的颜色图例应与room_type类别一致。 - 核对网格是在米制投影下生成的，面积接近 0.5 km²；平均值只使用与伦敦相交的网格。图 2 的颜色表示每格房源数区间，不是房源评分或按边缘有效面积校正的密度。
- 从新的 R 会话运行完整脚本。若更换数据快照，记录日期、有效房源数和过滤后的网格数，再用两三句话解释地图上的聚集区域及一项限制。
任务 2：小组实践与展示
小组从下列三种方案中选择一种，制作一份可在本地 Knit 的 R Markdown 报告。方案 A 比较房源的点分布与网格高评分占比；方案 B 比较同一批道路街段的昼夜噪声；方案 C 将伦敦绿色基础设施的需求与供给等级放进同一张地图。三种方案都要核对数据口径、坐标系和图例，并解释地图能够支持的观察。地图与分析代码由小组完成。
选择研究案例
方案 A｜城市房源的高评分占比
研究问题
评分至少 4.9 分的房源在城市内如何分布？整套住房（Entire home/apt）与独立房间（Private room）的 500 m
网格高评分占比有何差异？
4.9 分是本练习的分析阈值，并非平台官方等级。
数据入口
从 Inside Airbnb
数据页选择同时提供 listings.csv.gz 和
neighbourhoods.geojson 的城市；记录城市与快照日期。可从西雅图
2026 年 9 月 20 日房源表及同日社区边界起步。
两份文件须来自同一城市、同一快照。下载后存入报告项目的
data/，并以相对路径读取。
- 清理并定位。 检查
longitude、latitude、room_type和review_scores_rating；剔除坐标无效的记录，保留评分缺失的房源用于点图与房源总数。将经纬度声明为 WGS 84，核对边界 CRS，再把点与边界转换到适合所选城市的米制投影。记录落在研究边界外的点数及处理方式。 - 绘制点图。 只比较上述两类房源，并按房型分面。区分评分 ≥4.9、有效评分 <4.9 与评分缺失三种状态；两幅分面共用颜色规则和地图范围，以点径和透明度控制重叠。点位只用于讨论区域格局，不应解读为精确门牌位置。
- 建立并汇总网格。 在米制投影下生成边长 500
m 的完整正方形网格（每格 0.25
km²），设置网格
ID，并让每条房源只进入一个网格。按网格 ID ×
房型统计
listing_count（全部房源）、rated_count（有效评分）和high_rating_count（评分 ≥4.9）。计算high_rating_share = high_rating_count / rated_count；缺失评分不进入分母，仅在rated_count ≥ 5时着色。核对各格房源数之和与纳入点数。 - 绘制分面网格图。 合并两种房型的达标网格，计算一套 Fisher–Jenks 五级断点；两个分面共用断点、色盘和地图范围。分别标示无房源、无有效评分与有效评分不足 5 条的网格，不把它们当作低占比。占比取值不足以形成五个不同断点时，可减少级数并说明原因。沿海边缘只裁剪网格的显示范围，统计仍以完整方格为单位。
展开查看：方案 A 的西雅图参考图
下图采用西雅图 2026 年 9 月 20 日快照。它展示图例和分面网格的一种做法；选择其他城市时，应按该城市的数据重新计算分级断点，无需复现西雅图的区间数值。图中的方块是 500 m 统计网格，并非地籍地块；无方格的留白包括水域和研究范围外的邻近地区。独立房间达到每格五条有效评分门槛的网格较少，浅色状态格不能解释为低占比。
为便于定位，图中以 Westlake Park（市政府记录地址：401 Pine St）作为市中心参考点。虚线弧对应半径 2.5、5、7.5、10 km 的同心圆，仅显示在本图采用的社区边界内。它们表示直线距离，不能解读为步行距离或可达范围；Westlake Park 也不是城市边界的几何中心。
西雅图房源高评分占比：500 m 网格与距 Westlake Park 的直线距离弧。
方案 B｜苏黎世道路噪声的昼夜对照
研究问题
同一批苏黎世道路街段的模型噪声值在白天和夜间如何变化？哪些街段的昼夜差值较大？
lreday 与 lrenight
是道路声源的模型排放级（dB(A)），不等于路边实测声压或居民受到的噪声暴露；局部放大图标出具体街段，不代表整条道路的平均值。
数据入口
使用苏黎世市政府开放数据平台的道路噪声图层与同平台的 Stadtkreise 市区边界。
参考图从 streetemission 线要素中选取
ref_year == 2024 且 tram_noise == "nein" 的
4,489 条街段。原始图层有 7,295 条线；其余 2,806
条标记为
tram_noise == "ja"，本例为统一道路声源口径暂不纳入。排除记录不表示相应道路没有交通噪声。
下载官方数据：德语界面指引
- 打开上方道路噪声数据页，点击 Daten & Ressourcen，找到 GeoPackage（.gpkg） 资源，再点击同一行的 Herunterladen 按钮。边界数据从上方 Stadtkreise 页面以相同方式下载。
- 在打开的市政府下载配置页，将 Ausschnitt 设为 Gesamter Datensatz，核对 Format 为 GeoPackage（.gpkg），然后点击 Download。通知邮箱为可选项。
- 生成完成后，到 Meine Bestellungen
下载压缩包并解压。下载链接通常只保留五天；将解压后的数据放入报告项目的
data/文件夹，用相对路径读取。网站更新日期可能晚于数据的ref_year，请以字段核对分析年份。
| 德语页面用语 | English | 中文 |
|---|---|---|
| Daten & Ressourcen | Data & resources | 数据与资源 |
| Herunterladen | Download | 下载 |
| Ausschnitt · Gesamter Datensatz | Extent · Entire dataset | 范围 · 完整数据集 |
| Format · GeoPackage | Format · GeoPackage | 格式 · GeoPackage |
| Meine Bestellungen | My orders | 我的下载任务 |
展开查看：官网实际页面的三个下载步骤
① 找到官方数据集和资源入口。
苏黎世市政府道路噪声数据集页面，页面标题下方有
Daten & Ressourcen 入口。
② 在 GeoPackage 资源行点击 Herunterladen。 该行标明投影为 EPSG:2056。
GeoPackage 资源行标出 EPSG:2056，下载按钮文字为
Herunterladen。
③ 核对范围与格式，点击 Download。 截图显示完整数据集和 GeoPackage；邮箱栏可以留空。
市政府下载配置页显示 Ausschnitt 为 Gesamter
Datensatz、Format 为 Geopackage，右侧有 Download 按钮。
页面截图：苏黎世市政府开放数据平台，2026 年 10 月访问；界面可能更新，以页面现有字段为准。
推荐的 GeoPackage 为瑞士米制坐标系 EPSG:2056，便于检查比例尺与局部距离；若选 GeoJSON（EPSG:4326），进行米制距离操作前须先转换坐标系。下载后先查看文件内的图层名和字段，不预设压缩包内部文件名。
- 确定样本与范围。
读取道路和市区边界，检查几何、
st_crs()、ref_year与tram_noise。参考图采用上述 2024 年、tram_noise == "nein"的样本；若使用另一筛选口径，须说明并在三张图中保持一致。记录原始线数、筛选后街段数及白天或夜间原始值为 0 的数量。 - 绘制昼夜对照。 同一批街段分别按
lreday（06:00—22:00）与lrenight（22:00—06:00）着色。白天用浅色背景，夜间用深色背景；两图共用范围、5 dB 分级和噪声色盘。原始 0 值单列，暂不解释为“没有噪声”。两图均需市区边界、上北指北针、比例尺、图例及来源。 - 计算街段差值。
仅当同一街段的两个时段数值均大于 0 时计算
lreday − lrenight（单位为 dB）；原始 0 保留独立状态，不参与减法。差值图使用单独的顺序或发散色标，放大至少一处高差值街段，让主图与局部图可以相互定位。图注说明正负含义，不把分贝差写成百分比变化。 - 核对并解释。 列出可计算差值的街段数及差值范围，指出一处变化明显的区域或街段。说明颜色对应的指标、被排除或单列的记录，以及模型数据不能支持的判断。局部图做一处即可；参考图的三处放大仅作版式示例。
展开查看：方案 B 的昼夜对照与差值参考图
两图使用同一批 4,489 条道路街段。原始值 0 在白天有 6 条、夜间有 52 条，含义尚未从官方字段说明得到确认，因此参考图以蓝色单列；去除这 58 条后，4,431 条街段可以计算昼夜差值。5 dB 色阶与颜色参考荷兰 IPLO 的官方噪声地图色表；原色表用于环境受声地图，这里仅借用视觉编码，不表示本图指标的法规限值。
苏黎世道路声源模型值的昼夜对照：两图使用同一批街段、范围和 5 dB 色阶；原始值 0 单列。
同一街段的白天减夜间差值：右侧三张局部图分别标出一条高差值街段。
方案 C｜伦敦绿色基础设施的需求与供给
研究问题
伦敦哪些六边形同时呈现高需求、低供给？把两项官方评级放在同一张图后，这一组合在空间上如何分布？
这里的“需求”与“供给”均为伦敦绿色基础设施框架（LGIF）的综合等级，不能解释为绿地面积、覆盖率或政策效果。
数据入口
从 London
Datastore 的 LGIF 官方数据页下载
LGIF hex results (01.04.26).csv
修正版、Hex350_grid_GLA.zip、字段字典与色彩规范，核对等级含义。
如需叠加伦敦各区轮廓，可从 London
Boroughs 页面下载 GeoPackage；其约 2018
年的区界仅作位置参照。文件存入报告项目的
data/，以相对路径读取。
官网下载：修正版结果表与六边形网格
- 打开 LGIF 数据页，先确认说明中的 01.04.26 修正版。官网指出早期 CSV 有错误，已由此版本替代；不要沿用旧文件。
- 在 7 files 列表中下载 LGIF hex results (01.04.26).csv、Hex350_grid_GLA.zip、LGIF data dictionary.xlsx 与 LGIF webmap colour specification.pdf。CSV 存评级，ZIP 存几何；两份文件须配合使用。需要追溯指标构建时，再阅读技术方法报告。
- 若使用行政区轮廓，进入 London Boroughs 页面，在 3 files 中选择第三项 London Boroughs · Geopackage（.gpkg），不要误点同名的在线地图服务。该文件标注 British National Grid（EPSG:27700）；网页也说明其约 2018 年的来源。
| 页面文字 / 文件 | English | 中文用途 |
|---|---|---|
| LGIF hex results (01.04.26).csv | Corrected results table | 修正版六边形评级表 |
| Hex350_grid_GLA.zip | Hexagon geometry | 350 m 六边形网格 |
| LGIF data dictionary.xlsx | Data dictionary | 核对指标与字段含义 |
| London Boroughs · Geopackage | Borough boundary | 可选的伦敦区界参照 |
展开查看：官网真实页面与三个下载位置
① 在 LGIF 数据页核对修正版和开放许可。 右侧标注 Open Government Licence v3。
London Datastore 的 LGIF 数据页，正文标出
01.04.26 修正版替代早期 CSV，右侧显示开放许可。
② 在七个资源中下载结果表、网格、字段字典和色彩规范。
LGIF 页面七个下载资源的列表，包含修正版 CSV、350
m 六边形网格 ZIP、字段字典与网页色彩规范。
③ 如需区界，从另一个官方页面下载 GeoPackage。 图中第三项是文件；第二项是在线服务。
London Boroughs 页面，第三项 GeoPackage 显示
EPSG:27700，上方说明约 2018 年 Ordnance Survey 边界来源。
页面截图：London Datastore，2026 年 10 月访问；页面名称和排列可能更新，以官网现有说明为准。
两项数据页均标注 Open Government Licence v3。公开展示时注明 LGIF 的 Greater London Authority 来源；若使用区界，另注明 Contains OS data © Crown copyright and database right 2018。瓦片底图的署名须按所选服务补充。
- 核对并连接。
解压网格，检查几何、
Hex_ID、评级表的hex_id、重复 ID 与缺失值；按 ID 连接，不按行号拼接。参考版本应有 15,041 个六边形，连接后网格数不变、两项综合评级均有值。原始网格为英国国家格网（EPSG:27700）；CSV 另带 EPSG:3857 的 WKT，两套坐标不能直接混用。 - 整理双变量等级。 用
overall_need_txt与overall_provision_txt的官方低／中／高类别形成 3 × 3 组合，列出九类网格的数量和占比。保留官方分类，不重新计算 Fisher 断点，也不把类别当作连续分数。图例标明需求向上增加、供给向右减少。 - 绘制全域地图。 为九类设置含义一致的颜色，参考官网色彩规范；瓦片底图适度淡化，让区界、河流或道路提供背景而不盖住六边形。地图须有比例尺、上北指北针、九宫格图例与数据／底图署名。可自行放大一处区域，不必复刻参考图布局。
- 解释并核对。 报告高需求 × 低供给组合的网格数和占比，指出图上一处可辨认的分布区域；再自选另一类组合，报告其数量或占比并比较分布。用 4–6 句说明叠加分析多提供了什么信息，以及至少一项限制：综合评级的构建、350 m 网格尺度或约 2018 年区界仅作参照。不能据此认定地块绿地面积或干预效果。
展开查看：方案 C 的伦敦 GI 双变量地图参考图
参考图采用修正版的 15,041 个六边形、官方九色、淡化的无地名瓦片底图与约 2018 年区界。图例的两个方向和来源署名可供核对；地图中的 285 个高需求、低供给网格（1.9%）是该固定版本的结果，若使用后续版本应重新计算。参考图只展示一种成品形式，报告仍须自行完成数据核对、九类统计和解释。
伦敦绿色基础设施综合需求与供给：3 × 3 官方评级、区界与淡化瓦片底图。
共同交付要求
每份报告都要写明选用方案、数据来源与参考年或快照日期、研究范围、清理规则及 CRS，并提交该方案要求的地图与关键统计：
- 方案 A｜房源：按房型分面的点图与 500 m 网格图；有效评分数、高评分房源数、达标网格数及共用的 Fisher 断点。
- 方案 B｜噪声：昼夜对照图与街段差值图；原始街段数、原始 0 记录数、可计算差值的街段数及差值范围。
- 方案 C｜绿色基础设施：伦敦全域双变量图与九类统计；高需求、低供给及自选另一类组合的数量或占比。
用 4–6
句解释图中可核对的空间差异，并说明至少一项数据或方法限制。代码按读取、清理、空间处理和出图分段注释；下载文件放在项目
data/，以相对路径读取。
遇到分析单位或计算口径问题时，展开查看提示
- 方案 A：先把房源点分配给网格，再按“网格 ×
房型”统计。
listing_count包含评分缺失的房源；高评分占比只以有评分的房源数为分母。 - 方案
B：分析单位是同一街段上的两个时段值，无需另建网格。先用
st_crs()检查投影，再为昼夜图共用分级；dplyr::if_else()可帮助只对两个时段都大于 0 的街段计算差值。局部图可通过裁剪和高亮制作。 - 方案 C：官方六边形与评级表按共同 ID 连接。核对 ID
唯一性、连接后行数与
st_crs()；九格统计来自两个已有类别字段，无须重新分级或生成网格。底图置于专题网格之下，并保留署名。
上述提示只说明处理方向，地图和代码仍由小组自行完成。