r - 使用ggplot2绘制多边形shapefile和geom_points

Question

我一直在努力解决这个问题，希望能对您有所帮助。
我试图在我的geom_points上绘制一个多边形。到目前为止，这是我所做的:

> names(OT1)# my dataset
[1] "EID"       "latitude"  "longitude" "month"     "year"      "CPUE"      "TSUM" 
> dim(OT1)
[1] 2707    7
> head(OT1)
                EID latitude longitude month year CPUE TSUM
1   167-1-1996-1135 67.70000 -61.81667     9 1996    0    0
2  167-10-1996-1135 67.71667 -59.18333     9 1996    0    0
3 167-100-1996-1135 67.86667 -59.43333    10 1996    0    0
4 167-101-1996-1135 67.95000 -59.58333    10 1996    0    0
5 167-102-1996-1135 68.10000 -59.76667    10 1996    0    0
6 167-103-1996-1135 67.81667 -59.38333    10 1996    0    0

OTz<-OT1[with(OT1,OT1$TSUM=="0"),]#selecting only my zeros
OTc<-OT1[!with(OT1,OT1$TSUM=="0"),]

#plotting data with ggplot2 (see attached figure)
v<-ggplot() + geom_point(aes(longitude, latitude, size=TSUM),data= OTc, colour=alpha("red",0.2)) +  facet_wrap(~month, ncol=2)
v + geom_point(aes(longitude, latitude),size = 1,colour = alpha("black", 0.2), data = OTz) + opts(title="Otter trawl 1996-2011")

我想在每个图上绘制相同的多边形形状(请参见图2所示的多边形形状)。我遵循 R-help Re: another question on shapefiles and geom_point in ggplot2 和 Plotting polygon shapefiles上的说明。我可以绘制多边形，但是很难覆盖我的geom_points。

library(rgdal)
library(ggplot2)
library(sp)
library(maptools)
gpclibPermit()
div0A <- readOGR(dsn=".", layer="Projections")
> div0A <- readOGR(dsn=".", layer="Projections")
OGR data source with driver: ESRI Shapefile 
Source: ".", layer: "Projections"
with 1 features and 5 fields
Feature type: wkbPolygon with 2 dimensions
> names(div0A);dim(div0A)
[1] "Id"         "NPAzimutha" "UTM20"      "UTM19"      "AlberEA"   
[1] 1 5
> slotNames(div0A) # l
[1] "data"        "polygons"    "plotOrder"   "bbox"        "proj4string"

# add the 'id' variable to the shapefile and use it to merge both files
div0A@data$id = rownames(div0A@data)
div0A.df <- as.data.frame(div0A)# convert shapefile to dataframe
div0A_fort <- fortify(div0A,region="id")# fortify to plot with ggplot2 
head(div0A_fort)
> head(div0A_fort)
       long      lat order  hole piece group id
1 -73.50000 78.16666     1 FALSE     1   0.1  0
2 -73.50000 75.24043     2 FALSE     1   0.1  0
3 -73.38552 75.04169     3 FALSE     1   0.1  0
4 -72.95306 74.78239     4 FALSE     1   0.1  0
5 -70.11000 74.10167     5 FALSE     1   0.1  0
6 -68.62608 73.72649     6 FALSE     1   0.1  0
# Merge shapefile and the as.dataframe shapefile
div0A_merged <- join(div0A_fort,div0A.df, by ="id")
head(div0A_merged)
> head(div0A_merged)
       long      lat order  hole piece group id Id NPAzimutha    UTM20  UTM19  AlberEA
1 -73.50000 78.16666     1 FALSE     1   0.1  0  0     348877 349232.4 349162 348656.4
2 -73.50000 75.24043     2 FALSE     1   0.1  0  0     348877 349232.4 349162 348656.4
3 -73.38552 75.04169     3 FALSE     1   0.1  0  0     348877 349232.4 349162 348656.4
4 -72.95306 74.78239     4 FALSE     1   0.1  0  0     348877 349232.4 349162 348656.4
5 -70.11000 74.10167     5 FALSE     1   0.1  0  0     348877 349232.4 349162 348656.4
6 -68.62608 73.72649     6 FALSE     1   0.1  0  0     348877 349232.4 349162 348656.4
# Plot the shapefile
ggplot(div0A_merged, aes(long,lat,group=group)) +
  geom_polygon(aes(data=div0A_merged)) +
  geom_path(color="white") + theme_bw()

当我尝试使用以下代码作为测试时，我收到一条错误消息:“ [.data.frame中出现错误(plot $ data，setdiff(cond，names(df))，drop = FALSE):
未定义的列已选择” ...

p<-ggplot(div0A_merged, aes(long,lat,group=group)) +
  geom_polygon(aes(data=div0A_merged)) +
  geom_path(color="white") + theme_bw()

p + geom_point(aes(longitude, latitude, size=TSUM),data= OTc, colour=alpha("red",0.2)) +  facet_wrap(~month, ncol=2)

谢谢你!

Answer 1

好吧，我终于能够弄清楚我的问题!非常感谢ggplot2邮件列表上的Winston Chang和Felipe Carrillo。

这是在ggplot2版本0.8.9上执行此操作的一种方法。

library(ggplot2)

OT1 <- read.csv('OT1.csv')

OTz<-OT1[OT1$TSUM==0,]#selecting only my zeros
OTc<-OT1[OT1$TSUM!=0,]

# plotting data with ggplot2
library(scales)
v <- ggplot(OTc, aes(longitude, latitude, size=TSUM)) +
  geom_point(colour="red", alpha=0.1) + facet_wrap(~month, ncol=2)
v + geom_point(data = OTz,size = 1,colour = "black", alpha=0.2) +
  opts(title="Otter trawl 1996-2011")


library(rgdal)
library(sp)
library(maptools)
gpclibPermit()

div0A <- readOGR(dsn=".", layer="Projections")

names(div0A)
dim(div0A)

library(gpclib)
# add the 'id' variable to the shapefile and use it to merge both files
div0A@data$id = rownames(div0A@data)
div0A.df <- as.data.frame(div0A)# convert shapefile to dataframe
div0A_fort <- fortify(div0A, region="id")# fortify to plot with ggplot2 
head(div0A_fort)

# Merge shapefile and the as.dataframe shapefile
library(plyr)
div0A_merged <- join(div0A_fort, div0A.df, by="id")
head(div0A_merged)

# Get all the months used in OTc
monthdf <- data.frame(month = unique(OTc$month))

# Merge with div0A_merged 
# (replicate each row in div0A_merged for each month)
div0A_merged_month <- merge(div0A_merged, monthdf)

# Graph with the shapefile
ggplot(div0A_merged_month, aes(long, lat, group=group)) +
  geom_polygon() +
  geom_path(color="white") +
  geom_point(data=OTc, aes(x=longitude, y=latitude, size=TSUM),
             colour="red", alpha=0.2, inherit.aes=FALSE) +
  theme_bw() +
  facet_wrap(~ month, ncol=2)

希望这对其他人有帮助!