使用 GEOJSON 文件绘制县级和市级中国地图

Posted 2021-04-30 君泉计量

这篇教程的内容比较。。。丧心病狂。旨在演示灵活组合 geojson 文件绘制复杂的区域地图。我首先是将 34 个省级的 GEOJSON 文件组合起来绘制了市级中国地图，然后又将 344 个市级的 GEOJSON 文件组合起来绘制了县级中国地图。当然绘制县级和市级中国地图的最好办法是使用之前我提供的 shp 文件。

除此之外，本文还介绍了 here 包 和 fs 包的使用，here 包可以方便的管理多个文件夹；fs 包提供了一些方便文件索引的函数。

最后我还介绍了一个新的下载 GEOJSON 数据的网址。

绘制市级中国地图

首先是导入一些需要的 R 包和设置字体：

library(geojsonsf)library(sf)library(ggplot2)library(fs)library(here)library(hrbrthemes)library(cowplot)# devtools::install_github('czxa/hotkeys')# hotkeys 是我写的一个自用的 R 包，可以用于设定快捷键，这里使用这个包里面的 `%w/o%` 操作符library(hotkeys)cnfont = "STLibianTC-Regular"enfont = "CascadiaCode-Regular"

然后是将 34 个省级的 geojson 文件读取合并。

%w/o% 是放置在我编写的 hotkeys 包里面的一个函数，这个函数可以用于从某个向量里剔去指定元素，例如：

c(1, 2, 3, 4, 5) %w/o% 1
## [1] 2 3 4 5

dir_ls() 返回的是一个文件名向量：

setwd(here::here("d3js-geojson/china/geometryProvince"))city <- dir_ls()df <- st_read('71.json')
## Reading layer `71' from data source `/Users/czx/Desktop/使用 GEOJSON 文件绘制县级和市级中国地图/d3js-geojson/china/geometryProvince/71.json' using driver `GeoJSON'## Simple feature collection with 22 features and 2 fields## geometry type: MULTIPOLYGON## dimension: XY## bbox: xmin: 118.209 ymin: 21.90461 xmax: 122.0798 ymax: 26.38874## epsg (SRID): 4326## proj4string: +proj=longlat +datum=WGS84 +no_defs
for(i in city %w/o% '71.json'){ tempdf <- st_read(i, quiet = TRUE) %>%  select(id, name, geometry) df <- rbind(df, tempdf)}
ggplot() +  geom_sf(data = df, aes(geometry = geometry, fill = factor(name)), size = 0.05, color = "white") +  scale_fill_viridis_d() + theme_modern_rc(base_family = cnfont, subtitle_family = cnfont, caption_family = cnfont) +  worldtilegrid::theme_enhance_wtg() +  labs(title = "使用 GEOJSON 数据绘制市级中国地图", subtitle = "R Packages: geojsonsf + ggplot2 + sf", caption = "GEOJSON 数据来源：ufoe/d3js-geojson @ GitHub\n<https://github.com/ufoe/d3js-geojson/>") +  theme(legend.position = "none") +  coord_sf(crs = "+proj=laea +lat_0=23 +lon_0=113 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs")

绘制县级中国地图

绘制县级中国地图的方法是一样的，只不过需要合并 344 个 geojson 文件：

setwd(here::here('d3js-geojson/china/geometryCouties'))county <- dir_ls()# 绘制市级中国地图df <- st_read('710000.json')
## Reading layer `710000' from data source `/Users/czx/Desktop/使用 GEOJSON 文件绘制县级和市级中国地图/d3js-geojson/china/geometryCouties/710000.json' using driver `GeoJSON'## Simple feature collection with 20 features and 2 fields## geometry type: MULTIPOLYGON## dimension: XY## bbox: xmin: 119.3143 ymin: 21.88912 xmax: 123.6873 ymax: 25.93604## epsg (SRID): 4326## proj4string: +proj=longlat +datum=WGS84 +no_defs
for(i in county %w/o% '710000.json'){ tempdf <- st_read(i, quiet = TRUE) %>%  select(id, name, geometry) df <- rbind(df, tempdf)}
ggplot() +  geom_sf(data = df, aes(geometry = geometry, fill = factor(name)), size = 0.05, color = "white") +  scale_fill_viridis_d() + theme_modern_rc(base_family = cnfont, subtitle_family = cnfont, caption_family = cnfont) +  worldtilegrid::theme_enhance_wtg() +  theme(legend.position = "none") +  coord_sf(crs = "+proj=laea +lat_0=23 +lon_0=113 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs")

