search

ropensci / slopes

Package to calculate slopes of roads, rivers and trajectories
https://docs.ropensci.org/slopes/
GNU General Public License v3.0
70 stars 6 forks source link

slope_raster : difference between raster and terra #14

Closed temospena closed 3 years ago

temospena commented 3 years ago

I noticed some weird results if using terrain this function, when compared with the default raster.

### Weired results with terra
library(slopes)

#import files from package
demraster = raster::raster("https://github.com/ITSLeeds/slopes/releases/download/0.0.0/dem_lisbon.tif")
class(demraster)
#> [1] "RasterLayer"
#> attr(,"package")
#> [1] "raster"
demraster
#> class      : RasterLayer 
#> dimensions : 133, 200, 26600  (nrow, ncol, ncell)
#> resolution : 10, 10  (x, y)
#> extent     : -88285, -86285, -106485, -105155  (xmin, xmax, ymin, ymax)
#> crs        : NA 
#> source     : https://github.com/ITSLeeds/slopes/releases/download/0.0.0/dem_lisbon.tif 
#> names      : dem_lisbon 
#> values     : 0, 97.906  (min, max)

demterra = terra::rast("https://github.com/ITSLeeds/slopes/releases/download/0.0.0/dem_lisbon.tif")
class(demterra)
#> [1] "SpatRaster"
#> attr(,"package")
#> [1] "terra"
demterra
#> class       : SpatRaster 
#> dimensions  : 133, 200, 1  (nrow, ncol, nlyr)
#> resolution  : 10, 10  (x, y)
#> extent      : -88285, -86285, -106485, -105155  (xmin, xmax, ymin, ymax)
#> coord. ref. :  
#> data source : dem_lisbon.tif 
#> names       : dem_lisbon 
#> min values  :          0 
#> max values  :     97.906

segments = lisbon_road_segments #from package
class(segments)
#> [1] "sf"         "tbl_df"     "tbl"        "data.frame"

#plot them to check
raster::plot(demraster)
plot(sf::st_geometry(segments), add = TRUE)

raster::plot(demterra)
plot(sf::st_geometry(segments), add = TRUE)


#process both
segments$sloperaster = slope_raster(segments, e = demraster)
#> [1] TRUE
segments$slopeterra = slope_raster(segments, e = demterra, terra = T)
#> Warning in cbind(m, z): number of rows of result is not a multiple of vector
#> length (arg 2)

#compare results
summary(segments$sloperaster)
#>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
#> 0.00000 0.01246 0.03534 0.05462 0.08251 0.27583
summary(segments$slopeterra)
#>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
#>  0.1007  0.1199  0.1294  0.1547  0.1540  2.1752

print(segments[1:20,c(24,25)])
Simple feature collection with 20 features and 2 fields
geometry type:  LINESTRING
dimension:      XY
bbox:           xmin: -88274.94 ymin: -106371.5 xmax: -86363.56 ymax: -105162.1
projected CRS:  ETRS89 / Portugal TM06
   sloperaster      slopeterra          
1  0.003074005  0.1134249 
2  0.010870459  0.1346804 
3  0.004998315  0.1081541
4  0.014255667  0.1293899 
5  0.005333236  0.1281508
6  0.013103029  0.1135203 
7  0.013103029  0.1135203 
8  0.025857474  0.1429152
9  0.039232485  0.1247636 
10 0.002259499  0.1234486 
11 0.003420071  0.1321447 
12 0.063387649  0.1179675 
13 0.014752744  0.1204235 
14 0.023007322  0.1247678
15 0.007431299  0.1374304 
16 0.035662001  0.1233757
17 0.013103029  0.1135203
18 0.054831718  0.1263556
19 0.089791758  0.1536468
20 0.140735333  0.1137517

Created on 2020-10-27 by the reprex package (v0.3.0)

