title: extract value from raster stack from spatialpolygondataframe
body: <p>I have a raster stack with 27 rasters in it.  I have 27 corresponding polygons in a spatial polygon data frame.  I want to take polygon[i] overlay it on raster[i], extract and sum the values from raster [i], get a count of the number of cells within the polygon[i] and then divide the sum value by the # of cells.  In other words, the raster is a utilization distribution or a kernel density of use. I want to know much use is occurring in the area of the polygon where it is overlapping the raster.  I want to divide by the number of cells in the polygon to take into account the size of the polygon.  </p>
  
  <p>I have  a script that was given to me that does this, only it was written with the intention of extracting data from 1 raster only by any number of spatial polygons in the data frame.  It works, its ugly, and I now would like to convert it to something more stream line.  I only wish I had someone around me who could help because this might take a while?</p>
  
  <p>This is code Ive been given and my summary of what I think is going on:</p>
  
  <pre><code>msum99Kern07 = SpatialPolygonDataFrame (many polygons)
  KERNWolfPIX07m = Raster (this is a single raster, I have 27 rasters I put into a stack
  </code></pre>
  
  <p>)</p>
  
  <pre><code>#Extracting value from raster to many polygons 
  sRISK_Moose07m&lt;- extract(KERNWolfPIX07m, msum99Kern07,df=FALSE,method='bilinear')
  
  #Calculate THE SUM FOR EACH polygon#
  sRISK_Moose07m&lt;-unlist(lapply(sRISK_Moose07m, function(x) if (!is.null(x)) sum(x, na.rm=TRUE) else NA ))
  sRISK_Moose07m&lt;-as.data.frame(sRISK_Moose07m)
  
  #Im not sure why these next commands are needed Im only guessing
  #data.frame(levels) as there are many polygons creating a dataframe to put the info into
  ID_SUM_07&lt;-as.data.frame(levels(as.factor(msum07locs$ID2)))
  #ADD ID TO THE risk data frame 
  sRISK_Moose07m$ID&lt;-ID_SUM_07[,1] 
  
  #NUMBER OF CELLS WITHIN POLYGON EXTRACT CELLS/ POLYGON
  NB_SUM2007m&lt;-cellFromPolygon(KERNWolfPIX07m, msum99Kern07)
  NB_SUM07m&lt;-unlist(lapply(NB_SUM2007m, function(x) if (!is.null(x)) length(x) else NA ))
  
  #####CONVERT TO DATA FRAME
  NB_SUM07m&lt;-as.data.frame(NB_SUM07m)
  
  ###ADD THE NB OF CELLS TO THE RISK_SUM FILE###
  sRISK_Moose07m$NB_CELLS&lt;-NB_SUM07m[,1]
  
  ###DIVIDING VALUE by NB CELLS##
  sRISK_Moose07m$DIVID&lt;-sRISK_Moose07m$sRISK_Moose07m/sRISK_Moose07m$NB_CELLS 
  </code></pre>
  
  <p>Now, I have my spatial polygon data frame with 27 polygons and my raster stack with 27 rasters.  I want to select the raster[i] and polygon[i] and extract, sum, and calculate the kernel density of the overlapping area.  One side thing to keep in mind, I may get an error because it is possible that the polygon and raster do not overlap...I don't know how to check for this in R at all.  </p>
  
  <p>My script I have started:</p>
  
  <pre><code>moose99kern = spatial polygon data frame 27 moose
  Rastwtrial = stack of 27 rasters having the same unique name as the ID in moose99kern
  
  mkernID=unique(moose99kern$id)
  
  for (i in length(mkernID)){
             r = Rastwtrial[Rastwtrial[[i]]== mkernID[i]] #pick frm Rasterstack the raster that has the same name
              mp = moose99kern[moose99kern$id == mkernID[i]] #pick from spatialpolygondataframe the polygon that has the same name 
  
              RISK_MooseTrial&lt;- extract(r, mp, df=T, method'bilinear')
              risksum = (RISK_MooseTrial, function(x) if (!is.null(x)) sum(x, na.rm=TRUE) else NA )#sum all the values that were extracted from the raster
  </code></pre>
  
  <p>My script doesn't even start to work because I don't know how to index a raster stack.  But even still, going through 1 raster/1polygon at a time, Im not sure what to do next in the code.  If this is too much for StackOverflow I apologize.  Im just seriously stuck and have no where to turn.<br>
  Here is test data with 2 individuals for polygons</p>
  
  <pre><code> dput(mtestpoly) 
      new("SpatialPolygonsDataFrame"
          , data = structure(list(id = structure(1:2, .Label = c("F01001_1", "F07002_1"
      ), class = "factor"), area = c(1259.93082578125, 966.364499511719
      )), .Names = c("id", "area"), row.names = c("F01001_1", "F07002_1"
      ), class = "data.frame")
          , polygons = list(&lt;S4 object of class structure("Polygons", package = "sp")&gt;, 
          &lt;S4 object of class structure("Polygons", package = "sp")&gt;)
          , plotOrder = 1:2
          , bbox = structure(c(6619693.77161797, 1480549.31292137, 6625570.48348294, 
      1485861.5586371), .Dim = c(2L, 2L), .Dimnames = list(c("x", "y"
      ), c("min", "max")))
          , proj4string = new("CRS"
          , projargs = NA_character_
  </code></pre>
  
  <p>dput(Rastwtest)</p>
  
  <pre><code>new("RasterStack"
      , filename = ""
      , layers = list(&lt;S4 object of class structure("RasterLayer", package = "raster")&gt;, 
      &lt;S4 object of class structure("RasterLayer", package = "raster")&gt;)
      , title = character(0)
      , extent = new("Extent"
      , xmin = 1452505.6959799
      , xmax = 1515444.7110552
      , ymin = 6575235.1959799
      , ymax = 6646756.8040201
  )
      , rotated = FALSE
      , rotation = new(".Rotation"
      , geotrans = numeric(0)
      , transfun = function () 
  NULL
  )
      , ncols = 176L
      , nrows = 200L
      , crs = new("CRS"
      , projargs = NA_character_
  )
      , z = list()
      , layernames = "Do not use the layernames slot (it is obsolete and will be removed)\nUse function 'names'"
  )
  </code></pre>