I'm a relatively new to working with spatial data in R.

I'm working with Deforestation data (EPSG: 5880) for the Brazilian Amazon - consisting of 70,000 polygons of deforested land that covers all of Brazil. I'm only interested in polygons within 30 municipalities (which I have unioned and buffered - see below), so I've subsetted the deforestation data to these municipalities only.

The next step is to rasterize this data and calculate distance of each centroid to the municipality border (picture below) - where the value for each raster cell corresponds to the Total Area Deforested within that cell (i.e Total Area of the cell covered by deforestation polygons). Each black dot in the picture is a polygon of deforested land.

Unfortunately I'm struggling to rasterize this. I've tried st_rasterize, without success. I've also tried rasterizing the Municipality area through:

Municipalities_sf <- st_read("Munic.shp")
Deforest_sf <- st_read("Def.shp")

Municipalities <- as(Municipalities_sf, "Spatial")

r <- raster(nrow = 1000, ncol = 1000)
extent(r) <- extent(Municipalities)
Mun_Raster <- rasterize(Municipalities, r)

This rasterizes the Municipality polygon, but I can't figure out how to assign a value to each raster cell based on the Total Area Deforested - possibly by combining the above with st_intersection, but I'm not too sure.

How do I rasterize this Deforestation data within the Municipality boundary?

Deforestation within the Unioned Municipalities

1 Answer 1


We will use the terra library to rasterize the data. You have a few questions nested in your post so, it would be helpful to clarify your thinking and focus on exactly what you are after here. I will illustrate how you can rasterize your data and leave your distance question for another time.

Add libraries.


Read "deforest" data as an sf object. We will then create an empty raster "r" to act as our reference for rasterizing the deforestation data.

deforest <- st_read("Def.shp")
  r <- rast(ext(deforests), nrow = 1000, ncol = 1000)

Assuming that you have an attribute indicating area (eg., hectares) for each polygon, we can use it to have our raster values indicate area. Or, you can just output [0,1] and tally the "1" cells in another analysis. I an using "deforest_area" as a place holder for your attribute. Note; you have to coerce your sf object into a terra vect object.

Using an attribute file to rasterize

deforest_rast <- rasterize(vect(deforest), r, field="deforest_area")

Using uniform output "1" with a background value "0". The advantage of using a uniform value is that you can simply mask the data and return a cell count (see below).

deforest_rast <- rasterize(vect(deforest), r, values=1, background=0)

Let's perform a dummy analysis to illustrate a workflow

Here we just pretend that the nc polygons represent deforested areas and the buffer that we create are the municipalities, or any other aggregate unit.

Read the polygons (here we grab a random sample for partial coverage of the study area).

deforest <- st_read(system.file("shape/nc.shp", package="sf"))
  deforest <- deforest[sample(1:nrow(deforest), 20),]

Now, create our reference raster. We will rasterize the "AREA" field of the polygons, but just for illustration. We will be using a binary representation for analysis.

r <- rast(ext(deforest), nrow=1000, ncol=1000)
  plot(rdf_area <- rasterize(vect(deforest), r, field="AREA"))

We create a raster where deforested areas are one else zero. We will also create some buffers (pretend municipalities) to act as aggregate polygons to get cell frequencies.

rdf <- rasterize(vect(deforest), r, value=1, background=0)
  b <- st_cast(st_as_sf(st_buffer(st_sample(deforest, 5), 40000)),
    plot(st_geometry(b), add=TRUE)


Here,we can take a quick look at a subset for the first aggregate polygon.

plot(mask(crop(rdf, ext(b[1,])), vect(b[1,])))  

Now, things get a bit more complex. We will use lapply (to act as a loop) and iterate through the polygons to get cell counts for each [0,,1] class.

( f <- lapply(1:nrow(b), function(i){   
             freq(mask(crop(rdf, ext(b[i,])), 
             vect(b[i,]))) }) )

So, we can see our deforested cell counts for each polygon. We can now calculate percentages and add the results to our polygons.

b$deforest <- unlist(lapply(f, function(x) x[,3][2] / (x[,3][1] + x[,3][2])))
  b$deforest <- unlist(lapply(f, function(x) x[,3][2] / (x[,3][1] + x[,3][2])))

deforest fraction

  • Many thanks for your detailed answer (and yes, I do have a field for deforested area) - unfortunately all is good until deforest_rast, which returns an error: "unable to find inherited method for function 'rasterize' for signature "SpatVector", "RasterLayer". Could you please advise?
    – cstock
    Apr 18, 2022 at 9:58
  • First, make sure that the terra library is added then, you may have to explicitly call the function from the environment. Occasionally, when you call a function that also exists in the raster library, you get an object error. Not sure why but, I have seen this behavior with a few terra functions. The workaround is to just use terra::rasterize Always look at your the class of your object class(x)` For terra you need to read the data using terra::rast and not the raster functions eg., raster, stack, brick. The SpatVector class accounts for single and multi band. Apr 18, 2022 at 16:00

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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