Manuel--
The best way to proceed depends on the scale of your raster pixels versus
the detail of the polygons, and how 'precise' you need the estimate of
percent forest to be.

In my use case, I have small enough pixels that many are 100% or 0%, and I
mostly need to get the many edge pixels to perhaps within +/-5% of their
true cover.  My somewhat simpler R-only approach:

If your shapefile has polygons of various cover types, make a forest only
SpatialPolygonsDataFrame.
Next, create your desired raster, paying attention to the half-cell offset
(pixel centers v. pixel outer edges for the entire raster).  Then use the
getCover parameter in the raster::rasterize() function, which will compute
how much of each pixel in the raster is covered by the polygons in forest
(approximated by creating a 10x10 grid of subcells in each pixel and
testing whether the center of each subcell in in the polygons).

    forest <- nlcd[nlcd$Cover=='FOREST',]  # make SPDF of only forest
polygons
    bareRaster <- raster(...)  # create desired raster, paying attention to
half-cell offsets
    forestRaster <-
rasterize(forest,bareRaster,getCover=TRUE,filename='forestraster.grd')

See ?rasterize for details.
If you need a rasterBrick with layers for each level1 cover type, you can
roll the spdf subsetting and raster(...getCover) into a function, and then
use lappy or something from plyr.

In the alternative of many polygons smaller than each pixel (e.g., you are
computing global forest cover by 1deg cells), you may need to define a
raster with much smaller pixels (e.g., 100x100 of these smaller pixels for
each 1deg cell final pixel), use rasterize with getCover=FALSE to test the
center of each of the 10000 pixels in each large cell, then aggregate them.
 I haven't needed to do this.

For other folks who might find this thread while searching about forests,
note that most land cover data comes from RS imagery and starts out raster,
so there is some error introduced by converting to polygons and then more
error converting back to raster.  For the US, NLCD is available for several
dates at 30m pixels; NALC (North America Land Cover) is available for
MX/US/CA at 250m (and other) resolution.  I believe that there are global
(free) products with percent forest cover available, but those links are on
my office computer, and I'm locked out of work.

I hope that this helps...
Tom


On Tue, Oct 15, 2013 at 10:33 PM, Michael Treglia <mtreglia@gmail.com>wrote:

> Hi Manuel,
>
> My response may or may not be that helpful... I recently had to do
> something similar, for similar purposes, so I'll give you my general
> workflow. I used a few different software packages for this, thought it
> might be completely do-able in R (haven't tried all of this in R). I used
> SAGA GIS (freeware) a lot, and some Manifold GIS (paid software). But, I
> think this is probably all doable in freeware. I would definitely be
> curious to hear other ideas though... Steps 3, 3a, and 3b - Alternative are
> most relevant to your direct question. I apologize that I don't stick to R
> here, but figured any help might be useful.
>
> Getting the data from vector formats as polygons to rasters in the way that
> you describe can be difficult, and in my experience, needs intermediate
> processing. Again, I'll definitely be curious to hear of alternative
> strategies that anybody has.
>
> 1) I set up a vector file with grid cells, with the same grain size and
> extent as the desired raster. I actually started with a Raster, and
> converted it to polygons, using SAGA GIS. The attribute table from the
> polygons was useful in some additional data prep I needed to do, so I kind
> of wanted this vector GIS format anyway.
> 2) To transfer Raster data to these polygons, given that Raster data were
> frequently at different grain sizes than I needed, I used a function in
> SAGA GIS to Calculate Grid Statistics for Polygons. I calculated whatever
> features I needed (e.g., median Elevation...)  This automatically populated
> a new attribute (column)  in the attribute table for the Polygons. I tried
> this in R, using the extract() function, but my datasets were too large it
> seemed (>40,000 polygons).
>
> 3) For vector layers of environmental data, I clipped them to my study
> area, and if necessary I made dummy-layers. In other words, if there was a
> column of Land Cover, I would create separate vector layers for each land
> cover type that I was interested in. So, you would have a layer of Forest,
> where only polygons of Forest are included.
> 3a) I split the vector (environmental data) layers with the boundaries of
> the polygons. Then, did a spatial overlay to get the area of each
> environmental layer within each polygon as a new attribute. Then,
> calculated the percent of each polygon occupied by the given layer (e.g.,
> forest). I did this in Manifold GIS, but it might be do-able in R, or QGIS,
> or something else free.
> 3b - Alternative) Convert the Forest vector dataset to a fine grain Raster.
> Then either use a function in extract() (for R), or Grid Statistics for
> Polygons in SAGA, to get the count of pixels for Forest within each
> polygon, as a new attribute in the data table for the polygons, and convert
> that to the % of the polygon. *Naturally, you will lose information here,
> but the finer grain you make the Forest raster, the better off you'll be.
>
> 4) Either do the analyses using the tabular data from the Polygon dataset,
> or convert all of the layers to Rasters, and do analyses on those.
>
>
>
> To convert the grid-cell polygons to rasters (one raster per attribute), I
> used the rasterize() function, kind of like below. I'm sure the code could
> be cleaned up... this got me by and it works, but somebody else might be
> able to offer some improvements.
> library(rgdal)
> library(raster)
> #Import Polygon dataset
> polygons <- readOGR(... ) #enter appropriate info here
>
> #Import a raster that you want to use as a template... or create one from
> scratch
> raster.template<-raster("[file name here]")
>
> #Specify attributes that you want to convert to new raster layers - each
> polygon has a value for every attribute
> fields<-c("elevation", "forest cover")
>
> #make the raster template into a brick
> x<-brick(raster.template)
>
> #Create raster layer for each specified field
> for (i in fields){
>  x[[i]]<-rasterize(polygons, raster.template, field=i)
>  }
>
> #Set all of these up as a rasterStack... you can then plot these to make
> sure all layers look reasonable, or write them to a new file.
> y<-stack(x[[2]],x[[3]],x[[4]],x[[5]],x[[6]],x[[7]],x[[8]], x[[9]], x[[10]],
> x[[11]])
>
>
>
> Hope this helps at least a bit.. If nothing else, maybe it can help you
> think up alternative or simpler workflows...
> Mike
>
>         [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

	[[alternative HTML version deleted]]

