g

gtrobot

V1

2022/04/13阅读:12主题:全栈蓝

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. 反距离加权法(IDW)空间插值

chengdu_district_sp <- chengdu_district
coordinates(chengdu_district_sp) <- ~long+lat
grid <- expand.grid(x=seq(from=st_bbox(chengdu_map)[1],to=st_bbox(chengdu_map)[3],length.out=400),
                    y=seq(from=st_bbox(chengdu_map)[2],to=st_bbox(chengdu_map)[4],length.out=400))
coordinates(grid) <- ~x+y
gridded(grid) <- TRUE
idw <- idw(formula = Value ~ 1,locations = chengdu_district_sp,newdata = grid)
idw_output <- as.data.frame(idw)
colnames(idw_output)[1:3] <- c('long','lat','IDW Result')

5. IDW地图

chengdu_district_sf <- st_as_sf(chengdu_district,coords=c('long','lat'),crs=4326)

ggplot()+
  geom_tile(data = idw_output,mapping = aes(x=long,y=lat,fill=`IDW Result`))+
  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')
IDW插值地图
IDW插值地图

6. 裁剪

idw_output_sf <- st_as_sf(idw_output,coords = c('long','lat'),crs=4326)
mask_result <- st_intersection(idw_output_sf,st_make_valid(chengdu_map))

ggplot()+
  geom_sf(data = mask_result,mapping = aes(colour=IDW.Result))+
  geom_sf(data = chengdu_map,fill='NA',colour='black')+
  scale_colour_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')

分类:

数学

标签:

数学编程

作者介绍

g
gtrobot
V1