g

gtrobot

V1

2022/04/14阅读:25主题:全栈蓝

R语言画上海-核密度插值

1. 我的本期总结

# 使用sf包中的st_read函数读取【阿里云数据可视化平台】共享的上海市json格式地图数据
# 使用read_csv函数读取csv格式数据【这部分数据也是网上下载的】
# 数据是上海市的银行位置,当然不全
# 网址:http://guihuayun.com/poi/
# 新重点:发现了一个小工具【百度地图poi数据查询工具】
# 使用这个工具获取上海市全部银行位置信息
# 工具地址:https://www.bilibili.com/video/BV18U4y187MV/
# 使用choose.files函数以交互方式选择文件
# 使用geom_sf函数画地图,fill属性设置填充色
# 使用geom_point函数画点图
# 使用labs函数设置坐标轴和caption
# 使用theme函数设置主题细节
# 使用axis.text设置坐标轴标签属性
# 使用panel.background设置背景属性
# 使用panel.grid设置网格属性
# 使用ggtext包的plot.caption设置caption属性
# 使用ggspatial包中的annotation_north_arrow函数设置指北针,样式north_arrow_nautical,位置location
# 使用ggspatial包中的annotation_scale函数设置比例尺,位置location
# 使用sm包中的sm_density函数计算核密度估计
# 使用expand_grid函数生成网格
# 使用colorRampPalette生成新的调色板
# 使用rev反转向量
# 使用brewer.pal获取颜色

注:总结的包或函数不一定在下文代码出现,试错过程中遇到的。

2. 加载在线数据

library(tidyverse)
library(amapGeocode)
library(sf)
library(sp)
library(sm)
library(RColorBrewer)
library(ggtext)
library(ggspatial)

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

# 这部分数据是通过网站获取的
# 网址http://guihuayun.com/poi/
# 上海市,银行
shanghai_poi <- read_csv(file = choose.files())

3. 画上海地图

ggplot()+
  geom_sf(data = shanghai_map,fill='NA')+
  geom_point(data = shanghai_poi,mapping = aes(x=long,y=lat),color='blue',alpha=0.2,size=4)+
  labs(x='',y='',
       caption = "Visualization by <span style='color:#DD6449'>gtrobot</span>")+
  theme(axis.text = element_text(color = 'black'),
        panel.background = element_blank(),
        panel.grid = element_line(color = 'grey80',size = 0.2),
        plot.caption = element_markdown(face = 'bold'))+
  annotation_north_arrow(style = north_arrow_nautical(),location='tr')+
  annotation_scale(location='bl')

4. 核密度插值

point_dens <- sm.density(x=shanghai_poi[c('long','lat')],
                         display='image',
                         ngrid=400,
                         ylim=st_bbox(shanghai_map)[c(2,4)],
                         xlim=st_bbox(shanghai_map)[c(1,3)])
# 提取坐标
density_coords <- expand.grid(x=point_dens$eval.points[,1],y=point_dens$eval.points[,2])
# 提取估计值
density_estimate <- data.frame(kde_value=array(point_dens$estimate,length(point_dens$estimate)))
# 合并成一个数据框
density_result <- data.frame(density_coords,density_estimate)
# 转换成sf格式
density_result_sf <- st_as_sf(x=density_result,coords = c('x','y'),crs=4326)

5. 核密度上海地图

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

ggplot()+
  geom_tile(data = density_result,mapping = aes(x=x,y=y,fill=kde_value))+
  scale_fill_gradientn(colors = map_color)+
  geom_sf(data = shanghai_map,fill='NA',color='grey40',size=0.7)+
  labs(x='',y='',
       caption = "Visualization by <span style='color:#DD6449'>gtrobot</span>")+
  theme(axis.text = element_text(color = 'black'),
        panel.background = element_blank(),
        panel.grid = element_line(color = 'grey80',size = 0.2),
        plot.caption = element_markdown(face = 'bold'))+
  annotation_north_arrow(style = north_arrow_nautical(),location='tr')+
  annotation_scale(location='bl')

