How to make a raster (from a shapefile) with cells with forest cover percent.

classic Classic list List threaded Threaded
5 messages Options
Reply | Threaded
Open this post in threaded view
|

How to make a raster (from a shapefile) with cells with forest cover percent.

Manuel Spínola
Dear list members,

I have a shapefile with land cover (including forest cover) and I want to
make a raster that give the percent of forest cover in each cell of
specific size.  Is that possible? I was exploring the function "focal" in
the pacakge raster but I do not fully understand how it works?

I did the raster from the shapefile already but each cell take an integer
number according to the land cover type (forest = 1, urban areas = 2, etc.)

The final objective is to map the ocurrence of a species given some
covariates, including forest cover, and I am preparing the environmental
data.

Best regards,

Manuel


--
*Manuel Spínola, Ph.D.*
Instituto Internacional en Conservación y Manejo de Vida Silvestre
Universidad Nacional
Apartado 1350-3000
Heredia
COSTA RICA
[hidden email]
[hidden email]
Teléfono: (506) 2277-3598
Fax: (506) 2237-7036
Personal website: Lobito de río <https://sites.google.com/site/lobitoderio/>
Institutional website: ICOMVIS <http://www.icomvis.una.ac.cr/>

        [[alternative HTML version deleted]]


_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

Re: How to make a raster (from a shapefile) with cells with forest cover percent.

Michael Treglia
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
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

Re: How to make a raster (from a shapefile) with cells with forest cover percent.

Tom Philippi
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 <[hidden email]>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
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

Re: How to make a raster (from a shapefile) with cells with forest cover percent.

Manuel Spínola
In reply to this post by Michael Treglia
Thank you Michael.

I tried Tom Philippi suggestion and it works.

Manuel


2013/10/15 Michael Treglia <[hidden email]>

> 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
>


--
*Manuel Spínola, Ph.D.*
Instituto Internacional en Conservación y Manejo de Vida Silvestre
Universidad Nacional
Apartado 1350-3000
Heredia
COSTA RICA
[hidden email]
[hidden email]
Teléfono: (506) 2277-3598
Fax: (506) 2237-7036
Personal website: Lobito de río <https://sites.google.com/site/lobitoderio/>
Institutional website: ICOMVIS <http://www.icomvis.una.ac.cr/>

        [[alternative HTML version deleted]]


_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

Re: How to make a raster (from a shapefile) with cells with forest cover percent.

Manuel Spínola
In reply to this post by Tom Philippi
Thank you very much Tom.

It works.

Manuel


2013/10/16 Tom Philippi <[hidden email]>

> 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 <[hidden email]>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
>> [hidden email]
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>
>
>

--
*Manuel Spínola, Ph.D.*
Instituto Internacional en Conservación y Manejo de Vida Silvestre
Universidad Nacional
Apartado 1350-3000
Heredia
COSTA RICA
[hidden email]
[hidden email]
Teléfono: (506) 2277-3598
Fax: (506) 2237-7036
Personal website: Lobito de río <https://sites.google.com/site/lobitoderio/>
Institutional website: ICOMVIS <http://www.icomvis.una.ac.cr/>

        [[alternative HTML version deleted]]


_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo