rspatial / terra

R package for spatial data handling https://rspatial.github.io/terra/reference/terra-package.html
GNU General Public License v3.0
532 stars 86 forks source link

classify can not work well #1552

Closed SpatLyu closed 2 months ago

SpatLyu commented 2 months ago

I found that classify() did not work properly when testing a small demo with terra. I need to classify values greater than 1.96 as 1; values greater than or equal to -1.96 and less than or equal to 1.96 are classified as 2; values less than -1.96 are classified as 3, here are a few examples:

library(terra)
#> terra 1.7.78

r = rast(ncols=3, nrows=3)
values(r) = c(-2, -3, -1, -1.96, 0, 1, 2, 3, 1.96)

z = terra::classify(r, 
                    matrix(c(-1.96, 1.96, 2,
                              1.96, Inf, 1,
                             -Inf, -1.96, 3), byrow = T, nrow = 3),
                    include.lowest = TRUE, right = TRUE)
plot(z, col = c("white", "green", "red"))
coords = xyFromCell(r, 1:ncell(r))
original_vals = values(r)
for (i in 1:nrow(coords)) {
  text(coords[i, 1], coords[i, 2], labels=round(original_vals[i], 2), 
       cex=0.8, col="black")
}


z1 = terra::classify(r, 
                    matrix(c(-1.96, 1.96, 2,
                             1.96, Inf, 1,
                             -Inf, -1.96, 3), byrow = T, nrow = 3))
plot(z1, col = c("white", "green", "red"))
coords = xyFromCell(r, 1:ncell(r))
original_vals = values(r)
for (i in 1:nrow(coords)) {
  text(coords[i, 1], coords[i, 2], labels=round(original_vals[i], 2), 
       cex=0.8, col="black")
}


z2 = terra::classify(r, 
                     matrix(c(-1.96, 1.96, 2,
                              1.96, Inf, 1,
                              -Inf, -1.96, 3), byrow = T, nrow = 3),
                     include.lowest = FALSE, right = FALSE)
plot(z2, col = c("white", "green", "red"))
coords = xyFromCell(r, 1:ncell(r))
original_vals = values(r)
for (i in 1:nrow(coords)) {
  text(coords[i, 1], coords[i, 2], labels=round(original_vals[i], 2), 
       cex=0.8, col="black")
}


z3 = terra::classify(r, 
                     matrix(c(-1.96, 1.96, 2,
                              1.96, Inf, 1,
                              -Inf, -1.96, 3), byrow = T, nrow = 3),
                     include.lowest = FALSE, right = FALSE, others = 2)
plot(z3, col = c("white", "green", "red"))
coords = xyFromCell(r, 1:ncell(r))
original_vals = values(r)
for (i in 1:nrow(coords)) {
  text(coords[i, 1], coords[i, 2], labels=round(original_vals[i], 2), 
       cex=0.8, col="black")
}

Created on 2024-07-02 with reprex v2.1.0

As you can see, setting include.lowest and right to TRUE does not work. Also, when include.lowest and right are set to FALSE, the others parameter also not take effect.I don't know if I misunderstood

My terra package configuration information is as follows:

packageVersion("terra")
#> [1] '1.7.78'
terra::gdal(lib="all")
#>     gdal     proj     geos 
#>  "3.8.4"  "9.3.1" "3.12.1"

Created on 2024-07-02 with reprex v2.1.0

SpatLyu commented 2 months ago

I have tested the terra 1.7.80, the result is also wrong:

library(terra)
#> terra 1.7.80

r = rast(ncols=3, nrows=3)
values(r) = c(-2, -3, -1, -1.96, 0, 1, 2, 3, 1.96)

z = terra::classify(r, 
                    matrix(c(-1.96, 1.96, 2,
                             1.96, Inf, 1,
                             -Inf, -1.96, 3), byrow = T, nrow = 3),
                    include.lowest = TRUE, right = TRUE)