使用 shp 文件绘制市级和县级中国地图：

当然最好的办法是用我之前提供的 shp 文件：

# 市级地图setwd(here())mapdata <- read_sf(here("chinamap/中国市界.shp"))mapborder <- read_sf(here("chinaboundary/china_official_boundary.shp"))ggplot() + geom_sf(data = mapdata, aes(geometry = geometry, fill = NAME), size = 0.05, color = "white") +  geom_sf(data = mapborder, aes(geometry = geometry), size = 0.05, color = "white") +  scale_fill_viridis_d() +  hrbrthemes::theme_ft_rc(base_family = enfont) +  worldtilegrid::theme_enhance_wtg() +  theme(legend.position = "none")

# 县级地图mapdata <- read_sf(here("chinamap/中国县界.shp"))ggplot() + geom_sf(data = mapdata, aes(geometry = geometry, fill = NAME), size = 0.05, color = "white") +  geom_sf(data = mapborder, aes(geometry = geometry), size = 0.05, color = "white") +  scale_fill_viridis_d() +  hrbrthemes::theme_ft_rc(base_family = enfont) +  worldtilegrid::theme_enhance_wtg() +  theme(legend.position = "none")

从 Highcharts 网站获取 geojson 文件

wdf <- st_read('world-palestine-highres.geo.json')
## Reading layer `world-palestine-highres.geo' from data source `/Users/czx/Desktop/使用 GEOJSON 文件绘制县级和市级中国地图/world-palestine-highres.geo.json' using driver `GeoJSON'## Simple feature collection with 215 features and 15 fields## geometry type: MULTIPOLYGON## dimension: XY## bbox: xmin: -999 ymin: 4669 xmax: 9851 ymax: 9851## epsg (SRID): 54003## proj4string: +proj=mill +lat_0=0 +lon_0=0 +x_0=0 +y_0=0 +R_A +datum=WGS84 +units=m +no_defs
ggplot() + geom_sf(data = wdf, aes(geometry = geometry, fill = name), size = 0.05, color = "white") +  scale_fill_viridis_d() +  hrbrthemes::theme_ft_rc(base_family = enfont) +  worldtilegrid::theme_enhance_wtg() +  theme(legend.position = "none")

再例如中国地图：

cndf <- st_read('cn-all-sar-taiwan.geo.json')
## Reading layer `cn-all-sar-taiwan.geo' from data source `/Users/czx/Desktop/使用 GEOJSON 文件绘制县级和市级中国地图/cn-all-sar-taiwan.geo.json' using driver `GeoJSON'## Simple feature collection with 73 features and 22 fields## geometry type: MULTIPOLYGON## dimension: XY## bbox: xmin: -999 ymin: 927 xmax: 9851 ymax: 9851## epsg (SRID): 3415## proj4string: +proj=lcc +lat_1=18 +lat_2=24 +lat_0=21 +lon_0=114 +x_0=500000 +y_0=500000 +ellps=WGS72 +towgs84=0,0,1.9,0,0,0.814,-0.38 +units=m +no_defs
ggplot() + geom_sf(data = cndf, aes(geometry = geometry, fill = name), size = 0.05, color = "white") +  scale_fill_viridis_d() +  hrbrthemes::theme_ft_rc(base_family = enfont) +  worldtilegrid::theme_enhance_wtg() +  theme(legend.position = "none")

不过这个中国地图没有九段线，不是很完整（不过至今我还没有找到带有九段线的 GEOJSON 文件）。