來源 | CSDN 博客

作者 |?萬里寫入胸懷間

責編 |?Carol

出品 | CSDN云計算（ID：CSDNcloud）

疫情牽動大家，除了做好分內(nèi)工作，管好自己不給社會添亂，也就是只能持續(xù)關注疫情了。現(xiàn)在各大門戶平臺都上線了疫情實時地圖，但基本上數(shù)據(jù)都只展示到省這個級別。

當然現(xiàn)在各地存在感染病例的小區(qū)也都及時披露了，相信疫后應該會有高精度的傳染路徑、傳染擴散分布的分析。就目前而言，丁香園上的數(shù)據(jù)應該是目前聚合數(shù)據(jù)里面精度最高的了，可以獲取到設區(qū)市級別的數(shù)據(jù)。但丁香園的數(shù)據(jù)存在地名不規(guī)范的問題，以及省、市、縣三級行政區(qū)劃并存的問題，不進行清洗不便做圖。

底圖制備

一張合適的底圖是制作地圖的基礎。與行政區(qū)劃相關的底圖，最權威的數(shù)據(jù)源當然是從自然資源部的數(shù)據(jù)中心下載相關矢量數(shù)據(jù)再進行加工，最近微信公眾號數(shù)讀城市就在做這個浩大的工程。

我的方法是用高德API獲取行政區(qū)劃邊界數(shù)據(jù)然后自行制圖。高德數(shù)據(jù)的優(yōu)點在于可以獲取從全國、省、市、區(qū)（直轄市可以到鎮(zhèn)）各級別行政區(qū)劃邊界的數(shù)據(jù)，對于不同級別行政區(qū)數(shù)據(jù)并存的情況，可以按需靈活組合，快速生成需要的底圖。

丁香園公布的數(shù)據(jù)以設區(qū)市為主，臺灣、香港、澳門公布總體數(shù)字，縣（市）的數(shù)據(jù)較少。根據(jù)數(shù)據(jù)特點，我決定省級行政區(qū)劃就保留港、澳、臺三地，其余31個省、市、自治區(qū)保留設區(qū)市、自治州、省代管縣級市以及直轄市的區(qū)、縣。也就是高德API中，各省、自治區(qū)、直轄市子一級的行政單元。我是用sf包進行處理的，底圖對象叫china_city數(shù)據(jù)框結構如下。

adcode即城市編碼，下一步里從丁香園上下載的數(shù)據(jù)里有個locationId,和它是同一個意思。高德的adcode基于民政部的行政區(qū)劃代碼，并在民政部代碼的基礎上加了市轄區(qū)一級的代碼，會干擾后續(xù)的處理，需要把它清理掉。

Amap_adcode=readxl::read_xlsx('AMap_adcode_citycode.xlsx') %>% select(name=中文名,adcode=adcode) %>% mutate_at('adcode',as.numeric) %>% filter(!grepl('轄區(qū)',name))

此外，坐標我也轉換成wgs84。

疫情數(shù)據(jù)讀取

方法來自于《R語言獲取丁香園疫情數(shù)據(jù)動態(tài)網(wǎng)站》這篇帖子。這里就簡單展示一下代碼。

library(pacman) p_load(sf,tidyverse,rvest,jsonlite) setwd() #這里要自行設定工作路徑 #read web system('~/2019_nConv/phantomjs-2.1.1-windows/phantomjs-2.1.1-windows/bin/phantomjs.exe my_js_dxy.js') url='dxy_web.html' withjs=read_html(url) %>% html_nodes('script#getAreaStat') %>% html_text() mydata=str_remove_all(withjs,'try?\\{?window.getAreaStat?=?')?%>%?str_remove_all('\\}catch\\(e\\)\\{\\}')%>%?fromJSON()

《R語言獲取丁香園疫情數(shù)據(jù)動態(tài)網(wǎng)站》：

https://blog.csdn.net/yuanzhoulvpi/article/details/104099382

以上和鏈接帖子里完全一樣，但我要獲取的是城市的數(shù)據(jù)，所以還要再加一行：

city_data=mydata$cities?%>%?rlist::list.rbind()

清洗數(shù)據(jù)

清洗數(shù)據(jù)目的只有一個，那就是把citydata的locationId轉換成底圖數(shù)據(jù)adcode中包含的。這樣就可以鏈接兩張表并進行作圖了。主要有兩個方面。

1.手動糾正locationId值為（0，-1）地區(qū)的locationId值。

這部分數(shù)據(jù)的產(chǎn)生有幾個類別。