6. 裁剪后

# 注意都是sf格式数据
intersection <- density_result_sf[st_make_valid(shanghai_map),]

ggplot()+
  geom_sf(data = intersection,mapping = aes(color=kde_value))+
  scale_color_gradientn(colors = map_color)+
  geom_sf(data = shanghai_map,fill='NA',color='grey40',size=0.7)+
  labs(x='',y='',
       caption = "Visualization by <span style='color:#DD6449'>gtrobot</span>")+
  theme(axis.text = element_text(color = 'black'),
        panel.background = element_blank(),
        panel.grid = element_line(color = 'grey80',size = 0.2),
        plot.caption = element_markdown(face = 'bold'))+
  annotation_north_arrow(style = north_arrow_nautical(),location='tr')+
  annotation_scale(location='bl')

重复1. 加载详细poi数据

# 这部分数据是通过【百度地图poi数据查询工具】下载的
# 工具地址:https://www.bilibili.com/video/BV18U4y187MV/
shanghai_detail_poi <- read_csv(file = choose.files())

重复2. 画上海地图

ggplot()+
  geom_sf(data = shanghai_map,fill=NA)+
  geom_point(data = shanghai_detail_poi,mapping = aes(x=wgs84_lng,y=wgs84_lat),color='blue',alpha=0.1)+
  labs(x='',y='',
       caption = "Visualization by <span style='color:#DD6449'>gtrobot</span>")+
  theme(axis.text = element_text(color = 'black'),
        panel.background = element_blank(),
        panel.grid = element_line(color = 'grey80',size=0.2),
        plot.caption = element_markdown(face = 'bold'))

重复3. 核密度插值

point_dens <- sm.density(x=shanghai_detail_poi[c('wgs84_lng','wgs84_lat')],
                         display='image',
                         ngrid=400,
                         ylim=st_bbox(shanghai_map)[c(2,4)],
                         xlim=st_bbox(shanghai_map)[c(1,3)])
# 提取坐标
density_coords <- expand.grid(x=point_dens$eval.points[,1],y=point_dens$eval.points[,2])
# 提取估计值
density_estimate <- data.frame(kde_value=array(point_dens$estimate,length(point_dens$estimate)))
# 合并成一个数据框
density_result <- data.frame(density_coords,density_estimate)
# 转换成sf格式
density_result_sf <- st_as_sf(x=density_result,coords = c('x','y'),crs=4326)

重复4. 核密度上海地图

ggplot()+
  geom_tile(data = density_result,mapping = aes(x=x,y=y,fill=kde_value))+
  scale_fill_gradientn(colors = map_color)+
  geom_sf(data = shanghai_map,fill='NA',color='grey40',size=0.7)+
  labs(x='',y='',
       caption = "Visualization by <span style='color:#DD6449'>gtrobot</span>")+
  theme(axis.text = element_text(color = 'black'),
        panel.background = element_blank(),
        panel.grid = element_line(color = 'grey80',size = 0.2),
        plot.caption = element_markdown(face = 'bold'))+
  annotation_north_arrow(style = north_arrow_nautical(),location='tr')+
  annotation_scale(location='bl')

重复5. 裁剪后

intersection <- density_result_sf[st_make_valid(shanghai_map),]

ggplot()+
  geom_sf(data = intersection,mapping = aes(color=kde_value))+
  scale_color_gradientn(colors = map_color)+
  geom_sf(data = shanghai_map,fill='NA',color='grey40',size=0.7)+
  labs(x='',y='',
       caption = "Visualization by <span style='color:#DD6449'>gtrobot</span>")+
  theme(axis.text = element_text(color = 'black'),
        panel.background = element_blank(),
        panel.grid = element_line(color = 'grey80',size = 0.2),
        plot.caption = element_markdown(face = 'bold'))+
  annotation_north_arrow(style = north_arrow_nautical(),location='tr')+
  annotation_scale(location='bl')

对比图

数据少:151

数据多:6078

分类:

数学

标签:

数学编程

作者介绍

g
gtrobot
V1