[R-sig-Geo] Joining polygon and raster attributes

Empty Empty phytophthorasb at yahoo.com
Mon Aug 5 22:38:38 CEST 2013


[I apologize for reposting this – it showed up in the archive under the wrong subject, which would make it less useful as an archive.]


Hi. 
I’m trying to develop a repeatable process in R that I can use with a variety of data. Basically I'm trying to join/merge/combine/intersect (I’m not sure technically which word is right) a polygon and a raster so that I can do summaries on attributes from both. I have A) a polygon shape file of terrestrial ecoregions (contains a number of attributes that I am interested in) and B) a raster that is categories (0-2), but later I might have a continuous raster. Specifically, I want to ask things like what is the total area or proportion of area of biomes (from polygon A) that are >0 (from raster B). [Note that zonal statistics does Not seem to do what I want]. 

Below is some of the R code I have  tried and below that is the really cumbersome way I almost got it to work in  ArcGIS in case that clarifies what I’m trying to do – but I’m keen to do it in a scripted way since I will be changing the data later. 

I’m grateful for any suggestions. 
Thanks! 
  Juliann 


In R the closest I got is this, but I'm not sure I'm even on the right track. Perhaps I should be going about it in a completely different way. 

#(I probably don't need all these libraries. I was trying a lot of different things.) 

library(sp)            
library(raster)         
library(rgdal) 
library(maptools) 
library(rgeos) 

dryras <- raster("~overlay1.tif") #cells have values 0-2. 
teow <-readOGR(dsn=".", layer="teow_africa") 
#Terrestrial Ecoregions of the World 
pdf("quickplots.pdf") 
plot(dryras) 
plot(teow) 
dev.off()          
#Yup looks about the same as in arcgis 

# overlay points on raster and create new shapefile with extracted values in attribute table 

samp = extract(dryras, teow) 
#sampsm<-extract(dryras, teow, small=TRUE)  #maybe this would be better, but let’s see if it works first. 
#So far so good, it looks as though each polygon has different numbers of 0,1,2 associated with it (not sure how I'm going to get area from this, which is what I really want), but then it falls apart. I don't have column 
names on samp, for one thing (>names(samp) yields NULL) 

# combine with polygon data 
#One tangential problem is that I'm trying to use merge from the base library and raster also has merge, so I first removed raster before this step. I'm sure there's a simple and obvious solution!# 
test1 <- cbind(as.data.frame(teow), 
samp) 
test2<-merge(teow, samp) 
test3<-merge(teow, samp, by.x=[,1], by.y="OBJECTID", all.x=TRUE, 
all.y=TRUE, sort=FALSE) 


for test1& test2 the error is "Error in data.frame(c(0, 0, 0, 0, 0, 0), c(0, 0, 0, 0, 0, 0, 0, 0, 0,  :   arguments imply differing number of rows: . . ." 
for test 3 the error is "Error: unexpected '[' in "test2<-merge(teow, samp, by.x=["" 

Other things I've played with (over, gIntersection, rasterize) 

What I did in ArcGIS and JMP 
I came up with a really cumbersome and ugly brute force method (that I know for certain does not come out quite right): In ARCGIS, I converted the raster to a polygon, then used a union to combine this with the ecoregions polygon so that there were new polygons for each polygon-rastervalue combination that have all the attributes together. Then I added an area field to this new polygon table. This process didn't work 
very well as it left a few hundred weird (v. small) polygons that don't make sense and are missing attributes. Also, I don't think I used the right area calculator. Nevertheless, I took the table into JMP (just for the fun of 
looking at data) and was able to very easily do the sorts of summaries I was interested in. But I’d really like to script this and then I can do the summaries in R if I get the first part working.                   



More information about the R-sig-Geo mailing list