• 新疆生產(chǎn)建設兵團的數(shù)據(jù)。新疆生產(chǎn)建設兵團數(shù)據(jù)是以師團為單位上報的。由于師團往往分布在好幾個省下一級行政區(qū)中，我就以師部所在地所屬的行政區(qū)為準。比如第八師師部所在地石河子市，屬于阿勒泰地區(qū)，則將這一行的locationId改為阿勒泰地區(qū)的行政區(qū)劃代碼654300。

• 無行政區(qū)劃代碼的開發(fā)區(qū)、高新區(qū)等。如重慶的兩江新區(qū)、南昌的贛江新區(qū)等。部分高新區(qū)涉及到不同的行政區(qū)，這里就不詳加考證，劃入組成排序第一的城區(qū)。

• 因地名不規(guī)范而沒有提供區(qū)劃代碼的地區(qū)。如云南大理州、楚雄州等，也需要進行手動糾正。

要注意的是，這部分數(shù)據(jù)隨著每天新增感染病例在空間上的分散而不斷變化，以下代碼僅針對2月13日中午丁香園公布的數(shù)據(jù)。此外，還剩下一些數(shù)據(jù)因為地區(qū)有待明確，或被列為“外地來京”“外地來滬”等類型的病例，無法明確空間上的歸屬，出于嚴謹?shù)目紤]，這部分數(shù)據(jù)在制作分布圖的時候就不采用了。用city_data=city_data %>% filter(locationId>100)剔除掉。

# #change adcode manually { city_data$locationId[grepl('第四師',city_data$cityName)]=654300city_data$locationId[grepl('第六師',city_data$cityName)]=654300city_data$locationId[grepl('第七師',city_data$cityName)]=654000city_data$locationId[grepl('第八師',city_data$cityName)]=654300city_data$locationId[grepl('第九師',city_data$cityName)]=654200city_data$locationId[grepl('第十二師',city_data$cityName)]=650100 city_data$locationId[grepl('寧東',city_data$cityName)]=640100city_data$locationId[grepl('楊凌',city_data$cityName)]=610400city_data$locationId[grepl('兩江新區(qū)',city_data$cityName)]=500105city_data$locationId[grepl('高新區(qū)',city_data$cityName)]=500106city_data$locationId[grepl('保亭',city_data$cityName)]=469029city_data$locationId[grepl('昌吉',city_data$cityName)]=652300city_data$locationId[grepl('酉陽',city_data$cityName)]=500242city_data$locationId[grepl('楚雄',city_data$cityName)]=532300city_data$locationId[grepl('文山',city_data$cityName)]=532600city_data$locationId[grepl('贛江',city_data$cityName)]=360100city_data$locationId[grepl('萬盛',city_data$cityName)]=500110city_data$locationId[grepl('德宏',city_data$cityName)]=533100city_data$locationId[grepl('大理',city_data$cityName)]=532900city_data$locationId[grepl('紅河',city_data$cityName)]=532500city_data=city_data %>% filter(locationId>100) }

2.把縣（市）locationId轉換為上級行政區(qū)的locationId

這里的縣（市）指的是受設區(qū)市或自治區(qū)管轄的縣（市），省直管縣級市不在此列。經(jīng)過轉換之后，就可以在下一步的分類匯總中，將縣（市）數(shù)據(jù)匯總到上級行政區(qū)中。他們的轉換需要兩個步驟。

• 怎么識別縣（市）？經(jīng)過上一步驟，所有的地區(qū)都有了合法的行政區(qū)劃代碼。那么，代碼不在底圖中的地區(qū)，就是我們需要轉換的縣（市）了。于是，代碼就是row.no=which(!city_data$locationId %in% china_city$adcode)

• 其次，怎么轉換。比如我們得到的數(shù)據(jù)是石家莊市橋西區(qū)，其行政區(qū)劃代碼是130104，我們希望把它換成上級行政區(qū)也就是石家莊市的代碼130100。那我們就要在Amap_adcode這張表里，以橋西區(qū)所在行為基準，往上查找距其最近的可以被100整除的（如果用字符型的話，那就是距其最近的代碼最后兩位是’00’）的代碼。查找和替換的核心語句如下:

posi=grep(city_df$cityName[i],Amap_adcode$name) while(Amap_adcode$adcode[posi]%%100!=0) posi=posi-1 city_df$locationId[i]=Amap_adcode$adcode[posi]

該部分完整的代碼如下：

