使用 Stata 绘制长三角城市群高程地图

在之前的课程中我分享过使用 Stata 绘制长三角 26 个城市的地图的数据和代码。

使用 Stata 绘制长三角地图城市地图（带指北针和比例尺）：https://rstata.duanshu.com/#/brief/course/04578a56fede4d2d91a74ae8b92d442b

在使用本数据前请确保你已经学习过下面两个课程：

使用 Stata 绘制历年中国省级行政区划（小地图版本 + 长版）：https://rstata.duanshu.com/#/brief/course/e194d75fe1674f7c8921daec1ece7f7d

使用 Stata 绘制历年中国市级行政区划（小地图版本 + 长版）：https://rstata.duanshu.com/#/brief/course/0e9d78633ae64fefa684d6aad89633bc

在这些课程的基础上，今天我们继续讲解如何使用 Stata 绘制长三角城市群高程地图。

由于 Stata 没法处理栅格数据，所以我们还是需要先使用 R 语言初步处理下。处理的目标就是得到每个坐标点的高程数据。

附件中的 dem_500m 文件夹里面存放了 500m 分辨率的全国 DEM 数据。本来想把这个数据直接转换成 dta 格式的，不过转换结果太大了，反而没办法使用了。

下面我们使用 R 语言提取长三角城市群的高程数据并转换成 dta 格式的。

如果没办法理解这部分代码可以联系李老师帮忙提取。

加在所需 R 包：

library(tidyverse)
library(terra)

读取 DEM 栅格数据：

rast("dem_500m/w001001.adf") -> rst
rst

#> class       : SpatRaster
#> dimensions  : 10085, 17763, 1  (nrow, ncol, nlyr)
#> resolution  : 0.004031287, 0.004031287  (x, y)
#> extent      : 66.39632, 138.0041, 14.90839, 55.56392  (xmin, xmax, ymin, ymax)
#> coord. ref. : lon/lat WGS 84 (EPSG:4326)
#> source      : w001001.adf
#> name        : w001001
#> min value   :    -268
#> max value   :    8405

为了提取长三角地区的，我们需要准备一份长三角矢量数据：

library(sf)
read_sf("2021行政区划/市.shp") -> city

# 之前提供的地图数据使用的是这个 crs
mycrs <- "+proj=aea +lat_0=0 +lon_0=105 +lat_1=25 +lat_2=47 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs"

city %>%
  st_transform(mycrs) -> cityaea

readxl::read_xlsx("长三角.xlsx", sheet = 2) %>%
  pull(城市) -> csjcitylist

csjcitylist

#>  [1] "上海市"   "南京市"   "南通市"   "台州市"   "合肥市"   "嘉兴市"
#>  [7] "宁波市"   "安庆市"   "宣城市"   "常州市"   "扬州市"   "无锡市"
#> [13] "杭州市"   "池州市"   "泰州市"   "湖州市"   "滁州市"   "盐城市"
#> [19] "绍兴市"   "舟山市"   "芜湖市"   "苏州市"   "金华市"   "铜陵市"
#> [25] "镇江市"   "马鞍山市"

cityaea %>%
  filter(市 %in% csjcitylist) -> csjcity

csjcity

#> Simple feature collection with 26 features and 6 fields
#> Geometry type: MULTIPOLYGON
#> Dimension:     XY
#> Bounding box:  xmin: 1015102 ymin: 3077209 xmax: 1698065 ymax: 3776110
#> CRS:           +proj=aea +lat_0=0 +lon_0=105 +lat_1=25 +lat_2=47 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs
#> # A tibble: 26 × 7
#>    省     省代码 省类型 市       市代码 市类型                          geometry
#>  * <chr>   <dbl> <chr>  <chr>     <dbl> <chr>                 <MULTIPOLYGON [m]>
#>  1 安徽省 340000 省     安庆市   340800 地级市 (((1128423 3379673, 1128869 3379…
#>  2 安徽省 340000 省     池州市   341700 地级市 (((1221985 3344566, 1221675 3343…
#>  3 安徽省 340000 省     滁州市   341100 地级市 (((1209296 3606762, 1209300 3605…
#>  4 安徽省 340000 省     合肥市   340100 地级市 (((1122388 3525961, 1122728 3526…
#>  5 安徽省 340000 省     马鞍山市 340500 地级市 (((1242664 3482712, 1242619 3482…
#>  6 安徽省 340000 省     铜陵市   340700 地级市 (((1163089 3291762, 1162813 3291…
#>  7 安徽省 340000 省     芜湖市   340200 地级市 (((1231632 3425453, 1231657 3425…
#>  8 安徽省 340000 省     宣城市   341800 地级市 (((1283106 3407539, 1283108 3407…
#>  9 江苏省 320000 省     常州市   320400 地级市 (((1375461 3506300, 1376359 3505…
#> 10 江苏省 320000 省     南京市   320100 地级市 (((1292604 3431248, 1292348 3430…
#> # ℹ 16 more rows

# 合并
csjcity %>%
  st_union() -> csjcityall

然后我们就可以把刚刚读取的 rst 转换成 mycrs 坐标系，然后再裁减出 csjcityall 区域的部分：

rst %>%
  project(mycrs) %>%
  terra::crop(vect(csjcityall)) %>%
  terra::mask(vect(csjcityall)) -> rst2

plot(rst2)

最后我们再把这个数据转换成 dta 格式的：

library(raster)
raster(rst2) %>%
  rasterToPoints(spatial = T) %>%
  st_as_sf() -> rst2sf

bind_cols(
  rst2sf %>% st_drop_geometry(),
  rst2sf %>%
    st_coordinates() %>%
    as_tibble() %>%
    set_names("x", "y")
) %>%
  as_tibble() %>%
  rename(dem = w001001) -> rstdf

# 也就是这样的
rstdf

#> # A tibble: 1,507,495 × 3
#>        dem        x        y
#>      <dbl>    <dbl>    <dbl>
#>  1  1      1334296. 3774946.
#>  2  1.39   1334665. 3774946.
#>  3  3.81   1335035. 3774946.
#>  4  0.0920 1336512. 3774946.
#>  5 -2.15   1336881. 3774946.
#>  6  4.45   1339466. 3774946.
#>  7 -0.758  1332449. 3774577.
#>  8  0.314  1332819. 3774577.
#>  9  1.83   1333927. 3774577.
#> 10  1.01   1334296. 3774577.
#> # ℹ 1,507,485 more rows

rstdf %>%
  haven::write_dta("dempoint.dta", label = "数据处理：微信公众号 RStata")

下面我们就可以把这个数据绘制到地图上了。

为了让不同高程的点显示不同的颜色，我们需要对高程变量进行分组：

*- 添加散点
use dempoint.dta, clear
*- dem 分组，分成九组
egen group = cut(dem), group(9)
save pointdata, replace

对于文本标签，我们希望城市标签使用灰色、指北针和比例尺标签使用黑色，所以也需要生成一个分组变量：

*- 标签分组
use csjcity_label, clear
gen group = !inlist(cname, "N", "200km")
save csjcity_label2, replace

然后就可以绘制地图了：

use csjcity_db.dta, clear
grmap using csjcity_coord.dta, ///
 id(ID) osize(thin ...) ocolor(white ...) ///
 fcolor("white") ///
    leg(off) ///
 graphr(margin(medium)) ///
 polygon(data(csjpolygon) fcolor(black) ///
  osize(vvthin)) ///
 line(data(csjcity_line_coord.dta) ///
  size(*0.5) pattern(solid) ///
  color(black)) ///
 label(data(csjcity_label2) x(X) y(Y) by(group) ///
  label(cname) length(20) size(*0.7 *0.7) color(black gs12)) ///
 point(data(pointdata) x(x) y(y) by(group) ///
  size(vtiny ...) ///
  fcolor("252 251 253" "239 237 245" "218 218 235" ///
   "188 189 220" "158 154 200" "128 125 186" ///
   "106 81 163" "84 39 143" "63 0 125")) ///
 ti("使用 Stata 绘制长三角地区地形图") ///
 subti("绘制：微信公众号 RStata")
gr export pic1a.png, replace width(4800)

这种渐变色的颜色可以从这里获取：

https://tidyfriday.cn/colors/#type=sequential&scheme=BuGn&n=9

点击这里跳转到 RStata 短书平台获取附件：使用 Stata 绘制长三角城市群高程地图

相关文章

评论