Plotting shapefile using ggplot2

数据准备

演示的数据放置在data文件夹下。

library(rgdal)
library(dplyr)
options(stringsAsFactors = F)
mpa<-read.csv('data/mpa.csv')
admin.poly.dsn<-'data/省级行政区.shp'# windows下中文路径存在编码问题
admin.poly.layername<-ogrListLayers(admin.poly.dsn)[1]
admin.poly<-readOGR(admin.poly.dsn,layer = admin.poly.layername,verbose = F)

行政区边界数据的空间多边形如下图所示(点击多边形有提示)，

library(leaflet)
popup <- paste0("<strong>Name: </strong>", 
                        admin.poly$NAME)
leaflet()%>%addPolygons(data=admin.poly,popup=popup)%>%
    addTiles(urlTemplate = 'http://t0.tianditu.cn/DataServer?T=vec_w&X={x}&Y={y}&L={z}')%>%
    addMiniMap()

---

行政区边界的属性表如下所示，

DT::datatable(admin.poly@data,rownames = F, class = 'cell-border compact stripe')

---

专题数据如下所示，

DT::datatable(mpa,rownames = F, class = 'cell-border compact stripe')

---

地图可视化

基本出图

下面按pm2.5和rr列数值作图，

library(ggplot2)
p <- ggplot()+
geom_point(data=mpa,mapping=aes(x=long,y=lat,colour=rr,size=PM25))+theme_bw()+
    scale_colour_gradient(low = 'green',high = 'red')
p

南海诸岛

ggplot2使用多边形数据时，需要先清洗多边形数据的结构，提取分区所需的干净信息。

cn.df <- fortify(admin.poly)

## Regions defined for each Polygons

head(cn.df)

##       long      lat order  hole piece id group
## 1 112.0392 3.834796     1 FALSE     1  0   0.1
## 2 112.0327 3.833843     2 FALSE     1  0   0.1
## 3 112.0208 3.835330     3 FALSE     1  0   0.1
## 4 112.0104 3.837053     4 FALSE     1  0   0.1
## 5 112.0077 3.842897     5 FALSE     1  0   0.1
## 6 112.0067 3.850786     6 FALSE     1  0   0.1

为了不缺少南海诸岛，制图时只有把整个main.plot包含在内，使得地图比例显得不美观，地图排版显得不紧凑。

main.plot<-ggplot()+
    geom_path(data = cn.df, aes(long, lat, group = group),
              colour = "black")+
    geom_point(data=mpa,mapping=
                   aes(x=long,y=lat,colour=rr,size=PM25))+
    theme_bw()+ scale_colour_gradient(low = 'green',high = 'red')
main.plot

通常的制图方法是嵌入南海诸岛区域子地图，在R中可以按照下面的方法做。

sub.plot<-main.plot+
    theme_bw() + labs(x=NULL, y=NULL) +
    theme(axis.line=element_blank(),axis.text.x=element_blank(),
          axis.text.y=element_blank(),axis.ticks=element_blank(),
          axis.title.x=element_blank(),
          axis.title.y=element_blank(),
          legend.position = "none", 
          legend.text=element_text(size=12),
          legend.background =element_blank(),
          panel.background=element_blank(),
          panel.grid.major=element_blank(),panel.grid.minor=element_blank(),
          plot.background=element_blank())+
    coord_cartesian(xlim = c(106.8,122.0),ylim=c(3.307 ,25.25))
subplot<-ggplotGrob(sub.plot)
natioanl.plot<-main.plot+
    coord_cartesian(xlim = c(72.2,135.42),ylim=c(14.7 ,54.45))+
    annotation_custom(grob= subplot,xmin=123.5,xmax = Inf,ymin=-Inf,ymax=30)

全国尺度的最终效果

natioanl.plot

由于图例的占位影响了子区域地图的位置，南海诸岛部分没有严格对齐地图边缘。这个问题有待解决。

区域尺度制图

• 利用空间数据框的索引选择该城市的多边形作为底图

• 利用coord_cartesian(xlim = c(72.2,135.42),ylim=c(14.7 ,54.45))函数限定研究区空间范围

bj.idx<-!is.na(admin.poly$NAME)&admin.poly$NAME=='北京市'
bj.poly<-admin.poly[bj.idx,]
bj.df<-fortify(bj.poly)
library(raster)
bj.extent<-extent(bj.poly)
ggplot()+
    coord_cartesian(xlim = c(bj.extent@xmin,bj.extent@xmax),
                    ylim=c(bj.extent@ymin ,bj.extent@ymax))+
    geom_path(data = bj.df, aes(long, lat, group = group),
              colour = "black")+
    geom_point(data=mpa,mapping=
                   aes(x=long,y=lat,colour=rr,size=PM25))+
    theme_bw()+ scale_colour_gradient(low = 'green',high = 'red')