[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