多源数据分析与实践 · Week 5

1 / 30

第四章重点解析

主讲人：刘贇喆副教授 | 华东师范大学社会发展学院

空间数据导论 · 空间数据读写 · 空间数据处理 · 空间数据可视化 · 空间制图要求与美化

2 / 30

目录 CONTENTS

1.0 空间数据导论

3 / 30

SECTION01 / 05

1.0空间数据导论

Introduction to Spatial Data

1.0 空间数据导论

4 / 30

空间数据的两种形态

Vector Features & Raster Cells

空间要素 · 教学示意

行 × 列一格，一个值像元（栅格像素）

1.0 空间数据导论

5 / 30

空间点数据

Spatial Point Data

VECTOR LAYER / POINT点图层

SF FEATURE TABLE

点要素表

ID	属性	坐标对（x, y）

尚无记录—

1.0 空间数据导论

6 / 30

空间线数据

Spatial Line Data

VECTOR LAYER / LINE线图层

SF FEATURE TABLE

线要素表

ID	属性	顶点概览

尚无记录—

1.0 空间数据导论

7 / 30

空间面数据

Spatial Polygon Data

VECTOR LAYER / POLYGON面图层

SF FEATURE TABLE

面要素表

ID	属性	闭合顺序

尚无记录—

1.0 空间数据导论

8 / 30

sf：一个图层，一张要素表

One Layer, One Feature Table

library(sf)

三个独立的 sf 图层

st_read() → sf

一行一个要素，一列一个变量

非空间属性 · ATTRIBUTESID · 名称 · 类别 · 数值 · …

空间／几何属性geometrysfc 几何列：位置与形状

geometry 预览点击表格行查看坐标

1.0 空间数据导论

9 / 30

坐标参考系：位置与投影

Coordinate Reference Systems

地球表面与平面地图同一组位置 · 三种地图表达

地球表面点击球体查看经纬网

WGS 84 · EPSG:4326

平面地图经纬度直绘 · 高纬也会变形

2.0 空间数据读写

10 / 30

SECTION02 / 05

2.0空间数据读写

Reading and Writing Spatial Data

2.0 空间数据读写

11 / 30

空间数据存储方式

How Spatial Data Are Stored

真实文件 · sf 包

北卡罗来纳县界 Shapefile

100 个县级面 · EPSG:4267

同一目录中的配套文件

nc.shp几何形状

nc.shx几何索引

nc.dbf属性表

nc.prj坐标系定义

实际路径：sf/shape/nc.* · 四个文件同名并位于同一目录

读入 sf

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

读取时指定 .shp；复制或分享数据时一并保留配套文件。.prj 提供 CRS 定义。

真实文件 · sf 包

同一县界的单文件版本 GeoPackage

100 个县级面 · EPSG:4267

文件与实际图层

文件gpkg/nc.gpkg

图层nc.gpkg100 个面要素

GeoPackage 可存多个图层；这个样本经 st_layers() 核对，只有 1 个图层。

先查看，再读入

gpkg_path <- system.file("gpkg/nc.gpkg", package = "sf")
st_layers(gpkg_path)
nc_sf <- st_read(gpkg_path,
                 layer = "nc.gpkg", quiet = TRUE)

layer 使用实际图层名，不能把示意名称直接写进代码。

真实文件 · spData 包

伦敦共享单车站点 GeoJSON

532 个点 · sf 读入后 EPSG:4326

首个 Feature 的实际字段

type / geometryFeature / Point

osm_id"108539"

name · capacityWindsor Terrace · "24"

coordinates [lon, lat][-0.093387797474861, 51.529125213623047]

顶层为 FeatureCollection · 原文声明 CRS84 · spData/shapes/cycle_hire_osm.geojson

读入 sf

geojson_path <- system.file(
  "shapes/cycle_hire_osm.geojson", package = "spData")
cycle_sf <- st_read(geojson_path, quiet = TRUE)

文件保存要素属性和点坐标；普通 .json 文件未必符合 GeoJSON 结构。

