0

I have a land use raster (lu_raster) and several thousand shapefiles (counties). Within each county, I'd like to calculate the proportion of agricultural land use. So, I think I need to find a way to 1). sum all the raster pixels that are either value 81 or 82 (see here for reference) within each of the shapefiles, and 2). divide that by the total number of pixels in each of the county shapefiles.

I've been using exact_extract to do this kind of stuff. So I think the second part of my code will look something like this:

Results <- exact_extract(lu_raster, counties, 'count')

But I can't figure out how to do those counts conditionally (only counting the 81s and 82s).

0

Figured it out.

First, reclassify raster making agricultural values 1 and everything else 0:

agri_land_use_ras <- raster::reclassify(lu_raster, lu_matrix)

See here for how to set up a reclassification matrix.

Then, can just take the mean of the resulting raster:

Results <- exact_extract(agri_land_use_ras, counties, 'mean')
3
  • Classifying the raster first is conceptually easy but is also an unnecessary step. – Jeffrey Evans Apr 26 at 15:37
  • Why is it unnecessary? The only way that mean will give a land cover proportion is if the land cover of interest is classified as 1 and everything else is classified as 0. – pfadenhw Apr 26 at 23:35
  • Look at my answer, you can pull the number of pixels associated with given value(s), on the fly, in a function that can be passed to extract or operate on the extracted pixel values object. So, classifying the raster first is not necessary and, if large, could take sometime. – Jeffrey Evans Apr 27 at 0:29
0

It is just a matter of writing a function that returns the proportions and using lapply to apply it to the list returned from exact_extract.

Add libraries and create example data

library(sf)
library(raster)
library(exactextractr)

nc <- st_read(system.file("shape/nc.shp", package="sf"))
  nc <- st_cast(nc, "POLYGON")


r <- raster::raster(matrix(sample(c(81, 82,11, 21, 22, 41), 
                    500^2, replace=TRUE), 500, 500))
    extent(r) <- extent(nc)
    
plot(r)
  plot(st_geometry(nc), add=TRUE)

Now, extract the raster landcover data. We can then use lapply to operate on the resulting list object. The use of which and %in% allows us an easy way to apply a query.

e <- exact_extract(r, nc)

( ag <- lapply(e, function(x) { length(which(x$value %in% c(81,82))) / nrow(x) } ) )

nc$ag <- unlist(ag)
  plot( nc["ag"] ) 

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.