R-ggplot2+sf 核密度空间插值可视化绘制
R-ggplot2+sf 核密度空间插值可视化绘制
上篇推文我们介绍了使用Python的plotnine、Basemap包对空间kde插值结果进行了可视化绘制，当然也包括了具体的插值过程，详细内容大家可以点击下方链接查看：Python-plotnine 核密度空间插值可视化绘制 、Python-Basemap核密度空间插值可视化绘制。
本期推文我们就介绍下使用R进行核密度估计、空间插值计算以及ggplot2+sf的可视化绘制操作。涉及的主要知识点如下：
- R-sm包计算核密度估计结果
- R-SP包转换网格插值结果
- R-ggplot2+sf包绘制网格插值结果
- R-sf包实现完美“裁剪”
R-sm包计算核密度估计结果
sf包散点位置可视化
首先，我们需要读取的数据还是地图文件(jiangsu)和点数据(scatter_df),散点数据预览如下：
在计算核密度估计之前，我们先使用sf包进行散点的可视化绘制。首先我们使用sf包将其转换成sf地图对象，操作代码如下：
scatter_df_tro <- st_as_sf(scatter_df,coords = c("经度", "纬度"),crs = 4326)
结果如下(部分)：
接下来，我们将使用sf包进行点位置的可视化绘制，代码如下：
library(sf)
library(tidyverse)
library(ggspatial)
library(RColorBrewer)
library(ggtext)
library(hrbrthemes)
#自定义颜色
my_colormap <- colorRampPalette(rev(brewer.pal(11,'Spectral')))(32)
Map_point <- ggplot() +
geom_sf(data = jiangsu,fill="NA",size=.6,color="black") +
#coord_sf(crs = "+proj=laea +lat_0=40 +lon_0=104")+ #中国地图常用投影
geom_sf(data = scatter_df_tro,aes(fill=PM2.5),shape=21,size=5) +
scale_fill_gradientn(colours = my_colormap,)+
annotation_scale(location = "bl") +
# spatial-aware automagic north arrow
annotation_north_arrow(location = "tr", which_north = "false",
style = north_arrow_fancy_orienteering) +
labs(
title = "Map Charts in R Exercise 01: <span style='color:#D20F26'>Map point</span>",
subtitle = "processed map charts with <span style='color:#1A73E8'>geom_sf()</span>",
caption = "Visualization by <span style='color:#DD6449'>DataCharm</span>") +
theme_ipsum(base_family = "Roboto Condensed") +
theme(
plot.title = element_markdown(hjust = 0.5,vjust = .5,color = "black",
size = 20, margin = margin(t = 1, b = 12)),
plot.subtitle = element_markdown(hjust = 0,vjust = .5,size=15),
plot.caption = element_markdown(face = 'bold',size = 12),
)
可视化结果如下：
sm包计算核密度估计结果
在上述可视化结果之后，我们需要根据已有的点进行核密度估计，在R中，ks、gss、KernSmooth以及sm包都可以实现核密度估计操作，在考虑定制化设置上，我们最终选择sm包进行空间核密度计算，具体代码如下：
point_dens<- sm.density(data.frame(scatter_df$`经度`, scatter_df$`纬度`),
display= "image", ngrid=400,
ylim=st_bbox(jiangsu)[c(2,4)],
xlim=st_bbox(jiangsu)[c(1,3)])
library(sp)
Density_result<-SpatialPoints(expand.grid(x=point_dens$eval.points[,1],
y=point_dens$eval.points[,2]))
Density_result<-SpatialPixelsDataFrame(Density_result, data.frame(kde_value = array(point_dens$estimate,
length(point_dens$estimate))))
#准换成df格式
df_density<-data.frame(Density_result)
head(df_density)
结果如下(部分)：
知识点：
- sm.density()根据散点经纬度计算核密度估计结果
point_dens<- sm.density(data.frame(scatter_df$`经度`, scatter_df$`纬度`),
display= "image", ngrid=400,
ylim=st_bbox(jiangsu)[c(2,4)],
xlim=st_bbox(jiangsu)[c(1,3)])
其中：
2. sp包SpatialPoints、SpatialPixelsDataFrame方法将插值结果转网格再转空间df数据。
- 使用data.frame()将结果转成data.frame()类型便于ggplot2和sf包绘制。
R-ggplot2+sf包绘制网格插值结果
接下来，我们将上方的核密度估计结果进行可视化绘制，首先，我们绘制插值的网格结果：
#自定义颜色
my_colormap <- colorRampPalette(rev(brewer.pal(11,'Spectral')))(32)
Map_point_kde_nomask <- ggplot() +
geom_sf(data = jiangsu,fill="NA",size=.8,color="black") +
geom_tile(data = df_density,aes(x=x,y=y,fill=kde_value))+
scale_fill_gradientn(colours = my_colormap)+
annotation_scale(location = "bl") +
# spatial-aware automagic north arrow
annotation_north_arrow(location = "tr", which_north = "false",
style = north_arrow_fancy_orienteering) +
labs(
title = "Map Charts in R Exercise 01: <span style='color:#D20F26'>Map kde nomask</span>",
subtitle = "processed map charts with <span style='color:#1A73E8'>geom_title()</span>",
caption = "Visualization by <span style='color:#DD6449'>DataCharm</span>") +
theme_ipsum(base_family = "Roboto Condensed") +
theme(
plot.title = element_markdown(hjust = 0.5,vjust = .5,color = "black",
size = 20, margin = margin(t = 1, b = 5)),
plot.subtitle = element_markdown(hjust = 0,vjust = .5,size=15),
plot.caption = element_markdown(face = 'bold',size = 12),
)
可视化结果如下：
可以看到，还是出现了和我们使用Python绘制的结果一样，都是没有对感兴趣区域(地图文件)进行裁剪出来，接下里我们使用sf包进行 “裁剪” 操作。
R-sf包实现完美“裁剪”
st_intersection() 裁剪操作
sf包提供了多种方便、灵活、简单的空间数据操作函数(抽时间可以系统进行主要函数的详细讲解)，在这里，我们就使用了sf::st_intersection() 将我们感兴趣的区域进行获取，首先，我们将kde插值结果转换成sf对象类型，代码如下：
df_density_df <- st_as_sf(df_density,coords = c("x", "y"),crs = 4326)
head(df_density_df)
结果如下(部分)：
再结合之前的地图文件(jiangsu,sf包读取的geojson文件格式)，使用*st_intersection()*方法即可，如下：
mask_result <- sf::st_intersection(df_density_df,jiangsu)
geom_sf()裁剪结果可视化绘制
在获取裁剪的结果之后，我们就可以使用geom_sf()方法进行绘制了，这里注意:aes(color=kde_value)操作，具体代码如下：
Map_point_kde <- ggplot() +
geom_sf(data = mask_result,aes(color=kde_value)) +
geom_sf(data = jiangsu,fill="NA",size=.5,color="gray40") +
scale_color_gradientn(colours = my_colormap)+
annotation_scale(location = "bl") +
# spatial-aware automagic north arrow
annotation_north_arrow(location = "tr", which_north = "false",
style = north_arrow_fancy_orienteering) +
labs(
title = "Map Charts in R Exercise 01: <span style='color:#D20F26'>Map point kde</span>",
subtitle = "processed map charts with <span style='color:#1A73E8'>geom_sf()</span>",
caption = "Visualization by <span style='color:#DD6449'>DataCharm</span>") +
theme_ipsum(base_family = "Roboto Condensed") +
theme(
plot.title = element_markdown(hjust = 0.5,vjust = .5,color = "black",
size = 20, margin = margin(t = 1, b = 12)),
plot.subtitle = element_markdown(hjust = 0,vjust = .5,size=15),
plot.caption = element_markdown(face = 'bold',size = 12),
)
最终的可视化结果如下：
总结
这一篇推文我们详细介绍了R核密度估计、空间网格数据以及裁剪之后的可视化绘制结果，我们可以看出，R在操作空间数据上较Python 还是灵活下，特别是功能较为强大的sf包，此外，R在绘制地图可视化作品时，其元素添加也较为丰富(比例尺、指北针等)。接下里，我将继续使用R和Python(两个版本)， 探索空间插值应用较为广泛的方法及对应的可视化结果，感受空间可视化带给我们的视觉盛宴！希望小伙伴们能够喜欢
腾讯云开发者