Session info ``` r devtools::session_info() #> - Session info --------------------------------------------------------------- #> setting value #> version R version 3.6.2 (2019-12-12) #> os Windows 7 x64 SP 1 #> system x86_64, mingw32 #> ui RTerm #> language (EN) #> collate Portuguese_Portugal.1252 #> ctype Portuguese_Portugal.1252 #> tz Europe/London #> date 2020-10-27 #> #> - Packages ------------------------------------------------------------------- #> package * version date lib source #> assertthat 0.2.1 2019-03-21 [1] CRAN (R 3.5.3) #> backports 1.1.4 2019-04-10 [1] CRAN (R 3.5.3) #> callr 3.4.3 2020-03-28 [1] CRAN (R 3.6.3) #> class 7.3-15 2019-01-01 [2] CRAN (R 3.6.2) #> classInt 0.4-1 2019-08-06 [1] CRAN (R 3.5.3) #> cli 2.0.2 2020-02-28 [1] CRAN (R 3.6.3) #> codetools 0.2-16 2018-12-24 [2] CRAN (R 3.6.2) #> crayon 1.3.4 2017-09-16 [1] CRAN (R 3.5.1) #> curl 4.3 2019-12-02 [1] CRAN (R 3.6.3) #> DBI 1.0.0 2018-05-02 [1] CRAN (R 3.5.1) #> desc 1.2.0 2018-05-01 [1] CRAN (R 3.5.1) #> devtools 2.3.0 2020-04-10 [1] CRAN (R 3.6.3) #> digest 0.6.18 2018-10-10 [1] CRAN (R 3.5.3) #> e1071 1.7-1 2019-03-19 [1] CRAN (R 3.5.3) #> ellipsis 0.3.0 2019-09-20 [1] CRAN (R 3.5.3) #> evaluate 0.13 2019-02-12 [1] CRAN (R 3.5.3) #> fansi 0.4.0 2018-10-05 [1] CRAN (R 3.5.3) #> fs 1.4.1 2020-04-04 [1] CRAN (R 3.6.3) #> glue 1.3.1 2019-03-12 [1] CRAN (R 3.5.3) #> highr 0.8 2019-03-20 [1] CRAN (R 3.5.3) #> htmltools 0.4.0 2019-10-04 [1] CRAN (R 3.5.3) #> httr 1.4.1 2019-08-05 [1] CRAN (R 3.5.3) #> KernSmooth 2.23-16 2019-10-15 [2] CRAN (R 3.6.2) #> knitr 1.28 2020-02-06 [1] CRAN (R 3.5.3) #> lattice 0.20-38 2018-11-04 [2] CRAN (R 3.6.2) #> magrittr 1.5 2014-11-22 [1] CRAN (R 3.5.1) #> memoise 1.1.0 2017-04-21 [1] CRAN (R 3.5.1) #> mime 0.6 2018-10-05 [1] CRAN (R 3.5.2) #> pbapply 1.4-0 2019-02-05 [1] CRAN (R 3.5.2) #> pkgbuild 1.0.8 2020-05-07 [1] CRAN (R 3.6.3) #> pkgload 1.0.2 2018-10-29 [1] CRAN (R 3.5.1) #> prettyunits 1.0.2 2015-07-13 [1] CRAN (R 3.5.1) #> processx 3.4.4 2020-09-03 [1] CRAN (R 3.6.3) #> ps 1.3.0 2018-12-21 [1] CRAN (R 3.5.3) #> R6 2.4.0 2019-02-14 [1] CRAN (R 3.5.3) #> raster 3.3-13 2020-07-17 [1] CRAN (R 3.6.3) #> Rcpp 1.0.1 2019-03-17 [1] CRAN (R 3.5.3) #> remotes 2.1.1 2020-02-15 [1] CRAN (R 3.6.3) #> rgdal 1.4-3 2019-03-14 [1] CRAN (R 3.5.3) #> rlang 0.4.6 2020-05-02 [1] CRAN (R 3.6.3) #> rmarkdown 2.1 2020-01-20 [1] CRAN (R 3.5.3) #> rprojroot 1.3-2 2018-01-03 [1] CRAN (R 3.5.1) #> sessioninfo 1.1.1 2018-11-05 [1] CRAN (R 3.5.3) #> sf 0.9-3 2020-05-04 [1] CRAN (R 3.6.3) #> slopes * 0.0.0.9000 2020-05-11 [1] Github (itsleeds/slopes@ecc77ae) #> sp 1.4-1 2020-02-28 [1] CRAN (R 3.6.3) #> stringi 1.4.3 2019-03-12 [1] CRAN (R 3.5.3) #> stringr 1.4.0 2019-02-10 [1] CRAN (R 3.5.3) #> terra 0.8-6 2020-08-01 [1] CRAN (R 3.6.3) #> testthat 2.2.1 2019-07-25 [1] CRAN (R 3.5.3) #> units 0.6-2 2018-12-05 [1] CRAN (R 3.5.3) #> usethis 1.6.1 2020-04-29 [1] CRAN (R 3.6.3) #> withr 2.1.2 2018-03-15 [1] CRAN (R 3.5.1) #> xfun 0.6 2019-04-02 [1] CRAN (R 3.5.3) #> xml2 1.2.2 2019-08-09 [1] CRAN (R 3.5.3) #> yaml 2.2.1 2020-02-01 [1] CRAN (R 3.5.3) #> #> [1] D:/rosa/Documents/R/win-library/3.6 #> [2] C:/Program Files/R/R-3.6.2/library ```

Even for a large dataset, with 23k segments, terraretrieves this results:

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
0.09991 0.10951 0.12062 0.13564 0.14085 1.38011
Robinlovelace commented 3 years ago

That looks like a bug to me. Thanks for reporting, will check the code.

