Skip to content

overlay polygon with raster

5 messages · Frede Aakmann Tøgersen, José Hidasi, Nahm Lee +1 more

Original
Simone Ruzza
# 
Dear R-sig-geo users,

apologies for the simplicity of my question, but I have been
struggling with this for a bit I have not managed to find a solution.
Basically, I would like to overlay a polygon onto raster grid and
extract some statistics for the polygon in question i.e. the area of
every raster cell covered by the polygon. Here is a figure that might
make it clearer
<https://dl.dropboxusercontent.com/u/33966347/example.JPG> . I have
been fiddling around with the extract function from the raster
package, which works to some extent, as I am able to retrieve
information on the categories of each of the raster cells covered by
the polygon. However, as I said, I am interested in the area covered
by the polygon within each raster cell. I know that perhaps I should
try and convert my polygon into another raster and then overlay the
two rasters, but this perhaps would mean that the polygon should be
converted into a raster at a high resolution? If this was case, it is
not computationally feasible for me as some polygons I have cover a
huge area and it would take forever to overlay all those rasters. Any
help, suggestions would be highly appreciated!

thanks in advance!

Simone
Frede Aakmann Tøgersen
# 
Hi Simone

If I understand you correctly, then have a look at the man page of raster(). Here is the relevant sections:

     ## S4 method for signature 'Raster,SpatialPolygons'
     extract(x, y, fun=NULL, na.rm=FALSE, weights=FALSE, cellnumbers=FALSE,
          small=FALSE, df=FALSE, layer, nl, factors=FALSE, sp=FALSE, ...)

Where
       x: Raster* object

       y: points represented by a two-column matrix or data.frame, or
          'SpatialPoints*'; 'SpatialPolygons*'; 'SpatialLines';
          'Extent'; or a numeric vector representing cell numbers

and

weights: logical. If 'TRUE', the function returns, for each polygon, a
          matrix with the cell values and the approximate fraction of
          each cell that is covered by the polygon(rounded to 1/100).
          The weights can be used for averaging; see examples. This
          option can be useful (but slow) if the polygons are small
          relative to the cells size of the Raster* object.

I just googled for ' R area of raster cells within polygon' and got this hit with a reproducible small example:

http://stackoverflow.com/questions/17766989/extract-data-from-raster-with-small-polygons-rounded-weights-too-small 

See that page for a solution to some rounding errors in the case of the polygon being very small compared to raster cell size.

Yours sincerely / Med venlig hilsen


Frede Aakmann T?gersen
Specialist, M.Sc., Ph.D.
Plant Performance & Modeling

Technology & Service Solutions
T +45 9730 5135
M +45 2547 6050
frtog at vestas.com
http://www.vestas.com

Company reg. name: Vestas Wind Systems A/S
This e-mail is subject to our e-mail disclaimer statement.
Please refer to www.vestas.com/legal/notice
If you have received this e-mail in error please contact the sender.
José Hidasi
# 
Hello Simone,

You can also intersect the raster with the polygon and then calculate the
area covered by each raster cell. If that's what you want, you might need
to:

1) use "rasterToPolygons()" function to transform the raster into a polygon.

2) apply an equal area projection (e.g. lambert azimuthal) to both polygons
using "spTransform()" (raster/polygon and the other polygon).

3) create a new attribute in your "raster polygon" to identify each cell
(i.e. give a number to each cell). After creating the attribute, you can
use something like:
rasterpolygon$new.attribute <- c(1:length(rasterpolygon$new.attribute))

4) intersect both polygons using "intersect()".

5) use "gArea()" function to calculate the area of each "intersected cell",
but remember to use the "byid=T" argument. It should look like "area.cells
<- gArea(IntersectedPolygon, byid=T)". It will return a vector with the
area sizes in the units of the projection you used in step 2.

