5.4.1 震级的空间分布
地震发生在人口密度越大，经济越发达的地区，造成的破坏也越大，反之，若发生在荒无人烟的地区，地震再大，也不会造成伤亡或损失。因此，接下来要看看地震震源的空间分布，为了方便后续操作，首先用 sf 包 (Pebesma 2018) 将数据集 dat 做类型转换。这是一个空间点数据，10441 次地震观测，8 个字段，数据集 dat 前 10 次观测数据，如下：
library(sf)
# data.frame 转为 sf
dat_sf <- st_as_sf(dat, coords = c("经度(°)", "纬度(°)"), crs = 4326)
# 查看数据
dat_sf
## Simple feature collection with 10884 features and 8 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: -180 ymin: -69.9 xmax: 179.8 ymax: 86.94
## Geodetic CRS: WGS 84
## First 10 features:
## 发震时刻 震级(M) 深度(千米) 参考位置 时段 月份 年份
## 1 2022-12-30 07:26:47 5.2 10 印尼苏门答腊岛北部 07 12 2022
## 2 2022-12-29 21:13:36 5.2 120 菲律宾棉兰老岛 21 12 2022
## 3 2022-12-29 01:16:36 5.4 10 吉尔吉斯斯坦 01 12 2022
## 4 2022-12-29 00:34:20 6.0 10 洛亚蒂群岛东南 00 12 2022
## 5 2022-12-28 21:19:22 3.6 100 新疆阿克苏地区沙雅县 21 12 2022
## 6 2022-12-28 16:45:12 3.0 35 新疆克孜勒苏州阿克陶县 16 12 2022
## 7 2022-12-27 20:28:22 3.0 111 新疆喀什地区叶城县 20 12 2022
## 8 2022-12-27 11:30:10 4.1 10 青海海北州门源县 11 12 2022
## 9 2022-12-27 06:56:20 4.0 10 台湾花莲县海域 06 12 2022
## 10 2022-12-26 20:21:09 5.2 10 埃塞俄比亚 20 12 2022
## 日期 geometry
## 1 2022-12-30 POINT (97.45 -0.15)
## 2 2022-12-29 POINT (126.2 6.1)
## 3 2022-12-29 POINT (79.56 41.88)
## 4 2022-12-29 POINT (171.4 -21.3)
## 5 2022-12-28 POINT (83.12 40.16)
## 6 2022-12-28 POINT (75.89 38.88)
## 7 2022-12-27 POINT (76.27 37.32)
## 8 2022-12-27 POINT (101.3 37.74)
## 9 2022-12-27 POINT (122.5 23.84)
## 10 2022-12-26 POINT (40.1 14.5)
接下来，要将所有的点画在图上，单凭肉眼是无法从 10000 多个数据点看出什么空间分布的。
ggplot() +
geom_sf(data = dat_sf, aes(color = `震级(M)`), cex = .1) +
scale_color_viridis_c() +
theme_minimal() +
labs(
title = "2012-04-26 至 2022-08-12 的地震分布",
caption = "数据源：中国地震台网"
)
图 5.8: 世界近 10 年的地震分布
中国地震台网发布的数据除了中国境内的地震，还包含不少境外的地震，可能是仪器位于中国境内，在距离中国比较远的地方，只有较大的地震活动才能被监测到。本节主要关心中国的情况，先根据大致的经纬度范围，粗略地将中国及周边地区圈出来，观察地震活动的空间分布。
ggplot() +
geom_sf(data = dat_sf, aes(color = `震级(M)`), cex = .1) +
scale_color_viridis_c() +
coord_sf(xlim = c(73.68, 135.2), ylim = c(3.984, 53.65)) +
theme_minimal() +
labs(
title = "2012-04-26 至 2022-08-12 的地震分布",
caption = "数据源：中国地震台网"
)
图 5.9: 中国及周边近 10 年的地震分布
不看不知道，一看要吓一跳，似乎哪里都有地震，只是西北、西南比较多，小地震比较多。下面再将比较关心的 4.5 级以上的地震圈出来，如图。
ggplot() +
geom_sf(
data = subset(dat_sf, subset = `震级(M)` > 4.5),
aes(color = `震级(M)`), cex = .5
) +
scale_color_viridis_c() +
coord_sf(xlim = c(73.68, 135.2), ylim = c(3.984, 53.65)) +
theme_minimal() +
labs(
title = "2012-04-26 至 2022-08-12 的地震分布",
caption = "数据源：中国地震台网"
)
图 5.10: 中国及周边近 10 年 4.5 级及以上的地震分布
粗略来看，去掉没啥破坏力的小地震后，分布依然广泛，尤其是新疆、西藏、台湾的地震活动比较多。为了更加清晰地观察到地震在各个省或自治区的分布，从阿里云可视化平台获取来自高德开放平台的中国地图数据，在上图添加中国地理边界和各省的边界。
# 读取 GeoJSON 格式的多边形矢量地图数据
china_map <- read_sf("data/中华人民共和国.json", as_tibble = FALSE)
china_map
## Simple feature collection with 35 features and 10 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 73.5 ymin: 3.397 xmax: 135.1 ymax: 53.56
## Geodetic CRS: WGS 84
## First 10 features:
## adcode name adchar childrenNum level parent
## 1 110000 北京市
## 2 120000 天津市
## 3 130000 河北省
## 4 140000 山西省
## 5 150000 内蒙古自治区
## 6 210000 辽宁省
## 7 220000 吉林省
## 8 230000 黑龙江省
## 9 310000 上海市
## 10 320000 江苏省
## subFeatureIndex center centroid acroutes
## 1 0 116.4, 39.9 116.42, 40.19 100000
## 2 1 117.19, 39.13 117.35, 39.29 100000
## 3 2 114.50, 38.05 100000
## 4 3 112.55, 37.86 112.30, 37.62 100000
## 5 4 111.67, 40.82 114.08, 44.33 100000
## 6 5 123.4, 41.8 122.6, 41.3 100000
## 7 6 125.32, 43.89 126.2, 43.7 100000
## 8 7 126.64, 45.76 127.69, 48.04 100000
## 9 8 121.47, 31.23 121.44, 31.07 100000
## 10 9 118.77, 32.04 119.49, 32.98 100000
## geometry
## 1 MULTIPOLYGON (((117.3 40.58...
## 2 MULTIPOLYGON (((117.8 39.4,...
## 3 MULTIPOLYGON (((117.5 40.65...
## 4 MULTIPOLYGON (((110.4 34.6,...
## 5 MULTIPOLYGON (((97.17 42.8,...
## 6 MULTIPOLYGON (((123.5 39.79...
## 7 MULTIPOLYGON (((129.6 42.42...
## 8 MULTIPOLYGON (((123.6 46.22...
## 9 MULTIPOLYGON (((120.9 31.02...
## 10 MULTIPOLYGON (((117.3 34.56...
高德地图是具有甲级测绘资质的单位，一般来说，地图数据是比较权威可靠的。数据集 china_map 包含 35 个省或自治区，10 个字段，有行政编码、行政区名、市级行政区数量、区域中心、区域边界等信息。有的地理区域含有岛屿，由多个多边形组成，整个数据集是 MULTIPOLYGON 类型。
ggplot() +
geom_sf(data = china_map, color = "gray", fill = NA) +
geom_sf(
data = subset(dat_sf, subset = `震级(M)` >= 4.5),
aes(color = `震级(M)`), cex = .5
) +
scale_color_viridis_c() +
coord_sf(xlim = c(73.68, 135.2), ylim = c(3.984, 53.65)) +
theme_minimal() +
labs(
title = "2012-04-26 至 2022-08-12 的地震分布",
caption = "数据源：中国地震台网"
)
图 5.11: 中国及周边近 10 年 4.5 级以上的地震分布
添加了国界、省级地理边界后，可以清楚地看出新疆、西藏、青海、甘肃、四川、云南、台湾地震活动比较多。下面想看看能够造成较大破坏的大地震的分布情况，即 6.5 级及以上的地震分布。
ggplot() +
geom_sf(data = china_map, color = "gray", fill = NA) +
geom_sf(data = subset(dat_sf, subset = `震级(M)` >= 6.5), aes(color = `震级(M)`)) +
scale_color_viridis_c() +
coord_sf(xlim = c(73.68, 135.2), ylim = c(3.984, 53.65)) +
theme_minimal() +
labs(
title = "2012-04-26 至 2022-08-12 的地震分布",
caption = "数据源：中国地震台网"
)
图 5.12: 中国及周边近 10 年 6.5 级及以上的地震分布
过去这 10 多年里，分布在中国境内的大地震还是比较少，主要分布在西北地区。大地震的危害往往比较大，10 年发生 10 次可不能说算少了，更何况还不止 10 次呢，幸好都在人口相对较少的西部地区，但愿都是发生在人口稀少的地区，人们也能尽量避开有危险地区，比如板块交界的活跃地带。
得益于 sf 包强大的空间数据操作能力，可以圈选出经纬度位于中国地理边界内，震级在 6.5 级及以上的地震。每个省级地理区域是 MULTIPOLYGON，首先获得 MULTIPOLYGON 的凸包，
china_map_ch <- st_convex_hull(china_map)
再将每个省或自治区的凸包合并起来，组成一个更大的凸包。
china_map_union <- st_union(china_map_ch)
接着，将地震点在中国境内的数据过滤出来。
sub_dat_sf <- subset(dat_sf[china_map_union, op = st_within], subset = `震级(M)` >= 6.5)
sub_dat_sf
## Simple feature collection with 26 features and 8 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 74.04 ymin: 14.1 xmax: 121.7 ymax: 44.27
## Geodetic CRS: WGS 84
## First 10 features:
## 发震时刻 震级(M) 深度(千米) 参考位置 时段 月份 年份
## 313 2022-09-18 14:44:14 6.9 10 台湾花莲县 14 09 2022
## 324 2022-09-17 21:41:18 6.5 10 台湾台东县 21 09 2022
## 385 2022-09-05 12:52:18 6.8 16 四川甘孜州泸定县 12 09 2022
## 481 2022-07-27 08:43:27 7.0 10 菲律宾 08 07 2022
## 900 2022-03-23 01:41:38 6.6 20 台湾台东县海域 01 03 2022
## 943 2022-03-14 05:05:50 6.5 10 菲律宾群岛 05 03 2022
## 1165 2022-01-08 01:45:27 6.9 10 青海海北州门源县 01 01 2022
## 1969 2021-05-22 02:04:11 7.4 17 青海果洛州玛多县 02 05 2021
## 2813 2020-07-23 04:07:20 6.6 10 西藏那曲市尼玛县 04 07 2020
## 4257 2019-04-18 13:01:05 6.7 24 台湾花莲县海域 13 04 2019
## 日期 geometry
## 313 2022-09-18 POINT (121.3 23.15)
## 324 2022-09-17 POINT (121.2 23.05)
## 385 2022-09-05 POINT (102.1 29.59)
## 481 2022-07-27 POINT (120.5 17.7)
## 900 2022-03-23 POINT (121.5 23.45)
## 943 2022-03-14 POINT (119.4 14.1)
## 1165 2022-01-08 POINT (101.3 37.77)
## 1969 2021-05-22 POINT (98.34 34.59)
## 2813 2020-07-23 POINT (86.81 33.19)
## 4257 2019-04-18 POINT (121.7 24.02)
最后，将过滤出来的地震位置画在中国地图上，见下图。
ggplot() +
geom_sf(data = china_map, color = "gray", fill = NA) +
geom_sf(data = sub_dat_sf, aes(color = `震级(M)`)) +
scale_color_viridis_c() +
theme_minimal() +
labs(
title = "2012-04-26 至 2022-08-12 的地震分布",
caption = "数据源：中国地震台网"
)
图 5.13: 中国近 10 年来的地震分布
这样就比较清晰了，不过，有的地震正好在中国与其他国家的边界周围，凸包操作有个别误伤。