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!
