g
gtrobot
V1
2022/04/13阅读:40主题:全栈蓝
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')

作者介绍
g
gtrobot
V1