Heidi Rodenhizer
Heidi Rodenhizer

Reputation: 119

Layer selection is ignored using terra::as.polygons

I am using terra version 1.7-29 (previously had the same issue with 1.6-17).

My Data

I have a two-layer raster dataset (download it here):

prediction <- terra::rast('path/to/raster.tif')
names(prediction) <- c('prediction', 'validation')

It looks like this:

class       : SpatRaster 
dimensions  : 202, 294, 2  (nrow, ncol, nlyr)
resolution  : 2, 2  (x, y)
extent      : 1978598, 1979186, 889666, 890070  (xmin, xmax, ymin, ymax)
coord. ref. : WGS 84 / NSIDC Sea Ice Polar Stereographic North (EPSG:3413) 
source      : 000.tif 
names       : prediction, validation 

The original raster dataset

The Problem

I would like to create a vector polygon dataset from the 2nd layer, 'validation', but as.polygons() seems to be working from the 1st layer, 'prediction'.

val_poly <- as.polygons(subset(prediction, 'validation'))
plot(val_poly)

Why didn't the polygon get made from the 'validation' layer? This looks more like prediction > 0.5.

Other things I've tried

I have also tried subsetting the validation layer first and then providing that to as.polygons() (allowing me to check that the subset is working as expected) but it made no difference in the polygon output.

val_layer <- subset(prediction, 'validation')
val_poly <- as.polygons(val_layer)

subset(prediction, 'validation')
plot(subset(prediction, 'validation'))
class       : SpatRaster 
dimensions  : 202, 294, 1  (nrow, ncol, nlyr)
resolution  : 2, 2  (x, y)
extent      : 1978598, 1979186, 889666, 890070  (xmin, xmax, ymin, ymax)
coord. ref. : WGS 84 / NSIDC Sea Ice Polar Stereographic North (EPSG:3413) 
source      : 000.tif 
name        : validation 

This is the right layer

It's Not Reproducible

(with a handmade SpatRaster, at least)

When I make a 2-layer SpatRaster by hand and convert the 2nd layer to polygon, it works as expected. This has me considering two options:

  1. There may be an issue with my dataset, or
  2. there is something going wrong related to how terra handles files on disk (but this is way outside of my coding knowledge comfort-zone).
r1 <- rast(ncol=4, nrow=4)
names(r1) <- c('r1')
r2 <- deepcopy(r1)
names(r2) <- c('r2')
set.seed(0)
values(r1) <- c(0, 0, 0, 0,
                0, 1, 0, 0,
                0, 1, 1, 0,
                0, 0, 0, 0)
set.seed(1)
values(r2) <- c(0, 0, 0, 1,
                0, 0, 1, 1,
                0, 0, 1, 1,
                0, 0, 0, 0)
r <- c(r1, r2)
plot(r)

The test dataset

test <- as.polygons(r[['r2']])
plot(test)

enter image description here

Upvotes: 1

Views: 108

Answers (1)

Robert Hijmans
Robert Hijmans

Reputation: 47481

The below works with terra 1.7-32 (currently the development version); but with the current CRAN version there is an issue that needs a simple work-around.

library(terra)
# terra 1.7.29

x <- rast("000.tif")
x
#class       : SpatRaster 
#dimensions  : 202, 294, 2  (nrow, ncol, nlyr)
#resolution  : 2, 2  (x, y)
#extent      : 1978598, 1979186, 889666, 890070  (xmin, xmax, ymin, ymax)
#coord. ref. : WGS 84 / NSIDC Sea Ice Polar Stereographic North (EPSG:3413) 
#source      : 000.tif 
#names       : 000_1, 000_2 

x2 <- x[[2]] 
#x2 <- x[["000_2"]]
#x2 <- subset(x, "000_2")

## Work around for current version
x2 <- x2 * 1

p <- as.polygons(x2)
plot(x2); lines(p)

enter image description here You probably want a polygon for the cell values that are one, but not for the cells that are zero. In that case, you can do

x2 <- subst(x2, 0, NA)
p <- as.polygons(x2)

The problem is related to this being a file with a particular data type. That is why you could not reproduce it with a SpatRaster that has values in memory. Thank you for revealing this bug!

The development version is available from

install.packages('terra', repos='https://rspatial.r-universe.dev')

Upvotes: 1

Related Questions