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:
- 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.
Best Answer
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:
Now start GRASS, create a LOCATION (based on the above) and a MAPSET.
And now import the shapefile into GRASS by running:
Add a column for the area of each feature, and calculate its area:Now convert to raster, using the area_sqm column as the raster value.But first set the GRASS computational region:Now create a vector grid of polygons at the raster resolution:
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.db.addcolumn map=forest_grid_intersect columns="area_sqm INT" v.to.db forest_grid_intersect option=area column=area_sqm unit=metersand convert to raster, using the area attribute as the raster value:
And you're done. Hope I have it right this time :-) .