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:
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:
- First of all - is my solution and the result even correct?
- 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.