匹配 dbf 表中的栅格值并基于栅格进行过滤

问题描述 投票：0回答：1

我有一个覆盖澳大利亚的大栅格，我将 

crop()
和 
mask()
应用于使用多多边形 sf （嵌套集水区）。我使用 
foreign
包读取 dbf 表，使用 
terra::rast()
读取 tif。然后我确保每个文件都有相同的 CRS。

当我尝试将 dbf 与栅格中的单元格匹配时，它没有给出我所期望的结果，因为我的栅格没有价值。它只包含“森林分类（本地和非森林）。我需要将 dbf 表中的数据与我的栅格链接起来，以便我可以过滤单元格。

例如，我只想从 Forest_CATEGO 列中选择“原生森林”单元格，以及“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)

链接到我的数据

r raster dbf unique-values multipolygons
1个回答
0
投票

您的 SpatRaster 有一层，但您使用 

list(dbf_table)
添加了多个级别。根据您的描述，最好使用两列数据框简单地分配级别。

这是一个假设的表示：

  • dbf_table 是 SpatRaster 中每个值组合的摘要，例如SpatRaster 和 dbf_table 中有 10 个不同的类别（或者实际数据中存在的许多类别）
  • SpatRaster 中的值 == dbf_table 中的值

由于

head()
仅返回dbf_table的前六行，因此下面的解决方案将仅使用问题中所示的值（n = 6）。

首先，加载所需的包并创建示例数据：

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())

1

如果您的输出需要更多细微差别，您可以执行以下操作：

# 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 范围内： 2

最新问题
© www.soinside.com 2019 - 2024. All rights reserved.