會用到的包

library(magrittr) #管道函數 library(stringr) #字符處理 library(ggplot2) #繪圖 library(rgdal) #讀取shp文件 library(plyr) #合并數據 library(legendMap) #加比例尺 library(ggspatial) #加指北針 library(tidyfst) # count_dt library(dplyr) # left_join library(patchwork) #拼圖 library(RColorBrewer) # 挑顏色的##我不記得哪些是GitHub的包了，反正這個是 devtools::install_github("3wen/legendMap")

四處找了很多文章，但別人做的跟自己的需求總是不太一樣，東拼西湊加上自己碼代碼搞出了這篇文章。

導入地圖數據

setwd("E:/數據分析") shp_city <- readOGR("C:/Users/地圖/quxian",encoding = "UTF-8") shp_city <- spTransform(shp_city, CRS("+proj=longlat +datum=WGS84")) dtframe <- fortify(shp_city) # 轉化為數據框 data1 <- shp_city@data #提取行政區劃信息 data2 <- data.frame(id=row.names(data1),data1) #添加序號 map_data <- join(dtframe,data2, type = "full") #合并地理信息數據與經緯度數據 map_data$sheng_code <- str_sub(map_data$ZONECODE,1,2) map_data <-subset(map_data,sheng_code==13) #save(map_data,dt_qc,file = "地圖及數據.Rdata") load("地圖及數據.Rdata") #? fix(annotation_north_arrow) #設置為1cm

嘮叨幾句

• 我習慣跑代碼之前先設置一下工作路徑，這樣讀數據存數據的時候就不用輸入一長串的路徑也不用替換里面的“/”。

• 這里把我導數據過程遇到的問題提一下，我默認會點進來看的各位都是第一次用ggplot2畫地圖的啦。我用到的是整個中國的地圖，如果只用到某市某省的話，可以在R里按條件subset提取，也可以在GIS里提取完直接用，都ok。

• 劃重點：導入shp格式的數據最重要的是shp文件，但是如果路徑文件夾下只有shp文件的話是讀不了的，其他的文件也要有，所以你的文件夾里面大概是下面這個樣子。我用的是readOGR讀shp文件，跟讀csv之類的文件不一樣，你可以：
  （1）直接指定文件里的shp文件，也可以
  （2）直接讀有這個shp文件的文件夾
  

• 其實數據讀對了后面就都簡單了，我導入數據的時候，第一個問題就是中文是亂碼，這里你可以從下面兩句挑一種試一試，如果還有問題，可以再試試別的encoding格式或者換個數據（不然還能怎么辦呢…

readOGR("C:/Users/地圖/quxian", encoding = "UTF-8") 或者 readOGR("C:/Users/地圖/quxian", use_iconv = T , encoding = "UTF-8")

• 還有個問題就是，你的shp文件里的經緯度可能特別大（如果正常可以忽略這條），有666666這么大的那種，這樣的話你如果不用投影坐標系倒是也能畫，但是你得用呀，用下面這句試試：

shp_city <- spTransform(shp_city, CRS("+proj=longlat +datum=WGS84"))

• 接下來按照導數據的代碼跑就ok了，為啥那么跑我不知道，大家都是那么跑的，整理完最后用到的會是map_data，因為懶得每次導數據都整理一通，所以干脆把它存成.Rdata，以后直接讀.Rdata。

準備畫圖的病例數據

names_coun <- map_data[!duplicated(map_data$CNTY_CODE),] # 去重 names_coun <- names_coun[,c(8,9)];names_coun$CNTY_CODE <- as.character(names_coun$CNTY_CODE) sd_perm <- table(dt_qc$code)%>%data.frame() names(sd_perm)[1:2] <- c("CNTY_CODE","頻數") sd_perm$CNTY_CODE <- as.character(sd_perm$CNTY_CODE) sd_perm <- left_join(sd_perm,names_coun,by = c('CNTY_CODE')) sd_map_data <- join(map_data,sd_perm, type="full")#將業務數據與地理信息數據再次合并 # 區縣坐標 midpos <- function(x) mean(range(x,na.rm=TRUE))#取形狀內的平均坐標 centres <- ddply(sd_map_data,.(CNTY_CODE),colwise(midpos,.(long,lat))) centres <- join(centres,sd_perm, by="CNTY_CODE")%>%na.omit() #僅標記有病例的區縣#?

解釋一下代碼

• 第一行：去掉map_data中 CNTY_CODE 的重復數據，這里我是想取地圖數據中的行政區劃編碼和對應區縣中文名稱
• 第二行：提取行政區劃編碼和對應區縣中文名稱這兩列
• 第三行：統計一下每個行政區劃都有幾個病例，轉為數據框
• 第四行：改下數據框列名，為了下一步匹配合并
• 第五行：確保兩個需要合并數據框的數據格式一致
• 第六行：以sd_perm為主（合并完的行數還是sd_perm的行數）通過"CNTY_CODE"匹配進行合并
• 第七行：合并得到最終繪圖數據
• 第八-十行：得到想標記區縣的經緯度，取的是每個區縣整體經緯度均值

繪圖

out <- ggplot(sd_map_data, aes(long,lat)) +geom_polygon(aes(group = group,fill=頻數),colour="grey40",size=0.4) +scale_fill_gradient(low="#F7FCFD",high="#41AE76",na.value = "white") +coord_map() + # 選擇投影坐標系ggtitle("發病分布（區縣）")+ #添加標題geom_text(aes(label=NAME),data=centres,size=3)+ #標記市annotation_north_arrow(location = "tr")+scale_bar(lon =120.2,lat = 34.4, #比例尺位置distance_lon =100,distance_lat =16, #比例尺長度;比例尺寬度distance_legend =30,dist_unit ="km",orientation=F)+ #文字與比例尺的距離;比例尺單位 theme_bw()+theme(panel.grid = element_blank(),legend.position = c(0.90,0.35))+ #有順序annotate("text", x = 115, y = 38.2,label = "b",size=5) ggsave("E:/區縣.png",out,dpi=300)

解釋一下代碼

scale_fill_gradient() 設置填充顏色，從低到高梯度變化，并設置了na值的顏色

coord_map() 墨卡托投影（等角正切圓柱投影）; 也可以選擇coord_map(“polyconic”)多圓錐投影

annotation_north_arrow() 在topright放入指北針，這里提一下，要在fix(annotation_north_arrow)中修改源代碼中的圖例大小，默認1.5cm，可根據自己喜好調整

scale_bar()設置比例尺，orientation=F為不添加這個函數里的指北針【太丑

放個半成品【因為保密問題把標記和標題都去掉了

先這樣啦，如果有需要可以再出個批量畫圖的文章。

總結

以上是生活随笔為你收集整理的应用ggplot2绘制病例分级统计地图的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。