I am trying to find the area of forest (conserved patches) around crop fields. For example, I would like to know the area of forest within a 500m or 5km radius from fields within a county. My data is downloaded as a raster. Below, I describe what I have done so far. However, I am mostly curious about other approaches to do so (possibly using only raster/terra based methods) and that may not require intersect(). It has been suggested to me elsewhere that not going for intersections would be a better way to do what I am trying. I am interested in solutions that are not memory hungry, if possible, since most of the issues I have encountered so far are how to run my code in my notebook given that doing sf::st_buffer() and sf::st_intersection() is very memory intensive. The code I provide above is simplified and should not be very memory hungry. In sum, I am hoping for advice on how to get the area of forests around blueberry fields even (or specially) if it does not employ solutions similar to the code I already have.
CODE TO GET THE DATA I AM USING:
library(CropScapeR)
library(httr)
library(raster)
library(terra)
library(sf)
#THIS IS JUST TO GET THE DATA I AM WORKING WITH, THE CROPLAND DATA LAYER BY USDA
#To download CDL data on Linux using CropscapeR i need to do the below process as per
#the instructions in the package page (in Windows you can basically go directly to)
#GetCDLData(.) part
# Skip the SSL check
httr::set_config(httr::config(ssl_verifypeer = 0L))
# Automatically generate a temporary path to save the data
tif_file <- tempfile(fileext = '.tif')
ST_CDL <- GetCDLData(aoi = '34001', year = 2021, type = 'f', save_path = tif_file)
#I can write my data using terra or raster. Due to memory constraints I have been
#using the terra version (below) lately
terra::writeRaster(ST_CDL, "ST_CDL.tif", overwrite=TRUE)
ST_CDL = terra::rast("ST_CDL.tif")
#Creates SpatRasters for conserved areas and for blueberry fields with NA for all
#values except the ones I am interested in
replacements <- rbind(c(63, 63), c(64, 64), c(141, 141), c(142, 142), c(143, 143), c(152, 152))
ST_CDL_conserved <- terra::classify(ST_CDL, replacements, others=NA)
ST_CDL_blueberries <- terra::classify(ST_CDL, cbind(242, 242), others=NA))
Below, is what I have done so far to get at the areas. But as I said before, I am curious for solutions that don't require intersections or the use of sf. My code below is here to help better understand what I want to get at (and also because I am looking for other ways of coding this, but maybe this is indeed the best way to do so). In a nutshell, I polygonize the SpatRasters with conserved patches and blueberry fields with terra::as.polygons(, dissolve=TRUE) and then convert them to sf objects with sf::st_as_sf(). Then, I create buffers around the fields using sf::st_buffer(). I then intersect these buffers with the conserved areas polygon using sf::st_intersection() and calculate their areas.
I have used dissolve=TRUE in the terra::as.polygons(, dissolve=TRUE) step because I want to aggregate all the fields/gridpoints together. If I were to do one grid point at a time, I would get areas that are close to more than one gridpoint showing up in the area calculation more than once. That's also what have kept me from using terra::buffer and terra::intersect, because I could not figure out how to create the buffers and intersect them to forest without double counting areas. I also feel that dissolve=TRUE would make my code run faster and use less memory in the sf::st_buffer() and sf::st_intersection() steps, but I don't know if this is actually true.
#Creates SpatRasters with NA for all other values except the ones I am interested in
replacements <- rbind(c(63, 63), c(64, 64), c(141, 141), c(142, 142), c(143, 143), c(152, 152))
ST_CDL_conserved <- terra::classify(ST_CDL, replacements, others=NA)
ST_CDL_blueberries <- terra::classify(ST_CDL, cbind(242, 242), others=NA))
#Polygonizes the rasters and creates the sf objects
ST_CDL_conserved <- terra::as.polygons(ST_CDL_conserved, dissolve = TRUE)
ST_CDL_conserved <- sf::st_as_sf(ST_CDL_conserved)
ST_CDL_blueberries <- terra::as.polygons(ST_CDL_blueberries, dissolve = TRUE)
ST_CDL_blueberries <- sf::st_as_sf(ST_CDL_blueberries)
#Now I apply sf based methods to create a buffer and intersect it with forests
ST_CDL_blueberries_buffer <- sf::st_buffer(ST_CDL_blueberries, 500)
#Now I intersect it with forests
intersection <- sf::st_intersection(ST_CDL_blueberries_buffer,ST_CDL_conserved)
#Calculate the areas of the intersections. I still need to sum these areas
#together to get at the total area of intersection at the county.
area_intersection <- st_area(intersection)
Thanks for your help!
Here is how I show how you can do this. You need terra > 1.6.41 for this to work well.
Get the data
library(CropScapeR)
library(httr)
library(terra)
#terra 1.6.41
httr::set_config(httr::config(ssl_verifypeer = 0L))
tif_file <- tempfile(fileext = '.tif')
ST_CDL <- GetCDLData(aoi = '34001', year = 2021, type = 'f', save_path = tif_file)
terra::writeRaster(ST_CDL, "ST_CDL.tif", overwrite=TRUE)
ST_CDL = terra::rast("ST_CDL.tif")
replacements <- rbind(c(63, 63), c(64, 64), c(141, 141), c(142, 142), c(143, 143), c(152, 152))
r_conserved <- terra::classify(ST_CDL, replacements, others=NA)
r_blueberries <- terra::classify(ST_CDL, cbind(242, 242), others=NA)
Vector approach (here with "terra")
p_conserved <- terra::as.polygons(r_conserved)
p_blueberries <- terra::as.polygons(r_blueberries)
p_blueberries_buffer <- terra::buffer(p_blueberries, 500)
intersection <- terra::intersect(p_blueberries_buffer, p_conserved)
p_area <- terra::expanse(intersection)
sum(p_area)
#[1] 242535343
Raster approach
r_blueberries_buffer <- terra::buffer(r_blueberries, 500, background=NA)
m <- mask(r_conserved, r_blueberries_buffer)
expanse(m, transform=FALSE)
#[1] 235463400
The results are not exactly the same. I assume that the reason for that is that the buffer on a raster considers entire cells only, they are either in or out, whereas the polygon raster does not do that. The latter may seem more precise, but that is debatable given that your original data are raster data.