Robinlovelace commented 3 years ago 
library(slopes)
class(lisbon_road_segments)
#> [1] "sf"         "tbl_df"     "tbl"        "data.frame"
raster::plot(dem_lisbon_raster)
plot(sf::st_geometry(lisbon_road_segments), add = TRUE)
demterra = terra::rast(dem_lisbon_raster)
plot(sf::st_geometry(lisbon_road_segments), add = TRUE)

lisbon_road_segments$sloperaster = slope_raster(lisbon_road_segments, e = dem_lisbon_raster)
#> [1] TRUE
library(terra)
#> terra version 0.8.6 (beta-release)
plot(demterra)

lisbon_road_segments$slopeterra = slope_raster(lisbon_road_segments, e = demterra, terra = T)
#> Warning in cbind(m, z): number of rows of result is not a multiple of vector
#> length (arg 2)
summary(lisbon_road_segments$sloperaster)
#>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
#> 0.00000 0.01246 0.03534 0.05462 0.08251 0.27583
summary(lisbon_road_segments$slopeterra)
#>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
#>  0.1007  0.1199  0.1294  0.1547  0.1540  2.1752

Created on 2020-10-29 by the reprex package (v0.3.0)

Robinlovelace commented 3 years ago

This updated example shows what's going on: terra::extract() results in a matrix. I'm not sure when this was added but I think it was a change 'upstream' that caused the issue and we've only just noticed. Many thanks for identifying the bug.

library(slopes)
class(lisbon_road_segments)
#> [1] "sf"         "tbl_df"     "tbl"        "data.frame"
raster::plot(dem_lisbon_raster)
plot(sf::st_geometry(lisbon_road_segments), add = TRUE)
demterra = terra::rast(dem_lisbon_raster)
plot(sf::st_geometry(lisbon_road_segments), add = TRUE)

lisbon_road_segments$sloperaster = slope_raster(lisbon_road_segments, e = dem_lisbon_raster)
#> [1] TRUE
library(terra)
#> terra version 0.8.6 (beta-release)
plot(demterra)

lisbon_road_segments$slopeterra = slope_raster(lisbon_road_segments, e = demterra, terra = T)
#> Warning in cbind(m, z): number of rows of result is not a multiple of vector
#> length (arg 2)
summary(lisbon_road_segments$sloperaster)
#>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
#> 0.00000 0.01246 0.03534 0.05462 0.08251 0.27583
summary(lisbon_road_segments$slopeterra)
#>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
#>  0.1007  0.1199  0.1294  0.1547  0.1540  2.1752

# Looking at the code, this seems to be the culprit:

# res = as.numeric(terra::extract(e, m[, 1:2], method = method))
# vs
# res = as.numeric(raster::extract(e, m[, 1:2], method = method))
# lets try it...
method = "bilinear"
m = sf::st_coordinates(lisbon_road_segment)[1:2, ]
e = demterra
as.numeric(terra::extract(e, m[, 1:2], method = method))
#> [1]  1.00000  2.00000 92.31126 91.93055
as.numeric(terra::extract(dem_lisbon_raster, m[, 1:2], method = method))
#> [1] 92.31126 91.93055
res_terra = terra::extract(e, m[, 1:2], method = method)
class(res_terra)
#> [1] "matrix" "array"
res_terra
#>      ID       r1
#> [1,]  1 92.31126
#> [2,]  2 91.93055
as.numeric(res_terra[, "r1"])
#> [1] 92.31126 91.93055

Created on 2020-10-29 by the reprex package (v0.3.0)

Robinlovelace commented 3 years ago

Good news @temospena, this issue seems to have resulted in an important bug fix. Many thanks!

remotes::install_github("itsleeds/slopes")
#> Using github PAT from envvar GITHUB_PAT
#> Skipping install of 'slopes' from a github remote, the SHA1 (5c9120f4) has not changed since last install.
#>   Use `force = TRUE` to force installation
library(slopes)
class(lisbon_road_segments)
#> [1] "sf"         "tbl_df"     "tbl"        "data.frame"
raster::plot(dem_lisbon_raster)
plot(sf::st_geometry(lisbon_road_segments), add = TRUE)
demterra = terra::rast(dem_lisbon_raster)
plot(sf::st_geometry(lisbon_road_segments), add = TRUE)

lisbon_road_segments$sloperaster = slope_raster(lisbon_road_segments, e = dem_lisbon_raster)
#> [1] TRUE
library(terra)
#> terra version 0.8.6 (beta-release)
plot(demterra)

lisbon_road_segments$slopeterra = slope_raster(lisbon_road_segments, e = demterra, terra = T)
summary(lisbon_road_segments$sloperaster)
#>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
#> 0.00000 0.01246 0.03534 0.05462 0.08251 0.27583
summary(lisbon_road_segments$slopeterra)
#>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
#> 0.00000 0.01205 0.03534 0.05458 0.08251 0.27583

Created on 2020-10-29 by the reprex package (v0.3.0)