真实节选 · 本课件文件

伦敦共享单车站点 CSV + 坐标列

12 条记录 · 5 个字段

CSV 原文节选

osm_id,name,capacity,lon,lat
108539,Windsor Terrace,14,-0.09338779747486115,51.52912521362305

lon → X / 经度lat → Y / 纬度

实际文件：media/cycle_hire_osm_excerpt.csv · 与上面的 GeoJSON 来自不同历史快照

由坐标列生成 sf

cycle_tbl <- read_csv("cycle_hire_osm_excerpt.csv",
  col_types = cols(osm_id = col_character()))
cycle_sf <- st_as_sf(cycle_tbl,
  coords = c("lon", "lat"), crs = 4326)

CSV 本身不保存 CRS；这里的 EPSG:4326 已根据来源对象核对。运行代码时须先取得该 CSV 文件。

2.0 空间数据读写

12 / 30

三种格式读入，GeoJSON 写出

Read with sf · Export as GeoJSON

SOURCE FILES

选择一种输入格式

R / sf · 运行前加载 sf + tidyverse

读入 sf → 写出新 GeoJSON

# 01 · 读入 Shapefile
shp_path <- system.file("shape/nc.shp", package = "sf")
nc_sf <- shp_path %>% st_read(quiet = TRUE)

# 02 · 转为 WGS 84，写出 GeoJSON
out_geojson <- tempfile(fileext = ".geojson")
nc_sf %>%
  st_transform(4326) %>%
  st_write(
    dsn = out_geojson, driver = "GeoJSON",
    layer_options = "RFC7946=YES", quiet = TRUE)

# 03 · 重读并核对
nc_out <- out_geojson %>% st_read(quiet = TRUE)
c(rows = nrow(nc_out), epsg = st_crs(nc_out)$epsg)

# 01 · 查看图层，读入 GeoPackage
gpkg_path <- system.file("gpkg/nc.gpkg", package = "sf")
gpkg_path %>% st_layers()
nc_gpkg_sf <- gpkg_path %>% st_read(layer = "nc.gpkg", quiet = TRUE)

# 02 · 转为 WGS 84，写出 GeoJSON
out_geojson <- tempfile(fileext = ".geojson")
nc_gpkg_sf %>%
  st_transform(4326) %>%
  st_write(
    dsn = out_geojson, driver = "GeoJSON",
    layer_options = "RFC7946=YES", quiet = TRUE)

# 03 · 重读并核对
nc_out <- out_geojson %>% st_read(quiet = TRUE)
c(rows = nrow(nc_out), epsg = st_crs(nc_out)$epsg)

# 01 · 读入 GeoJSON
geojson_path <- system.file(
  "shapes/cycle_hire_osm.geojson", package = "spData")
cycle_sf <- geojson_path %>%
  st_read(quiet = TRUE)

# 02 · 另存为新的 GeoJSON 文件
out_geojson <- tempfile(fileext = ".geojson")
cycle_sf %>%
  st_write(
    dsn = out_geojson, driver = "GeoJSON",
    layer_options = "RFC7946=YES", quiet = TRUE)

# 03 · 重读并核对
cycle_out <- out_geojson %>% st_read(quiet = TRUE)
c(rows = nrow(cycle_out), epsg = st_crs(cycle_out)$epsg)

2.0 空间数据读写

13 / 30

从 CSV 坐标表到 sf

From a CSV Table to Spatial Features

伦敦共享单车 · spData / © OpenStreetMap contributors

`cycle_tbl`

# A tibble: 12 × 5

osm_id <chr>	name <chr>	capacity <dbl>	lon <dbl>	lat <dbl>

TIBBLE COLUMN ROLES

识别属性与坐标列

CSV 中 lon / lat 是位置数值列；生成 sf 后才有 geometry。

OBSERVATION 01Windsor Terrace

X / lon −0.093388° · Y / lat 51.529125°

伦敦共享单车 · spData / © OpenStreetMap contributors

