Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

inverting values of rast object instead of subtraction #1005

Closed
atsyplenkov opened this issue Feb 8, 2023 · 2 comments
Closed

inverting values of rast object instead of subtraction #1005

atsyplenkov opened this issue Feb 8, 2023 · 2 comments

Comments

@atsyplenkov
Copy link

I have faced up with strange behaviour of the terra package. I was trying to scale DEM values to the maximum height. While the raster package gives the expected results, the terra package gives me an inverted one. The problem seems to be associated with a subtraction syntax. I solved it by multiplying the rast object by (-1). Please take a look at the example below. Is it a bug or a feature? :-)

library(raster) # CRAN v3.6-14
library(terra)  # CRAN v1.6-53

# load datasets -----------------------------------------------------------
dem <- 
  system.file("ex/elev.tif", package="terra")

dem_terra <- 
  rast(dem)

dem_raster <- 
  raster(dem)

# helper functions --------------------------------------------------------
# function to compare two rast objects
compare_rasts <- 
  function(.rast1, .rast2, .rast3){
    
    rbind(
      .get_global_stats(.rast1),
      .get_global_stats(.rast2),
      .get_global_stats(.rast3)
    )
    
  }

.get_global_stats <- 
  function(.rast){
    
    tibble::tibble(
      rast = names(.rast),
      res = res(.rast)[1],
      mean = global(.rast, "mean", na.rm = T)[1,1],
      med = global(.rast, function(x) median(x, na.rm = T))[1,1],
      sd = global(.rast, function(x) sd(x, na.rm = T))[1,1],
      max = global(.rast, function(x) max(x, na.rm = T))[1,1],
      min = global(.rast, function(x) min(x, na.rm = T))[1,1],
      na = global(.rast, function(x) sum(is.na(x)))[1,1]
    )
    
  }

# function to get maximum values
.max_terra <- 
  function(x){
    
    global(x, function(i) max(i, na.rm = T))[1,1]
    
  }

# scale process -----------------------------------------------------------
# raster approach
norm_raster <- 1 - (dem_raster / (maxValue(dem_raster)))
names(norm_raster) <- "raster"
norm_raster <- rast(norm_raster)

# terra straightforward approach
norm_terra_inverted <- 1 - (dem_terra / (.max_terra(dem_terra)))
names(norm_terra_inverted) <- "terra inverted"

# terra solution...
norm_terra <- 1 + (-1) * (dem_terra / (.max_terra(dem_terra)))
names(norm_terra) <- "terra"

# compare -----------------------------------------------------------------
compare_rasts(norm_terra_inverted, norm_raster, norm_terra)

# A tibble: 3 × 8
  rast               res   mean    med    sd   max    min    na
  <chr>            <dbl>  <dbl>  <dbl> <dbl> <dbl>  <dbl> <int>
1 terra inverted 0.00833 -0.363 -0.391 0.147 0     -0.742  3942
2 raster         0.00833  0.363  0.391 0.147 0.742  0      3942
3 terra          0.00833  0.363  0.391 0.147 0.742  0      3942
@brownag
Copy link
Contributor

brownag commented Feb 9, 2023

I think this is a bug and was fixed in #978

@rhijmans
Copy link
Member

rhijmans commented Feb 9, 2023

With the current CRAN version I see no differences.

library(raster)
#Loading required package: sp
library(terra)
#terra 1.7.3

f <- system.file("ex/elev.tif", package="terra")
dem_terra <- rast(f)
dem_raster <- raster(f)

global_stats <- function(x){
    data.frame(
      rast = names(x),
      res = res(x)[1],
      mean = global(x, "mean", na.rm = T)[1,1],
      med = global(x, function(x) median(x, na.rm = T))[1,1],
      sd = global(x, function(x) sd(x, na.rm = T))[1,1],
      max = global(x, function(x) max(x, na.rm = T))[1,1],
      min = global(x, function(x) min(x, na.rm = T))[1,1],
      na = global(x, function(x) sum(is.na(x)))[1,1]
    )
}

compare_rasts <- function(x, y, z){   
    rbind(
      global_stats(x),
      global_stats(y),
      global_stats(z)
    )
}

norm_raster <- 1 - (dem_raster / (maxValue(dem_raster)))
names(norm_raster) <- "raster"
norm_raster <- rast(norm_raster)

norm_terra_inverted <- 1 - (dem_terra / ( minmax(dem_terra)[2]))
names(norm_terra_inverted) <- "terra inverted"

norm_terra <- 1 + (-1) * (dem_terra / (minmax(dem_terra)[2]))
names(norm_terra) <- "terra"

compare_rasts(norm_terra_inverted, norm_raster, norm_terra)

#            rast         res      mean       med        sd       max min   na
#1 terra inverted 0.008333333 0.3631872 0.3912249 0.1466524 0.7422303   0 3942
#2         raster 0.008333333 0.3631872 0.3912249 0.1466524 0.7422303   0 3942
#3          terra 0.008333333 0.3631872 0.3912249 0.1466524 0.7422303   0 3942

The version you are using was the CRAN version for about a week; unfortunately you got that one.

@rhijmans rhijmans closed this as completed Feb 9, 2023
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

3 participants