您可以使用terra轻松完成此操作，避免转换和潜在错误。请参阅https://gis.stackexchange.com/questions/421821/how-to-subset-a-spatraster-by-value-in-r，terra开发人员对此进行了进一步解释。
一个mock文件的例子：

library(terra)
#> terra 1.7.23

# Mock data
f <- system.file("extdata/asia.tif", package = "tidyterra")
xx <- rast(f)
# End of Mock data

# xx <- rast("MyRaster.asc")

xx
#> class       : SpatRaster 
#> dimensions  : 232, 432, 1  (nrow, ncol, nlyr)
#> resolution  : 22550.66, 22512.94  (x, y)
#> extent      : 7619120, 17361007, -1304745, 3918256  (xmin, xmax, ymin, ymax)
#> coord. ref. : WGS 84 / Pseudo-Mercator (EPSG:3857) 
#> source      : asia.tif 
#> name        : file44bc291153f2 
#> min value   :        -10071.50 
#> max value   :          6064.73
ncell(xx)
#> [1] 100224
plot(xx)

# Clamp it

# Water
xx_water <- clamp(xx, upper = 1, value = FALSE)
xx_water
#> class       : SpatRaster 
#> dimensions  : 232, 432, 1  (nrow, ncol, nlyr)
#> resolution  : 22550.66, 22512.94  (x, y)
#> extent      : 7619120, 17361007, -1304745, 3918256  (xmin, xmax, ymin, ymax)
#> coord. ref. : WGS 84 / Pseudo-Mercator (EPSG:3857) 
#> source(s)   : memory
#> name        : file44bc291153f2 
#> min value   :    -1.007150e+04 
#> max value   :     9.916311e-01

plot(xx_water, col = hcl.colors(10, "Blues2"))

# Land
xx_land <- clamp(xx, lower = 1, value = FALSE)
xx_land
#> class       : SpatRaster 
#> dimensions  : 232, 432, 1  (nrow, ncol, nlyr)
#> resolution  : 22550.66, 22512.94  (x, y)
#> extent      : 7619120, 17361007, -1304745, 3918256  (xmin, xmax, ymin, ymax)
#> coord. ref. : WGS 84 / Pseudo-Mercator (EPSG:3857) 
#> source(s)   : memory
#> name        : file44bc291153f2 
#> min value   :         1.009122 
#> max value   :      6064.730469

plot(xx_land, col = hcl.colors(10, "BuGn"))

创建于2023-04-21带有reprex v2.0.2
我尝试了你的文件（~1.3Gb），这个过程可能需要一点时间，但上面的代码应该仍然可以工作。只看水的结果：

library(terra)
#> terra 1.7.23

xx <- rast("del/sfbaydeltadem10m2016.asc")

xx

# class       : SpatRaster 
# dimensions  : 15795, 13588, 1  (nrow, ncol, nlyr)
# resolution  : 10, 10  (x, y)
# extent      : 520545, 656425, 4136705, 4294655  (xmin, xmax, ymin, ymax)
# coord. ref. : +proj=utm +zone=10 +datum=NAD83 +units=m +no_defs 
# source      : sfbaydeltadem10m2016.asc 
# name        : sfbaydeltadem10m2016 
# min value   :            -114.8981 
# max value   :             432.1640 
ncell(xx)
#> [1] 214622460

# Clamp it

# Water
xx_water <- clamp(xx, upper = 1, value = FALSE)
plot(xx_water, col = hcl.colors(10, "Blues2"))

读取数据：

library(terra)
xx = rast("sfbaydeltadem10m2016.asc")

创建二进制0/1栅格，其中陆地为1，海洋为0：

landsea = xx > 1

转换为多边形（需要一段时间）：

pol = as.polygons(landsea)

显示（也需要一段时间）：

plot(pol, col=c("red","blue")

请注意，这里的两个“多边形”是非常复杂的多多边形，每行pol表示该值的所有小区域。
要完成将栅格裁剪到水域范围的最后一步，您实际上不必执行此多边形化过程，只需将所有土地值设置为缺失值并从所有侧面切割完整的行和列-terra中可能有一个函数可以做到这一点-是的，它称为trim()，您可以向其传递一个值以计数为缺失。