`cycle_sf`

# sf: 12 × 4 · POINT · EPSG:4326

osm_id <chr>	name <chr>	capacity <dbl>	geometry <POINT>

st_as_sf()

生成点时发生了什么？

coords

lon → X，lat → Y列的顺序决定坐标顺序。

crs

4326 · WGS 84依据已核对的源对象声明 CRS；此步没有转换坐标。

remove

默认 TRUE · 不保留原列lon / lat 合成 geometry 后，不再作为属性列显示。

选中行的 geometry点 (−0.093388, 51.529125)

R · 已加载 sf + tidyverse

读取 CSV 数据

cycle_tbl <- read_csv("cycle_hire_osm_excerpt.csv",
                      col_types = cols(osm_id = col_character()))

cycle_sf <- cycle_tbl %>%
  st_as_sf(coords = c("lon", "lat"), crs = 4326)

cycle_bng <- cycle_sf %>% st_transform(27700)

2.0 空间数据读写

14 / 30

把指标接到边界上

Join a Data Table to Spatial Boundaries

边界 `boundary_sf`

geo_id	名称	geometry
	A 区	POLYGON
	B 区	POLYGON
	C 区	POLYGON

指标 `attr_tbl`

geo_id	指标值
`"01"`	8.4
`"03"`	13.1

结果 `joined_sf`示意数据

点击上方代码，查看连接结果

geo_id	名称	指标值	geometry
`"01"`	A 区	8.4	POLYGON
`"02"`	B 区	NA	POLYGON
`"03"`	C 区	13.1	POLYGON

3.0 空间数据处理

15 / 30

SECTION03 / 05

3.0空间数据处理

Spatial Data Processing

3.0 空间数据处理

16 / 30

从属性处理到空间处理

From Attributes to Spatial Relationships

ONE SF OBJECT

同一份 `sf` 数据

一行一个地理要素

geo_id	名称	人口	geometry
01	A 区	2.4 万	POLYGON
02	B 区	1.8 万	POLYGON
03	C 区	3.1 万	POLYGON

沿用第二章

sf：带 geometry 列的 data.frame
filter() 筛选要素 · mutate() 新增属性

GEOMETRY IN ACTION

让位置参与处理

由 geometry 判断要素之间的空间关系

01

判断空间关系

点落在哪个区域内？

st_within()

02

空间连接

按位置匹配，再把区域属性接给点

st_join()

03

几何运算

按设定距离生成缓冲范围

st_buffer()

3.0 空间数据处理

17 / 30

空间关系：按位置找到匹配

Spatial predicates · Which features match?

教学示意点是否位于面 A 内部？

边界上的 Q 与内部的 P，判断结果不同。

3.0 空间数据处理

18 / 30

空间连接与叠置

Spatial Join & Overlay

点落入街区，获得街区属性

教学示意 · TWO SF LAYERS

站点落在哪个街区？

points_sfregions_sf

A居住区人口 24,000

B商业区人口 18,000

C混合区人口 8,000

区域外

P1位于 A 区 · 居住区 · 人口 24,000连接前核对两图层 CRS

POINT IN POLYGON

街区属性写入站点表

点原有	从面带入的属性			点几何
point_id	region_id	area_type	population	geometry
P1	A	居住区	24000	POINT
P2	A	居住区	24000	POINT
P3	B	商业区	18000	POINT
P4	B	商业区	18000	POINT
P5	NA	NA	NA	POINT

点击代码，查看面的属性如何进入点表。

面积加权 · 教学示意

A、B 各覆盖哪些分区？

街区target_sf来源区source_sf

A总面积 10 km²；交叠 4 / 4 / 2 km² → 36%等积 CRS 测面积

OVERLAY → AREA-WEIGHTED MEAN

按交叠面积计算街区绿地率

目标面	面积加权	目标面几何
A	36%	POLYGON
B	60%	POLYGON

假设来源区内绿地率均匀；按交叠面积加权。

3.0 空间数据处理

19 / 30

