Create raster with grid cells values based on area covered by polygons

Question

I have a polygon data with area covered by the forests (data is here - https://www.dropbox.com/s/zgckliydalljw6a/sp_data.zip?dl=0). I want to convert polygons to raster. The value of each grid cell should be based on the area covered by polygon. For example, if grid cell size is 100m x 100m (10000m2) and polygons covered 5000m2 of this cell - the value should be 50 (or 0.5). Below, you can find my proposed solution:

Raw data:

library('rgdal')
library('raster')
library('rgeos')
library('maptools')

sp_data <- readOGR(".", "sp_data")
spplot(sp_data)

I've created new raster with 100 meters resolution and extent equal to the polygon bounding box:

r <- raster()
bb <- extent(290000, 300000, 500000, 510000)
extent(r) <- bb
res(r) <- 100

I've also created background polygon - with the values for the areas without forest equal to 0. Afterwards, I've merged these two polygon data:

sp_back <- as(extent(r), "SpatialPolygons")
proj4string(sp_back) <- proj4string(sp_data)
sp_back <- gDifference(sp_back, sp_data)
plot(sp_back, col='red')
plot(sp_data, col='green')
sp_back <- SpatialPolygonsDataFrame(sp_back, data=data.frame(value=rep(0, length(sp_back)), row.names=row.names(sp_back)))
sp_back <- spChFIDs(sp_back, "new_id")
sp_bind <- spRbind(sp_data, sp_back)
spplot(sp_bind)

Additionally, I've created new raster with reduced resolution (500 meters):

r2 <- r
res(r2) <- 500

I've rasterized polygon data to the first raster:

sp_raster <- rasterize(sp_bind, r, field="value", fun=max)
spplot(sp_raster, aspect='iso')

At the end, I've resampled the first raster values into the second raster:

sp_raster2 <- resample(sp_raster, r2, method='bilinear')
spplot(sp_raster2, aspect='iso')

Based on the results, I've have some questions:

1. First of all - is my solution and the result even correct?
2. Conversion between vector and raster somehow simplify my data (for example, loss of small polygons) and probably add some errors. The result of the resampling is affected by this simplification. Are there any alternative solution to calculate the share of polygon area in the grid cell?

Answers using R, Grass or SAGA will be welcome.

Answer 1

Since you put GRASS into the tags, here's a solution based on GRASS:

First, you need to know exactly what coordinate system the original data is in (as always with GRASS). I see that the *.prj file contains "TRANSVERSE MERCATOR" but it's not one of the standard UTM zones. Since you have not mentioned the EPSG code of this data, here's the proj4 string for creating a new matching GRASS Location:

+proj=tmerc +lat_0=0 +lon_0=19 +k=0.9993 +x_0=500000 +y_0=-5300000 +ellps=WGS84 +units=m +no_defs

Now start GRASS, create a LOCATION (based on the above) and a MAPSET.

And now import the shapefile into GRASS by running:

v.in.ogr input="sp_data.shp" output=sp_data

Add a column for the area of each feature, and calculate its area:

v.db.addcolumn map=sp_data columns="area_sqm INT"
v.to.db sp_data option=area column=area_sqm unit=meters

Now convert to raster, using the area_sqm column as the raster value. But first set the GRASS computational region:

g.region -p vector=sp_data res=100

Now create a vector grid of polygons at the raster resolution:

v.mkgrid map=grid

Intersect with the forest polygons, and calculate areas of these intersect polygons: Intersect with the forest polygons, dissolve the polygons inside the same grid cell, calculate areas of these intersect polygons, and move the area value to the grid cells:

v.overlay ain=grid bin=sp_data operator=and output=forest_grid_intersect
v.dissolve --overwrite input=forest_grid_intersect@jn column=a_cat output=forest_grid_intersect_dissolve
v.db.addtable map=forest_grid_intersect_dissolve columns="area_sqm INT"
v.to.db forest_grid_intersect_dissolve option=area column=area_sqm unit=meters 
v.db.join map=grid@jn column=cat other_table=forest_grid_intersect_dissolve other_column=cat

v.db.addcolumn map=forest_grid_intersect columns="area_sqm INT" v.to.db forest_grid_intersect option=area column=area_sqm unit=meters

and convert to raster, using the area attribute as the raster value:

v.to.rast input=grid type=area output=sp_data_rast use=attr attribute_column=area_sqm

And you're done. Hope I have it right this time :-) .

Answer 2

I was just trying to figure out how to do the same thing in r. Here is the approach I came up with:

-Create a raster that has the dimensions and cell size you want

-Use the function 'rasterize' from the raster package, using the "getCover" option

library(raster)
wetlands <- readShapePoly("data/wetlands.shp", proj4string = CRS("+init=EPSG:26907))

r <- raster(ncols = 200, nrows = 200, resolution = 5, crs = CRS("+init=EPSG:26907))

wetlandsAREA <- rasterize(wetlands, r, getCover = TRUE)

I think this essentially does the the same process as in the original post. The rasterize function divides the raster cell into 100 subunits, and the output is the proportion of those that were covered by the polygon.