上海疫情数据可视化展示-丁
丁思怡
2023-10-15
library(sp)## The legacy packages maptools, rgdal, and rgeos, underpinning the sp package,
## which was just loaded, were retired in October 2023.
## Please refer to R-spatial evolution reports for details, especially
## https://r-spatial.org/r/2023/05/15/evolution4.html.
## It may be desirable to make the sf package available;
## package maintainers should consider adding sf to Suggests:.library(raster)
library(sf)## Linking to GEOS 3.11.0, GDAL 3.5.3, PROJ 9.1.0; sf_use_s2() is TRUElibrary(ggplot2)
library(gganimate)##
## 载入程辑包：'gganimate'## The following object is masked from 'package:raster':
##
## animatelibrary(RColorBrewer)
library(magick)## Linking to ImageMagick 6.9.12.93
## Enabled features: cairo, fontconfig, freetype, heic, lcms, pango, raw, rsvg, webp
## Disabled features: fftw, ghostscript, x11Step1 - 在https://geo.datav.aliyun.com找到上海的地图坐标数据（json格式），然后plot出来，就得到一个上海的空地图
*顺便给大家再提供一个中国地图数据china_map = st_read(dsn = “http://xzqh.mca.gov.cn/data/quanguo_Line.geojson”, stringsAsFactors = FALSE) 是中华人民共和国民政部提供的， GeoJSON 格式的，属于权威数据源，推荐使用
shanghai_map = st_read(dsn = "https://geo.datav.aliyun.com/areas_v2/bound/310000_full.json", stringsAsFactors = FALSE)## Reading layer `310000_full' from data source
## `https://geo.datav.aliyun.com/areas_v2/bound/310000_full.json'
## using driver `GeoJSON'
## Simple feature collection with 16 features and 9 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 120.8568 ymin: 30.67559 xmax: 122.2471 ymax: 31.87272
## Geodetic CRS: WGS 84shanghai_map #展示数据格式## Simple feature collection with 16 features and 9 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 120.8568 ymin: 30.67559 xmax: 122.2471 ymax: 31.87272
## Geodetic CRS: WGS 84
## First 10 features:
## adcode name childrenNum level parent subFeatureIndex
## 1 310101 黄浦区 0 district { "adcode": 310000 } 0
## 2 310104 徐汇区 0 district { "adcode": 310000 } 1
## 3 310105 长宁区 0 district { "adcode": 310000 } 2
## 4 310106 静安区 0 district { "adcode": 310000 } 3
## 5 310107 普陀区 0 district { "adcode": 310000 } 4
## 6 310109 虹口区 0 district { "adcode": 310000 } 5
## 7 310110 杨浦区 0 district { "adcode": 310000 } 6
## 8 310112 闵行区 0 district { "adcode": 310000 } 7
## 9 310113 宝山区 0 district { "adcode": 310000 } 8
## 10 310114 嘉定区 0 district { "adcode": 310000 } 9
## center centroid acroutes
## 1 121.49032, 31.22277 121.48357, 31.21595 100000, 310000
## 2 121.43752, 31.17997 121.43940, 31.16299 100000, 310000
## 3 121.42220, 31.21812 121.38095, 31.20737 100000, 310000
## 4 121.4482, 31.2290 121.45066, 31.27082 100000, 310000
## 5 121.3925, 31.2417 121.39206, 31.25789 100000, 310000
## 6 121.49183, 31.26097 121.48544, 31.27665 100000, 310000
## 7 121.52280, 31.27076 121.52930, 31.29835 100000, 310000
## 8 121.37597, 31.11166 121.41890, 31.08721 100000, 310000
## 9 121.4899, 31.3989 121.40486, 31.39211 100000, 310000
## 10 121.25033, 31.38352 121.24439, 31.35814 100000, 310000
## geometry
## 1 MULTIPOLYGON (((121.476 31....
## 2 MULTIPOLYGON (((121.4577 31...
## 3 MULTIPOLYGON (((121.4395 31...
## 4 MULTIPOLYGON (((121.4577 31...
## 5 MULTIPOLYGON (((121.3542 31...
## 6 MULTIPOLYGON (((121.4854 31...
## 7 MULTIPOLYGON (((121.5165 31...
## 8 MULTIPOLYGON (((121.3587 30...
## 9 MULTIPOLYGON (((121.4253 31...
## 10 MULTIPOLYGON (((121.3365 31...#plot(shanghai_map)
plot(shanghai_map[1])
Step2 读取疫情数据（各区每日新增阳性患者数）
上海分区<-read.csv('/Users/sanhuai/Desktop/上海分区.csv')
上海分区#看一眼数据格式## area X2022.4.1 X2022.4.2 X2022.4.3 X2022.4.4 X2022.4.5 X2022.4.6
## 1 徐汇区 639 1042 499 1229 920 1107
## 2 松江区 493 567 263 559 796 7S1
## 3 黄浦区 260 659 824 970 658 1044
## 4 嘉定区 200 617 415 237 481 1408
## 5 虹口区 61 342 188 608 410 668
## 6 杨浦区 134 274 374 220 623 630
## 7 普陀区 245 388 321 254 4S3 1033
## 8 青浦区 87 231 232 314 3S4 470
## 9 静安区 189 324 338 50 302 545
## 10 金山区 49 60 103 52 77 79
## 11 奉贤区 161 157 119 64 144 278
## 12 崇明区 84 113 165 47 79 63
## 13 长宁区 37 100 104 33 84 350
## 14 宝山区 45 485 467 265 554 660
## 15 闵行区 1043 829 940 1381 2937 2409
## 16 浦东新区 2584 2038 3654 7071 8145 8457
## X2022.4.7 X2022.4.8 X2022.4.9 X2022.4.10
## 1 2076 1629 1150 3185
## 2 751 771 500 1825
## 3 1380 2607 548 1761
## 4 933 1505 369 1405
## 5 594 378 349 1238
## 6 601 701 1080 1877
## 7 957 1094 653 1001
## 8 493 342 557 877
## 9 381 665 602 603
## 10 129 90 68 57
## 11 92 97 124 38
## 12 262 170 171 55
## 13 852 614 757 405
## 14 414 2821 2261 1839
## 15 2257 2854 4624 3189
## 16 9050 7286 11130 6732names(上海分区)[1]<-c('name')
##这个数据的第一列是“area”但是为了和shanghi_map的列名“name”匹配，这里需要修改成“name”
#这里的一个知识点是names()函数可以直接输出表头Step3 准备等下图片的label要用到的日期数据
##从【上海分区】数据里面可以看到，第一份疫情数据的标题x2022.4.1，R不能把这样格式的数据自动识别为日期数据，使用我们手动来生成日期列表mytime
startTime <- as.Date("2022-4-01",'%Y-%m-%d')
mytime<-startTime+(0:9)（以”2022-4-01”的徐汇区为例，先画第一张图） Step4 把每日新增病例数分段，以便之后用不同的颜色填充代表新增病例数量/疫情发展的严重程度
上海分区$case = cut(as.numeric(上海分区[1:16,2]),c(1,10,100,200,500,1000,2000,100000),labels=c('<10',"10-100","100-200", "200-500","500-1000","1000-2000",'>2000'))Step5 合并疫情数据和地图数据
merge_data=merge(shanghai_map,上海分区, by="name", all.x=TRUE)**提问：为什么要使用merge函数而非直接把疫情数据新增在地图数据之后？
Step6 作图
##先产生一个基本地图
map_base <- ggplot(merge_data) + geom_sf() + xlab("longitude") + ylab("latitude")
map_base
#再在基本地图的基础上
map_full<-map_base + geom_sf(aes(fill = merge_data$case))+scale_fill_brewer(palette = "OrRd",direction = 1)+guides(fill=guide_legend(title = "Number of confirmed cases",reverse = 1))+labs(title = mytime[1])
map_full## Warning: Use of `merge_data$case` is discouraged.
## ℹ Use `case` instead.
#存储生成的疫情可视化地图
ggsave(filename = paste(mytime[1],".png",sep=''),plot = map_full,path = "/Users/sanhuai/Desktop",width = 20, height = 20, units = "cm")## Warning: Use of `merge_data$case` is discouraged.
## ℹ Use `case` instead.Step7 循环作图
for (i in (2:11))
#把case数分段
{
上海分区$case = cut(as.numeric(上海分区[1:16,2]),c(1,10,100,200,500,1000,2000,100000),labels=c('<10',"10-100","100-200", "200-500","500-1000","1000-2000",'>2000'))
###画图
map_base <- ggplot(merge_data) + geom_sf() + xlab("经度") + ylab("纬度") ##基本地图
map_full<-map_base + geom_sf(aes(fill = merge_data$case))+scale_fill_brewer(palette = "OrRd",direction = 1)+guides(fill=guide_legend(title = "Number of confirmed cases",reverse = 1))+labs(title = mytime[1])
ggsave(filename = paste(mytime[i-1],".png",sep=''),plot = map_full,path = "/Users/sanhuai/Desktop/maps",width = 20, height = 20, units = "cm")
}## Warning: Use of `merge_data$case` is discouraged.
## ℹ Use `case` instead.
## Use of `merge_data$case` is discouraged.
## ℹ Use `case` instead.
## Use of `merge_data$case` is discouraged.
## ℹ Use `case` instead.
## Use of `merge_data$case` is discouraged.
## ℹ Use `case` instead.
## Use of `merge_data$case` is discouraged.
## ℹ Use `case` instead.
## Use of `merge_data$case` is discouraged.
## ℹ Use `case` instead.
## Use of `merge_data$case` is discouraged.
## ℹ Use `case` instead.
## Use of `merge_data$case` is discouraged.
## ℹ Use `case` instead.
## Use of `merge_data$case` is discouraged.
## ℹ Use `case` instead.
## Use of `merge_data$case` is discouraged.
## ℹ Use `case` instead.Step8 制作动态地图
animate_p=image_animate(image=image_read(path=paste("/Users/sanhuai/Desktop/maps","/",mytime,".png",sep='')))
anim_save(filename = "疫情地图可视化动态图.gif",animation = animate_p,path="/Users/sanhuai/Desktop/maps")