# 加载运行环境

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(dem.sj)) %>%  #设置分辨率

  mask(bdf) %>% as.data.frame(xy = TRUE) 

SitesinRegion <- readRDS("02 DEM/Liaoning/pointsinRegions.rds")

pal8 <- c("#33A02C", "#B2DF8A", "#FDBF6F", "#1F78B4", "#999999", "#E31A1C", "#E6E6E6", "#A6CEE3")

G1 <- ggplot(dem.sj) + geom_raster(aes(x, y, fill = band1), na.rm = TRUE) +

  scale_fill_gradientn(colours = pal8, na.value = "transparent")+

  geom_point(data = SitesinRegion,aes(x=X,y=Y,size="a"))+

  scale_size_manual(

    values = c(1),

    breaks = c("a"),

    labels = c("Site")

  ) +

  geom_text(data = SitesinRegion, 

            aes(x = X, y = Y, label = EnStationN),

            size=2.5,check_overlap =TRUE,nudge_y=0.1) +

  north(sjpy,location ="topleft" ,0.1)+

  ylim(38.3,43.3)+

  coord_sf(crs = crs.geo)+

  theme(

    plot.background=element_rect(fill = "transparent"),

    panel.ontop = TRUE, 

    panel.background = element_rect(color = NA, fill = NA),

    panel.grid=element_line(color = "grey",size = 0.25),

    panel.border=element_rect(color = "black",fill = NA ),

    legend.key = element_rect(fill=NA,size=0.5,color = NA),

    legend.title = element_text(size=8),

    legend.position = c(1,0),

    legend.justification = c(1,0.08),

    legend.spacing.y = unit(-0.2,"cm"),

    legend.background = element_rect(fill = "transparent"),

    plot.margin = margin(0,0,-0.2,-0.3,"cm")

  )+

  labs(fill="Elevation(m)",x="",y="",size="")+

  guides(size = guide_legend(order=2,title=NULL),

         fill=guide_colorbar(order=1,title.vjust=10,barheight =5 )  )

# China mini map ------------------------------------------------------------------------------

China <- read_sf(

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

  options = "ENCODING=gb2312",

  stringsAsFactors = FALSE

)  %>%

  st_transform(crs = crs.geo)

G.China <- ggplot(China)+geom_sf(fill="white",color="black",size=0.1)+

  geom_sf(data=sjpy,fill="blue",color="black",size=0.1)+

  coord_sf(crs = crs.geo)+theme_void()

# insert mini China map into the plot ---------------------------------------------------------

Gs <- cowplot::ggdraw()+

  cowplot::draw_grob(ggplotGrob(G.China),x = 0,y = 0,

                     scale = 0.3,

                     hjust = 0.24, vjust = 0.3)+

  draw_plot(G1)

# draw and preview ----------------------------------------------------------------------------

Taotao::preview.ggplot(Gs,width = 12,height = 10,dpi = 600)