[R-sig-Geo] problem with extract() raster package using weight=TRUE

[Ariel Ortiz-Bobea]

That's much clever than my longer workaround using rasterize and zonal.

By the way, Umberto, my workaround could help you get more precision on the
extract weights. Hope this help others looking to do the same. See code
below:

---------------

# Create raster
   library(raster)
   r <- raster(ncol=36, nrow=18)
   r[] <- 1:ncell(r)
   cds1 <- rbind(c(-180,-20), c(-160,5), c(-60, 0), c(-160,-60),
c(-180,-20))
   cds2 <- rbind(c(80,0), c(100,60), c(120,0), c(120,-55), c(80,0))
   polys <- SpatialPolygons(list(Polygons(list(Polygon(cds1)), 1),
         
# visualize
   plot(r)
   plot(polys, add=T)
     
# Extract the cellnumbers and weights using extract()
   v <- extract(r, polys, weights=TRUE, cellnumbers=TRUE, small=TRUE)

# Extract the cellnumbers and weights using rasterize() and zonal()
   r.mod <- r
   values(r.mod) <- 1:ncell(r)
   fac <- sqrt(1000) # this is where you set the resolution factor
   r.mod <- disaggregate(r.mod, fact=c(fac, fac)) # increases raster
resolution by a factor of "fac"
   r.polys <- rasterize(x=polys, y=r.mod) # rasterize county polygons

   # Compute weights
   w <- lapply(unique(r.polys), function(x) {
      m <- r.polys
      m[getValues(m)!=x]  <- NA
      m[getValues(m)==x]  <- 1
      out <- zonal(x=m, z=r.mod, fun="sum")
      colnames(out) <- c("cell", "weight")
      out[,"weight"] <- out[,"weight"]/(fac^2)
      out[out[,"weight"]!=0,]
     })
   names(w) <- unique(r.polys)

# See weigths
 head(v[[1]])
 head(w[[1]]) 



-----
Ariel Ortiz-Bobea
Fellow at Resources for the Future
--
View this message in context: http://r-sig-geo.2731867.n2.nabble.com/problem-with-extract-raster-package-using-weight-TRUE-tp7584665p7584693.html
Sent from the R-sig-geo mailing list archive at Nabble.com.