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

focal resets min position to preferred position 4 depending on where min value is in window, i.e. unexpected #1364

Closed chris-english closed 9 months ago

chris-english commented 9 months ago

Context: line mapping a descent, as water, by steepest path, over dem.

This is with ubuntu-22.04 R-4.2.1

library(terra)
# terra 1.7.62  
elev <- c(78, 72, 69, 71, 58, 49,
          74, 67, 56, 49, 46, 50,
          69, 53, 44, 37, 38, 48,
          64, 58, 55, 22, 31, 24,
          68, 61, 47, 21, 16, 19,
          74, 53, 34, 12, 11, 12) |> matrix(ncol = 6, byrow = TRUE)
dem <- rast(elev) 
w <- matrix(c(1,1,NA,1,1,NA,NA,NA,NA), 3, 3)
fdem <- focal(dem, w =w, fun = 'min') 
dem_foc_min <- cells(fdem)
dem_foc_min
# [1]  8  9 10 11 12 14 15 16 17 18 20 21 22 23 24 26 27 28 29 30 32 33 34 35 36

Not all of these values will be used as cells(focal( %%6 == 0)) wrap RHS to LHS, but useful starting point. Indexing the focal window in this example

min_win_cells <- matrix(c(seq(1, ncell(dem)-7,1), seq(2, ncell(dem)-6,1), 
    seq(7, ncell(dem)-1,1), seq(8, ncell(dem), 1)),  
    nrow = length(seq(8, ncell(dem),1)), ncol = 4)

# original dem vs focal dem for function min.
kind1234 = vector(mode = 'integer', length = 29L)
for(i in 1:nrow(min_win_cells)) {
    kind1234[i] <- which.min(unname(unlist(dem[min_win_cells[i, ]])))
 }
kind1234
# [1] 4 4 4 4 3 1 4 4 4 3 3 3 2 2 4 3 4 3 2 4 4 4 3 3 4 4 4 4 3
#min_position_mtx = matrix(rep(rep(NA, 36), 29), nrow = 36)

for (i in 1:length(kind1234)) {
    min_position_mtx[c(min_win_cells[i, ]), i] <- kind1234[i]
 }

dem_foc_min_mtx = matrix(c(values(dem), values(dem_foc_min)), nrow = 36)
val_position_mtx =matrix(c(dem_foc_min_mtx, min_position_mtx), nrow = 36)
which(kind1234 == 3)
#[1]  5 10 11 12 16 18 23 24 29
dem[min_win_cells[5, ]]
#  lyr.1
#1    58
#2    49
#3    46
#4    50
dem_foc_min[min_win_cells[5, ]]
#  focal_min
#1        58
#2        49
#3        46
#4        46
which(kind1234 == 2)
#[1] 13 14 19
dem[min_win_cells[13, ]]
#  lyr.1
#1    69
#2    53
#3    64
#4    58
dem_foc_min[min_win_cells[13, ]]
#  focal_min
#1        69
#2        53
#3        64
#4        53
## kind1234 == 4 are uneffected
rhijmans commented 9 months ago

I have edited your question for legibility, but I do not know what the issue is. Can you please explain? What do you expect to find, and what do you find instead? I am guessing that you do not think that the min value returned by focal is always correct; and that you show that by computing them "manually" using the adjacent cells. Can you point to a particular cell in fdem that is wrong?

Here is how I would compare focal with an adjacency approach, suggesting that the results are correct.

w <- matrix(c(1,1,NA,1,1,NA,NA,NA,NA), 3, 3)
fdem <- focal(dem, w =w, fun = 'min', na.rm=T) 
w[is.na(w)] <- 0
a <- adjacent(dem, 1:ncell(dem), w)
da <- setValues(dem, apply(a, 1, \(i) min(dem[i], na.rm=T)))
all(values(da) == values(fdem), na.rm=T) 
# [1] TRUE
chris-english commented 9 months ago

I see I wasn't properly attending to NA via na.rm=T (and neglected to explore the difference) . Sorry for the distraction and thanks for the above code.

chris-english commented 9 months ago

Well, there's something going on here, as using your code above: all(values(da) == values(fdem), na.rm = T) [1] TRUE but all(values(dem) == values(fdem), na.rm = T) # [1] FALSE which(matrix(c(values(dem), values(fdem)), nrow = 36, ncol = 2)[, 1] != matrix(c(values(dem), values(fdem)), nrow = 36, ncol = 2)[, 2]) [1] 4 12 17 18 20 21 23 25 26 30 31 36 Isn't dem the ground truth of its values?

rhijmans commented 9 months ago

Given this

dem <- rast(elev) 
w <- matrix(c(1,1,NA,1,1,NA,NA,NA,NA), 3, 3)
fdem <- focal(dem, w =w, fun = 'min') 

It would be very surprising indeed if this were TRUE

all(values(dem) == values(fdem), na.rm = T) 
#[1] FALSE

But we should expect this to be TRUE

all(values(dem) >= values(fdem), na.rm = T) 
#[1] TRUE
chris-english commented 9 months ago

Noted, as >= is the case. I misconstrued that focal could be used for indexing.