洗衣机的博客

library(elevatr)  # Get rasterlay from AWS by `get_elev_raster` fucntion
library(raster)  # Manipulate RasterLayer object
library(tidyverse)  # Tidy is everything.

# Set the extent we want to plot
ext_sample <- extent(70, 105, 25, 45)

# Preparing for getting the elevation raster data, make a blank RasterLayer,
# becasue the first parameter of get_elev_raster is a target Rasterlayer.
bg_init <- raster(ext = ext_sample, resolution = 0.01)
# Get elevation raster with zoom 5, then only keep the extend we want to plot
# later.
bg_rst <- get_elev_raster(bg_init, z = 5) %>% crop(ext_sample)

# Let's check the detail of bg_rst, the Background RasterLayer
bg_rst

library(rgdal)
# Read the SpatialPolygon File from DBATP_Polygon.shp
tp_ext <- readOGR("DBATP/DBATP_Polygon.shp")

# Keep the shpae of the Tibetan Plateau
tp_rst <- mask(bg_rst, tp_ext)

# Create the list of landmarks which we want to mark
city_ls <- list(x = c(91.1, 86.925278, 101.7781), y = c(29.65, 27.988056, 36.6169), 
    label = c("Lhasa", "Qomolangma", "Xi'Ning"))
str(city_ls)

# Read point with latitude and longitude. This operation is not needed for
# everyone, actually it depends on the actual data structure.
hbt_coord <- read_rds("hbt_coord.rds") %>% mutate(Ecosystem = ifelse(hbt == "M", 
    "Alpine Meadow", "Alpine Steppe"))
str(hbt_coord)

# scales package provide rescale function which can convert the range of numbers
# list to another range.
library(scales)

# First convert RasterLayer as Data.Frame with xy coordinate system.  Then
# rescale the elevaion to alpha, as the background part, super high alpha value
# is not a good idea, which range is the best? It depends by the actual. Save the
# alpha value with colname 'alpha'
bg_rst_df <- as.data.frame(bg_rst, xy = TRUE) %>% mutate(alpha = rescale(layer, to = c(0.25, 
    0.75)))

# NA will be generated by the mask function, if use NA and evelation as the alpha
# of Topographic figure, it will be dizzy, and for color Topographic figure
# please don't use the greyscale and color for the evelation simultaneously. Just
# use alpha to control the shape of regional shpae.
tp_rst_df <- as.data.frame(tp_rst, xy = TRUE) %>% na.omit()

# sacle_parm as a parameter controls the sclae of the hole figure, it will be
# used to control the size of text, point, legend, etc. to fit the size of
# figure.
scale_parm <- 2
# Init ggplot
# plot the backgroun layer, set the alpha without color will make a grey
# background
gmap <- ggplot() + 
# plot the topographic layer, set alpha to keep the shape of Tibetan Plateau
# (TP). Color indicates the elevation.
geom_raster(data = bg_rst_df, aes(x = x, y = y, alpha = alpha)) + 
 # terrain.colors is an built-in function to generate a list of color palettes.
# set the legend title of evelation by name parameter.
geom_raster(data = tp_rst_df, aes(x = x, y = y, fill = layer, alpha = alpha)) +
 # As we said before, the alphas is used to determine the shpae of TP, we don't
# need to show them as legends.
scale_fill_gradientn(colours = terrain.colors(100), name = "Elevation (m)") +
 # Project this figure as a map but not a normal figure
scale_alpha(guide = "none") +
 # Set preset theme makes things easire
coord_quickmap() +
 # Set the limititions of axes. `expand` parameter will remove the gaps between
theme_minimal() +
 # Set the titles of axis
scale_x_continuous(limits = c(70, 105), expand = c(0, 0)) + scale_y_continuous(limits = c(25, 
    45), expand = c(0, 0)) +
labs(x = "Longtitude (E)", y = "Laitude (N)") +
 # remove the background color and background grid, you know the classical
# ggplot's grid, don't you?
theme(panel.grid = element_blank(), panel.background = element_blank()) + 
# Set the size of axis and legend
theme(axis.title = element_text(size = 7 * scale_parm), axis.text = element_text(size = 6 * 
    scale_parm)) + theme(legend.key.width = unit(0.2 * scale_parm, "cm"), legend.key.height = unit(0.5 * 
    scale_parm, "cm"), legend.text = element_text(size = 5 * scale_parm), legend.title = element_text(size = 6 * 
    scale_parm))

# Preview will slow down the process of operations, I highly recommand do not
# preveiw the ggplot and save it as a file directly.

gmap <- gmap +
  # Add the city_ls to the main plot as landmarks.
  geom_text(
    mapping = aes(x = x, y = y, label = label),
    # geom_text don't support the structure we used. 
    # convert the list into data.frame, every element is used as column here.
    data = bind_cols(city_ls),
    size = 2 * scale_parm
  )

gmap <- gmap +
  # Add the sample sites to the main plot as points.
  # Due to the copyright of my scholar article and data share policy, I won't point my sample sites to the picuture, the coordinations of point is 0,0 here.
  geom_point(
    mapping = aes(
      # x = lon,
      x = 0,
      y = 0,
      # y = lat,
      col = Ecosystem,
      shape = Ecosystem
    ),
    data = hbt_coord,
    # Don't set the size as 0 until you don't want to see anything here.
    # Don't set the size as 0 until you don't want to see anything here.
    # Don't set the size as 0 until you don't want to see anything here.
    size = 0 * scale_parm
  ) +
  # Convert the color as legend class, becasuse the shape legend is legend class.
  # If there is no class conversion, the shape and color will be showed as two legends.
  # Then select colour as guides and site a larger size to make it more readable.
  # Don't know what's these means? Commit below code will show you everything.
  scale_color_discrete(guide = "legend") +
  guides(colour = guide_legend(override.aes = list(size = .8 * scale_parm)))

gmap <- gmap +
  theme(plot.margin =
          # Set marigns of figure, the order of parameters is top, right, bottom, left
          unit(
            c(0 * scale_parm, 0 * scale_parm, -.2 * scale_parm, .2 * scale_parm),
            "cm"
          ))

ggsave(filename = "gmap.pdf", plot = gmap, width = 9 * scale_parm, height = 6.2 * 
    scale_parm, units = "cm", dpi = 600)