英文:
Extracting values from one raster based on condition in another raster and distinguished by polygons in a vector file
问题
我有一个关于土地覆盖的栅格文件(30x30米像素),一个关于土壤有机碳(SOC)的栅格文件(250x250米像素),以及一个包含不同大小连续多边形的矢量文件。
我需要的是每个多边形的平均SOC值。但只基于具有特定土地覆盖(草地)的像素。
首先,我确保这两个栅格文件具有相同的范围和分辨率,所以我将SOC地图缩放到30x30米。然后,我将矢量地图转换为SpatRaster(也许不是必要的),并意识到其范围与两个栅格文件略有不同(矢量文件加载到R时具有“边界框”而不是范围,也许与此相关)。我还叠加了土地覆盖和SOC文件。
我现在卡住了。首先,我不确定矢量文件是否具有不同的范围是否是问题,如果是问题,如何将其与栅格文件对齐。其次,我不确定是否首先将矢量文件转换为SpatRaster,还是保持为矢量文件就可以了。第三,我应该使用什么代码来提取每个多边形的草地的平均SOC值(即条件if子句)。
英文:
I have one raster file on land cover (30x30m pixels), one raster file on soil organic carbon (SOC) (250x250m pixels), and one vector file with contiguous polygons of different sizes.
What I require is a mean SOC value for each polygon. But only based on pixels that have a select land cover, here grassland.
I first made sure that the two raster files have the same extent and resolution, so I scaled the SOC map to 30x30m. I then turned the vector map into a SpatRaster (perhaps not needed), and realised that the extent is slightly different from the two raster files (the vector file when loaded into R has a 'bounding box' rather than extent, perhaps it relates to that). I also stacked the land cover and SOC files.
I am stuck here now. First, I am not sure whether it is an issue that the vector file has a different extent, and if so, how to align it with the raster files. Second, I am not sure whether I had to to convert the vector file to a SpatRaster in the first place, or if keeping it as a vector file is just fine. Third, what code could I use to now extract the mean SOC values for grassland for each polygon (i.e. a conditional if clause).
答案1
得分: 1
以下是翻译好的代码部分:
# 示例数据如下
library(terra)
pols <- vect(system.file("ex/lux.shp", package="terra"))
carbon <- rast(system.file("ex/elev.tif", package="terra"))
names(carbon) <- "carbon"
landuse <- round(carbon/100)
# 解决方案,假设草地在'landuse'中的值为3
grass <- subst(landuse, 3, 1, others=NA)
grass_carbon <- mask(carbon, grass)
e <- extract(grass_carbon, pols, mean, na.rm=TRUE, ID=FALSE)
data.frame(pols$NAME_2, e)
# pols.NAME_2 carbon
#1 Clervaux 344.6667
#2 Diekirch 305.2462
#3 Redange 300.0202
#4 Vianden 295.0000
#5 Wiltz 329.2917
#6 Echternach 312.9183
#7 Remich 281.1053
#8 Grevenmacher 294.5238
#9 Capellen 321.3462
#10 Esch-sur-Alzette 300.4255
#11 Luxembourg 299.0872
#12 Mersch 301.0226
# 如果多边形内部的单元格尺寸差异相对较小,这样就足够了。否则,您可以使用'zonal'函数,将单元格尺寸视为权重。如果地表覆盖类别表示为格网单元格的分数,也可以使用它。
# 要绘制结果,您可以执行以下操作
pols$carbon <- e
plot(pols, "carbon")
希望这有所帮助。如果您有其他问题,欢迎提出。
英文:
Here is some example data
library(terra)
pols <- vect(system.file("ex/lux.shp", package="terra"))
carbon <- rast(system.file("ex/elev.tif", package="terra"))
names(carbon) <- "carbon"
landuse <- round(carbon/100)
Here is a solution. Assuming that grassland has value "3" in landuse
grass <- subst(landuse, 3, 1, others=NA)
grass_carbon <- mask(carbon, grass)
e <- extract(grass_carbon, pols, mean, na.rm=TRUE, ID=FALSE)
data.frame(pols$NAME_2, e)
# pols.NAME_2 carbon
#1 Clervaux 344.6667
#2 Diekirch 305.2462
#3 Redange 300.0202
#4 Vianden 295.0000
#5 Wiltz 329.2917
#6 Echternach 312.9183
#7 Remich 281.1053
#8 Grevenmacher 294.5238
#9 Capellen 321.3462
#10 Esch-sur-Alzette 300.4255
#11 Luxembourg 299.0872
#12 Mersch 301.0226
This is good enough if differences in cell sizes within polygons are relatively small. Otherwise, you could use zonal
with cell sizes as weights. You could also use that if the land cover classes were expressed as fractions of the grid cells.
To plot the results, you can do
pols$carbon <- e
plot(pols, "carbon")
通过集体智慧和协作来改善编程学习和解决问题的方式。致力于成为全球开发者共同参与的知识库,让每个人都能够通过互相帮助和分享经验来进步。
评论