#R語言中行和列的待遇大不一樣，好像還沒有像select和mutate對列操作那樣，對行在選擇的同時進行操作的函數(shù)。所以這里就專門寫了一個函數(shù)對行進行操作。 row_change=function(dataframe,row.no,fun){df_unfilter=dataframe[-row.no,]df_filter=dataframe[row.no,] %>%fun() %>%rbind(df_unfilter)return(df_filter) } #定義一個對縣locationId進行向上替換的函數(shù) upper_class_adcode=function(city_df) { for (i in 1:nrow(city_df)){if(city_df$locationId[i]>100){posi=grep(city_df$cityName[i],Amap_adcode$name)while(Amap_adcode$adcode[posi]%%100!=0) posi=posi-1 city_df$locationId[i]=Amap_adcode$adcode[posi]}}return(city_df) } #得到縣（市）數(shù)據(jù)所在行 row.no=which(!city_data$locationId %in% china_city$adcode) #利用以上兩個自定義函數(shù)進行替換 city_data_new=city_data %>% row_change(row.no,upper_class_adcode) #從原始數(shù)據(jù)中單獨提取港澳臺數(shù)據(jù) GAT=mydata[lapply(mydata$cities,is_empty) %>% unlist,] %>% select(cityName=provinceName,currentConfirmedCount,confirmedCount,suspectedCount,curedCount,deadCount,locationId) #合并形成完整的城市疫情數(shù)據(jù) city_data_new?=city_data_new??%>%?rbind(GAT)

分類匯總、生成包含數(shù)據(jù)的sf對象

得到完整的城市疫情數(shù)據(jù)之后，以locationId為分類變量進行匯總，就可以把同一設區(qū)市的縣級數(shù)據(jù)相加，形成設區(qū)市的匯總數(shù)據(jù)。鑒于丁香園上的地名并不規(guī)范，我們通過與Amap_code表格進行聯(lián)結，用規(guī)范的地名替換原地名。

city_data_sum=city_data_new %>%select(-cityName)%>% group_by(locationId) %>% summarise_all(sum) %>% left_join(Amap_adcode,by=c('locationId'='adcode'))

再下一步，將底圖的sf對象與匯總后的城市疫情數(shù)據(jù)鏈接，就可以生成包含數(shù)據(jù)的sf對象。

city_data_sf=china_city %>%left_join(,city_data_sum,by=c('adcode'='locationId'))%>%st_as_sf() city_data_sf$confirmedCount[is.na(city_data_sf$confirmedCount)]=0??????????????#將未報告患病人數(shù)的城市確診數(shù)設為0

作圖

這個方面，大神太多，我也是用的其他帖子《2019-nCoV疫情地圖動態(tài)可視化》。基本上原封不動的抄過來。在我的電腦上，顯示中文字體時R語言會崩潰，需要使用showtext包中的showtext_auto這個函數(shù)。

mybreaks <- c(0, 1, 10, 50, 100, 500, 1000, 5000, 100000) mylabels <- c("0", "1-9", "10-49", "50-99", "100-499","500-999", "1000-4999", ">=5000") city_data_sf=mutate(city_data_sf,conf_level=cut(confirmedCount,breaks=mybreaks,labels=mylabels,include.lowest = T,right=F,ordered_result = T)) ggplot(city_data_sf) + geom_sf(aes(fill = conf_level)) + coord_sf() + scale_fill_brewer(palette = "YlOrRd",) + guides(fill = guide_legend(title = "確診病例", reverse = T)) + labs(title = "2019-ncov確診患者空間分布",subtitle = Sys.Date(),caption = "數(shù)據(jù)來自丁香園") + theme(# 標題plot.title = element_text(face = "bold", hjust = 0.5,color = "black"),plot.subtitle = element_text(face = "bold", hjust = 0.5, size = 20,color = "red"),plot.caption = element_text(face = "bold", hjust = 1,color = "blue"),# 圖例legend.title = element_text(face = "bold",color = "black"),legend.text = element_text(face = "bold",color = "black"),legend.background = element_rect(colour = "black"),legend.key = element_rect(fill = NA), # 圖???????薇???legend.position = c(0.85, 0.2),axis.ticks = element_blank(),axis.text = element_blank(),# 繪圖面板panel.background = element_blank(),panel.border = element_rect(color = "black", linetype = "solid", size = 1, fill = NA))

推薦閱讀：Python + ElasticSearch：有了這個超級武器，你也可以報名參加詩詞大會了！| 博文精選 如何使用5個Python庫管理大數(shù)據(jù)？ 金三銀四網(wǎng)絡面經(jīng)之 DNS 詳解！ 愿得一心人：硅谷億萬富豪們的婚姻怎樣？有人白首相守七十年 一文告訴你，如何使用Python構建一個“谷歌搜索”系統(tǒng) | 內(nèi)附代碼 2019年度區(qū)塊鏈安全復盤總結 真香，朕在看了！猛戳“閱讀原文”，填寫中國遠程辦公-調(diào)查問卷

總結

以上是生活随笔為你收集整理的利用丁香园数据生成疫情分布地图（R语言）| 博文精选的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯，歡迎將生活随笔推薦給好友。