6

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)

enter image description here

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)

enter image description here

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')

enter image description here

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')

enter image description here

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.

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 :-) .

  • Correct Me If I'm Wrong: your results are grid cells with value of overlapping polygon area. For example, we have big, irregular polygon (let say 300sq meters) and if this polygon is overlapped by seven grid cells, each of the grid cell will get a value of 300? – Jot eN Jul 10 '15 at 11:35
  • Based on the example above - The values in the grid cells should be related to the part of the polygon in this grid cell area. So, if the polygon area is 300 m2, and is overlapped by seven grid cells - each cell should have different value (for example, cell#1=100, cell#2=50, cell#3=50, cell#4=35, cell#5=30, cell#6= 20, cell#7=15. – Jot eN Jul 10 '15 at 11:40
  • Yes, you're right. Each cell will have the full area of the polygon that it intersects. I misunderstood your question. I think you can get the part of each raster cell covered by the polygon by first creating a vector of polygon grids at the raster resolution, then intersecting that with the forest polygons. See my edits above. – Micha Jul 11 '15 at 11:10
  • Thank you again for everything you’ve done. You're quite close, but it's still not an answer. The result (based on your last answer) is fine when there is only one polygon in one grid cell and this polygon covers the center of the grid cell. However, there are two unwanted situations: (1) Polygon is not in the center of the grid cell = grid cell value is 'No data' (2) There are more than one polygon in the grid cell = grid cell value is equal to the polygon covering the center of the grid cell. It should be equal to the sum of polygons area. – Jot eN Jul 21 '15 at 9:48
  • 1
    Great! and thanks for adding the edits, so the solution is complete for future. – Micha Jul 24 '15 at 10:49
1

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.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.