hunzikp / velox

https://hunzikp.github.io/velox/
119 stars 23 forks source link

extract: return raw "inclusion" values, not necesarily aggregating them with fun #8

Closed MatthieuStigler closed 6 years ago

MatthieuStigler commented 7 years ago

It would be great if the extract function could return also the raw values indicating on which cells a given polygon falls. This would emulate the functioning of raster::extract(), with fun=NULL.

This does not seem too difficult right now:

In pseudo-code (working only for nbands=1)

Before the p loop:

if(is.null(fun)) out <- vector(mode = "list", length=length(sp))
IDs <- sapply(sp@polygons, function(x) x@ID)

Within the p loop, line 87:

valmat <- hitmat[,3:ncol(hitmat),drop=FALSE]
if (!is.null(fun) & nrow(valmat) > 0) {
  for (k in 1:ncol(valmat)) {
    out[p,k] <- fun(valmat[,k])
  }
} else {
  out[[p]] <- data.frame(ID=IDs[p], valmat)
}

After the p loop:

if(is.null(fun)) out <- do.call( "rbind", out)
markwestcott34 commented 6 years ago

Yes, this would be really useful for me too; I want to run table() on the extracted points to get a summary of the values within given polygons.

In the meantime, I've added a extract_to_list() function which allows me to pass fun=table as an argument and get the tables back as elements of a list. (currently extract needs fun to return a single numeric, I think, as the results are stored in a matrix). However, this is probably not a proper solution as it only works for nbands = 1.

MatthieuStigler commented 6 years ago

Yes, I think indeed for multi-outputs functions, the best is to have extract() return a data-frame, and let user do the computation. Otherwise, would be quite complicated structure with the multi output/multi-band.

AronBoettcher commented 6 years ago

Upbump for visibility. Came here and posted the identical issue. This was the only reason I wanted to use this function. Either update it for multiple outputs, or allow for the raw values to be extracted. As it stands it is not particularly useful.

hunzikp commented 6 years ago

Hi everyone. Thanks for your comments.

The extract function now returns 'raw' raster values if argument fun is set to NULL. The return value is a list with one list element per polygon, and each list element consisting of a matrix with as many columns as there are raster bands.

## Make VeloxRaster with two bands
set.seed(0)
mat1 <- matrix(rnorm(100), 10, 10)
mat2 <- matrix(rnorm(100), 10, 10)
brk <- brick(raster(mat1), raster(mat2))
vx <- velox(list(mat1, mat2), extent=c(0,1,0,1), res=c(0.1,0.1),
            crs="+proj=longlat +datum=WGS84 +no_defs")

## Make SpatialPolygons
coord <- matrix(c(0,0,1,1), 2, 2, byrow = TRUE)
spoint <- SpatialPoints(coords=coord)
spols <- gBuffer(spgeom=spoint, width=0.2, byid = TRUE)

## Extract raw values as list of matrices
vx.elist <- vx$extract(sp=spols, fun = NULL)

## Print
print(vx.elist)
[[1]]
[,1]       [,2]
[1,] -0.005767173 -0.5006966
[2,]  2.404653389  1.6782972
[3,] -1.237538422  0.1380527

[[2]]
[,1]        [,2]
[1,] -0.7970895 -0.07041738
[2,]  1.2993123 -0.01674826
[3,] -0.8732621  0.16178863
MatthieuStigler commented 6 years ago

nice, thanks a lot for the great work!

AronBoettcher commented 6 years ago

Great! Can't wait to put this into action. Thanks!