plot(z, col = c("white", "green", "red"))
coords = xyFromCell(r, 1:ncell(r))
original_vals = values(r)
for (i in 1:nrow(coords)) {
  text(coords[i, 1], coords[i, 2], labels=round(original_vals[i], 2), 
       cex=0.8, col="black")
}


z1 = terra::classify(r, 
                     matrix(c(-1.96, 1.96, 2,
                              1.96, Inf, 1,
                              -Inf, -1.96, 3), byrow = T, nrow = 3))
plot(z1, col = c("white", "green", "red"))
coords = xyFromCell(r, 1:ncell(r))
original_vals = values(r)
for (i in 1:nrow(coords)) {
  text(coords[i, 1], coords[i, 2], labels=round(original_vals[i], 2), 
       cex=0.8, col="black")
}


z2 = terra::classify(r, 
                     matrix(c(-1.96, 1.96, 2,
                              1.96, Inf, 1,
                              -Inf, -1.96, 3), byrow = T, nrow = 3),
                     include.lowest = FALSE, right = FALSE)
plot(z2, col = c("white", "green", "red"))
coords = xyFromCell(r, 1:ncell(r))
original_vals = values(r)
for (i in 1:nrow(coords)) {
  text(coords[i, 1], coords[i, 2], labels=round(original_vals[i], 2), 
       cex=0.8, col="black")
}


z3 = terra::classify(r, 
                     matrix(c(-1.96, 1.96, 2,
                              1.96, Inf, 1,
                              -Inf, -1.96, 3), byrow = T, nrow = 3),
                     include.lowest = FALSE, right = FALSE, others = 2)
plot(z3, col = c("white", "green", "red"))
coords = xyFromCell(r, 1:ncell(r))
original_vals = values(r)
for (i in 1:nrow(coords)) {
  text(coords[i, 1], coords[i, 2], labels=round(original_vals[i], 2), 
       cex=0.8, col="black")
}

Created on 2024-07-03 with reprex v2.1.0

rhijmans commented 2 months ago

I think you want right=NA in this case

library(terra)
r <- rast(ncols=3, nrows=3)
values(r) <- c(-3, -2, -1.96, -1, 0, 1, 1.96, 2, 3)
m <- matrix(c(-1.96, 1.96, 2,     1.96, Inf, 1,    -Inf, -1.96, 3),  byrow = TRUE, nrow = 3)

z <- terra::classify(r, m, right = NA, include.lowest=TRUE)
values(c(r, z))
#      lyr.1 lyr.1
# [1,] -3.00     3
# [2,] -2.00     3
# [3,] -1.96     2
# [4,] -1.00     2
# [5,]  0.00     2
# [6,]  1.00     2
# [7,]  1.96     2
# [8,]  2.00     1
# [9,]  3.00     1

But note that you are using 1.96 as if it were an integer. That may not be appropriate because

1.96 == 1.96000000000000000001
#[1] TRUE
1.96 == 1.9599999999999999999
# [1] TRUE
rhijmans commented 2 months ago

This would probably be a better approach

m2 <- matrix(c(-1.96, 1.96, 2,
          1.96, Inf, 1), byrow = TRUE, ncol = 3)
zx <- terra::classify(r, m2, right = T, include.lowest=TRUE, others=3)

You can also use ifel. That has a clear syntac, but it is much less efficient

zy <- ifel(r < -1.96, 3, ifel(r > 1.96, 1, 2))

values(c(r, z, zx, zy))
#      lyr.1 lyr.1 lyr.1 lyr.1
# [1,] -3.00     3     3     3
# [2,] -2.00     3     3     3
# [3,] -1.96     2     2     2
# [4,] -1.00     2     2     2
# [5,]  0.00     2     2     2
# [6,]  1.00     2     2     2
# [7,]  1.96     2     2     2
# [8,]  2.00     1     1     1
# [9,]  3.00     1     1     1