[R-sig-Geo] Intersect polygon with grid

[Robert J. Hijmans]

Johannes,

Here is another approximation, based on the concept that Paul
suggested (first disaggregating, then aggregating the results):

library(raster)

# get some data
data(meuse.riv)
pol = SpatialPolygons(list(Polygons(list(Polygon(meuse.riv)),"meuse")))
data(meuse.grid)
coordinates(meuse.grid) = c("x", "y")
gridded(meuse.grid) = TRUE

# create a RasterLayer (in your case from file)
r <- raster(meuse.grid)

rr <- disaggregate(raster(r), 10)
pr <- polygonsToRaster(pol, rr, -1, progress='text')
ra <- aggregate(pr, 10, sum, progress='text')

plot(ra)
plot(pol, add=TRUE)

These should work on a raster of any size (and ample patience).

Robert



On Tue, Dec 22, 2009 at 4:02 AM, Paul Hiemstra <p.hiemstra at geo.uu.nl> wrote:
> Hi,
>
> A simple approach would be the following:
>
> 1. Create a point dataset that comprises of the center points of the grid
> 2. Overlay the pointset with the polygon set using the overlay() command
> 3. Convert the pointset into a grid using the gridded() command
>
> A problem could be that if your polygons are small in regard to the
> gridsize, only the center point of the gridcell determines the value of the
> gridcell. It does not look at how much area from several polygons is in a
> gridcell, and creating a weigted average by area. You could remedy this by
> making the pointset much denser, doing step 1-3 and then aggregating the
> dense grid into a coarser grid, calculating for example the mean. This can
> be done using the pointsToRaster() function from the raster package.
>
> cheers,
> Paul
>
> Johannes Signer wrote:
>>
>> Dear List,
>>
>> I am trying to do the following:
>> I have a shapefile, that I want to overlay with a raster and then assign
>> each raster cell a value that represents the percentage covered by the
>> shape
>> file. I have found a way to do that with small data sets, but I am
>> struggling with large datasets (1mio+ cells).
>>
>> My approach so far is:
>>
>> 1.) Create rectangular polygons (covering the whole area), that will be
>> used
>> for the grid afterwards.
>> 2.) Intersect the each polygon with the shapefile and calculate the area.
>> 3.) Create a grid from the rectangular polygons.
>>
>> Any suggestions are very welcome!
>>
>> Thanks in advance
>>
>> Johannes
>>
>>        [[alternative HTML version deleted]]
>>
>> _______________________________________________
>> R-sig-Geo mailing list
>> R-sig-Geo at stat.math.ethz.ch
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>
>
>
> --
> Drs. Paul Hiemstra
> Department of Physical Geography
> Faculty of Geosciences
> University of Utrecht
> Heidelberglaan 2
> P.O. Box 80.115
> 3508 TC Utrecht
> Phone:  +3130 274 3113 Mon-Tue
> Phone:  +3130 253 5773 Wed-Fri
> http://intamap.geo.uu.nl/~paul
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at stat.math.ethz.ch
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>