利用R绘制指定区域的DEM图

# 加载运行环境

library(tidyverse)

library(sf)

library(raster)

library(rgdal)

library(glue)

library(data.table)

library(ggsn)

library(cowplot)



## read dem file

myfile <- "00DataSet/dem_250m/dem_250m/w001001.adf"

dem <- new("GDALReadOnlyDataset", myfile) %>% asSGDF_GROD

crs.geo <- CRS("+proj=longlat +ellps=WGS84 +datum=WGS84")





## get the boundary of Three River Plain



sjpy <-

  read_sf("00DataSet/Province/CN-sheng-A.shp",

          options = "ENCODING=gb2312",

          stringsAsFactors = FALSE) %>%

  filter(name == "辽宁") %>%

  st_transform(crs = crs.geo)



bdf <-sjpy %>%   getElement("geometry")  %>%

  st_transform(crs = crs.geo) %>% as_Spatial()



## mask dem using the boundary 

dem.sj <- raster(dem)

dem.sj <- dem.sj %>% crop(bdf) %>%

  raster::aggregate(fact = 0.01 / res(d