I have a raster that is too large to post here but can be downloaded from here: https://login.filesanywhere.com/fs/v.aspx?v=8c6c688b586472bcab6a . When I plot it I see a mix of water and land. I would like to separate it with the following criteria: If the values are greater than 1 then is land and if values are less than 1 then is water.
Then, convert the water class to polygon feature. Then use it to clip the original raster to show only water. Any help is very appreciated.
Here is what I tried:
library(terra)
library(mapview)
library(raster)
xx <- rast("MyRaster.asc")
xx
ncell(xx)
plot(xx)
#Convert to raster first to vizualise with mapview
xxx <- raster(xx)
mapview(xxx)
You can do that easily with terra, avoiding conversions and potential errors. See https://gis.stackexchange.com/questions/421821/how-to-subset-a-spatraster-by-value-in-r where this is further explained by the terra developer.
An example with a mock file:
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"))
Created on 2023-04-21 with reprex v2.0.2
I tried with your file (~1.3Gb) and the process can take a bit of time, but the code above should still work. See the result for water only:
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"))