R语言画成都地图-克里金空间插值

1. 我的本期总结

# 裁剪地图一直失败，提示有无效的球面几何，成功解决了，就这么牛
# 使用st_make_valid函数将chengdu_map的球面几何变得有效
# 然后正常裁剪就行了，使用方法2更快，两个sf数据
# 使用sf包中的st_read函数读取阿里云数据平台的json格式地图文件
# 使用amapGeocode包的getAdmin函数读取成都市各区的经纬度
# 使用key属性读取高德地图数据
# 使用colnames函数更改列名
# 使用runif函数生成n个均匀分布的随机数
# 使用RColorBrewer包中的brewer.pal获取颜色
# 使用rev函数反转向量
# 使用colorRampalette函数生成新的色板
# 使用geom_sf函数画地图
# 使用fill属性设置填充色为NA
# 使用geom_point函数画点图
# 使用scale_fill_gradientn函数设置离散填充色
# 使用geom_tile函数画矩形

注：总结的东西不一定出现在下文的代码中，「试错过程」中出现的。

2. 读取在线地图

library(tidyverse)
library(sf)
library(ggspatial)
library(amapGeocode)
library(gstat)
library(RColorBrewer)
library(ggtext)

chengdu_map <- st_read('https://geo.datav.aliyun.com/areas_v3/bound/510100_full.json')[c('adcode','name','geometry')]

# 先申请高德开放平台地图API开发者，获得API Key，然后就能读取数据了
# 读取数据详细过程参考：https://zhuanlan.zhihu.com/p/264281505
chengdu_district <- getAdmin(keywords = '成都市',key = 'f8c0968a7636416a2f58587c307bf231')
colnames(chengdu_district)[1] <- 'long'
chengdu_district$Value <- runif(n=20,min = 10,max = 90)

3. 画成都地图

my_colourmap <- colorRampPalette(rev(brewer.pal(n=11,name = 'Spectral')))(32)

ggplot()+
  geom_sf(data = chengdu_map,fill='NA',colour='black')+
  geom_point(data = chengdu_district,mapping = aes(x=long,y=lat,fill=Value),shape=21,size=5)+
  scale_fill_gradientn(colours = my_colourmap)+
  labs(caption = "Visualization by <span style='color:#DD6449'>gtrobot</span>")+
  theme(axis.title = element_blank(),
        axis.text = element_text(colour = 'black'),
        panel.background = element_rect(fill = 'NA'),
        panel.border = element_blank(),
        panel.grid.major = element_line(colour = 'grey80'),
        plot.caption = element_markdown(face='bold'))+
  annotation_north_arrow(style = north_arrow_nautical(),location='tr')+
  annotation_scale(location='bl')

4. 克里金空间插值

# sf格式数据
semivariog <- variogram(object = Value~1,locations = chengdu_district_sf,data = chengdu_district_sf)
plot(semivariog)

model_variog <- vgm(psill = 600,model = 'Lin',nugget = 0,range = 20)
fit_variog <- fit.variogram(object = semivariog,model = model_variog)
plot(semivariog,fit_variog)

# sp格式数据
krige <- krige(formula = Value~1,locations = chengdu_district_sp,newdata=grid,model=model_variog)

krige_output <- as.data.frame(krige)
colnames(krige_output)[1:3] <- c('long','lat','OK pred')

5. 克里金成都地图

ggplot()+
  geom_tile(data = krige_output,mapping = aes(x=long,y=lat,fill=`OK pred`))+
  geom_sf(data = chengdu_map,fill='NA',colour='black')+
  geom_sf(data = chengdu_district_sf,mapping = aes(fill=Value),shape=21,size=5)+
  scale_fill_gradientn(colours = my_colourmap)+
  labs(caption = "Visualization by <span style='color:#DD6449'>gtrobot</span>")+
  theme(axis.title = element_blank(),
        axis.text = element_text(colour = 'black'),
        panel.background = element_rect(fill = 'NA'),
        panel.border = element_blank(),
        panel.grid.major = element_line(colour = 'grey80'),
        plot.caption = element_markdown(face='bold'))+
  annotation_north_arrow(style = north_arrow_nautical(),location='tr')+
  annotation_scale(location='bl')

6. 裁剪后

krige_output_sf <- st_as_sf(x=krige_output,coords = c('long','lat'),crs=4326)
mask_result_krige <- krige_output_sf[st_make_valid(chengdu_map),]

ggplot()+
  geom_sf(data = mask_result_krige,mapping = aes(colour=`OK pred`))+
  scale_colour_gradientn(colours = my_colourmap)+
  geom_sf(data = chengdu_map,fill='NA',colour='black')+
  geom_sf(data = chengdu_district_sf,mapping = aes(fill=Value),shape=21,size=5)+
  scale_fill_gradientn(colours = my_colourmap)+
  guides(fill='none')+
  labs(caption = "Visualization by <span style='color:#DD6449'>gtrobot</span>")+
  theme(axis.title = element_blank(),
        axis.text = element_text(colour = 'black'),
        panel.background = element_rect(fill = 'NA'),
        panel.border = element_blank(),
        panel.grid.major = element_line(colour = 'grey80'),
        plot.caption = element_markdown(face='bold'))+
  annotation_north_arrow(style = north_arrow_nautical(),location='tr')+
  annotation_scale(location='bl')