subset raster to separate two classes

Question

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)

What did your code do? Did you try comparing to "1" to create land/sea? — Spacedman, Apr 21 '23 at 09:24

dieghernan · Accepted Answer · 2023-04-21T10:58:52.840

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

Looks great, however, as I zoom in I see a lot of cells that appear to be land. How can I aggregate all cells within a "12 meter" distance of each other creating a single large SpatRaster representing continuous water that can be separated from the rest of the cells? — Salvador, Apr 21 '23 at 20:14
To be clear, I want a raster with continuous water. So, group all the cells that are less than 12 and make a final raster. Anything greater than 12 should be excluded. Can I use `clamp` again here? something like `clamp(xx_water, upper = 12, value = FALSE` ? — Salvador, Apr 22 '23 at 02:08
It is better to open a new question with all the details on this and also adding the code you try, the results and why that is not your expected result — dieghernan, Apr 22 '23 at 06:52

score 0 · Answer 2 · answered Apr 21 '23 at 11:34

Read the data:

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

Create a binary 0/1 raster where land is 1 and sea is 0:

landsea = xx > 1

Convert to polygons (takes a while):

pol = as.polygons(landsea)

Show (also takes a while):

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

Note your two "polygons" here are very complex multi-polygons with each row of pol representing all the little areas of that value.

To do this final step of cropping the raster to the water extent, you don't actually have to go through this polygonisation process, just set all the land values to the missing value and cut complete rows and column from all sides - there's probably a function in terra to do that - yes its called trim(), and you can pass it a value to count as missing.

subset raster to separate two classes

2 Answers2