我试图找到农作物田周围的森林面积(保护斑块)。例如，我想知道一个县内距离田野500米或5公里半径内的森林面积。我的数据以栅格格式下载。下面，我将描述到目前为止我所做的工作。然而，我最好奇的是其他方法(可能只使用基于栅格/terra的方法)，这可能不需要intersect()。有人向我建议，不走十字路口会是一个更好的方式来做我正在尝试的事情。如果可能的话，我对不需要内存的解决方案感兴趣，因为到目前为止我遇到的大多数问题都是如何在我的笔记本中运行我的代码，因为做sf::st_buffer()和sf::st_intersection()是非常内存密集的。我上面提供的代码是简化的，不应该占用太多内存。总之，我希望关于如何获得蓝莓田周围的森林面积的建议，即使(或特别)如果它不采用类似于我已经拥有的代码的解决方案。
获取我正在使用的数据的代码:
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))
下面，是我到目前为止所做的，以获得的区域。但正如我之前所说，我对不需要交集或使用sf的解决方案很好奇。我下面的代码是为了帮助更好地理解我想要得到什么(也因为我正在寻找其他的编码方式，但也许这确实是最好的方法)。简而言之，我使用terra::as.polygons(, dissolve=TRUE)将带有保守补丁和蓝莓地的spatraster多边形化，然后使用sf::st_as_sf()将它们转换为sf对象。然后，我使用sf::st_buffer()在字段周围创建缓冲区。然后，我使用sf::st_intersection()将这些缓冲区与保守区域多边形相交，并计算它们的面积。
我在terra::as.polygons(, dissolve=TRUE)步骤中使用了dissolve=TRUE，因为我想将所有字段/网格点聚合在一起。如果我一次处理一个网格点，我将得到与多个网格点接近的区域，这些区域在面积计算中不止一次出现。这也是让我不使用terra::buffer和terra::intersect的原因，因为我不知道如何创建缓冲区并在没有重复计算区域的情况下将它们相交到森林。我也觉得dissolve=TRUE会让我的代码运行得更快，在sf::st_buffer()和sf::st_intersection()步骤中使用更少的内存，但我不知道这是否真的。
#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)
谢谢你的帮助!
我将演示如何做到这一点。你需要土地。1.6.41，
获取数据
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)
矢量法(这里有"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
栅格方法
r_blueberries_buffer <- terra::buffer(r_blueberries, 500, background=NA)
m <- mask(r_conserved, r_blueberries_buffer)
expanse(m, transform=FALSE)
#[1] 235463400
结果并不完全相同。我认为这样做的原因是栅格上的缓冲区只考虑整个细胞，它们要么在里面要么在外面，而多边形栅格不这样做。后者可能看起来更精确，但考虑到您的原始数据是栅格数据，这是有争议的。