使用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") 

我想繪製相同的多邊形隨各圖中（見polygone形狀附圖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（情節$數據錯誤，setdiff（條件，名稱（DF）），降= 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

嗯，我終於能夠解決我的問題！非常感謝Winston Chang & Felipe Carrillo在ggplot2郵件列表上。

這是在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) 

希望這可以幫助別人！