缓冲区分析：从点生成面

Buffer Analysis · From Point to Polygon

伦敦共享单车站点 · 中心相距约 670 m

两座站点，同一缓冲半径

Waterloo Embankment 交叠

站点：spData 历史 OSM 快照 · © OpenStreetMap contributors　|　英国国家格网 EPSG:27700

R / sf

两个点，各生成一个面

EPSG:4326 经纬度EPSG:27700 米制投影

library(sf)
library(tidyverse)
data("cycle_hire_osm", package = "spData")
stations_m <- cycle_hire_osm %>% st_transform(27700)
centres_m <- stations_m %>%
  filter(name %in% c("Waterloo Bridge, South Bank",
                     "Embankment (Horse Guards)"))
radius_m <- 500
buffers_sf <- centres_m %>% st_buffer(dist = radius_m)

缓冲半径

2 POINTS 2 POLYGONS

站间距约 670 m

每站各得一个面，交叠处不会自动合并；500 m 是直线半径，并非步行距离。

3.0 空间数据处理

20 / 30

st_drop_geometry()：只留下属性

Remove Geometry · Keep the Attributes

输入 · SF

`points_joined`

point_id	region_id	area_type	population	geometry
P1	A	居住区	24000	POINT
P2	A	居住区	24000	POINT
P3	B	商业区	18000	POINT
P4	B	商业区	18000	POINT
P5	NA	NA	NA	POINT

运行函数

保留属性

只移除 geometry；原 sf 不变

输出 · 普通属性表

`point_tbl`

运行中间的函数，观察哪一列消失。

去掉几何后，继续处理属性

先运行 st_drop_geometry()，再选择一种属性运算。

筛出所在区域人口 ≥ 20000 的站点

selected_points <- point_tbl %>%
  filter(population >= 20000)

selected_points · 2 行

point_id	region_id	area_type	population
P1	A	居住区	24000
P2	A	居住区	24000

filter() 也可直接筛选 sf，并保留几何；这里先去几何，是为了得到纯属性表。

按区域统计站点数

point_counts <- point_tbl %>%
  filter(!is.na(region_id)) %>%
  group_by(region_id) %>%
  summarise(n_points = n(), .groups = "drop")

point_counts · 2 行

region_id	n_points
A	2
B	2

直接汇总 sf 会合并组内几何；先去几何只生成属性统计。P5 在计数前排除。

4.0 空间数据可视化

21 / 30

SECTION04 / 05

4.0空间数据可视化

Visualizing Spatial Data

4.0 空间数据可视化

22 / 30

ggplot + sf：用 geom_sf() 绘图

The Same Grammar · A Spatial Geometry Layer

county_map.R

2025 · 3,108 个县级单元

已整理的县界 sfclass_43_counties

class_43_counties %>%
  ggplot()

class_43_counties %>%
  ggplot() +
  geom_sf()

class_43_counties %>%
  ggplot() +
  geom_sf(fill = "#E8D7BD",
          colour = "#8B6557",
          linewidth = 0.15) +
  labs(title = "US counties · 2025") +
  theme_minimal()

ggplot() 接收 sf 数据；此时尚无几何图层。

DRAWN FROM sf

美国本土县级边界

01 / 03

class_43_counties %>% ggplot()

画布已建立，尚未加入绘图图层。

4.0 空间数据可视化

23 / 30

同一数据，三种分级断点

Classification Methods · Same Data, Different Breaks

2025 年美国本土 · 3,108 个县级单元人口密度（人/km²） → log10(1 + 密度) → 五级分色

分级方法各级县数 · 由浅到深四个内部断点 · 人/km²

Fisher–Jenks自然断点

寻找组内差异较小的分组

557 · 894 · 931 · 514 · 212

3.915.754.4269.6

Quantile分位数

每级县级单元数尽量接近

622 · 621 · 622 · 621 · 622

4.712.325.066.4

Equal interval对数等间距

在对数尺度划分等宽区间

