最近有朋友需要一張地形圖，碰巧自己沒接觸過，索性趁著周末琢磨了一番；下面就給各位觀眾老爺介紹一下如何如何通過ggplot2繪制一張地形圖，全部細節(jié)公開，

加載R包

由于數(shù)據(jù)處理涉及很多細節(jié)因此需要加載如下R包，別問為什么個個都有用

library(geoviz)
library(tidyverse)
library(sf)
library(terra)
library(rasterVis)
library(ggspatial)
library(rgdal)
library(rnaturalearth)
library(rnaturalearthdata)
library(raster)

作為一個陜西人，首先來介紹如何繪制一張陜西省的地形圖，下面讀入陜西省的經(jīng)緯度信息

lat=35.3940
long=109.1880
square_km=400

核心點1(DEM數(shù)據(jù)的獲取)

可以通過geoviz包的mapbox_dem函數(shù)獲取DEM數(shù)據(jù)，由于此函數(shù)需要API密鑰因此首先需要注冊賬號才能獲得密鑰，也可以使用下面我的密鑰

dem <- mapbox_dem(lat,long,square_km,
                  api="pk.eyJ1IjoiYmVueXNmIiwiYSI6ImNrczBtdWE0ajBwNjcydnBqMjRyZDdsOXkifQ.sUcMdooE7b9uQqzfrnWdSQ")

核心點2(獲取陜西省地圖)

可通過如下Datav獲取各個省的地圖

shp <- sf::read_sf("shanxi.json") # 讀入陜西省地圖
shanxi<- raster::mask(dem,shp) %>% # 將地圖與DEM數(shù)據(jù)結(jié)合
  crop(.,extent(shp))

df_shanxi <- as.data.frame(as(shanxi,"Raster"),xy=T) #格式轉(zhuǎn)換

核心點3(獲取河流信息)

通過ne_download下載世界河流信息，后續(xù)要將河流信息添加到ggplot2圖中，因此通過 st_as_sf將其轉(zhuǎn)換為sf

rivers <- ne_download(scale = 10, type = 'rivers_lake_centerlines',category = 'physical')

設(shè)置河流經(jīng)緯度信息

此處的xmin，ymin代表河流的經(jīng)緯度信息，應該根據(jù)后續(xù)繪制出來的地圖判讀，為了代碼的整潔性我將其提前放置了

rivers_cropped <- st_crop(st_as_sf(rivers), xmin = 107, xmax =111,
                          ymin = 32, ymax =39)

完成以上3大核心步驟可以說我們已經(jīng)就成功了，ggplot2繪制地形圖關(guān)鍵點就在于前3點，此前參考氣象水文科研貓的文檔，作者只介紹了繪圖上面三大核心點真是一筆帶過，在此感謝作者的簡單介紹

自定義顏色

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

ggplot2數(shù)據(jù)可視化

ggplot()+
  geom_sf(data=shp,aes(fill = NULL))+
  annotation_scale(location = "bl") + # 設(shè)置距離刻度尺
  annotation_north_arrow(location="tl",
                         style = north_arrow_nautical(
                           fill = c("grey40","white"),
                           line_col = "grey20"))+  # 添加指北針
  geom_tile(data=df_shanxi,aes(x=x,y=y,fill=layer),show.legend = F)+
  geom_sf(data = rivers_cropped,col='blue',size=1)+ # 添加河流信息
  scale_fill_gradientn(colours=colors,na.value="transparent")+ 
  labs(x=NULL,y=NULL)+
  geom_sf(data=shp,fill="NA",size=1,color="black")+ # 添加地圖邊界
  theme_bw()+
  theme(panel.grid.major=element_blank(),
        panel.grid.minor=element_blank(),
        panel.background = element_blank(),
        legend.title = element_blank())

喜歡的小伙伴歡迎關(guān)注我的公眾號

R語言數(shù)據(jù)分析指南，持續(xù)分享數(shù)據(jù)可視化的經(jīng)典案例及一些生信知識，希望對大家有所幫助