6) finally, get the area for each cell while identifying its number before
the intersection. To do that, you can bring up together "area.cells" (from
step 5) and the attribute you created in step 3 (which was maintained after
the intersection). For example:
info <- cbind(intersected.raster$new.attribute, area.cells)


Just an idea.


All the best,
Jose


On Wed, Dec 3, 2014 at 10:15 PM, Simone Ruzza <simone.ruzza12 at gmail.com>
wrote:

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Nahm Lee
    

    

      
      #
    

  


  

      
Yes, there is a round issues when you use "extract or mask", http://gis.stackexchange.com/questions/112274/is-this-a-known-issue-in-gaps-between-masked-raster-and-spatial-polygon-in-r-wit.

I use "disaggregate" to improve(?) raster resolution.
This is an example

yourraster.dis <- disaggregate("yourraster" , fact = 2)

Nahm Lee
nlee at valleywater.org
Santa Clara Valley Water District

-----Original Message-----
From: R-sig-Geo [mailto:r-sig-geo-bounces at r-project.org] On Behalf Of Frede Aakmann T?gersen
Sent: Wednesday, December 03, 2014 10:37 PM
To: Simone Ruzza; r-sig-geo at r-project.org
Subject: Re: [R-sig-Geo] overlay polygon with raster

Hi Simone

If I understand you correctly, then have a look at the man page of raster(). Here is the relevant sections:

     ## S4 method for signature 'Raster,SpatialPolygons'
     extract(x, y, fun=NULL, na.rm=FALSE, weights=FALSE, cellnumbers=FALSE,
          small=FALSE, df=FALSE, layer, nl, factors=FALSE, sp=FALSE, ...)

Where
       x: Raster* object

       y: points represented by a two-column matrix or data.frame, or
          'SpatialPoints*'; 'SpatialPolygons*'; 'SpatialLines';
          'Extent'; or a numeric vector representing cell numbers

and

weights: logical. If 'TRUE', the function returns, for each polygon, a
          matrix with the cell values and the approximate fraction of
          each cell that is covered by the polygon(rounded to 1/100).
          The weights can be used for averaging; see examples. This
          option can be useful (but slow) if the polygons are small
          relative to the cells size of the Raster* object.

I just googled for ' R area of raster cells within polygon' and got this hit with a reproducible small example:

http://stackoverflow.com/questions/17766989/extract-data-from-raster-with-small-polygons-rounded-weights-too-small 

See that page for a solution to some rounding errors in the case of the polygon being very small compared to raster cell size.

Yours sincerely / Med venlig hilsen


Frede Aakmann T?gersen
Specialist, M.Sc., Ph.D.
Plant Performance & Modeling

Technology & Service Solutions
T +45 9730 5135
M +45 2547 6050
frtog at vestas.com
http://www.vestas.com

Company reg. name: Vestas Wind Systems A/S This e-mail is subject to our e-mail disclaimer statement.
Please refer to www.vestas.com/legal/notice If you have received this e-mail in error please contact the sender.

            

      
      
      
    

        

      
      
    

        
_______________________________________________
R-sig-Geo mailing list
R-sig-Geo at r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

  
  


  


      

    

    

      
      

        


  

        

      Simone Ruzza
    

    

      
      #
    

  


  

      
Thank you for your prompt and kind replies! I think that maybe what
Jose Hidasi was suggesting is what I want. I will give it a go! Just
to clarify, what I would like to extract is the only the area (or
proportion) of a particular cell covered by the polygon and the
corresponding raster category. I am not sure whether this makes it
clearer but I have drawn a rough sketch of the information I am
interested in, highlighted in different colours...
<https://dl.dropboxusercontent.com/u/33966347/query1.JPG>Ideally, I
would like to extract the area or the proportion of each the sections
in different colours and the corresponding category number for the
raster below it...

thanks in advance!

Simone

        
On Thu, Dec 4, 2014 at 9:32 PM, Nahm Lee <nlee at valleywater.org> wrote: