将R中的栅格值与dbf表进行匹配，并根据栅格进行过滤
问题描述 投票：0回答：1
如果有人能指导我，我将不胜感激。 我在澳大利亚全国有一个大栅格。我用多个多边形（嵌套集水区）掩盖了栅格。我使用foreign包读取dbf表，我扫描tif，然后匹配crs，现在当我想将dbf与栅格中的单元格匹配时，它没有给我我所期望的，因为我的栅格没有价值。它只包含“森林分类（本地和非森林）。 我需要它的原因是我想将 dbf 表与我的栅格链接起来，然后使用条件要求 R 过滤具有特定属性的单元格，这些属性在包含很多列的 dbf 表中。我想说的是，仅从 Forest_Category 列、“HT_CODE”和“COVER_CODE”列 %in% c(1,2,3) 中选择“本地森林”单元格。然后我想要求R只过滤其区域内包含超过60％这种条件的集水区（这里我应该考虑集水区是嵌套的）我不知道这种过滤是否可能？这是一个视图嵌套集水区.
这是我的代码：
library(terra)
#> terra 1.7.78
library(foreign)
library(dplyr)
#>
#> Attaching package: 'dplyr'
#> The following objects are masked from 'package:terra':
#>
#> intersect, union
#> The following objects are masked from 'package:stats':
#>
#> filter, lag
#> The following objects are masked from 'package:base':
#>
#> intersect, setdiff, setequal, union
install.packages("foreign")
#> Warning: package 'foreign' is in use and will not be installed
# File paths
raster_file <- "C:/Users/for23.tif"
dbf_file <- "C:/Users/for23.tif.vat.dbf"
shapefile <- "C:/Users/boundry.shp"
# Load the raster, .dbf table, and shapefile
rast <- rast(raster_file)
rast
#> class : SpatRaster
#> dimensions : 38370, 40100, 1 (nrow, ncol, nlyr)
#> resolution : 100, 100 (x, y)
#> extent : -1888000, 2122000, -4847000, -1010000 (xmin, xmax, ymin, ymax)
#> coord. ref. : GDA_1994_Albers
#> source : aus_for23.tif
#> categories : FOR_CATEGO, FOR_TYPE, FOR_CODE, COVER_CODE, HT_CODE, FOREST, STATE, FOR_SOURCE
#> name : FOR_CATEGO
#> min value : Non-forest
#> max value : Native forest
level <- levels(rast)
head(level)
#> [[1]]
#> VALUE FOR_CATEGO
#> 1 1 Non-forest
#> 2 2 Non-forest
#> 3 3 Native forest
#> 4 4 Native forest
#> 5 5 Native forest
#> 6 6 Native forest
#> 7 7 Native forest
#> 8 8 Native forest
#> 9 9 Native forest
#> 10 10 Native forest
dbf_table <- read.dbf(dbf_file)
head(dbf_table)
#> VALUE COUNT FOR_CATEGO FOR_TYPE FOR_CODE COVER_CODE HT_CODE FOREST
#> 1 1 121022073 Non-forest Non-forest 0 6 6 0
#> 2 2 5875 Non-forest Non-forest 0 6 6 0
#> 3 3 17703 Native forest Mangrove 63 3 2 1
#> 4 4 335 Native forest Mangrove 61 1 2 1
#> 5 5 16724 Native forest Rainforest 98 3 2 1
#> 6 6 4555 Native forest Rainforest 96 3 1 1
#> STATE FOR_SOURCE
#> 1 Qld <NA>
#> 2 <NA> <NA>
#> 3 Qld Global Mangroves
#> 4 Qld Global Mangroves
#> 5 Qld NVIS 6.0
#> 6 Qld NVIS 6.0
catchments <- vect(shapefile)
shapefile1<- st_read(shapefile)
#> Error in st_read(shapefile): could not find function "st_read"
if (st_crs(shapefile1) != crs(rast)) {
shapefile_data <- st_transform(shapefile1, st_crs(rast))
}
#> Error in st_crs(shapefile1): could not find function "st_crs"
# Inspect the .dbf table
head(dbf_table)
#> VALUE COUNT FOR_CATEGO FOR_TYPE FOR_CODE COVER_CODE HT_CODE FOREST
#> 1 1 121022073 Non-forest Non-forest 0 6 6 0
#> 2 2 5875 Non-forest Non-forest 0 6 6 0
#> 3 3 17703 Native forest Mangrove 63 3 2 1
#> 4 4 335 Native forest Mangrove 61 1 2 1
#> 5 5 16724 Native forest Rainforest 98 3 2 1
#> 6 6 4555 Native forest Rainforest 96 3 1 1
#> STATE FOR_SOURCE
#> 1 Qld <NA>
#> 2 <NA> <NA>
#> 3 Qld Global Mangroves
#> 4 Qld Global Mangroves
#> 5 Qld NVIS 6.0
#> 6 Qld NVIS 6.0
croped<- crop(rast, vect(shapefile_data))
#> Error in h(simpleError(msg, call)): error in evaluating the argument 'x' in selecting a method for function 'vect': object 'shapefile_data' not found
masked_hrs<- mask(rast, vect(shapefile_data))
#> Error in h(simpleError(msg, call)): error in evaluating the argument 'mask' in selecting a method for function 'mask': error in evaluating the argument 'x' in selecting a method for function 'vect': object 'shapefile_data' not found
levels(masked_hrs)<- list(dbf_table)
#> Error: object 'masked_hrs' not found
# Define conditions
condition <- dbf_table$FOR_CATEGO == "Native forest" &
dbf_table$COVER_CODE %in% c(1, 2, 3) &
dbf_table$HT_CODE %in% c(1, 2, 3)
1个回答
0
投票
投票
您的 SpatRaster 有一层，但您使用
list(dbf_table)这是一个假设的表示：
- dbf_table 是 SpatRaster 中每个值组合的摘要，例如SpatRaster 和 dbf_table 中都有 10 个不同的类别
- SpatRaster 中的值 == dbf_table 中的值
由于
head()首先，加载所需的包并创建示例数据：
library(terra)
library(sf)
library(ggplot2)
library(tidyterra)
# Create example dbf file
dbf_table <- structure(list(VALUE = 1:6, COUNT = c(121022073L, 5875L, 17703L,
335L, 16724L, 4555L), FOR_CATEGO = c("Non-forest", "Non-forest",
"Native forest", "Native forest", "Native forest", "Native forest"
), FOR_TYPE = c("Non-forest", "Non-forest", "Mangrove", "Mangrove",
"Rainforest", "Rainforest"), FOR_CODE = c(0L, 0L, 63L, 61L, 98L,
96L), COVER_CODE = c(6L, 6L, 3L, 1L, 3L, 3L), HT_CODE = c(6L,
6L, 2L, 2L, 2L, 1L), FOREST = c(0L, 0L, 1L, 1L, 1L, 1L), STATE = c("Qld",
"<NA>", "Qld", "Qld", "Qld", "Qld"), FOR_SOURCE = c("<NA>", "<NA>",
"Global Mangroves", "Global Mangroves", "NVIS 6.0", "NVIS 6.0"
)), row.names = c("1", "2", "3", "4", "5", "6"), class = "data.frame")
# Create SpatRaster based on your metadata (lower resolution example)
r <- rast(nrows = 3837, ncols = 4010,
ext = ext(-1888000, 2122000, -4847000, -1010000),
crs = "EPSG:3577")
# Assign values to SpatRaster in order to replicate your actual data
set.seed(123)
values(r) <- sample(dbf_table$VALUE,
size = ncell(r),
replace = TRUE,
prob = 1 / sum(dbf_table$COUNT) * dbf_table$COUNT)
# Set the levels for the SpatRaster
levels(r) <- dbf_table[, c("VALUE", "FOR_CATEGO")]
head(levels(r))
# [[1]]
# VALUE FOR_CATEGO
# 1 1 Non-forest
# 2 2 Non-forest
# 3 3 Native forest
# 4 4 Native forest
# 5 5 Native forest
# 6 6 Native forest
# Create example mask sf
sf_poly <- st_as_sfc(st_bbox(c(xmin = -1.5e6,
ymin = -3.5e6,
xmax = -1.44e6,
ymax = -3.44e6))) %>%
st_as_sf(crs = 3577)
# Mask and crop r
r1 <- crop(r, vect(sf_poly))
r1 <- mask(r1, vect(sf_poly))
考虑到您想要对数据进行子集化（根据您定义的条件），将所有其他值转换为 NA 是有意义的。如果您需要保留所有值/级别，请跳过下一步（请注意，如果您这样做，则可能还需要修改下面的
ggplot()# Keep only wanted values, change others to NA
r1[!(r1 %in% dbf_table[dbf_table$FOR_CATEGO == "Native forest" &
dbf_table$COVER_CODE %in% 1:3 &
dbf_table$HT_CODE %in% 1:3, 1])] <- NA
# Change levels to desired variable e.g. FOR_TYPE
levels(r1) <- dbf_table[, c("VALUE", "FOR_TYPE")]
r1
# class : SpatRaster
# dimensions : 60, 60, 1 (nrow, ncol, nlyr)
# resolution : 1000, 1000 (x, y)
# extent : -1500000, -1440000, -3500000, -3440000 (xmin, xmax, ymin, ymax)
# coord. ref. : GDA94 / Australian Albers (EPSG:3577)
# source(s) : memory
# categories : FOR_TYPE
# name : FOR_TYPE
# min value : Mangrove
# max value : Rainforest
ggplot() +
geom_spatraster(data = r1, maxcell = Inf) +
scale_fill_discrete(na.translate = FALSE,
name = "Type") +
coord_sf(expand = FALSE) +
theme(panel.background = element_blank())
如果您的输出需要更多细微差别，您可以执行以下操作：
# Change levels to desired variable e.g. FOR_CODE
levels(r1) <- dbf_table[, c("VALUE", "FOR_CODE")]
r1
# class : SpatRaster
# dimensions : 60, 60, 1 (nrow, ncol, nlyr)
# resolution : 1000, 1000 (x, y)
# extent : -1500000, -1440000, -3500000, -3440000 (xmin, xmax, ymin, ymax)
# coord. ref. : GDA94 / Australian Albers (EPSG:3577)
# source(s) : memory
# categories : FOR_CODE
# name : FOR_CODE
# min value : 63
# max value : 96
# For legend (illustrative purposes only)
legend <- subset(dbf_table, dbf_table$VALUE %in% as.numeric(unique(values(r1))))
ggplot() +
geom_spatraster(data = r1, maxcell = Inf) +
scale_fill_manual(name = "Type/Code",
labels = paste(legend$FOR_TYPE, legend$FOR_CODE),
values = c("#61D04F" ,"#009292", "darkgreen"),
na.translate = FALSE) +
coord_sf(expand = FALSE) +
theme(panel.background = element_blank())
请注意，即使您定义的条件从本示例中使用的 dbf_table 返回四行，但只有三个不同级别位于 sf_polygon 范围内：
最新问题
- Android Room：如何从关系表中提取数据？
- 如何在Delphi中使用记录实现Builder模式而不创建副本？
- 使用 cv2 和 Streamlit lib 时发生 Streamlit 部署错误
- GDB 可以调试通过没有可用源代码的程序启动的劣质进程吗？
- Configuration.set可以在Mapper中使用吗？
- 每次我尝试运行我的应用程序时，都会收到“运行‘main.dart’时出错：入口点不在 Flutter pub root 内”
- 在UILabel中显示分数
- 通过命令行或 powershell 5.1 调用时，Powershell 7.2 操作不起作用
- 多处理嵌套 DictProxy 对象不允许使用 .items()
- 如何让hadoop put创建不存在的目录
- Hono Zod 类型“json”的参数不可分配给类型“never”的参数
- 如何在 Compose 中查看应用的位图旋转
- nginx 重定向到索引页面
- 如何在Python中计算math.sqrt(-1)？
- Angular 19 独立组件错误：尽管导入了 HttpClientModule，“没有 _HttpClient 的提供程序”
- 不同加数问题贪心算法
- 我无法在 Ajax 中执行 Perl 脚本
- AVX2 / gcc：通过使用不同的寄存器来提高CPU级并行性
- 我无法在ajax中执行perl脚本
- 实体框架。直接编辑集合导航属性时的奇怪行为
© www.soinside.com 2019 - 2024. All rights reserved.