Rasterise sf Polygon and Assign Cell Value based on Total Area Deforested

Question

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?

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.

library(terra)
library(sf)

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)),
               "POLYGON")
  
  plot(rdf)
    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])))
    plot(b["deforest"])