808 · 1,626 · 548 · 117 · 9

7.060.5474.23,670.7

4.0 空间数据可视化

24 / 30

三种分级方法的地图对比

Three Classification Methods · One County Map

美国县级人口密度 · 2025

Fisher–Jenks 五级分层设色图：美国本土 48 州的 3,108 个县级单元按人口密度着色；附真北箭头和 0–800 km 投影距离比例尺

classInt::classIntervals(
  ..., style = "fisher"
)

密度区间（约）· 人/km²县级单元数

0.03–3.9557
3.9–15.7894
15.7–54.4931
54.4–269.6514
269.6–28,370212

分级变量：log10(1 + 人口密度)

5.0 空间制图要求与美化

25 / 30

SECTION05 / 05

5.0空间制图要求与美化

Cartographic Standards and Design

5.0 空间制图要求与美化

26 / 30

中国地图：先过三道检查

Basemap Source · Map Content · Intended Use

01 BASEMAP

底图有来源

优先在自然资源部标准地图服务系统选择适合图幅的底图。

图名审图号下载日期

仅凭第三方文件名，不能确认地图内容已审核。

02 EXTENT

内容按图幅核对

中国全图应表达完整范围，并核对图幅和比例尺下应表示的重要岛屿。

大陆海南岛台湾岛南海诸岛

南海附图须清楚。西藏周边等国界关键部位应对照国家标准画法核对；裁切或叠图不得改变、遮挡界线。

03 USE

按用途核对

课堂练习、直接使用标准图、公开新增专题图，分别记录来源与要求。

课堂练习 注明来源与限制直接用标准图 保留原图与审图号公开新图 核对审核要求

底图审图号 ≠ 新图审图号

三条信息分开写 底图图名 · 审图号 专题数据 来源 · 年份新图用途 · 审核状态

5.0 空间制图要求与美化

27 / 30

用 R 制作中国专题地图

A County-Level Choropleth with a South China Sea Inset

Albers 等面积投影

2020 年中国县级常住人口暖色深色地图。主图聚焦大陆、海南岛与台湾岛，南海附图完整展示诸岛及十段线；主图右下海面标有当地真北箭头，两图各有局地比例尺。台湾省资料暂缺写于图下注释。

2020 年中国县级常住人口浅色地图。主图聚焦大陆、海南岛与台湾岛，南海附图完整展示诸岛及十段线；主图右下海面标有当地真北箭头，两图各有局地比例尺。台湾省资料暂缺写于图下注释。

5.0 空间制图要求与美化

28 / 30

四种设色方案的地图效果

Categorical · Diverging · Continuous · Classified

美国本土 48 州县级地图，按美国人口调查局四大普查区域以 Set2 四色填充，含经纬网、图例、比例尺与真北指示。

5.0 空间制图要求与美化

29 / 30

双变量专题地图

A + B → C · Bivariate Choropleth

A2025 · 人 / km²

人口密度

按县级数值分为四档

Q1 低Q4 高

B2020—2025 · %

人口变化率

同一批县级单元，另分四档

Q1 低Q4 高

A×B4 × 4

美国本土 48 州、3,108 个县级单元的双变量专题地图。每县按 2025 年人口密度和 2020—2025 年人口变化率各自的四分位等级组合填色；图中含经纬网、局地比例尺和真北指示，右侧二维图例说明 16 种颜色。

READ THE COLOUR

一格颜色
读取两个等级

密度分位 ↑

Q4Q3Q2Q1

Q1Q2Q3Q4

变化率分位 →

多源数据分析与实践 · Week 5

30 / 30

Thanks

第五周：第四章重点解析

主讲人：刘贇喆副教授 | 华东师范大学社会发展学院

第四章 重点解析

1.0空间数据导论

2.0空间数据读写

3.0空间数据处理

去掉几何后，继续处理属性

4.0空间数据可视化

5.0空间制图要求与美化

底图有来源

内容按图幅核对

按用途核对

人口密度

人口变化率

Thanks

第四章重点解析