2012-02-27 194 views
6

我一直在为这个情节挣扎,并希望得到任何帮助。 我想绘制一个多边形在我的geom_points上。这是我迄今所做的:使用ggplot2绘制多边形shapefile和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") 

enter image description here

我想绘制相同的多边形随各图中(见polygone形状附图2)。我遵循R-help Re: another question on shapefiles and geom_point in ggplot2Plotting 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 
[email protected]$id = rownames([email protected]) 
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() 

enter image description here

当我尝试这样做下面的代码作为测试,我得到一个错误信息:“在[.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) 

谢谢!

回答

12

嗯,我终于能够解决我的问题!非常感谢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 
[email protected]$id = rownames([email protected]) 
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) 

希望这可以帮助别人!

+0

Hi GodinA,你能不能分享你的数据,以便我们更容易理解并复制图表。谢谢。 – 2013-01-17 00:24:26

+0

伟大的帮助和伟大的例子!感谢分享。 – Chris 2013-04-30 11:34:12