raster size limit?

[olenaboiko303]

Is there a friction raster size limit to process with gdistance?
When I run the tool on a 4801x4381 raster, I produce the expected output. When I run it on a friction surface with 9295x12929 size, the output is still being produced (I use HPC resources to have enough memory, I used 256Gb) but it's completely wrong (for example, the pixel which fell into one of the targets is expected to have 0 for travel time output, but I'm seeing 8574.740234 value).

[DemetriosG]

Same issue as olenaboiko303. I write to provide more detailed information:

I have been using the gdistance package for years for cost path analysis. Thank you for the great work. Recently, I was tasked to support the computation of accumulated cost using a cost surface (raster) that represents travel time on a road network. The road network cost surface r is

class : RasterLayer
dimensions : 23906, 34740, 830494440 (nrow, ncol, ncell)
resolution : 50, 50 (x, y)
extent : -2003550, -266550, 2995750, 4191050 (xmin, xmax, ymin, ymax)
crs : +proj=aea +lat_0=23 +lon_0=-96 +lat_1=29.5 +lat_2=45.5 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
source(s) : memory
name : max
min value : 1721
max value : 7456

The transition layer is computed as:
tr <- transition( r, function(x) { ( 2 * res(r)[1] ) / ( x[1] + x[2] ) }, directions=8 ) # res = cell resolution
tr <- geoCorrection( tr )
cs <- accCost( tr , facility.xy )

The tr objects seems to compute normally:

str(summary(tr))
Classes ‘sparseSummary’ and 'data.frame': 6030184 obs. of 3 variables:
$ i: int 13458 48199 82940 117681 152422 187163 221904 256645 291386 326127 ...
$ j: int 48199 82940 117681 152422 187163 221904 256645 291386 326127 360868 ...
$ x: num 0.0291 0.0291 0.0291 0.0291 0.0291 ... - attr(*, "header")= chr "830494440 x 830494440 sparse Matrix of class "dsCMatrix", with 6030184 entries"

The output raster though (cs object) only contains 0 values.
class : SpatRaster
dimensions : 23906, 34740, 1 (nrow, ncol, nlyr)
resolution : 50, 50 (x, y)
extent : -2003550, -266550, 2995750, 4191050 (xmin, xmax, ymin, ymax)
coord. ref. : +proj=aea +lat_0=23 +lon_0=-96 +lat_1=29.5 +lat_2=45.5 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
source : r_tmp_2023-02-28_014238_1720_03797.grd
name : layer
min value : 0
max value : 0

I next started clipping the input raster: At 600x600km centered on facility.xy I get:

str(summary(tr))
Classes ‘sparseSummary’ and 'data.frame': 6512490 obs. of 3 variables:
$ i: int 727 2425 2452 2453 2454 2580 2581 2582 2583 2584 ...
$ j: int 728 2426 2453 2454 2455 2581 2582 2583 2584 2585 ...
$ x: num 0.0291 0.0201 0.0156 0.0156 0.0156 ... attr(*, "header")= chr "137208000 x 137208000 sparse Matrix of class "dsCMatrix", with 6512490 entries"

Note that the sparse matrix in tr now has more value cells, 6512490. With the entire raster, the sparse matrix had 6030184 value cells.

class : SpatRaster
dimensions : 12000, 12000, 1 (nrow, ncol, nlyr)
resolution : 50, 50 (x, y)
extent : -1768250, -1168250, 995750, 1567450 (xmin, xmax, ymin, ymax)
coord. ref. : +proj=aea +lat_0=23 +lon_0=-96 +lat_1=29.5 +lat_2=45.5 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
source(s) : memory
name : layer
min value : 1.797693e+308
max value : Inf

When the area is clipped down to 500x500km centered on facility.xy I get the expected result:

str(summary(tr))
Classes ‘sparseSummary’ and 'data.frame': 6205430 obs. of 3 variables:
$ i: int 1209 1210 1465 1466 1543 1549 1550 1551 1635 1636 ...
$ j: int 1210 1211 1466 1467 1544 1550 1551 1552 1636 1637 ...
$ x: num 0.0156 0.0156 0.0156 0.0156 0.0156 ... attr(*, "header")= chr "100000000 x 100000000 sparse Matrix of class "dsCMatrix", with 6205430 entries"

cs
class : SpatRaster
dimensions : 10000, 10000, 1 (nrow, ncol, nlyr)
resolution : 50, 50 (x, y)
extent : -1718250, -1218250, 3017450, 3517450 (xmin, xmax, ymin, ymax)
coord. ref. : +proj=aea +lat_0=23 +lon_0=-96 +lat_1=29.5 +lat_2=45.5 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
source(s) : memory
name : layer
min value : 0.00
max value : 21922.53

The workstation has 250GB of RAM and about 180GB remain available at runtime.

Any resolution on this issue would be greatly appreciated.

Demetrios

sessionInfo()
R version 4.2.2 (2022-10-31 ucrt)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 19045)
Matrix products: default
locale:
[1] LC_COLLATE=English_United States.utf8 LC_CTYPE=English_United States.utf8 LC_MONETARY=English_United States.utf8 LC_NUMERIC=C
[5] LC_TIME=English_United States.utf8

attached base packages: [1] stats graphics grDevices utils datasets methods base other attached packages: [1] doSNOW_1.0.20 snow_0.4-4 iterators_1.0.14 foreach_1.5.2 gdistance_1.6 Matrix_1.5-1 igraph_1.4.1 raster_3.6-14 sp_1.6-0 terra_1.7-15

loaded via a namespace (and not attached): [1] Rcpp_1.0.10 lattice_0.20-45 codetools_0.2-18 grid_4.2.2 magrittr_2.0.3 rlang_1.0.6 cli_3.6.0 tools_4.2.2 parallel_4.2.2 compiler_4.2.2 [11] pkgconfig_2.0.3

[eev38]

I may be having a similar issue. I need to create a conductance layer with a large extent(-18, 142, -32, 44), but the costDistance function does not work . I get the following error:

`Error in validObject(.Object): invalid class “dsparseVector” object: 'i' slot is not strictly increasing after truncation towards zero
Traceback:

1. costDistance(Conductance, A, B)
2. costDistance(Conductance, A, B)
3. .cd(x, fromCoords, toCoords)
4. isSymmetric(y)
5. isSymmetric(y)
6. .local(object, ...)
7. as(object, "sparseVector")
8. asMethod(object)
9. new(cl, i = ii, length = n, x = from@x)
10. initialize(value, ...)
11. initialize(value, ...)
12. .local(.Object, ...)
13. callNextMethod()
14. .nextMethod(.Object = .Object, i = i, x = x, ... = ...)
15. validObject(.Object)
16. stop(msg, ": ", errors, domain = NA)
    Raster layer: class : RasterLayer
    dimensions : 18000, 43200, 777600000 (nrow, ncol, ncell)
    resolution : 0.008333334, 0.008333334 (x, y)
    extent : -180, 180, -60, 90.00001 (xmin, xmax, ymin, ymax)
    crs : +proj=longlat +datum=WGS84 +no_defs
    source : w001001.adf
    names : w001001
    values : -431, 8233 (min, max)`

Conductance layer:
class : TransitionLayer dimensions : 9120, 19200, 175104000 (nrow, ncol, ncell) resolution : 0.008333334, 0.008333334 (x, y) extent : -17.99999, 142, -32, 44.00001 (xmin, xmax, ymin, ymax) crs : +proj=longlat +datum=WGS84 +no_defs values : conductance matrix class: dgCMatrix

This conductance layer with a smaller extent does work as it should:

class : TransitionLayer dimensions : 9000, 9000, 8.1e+07 (nrow, ncol, ncell) resolution : 0.008333334, 0.008333334 (x, y) extent : -19.99999, 55.00001, -35, 40.00001 (xmin, xmax, ymin, ymax) crs : +proj=longlat +datum=WGS84 +no_defs values : conductance matrix class: dgCMatrix

My R session info:
`R version 4.2.2 (2022-10-31)
Platform: x86_64-conda-linux-gnu (64-bit)
Running under: Red Hat Enterprise Linux 8.8 (Ootpa)

Matrix products: default
BLAS/LAPACK: /storage/home/x/.conda/envs/TDist/lib/libopenblasp-r0.3.21.so

locale:
[1] LC_CTYPE=C.UTF-8 LC_NUMERIC=C LC_TIME=C
[4] LC_COLLATE=C LC_MONETARY=C LC_MESSAGES=C
[7] LC_PAPER=C LC_NAME=C LC_ADDRESS=C
[10] LC_TELEPHONE=C LC_MEASUREMENT=C LC_IDENTIFICATION=C

attached base packages:
[1] stats graphics grDevices utils datasets methods base

other attached packages:
[1] gdistance_1.6 Matrix_1.5-3 igraph_1.3.1 raster_3.6-13 sp_1.5-1

loaded via a namespace (and not attached):
[1] Rcpp_1.0.10 magrittr_2.0.3 uuid_1.1-0 lattice_0.20-45
[5] rlang_1.0.6 fastmap_1.1.0 fansi_1.0.4 tools_4.2.2
[9] rgdal_1.5-32 grid_4.2.2 utf8_1.2.3 terra_1.7-3
[13] cli_3.6.0 htmltools_0.5.4 digest_0.6.31 lifecycle_1.0.3
[17] crayon_1.5.2 IRdisplay_1.1 repr_1.1.6 base64enc_0.1-3
[21] vctrs_0.5.2 codetools_0.2-19 IRkernel_1.3.2 glue_1.6.2
[25] evaluate_0.20 pbdZMQ_0.3-9 compiler_4.2.2 pillar_1.8.1
[29] jsonlite_1.8.4 pkgconfig_2.0.3 `

and the R script used to create the conductance layer:

`library(gdistance)

alt <-raster('/filepath/')
e <- extent(-18, 142, -32, 44)
alt.c <- crop(alt, e)

w <- alt.c
values(w) <- 0
alt.nona <- cover(alt.c, w)

altDiff <- function(x){x[2] - x[1]}#
hd <- transition(alt.nona, altDiff, 8, symm=FALSE)
slope <- geoCorrection(hd)

adj <- adjacent(alt.nona, cells=1:ncell(alt.nona), pairs=TRUE, directions=8)
speed <- slope
speed[adj] <- 6 * exp(-3.5 * abs(slope[adj] + 0.05))
Conductance <- geoCorrection(speed)

s <- calc(alt.c, fun=function(x){ x[x != "NA"] <- 1; return(x)} )
w <- alt.c
values(w) <- 1.1
alt.land1.ocean1.1 <- cover(s,w)
hd.2 <- transition(alt.land1.ocean1.1, mean, 8, symm=FALSE)
con_ocean1.1x <- Conductance*hd.2
`

Has anyone figured out a solution for this? If not, is it possible to splice two separate layers? I would appreciate any insight.

[andrewmarx]

I don't expect it to help, but I would start by updating gdistance, Matrix, raster, and terra. They've all had significant updates and bug fixes since the versions you have, and maybe we'll get lucky.

[eev38]

I think that did the trick - thanks Andrew. I do get a warning message, but fortunately the output is the same as when I tested the function with a smaller conductance layer in the old environment. Here is my R session info if anyone is interested:

R version 4.3.1 (2023-06-16)
Platform: x86_64-conda-linux-gnu (64-bit)
Running under: Red Hat Enterprise Linux 8.8 (Ootpa)

Matrix products: default
BLAS/LAPACK: /.conda/envs/UTDist/lib/libopenblasp-r0.3.24.so;  LAPACK version 3.11.0

locale:
 [1] LC_CTYPE=C.UTF-8    LC_NUMERIC=C        LC_TIME=C          
 [4] LC_COLLATE=C        LC_MONETARY=C       LC_MESSAGES=C      
 [7] LC_PAPER=C          LC_NAME=C           LC_ADDRESS=C       
[10] LC_TELEPHONE=C      LC_MEASUREMENT=C    LC_IDENTIFICATION=C

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] gdistance_1.6.4 Matrix_1.6-1.1  igraph_1.5.1    raster_3.6-23  
 [5] sp_2.1-0        lubridate_1.9.3 forcats_1.0.0   stringr_1.5.0  
 [9] dplyr_1.1.3     purrr_1.0.2     readr_2.1.4     tidyr_1.3.0    
[13] tibble_3.2.1    ggplot2_3.4.3   tidyverse_2.0.0

loaded via a namespace (and not attached):
 [1] gtable_0.3.4     jsonlite_1.8.7   compiler_4.3.1   crayon_1.5.2    
 [5] Rcpp_1.0.11      tidyselect_1.2.0 IRdisplay_1.1    scales_1.2.1    
 [9] uuid_1.1-1       fastmap_1.1.1    IRkernel_1.3.2   lattice_0.21-9  
[13] R6_2.5.1         generics_0.1.3   munsell_0.5.0    pillar_1.9.0    
[17] tzdb_0.4.0       rlang_1.1.1      utf8_1.2.3       terra_1.7-46    
[21] stringi_1.7.12   repr_1.1.6       timechange_0.2.0 cli_3.6.1       
[25] withr_2.5.1      magrittr_2.0.3   digest_0.6.33    grid_4.3.1      
[29] base64enc_0.1-3  hms_1.1.3        pbdZMQ_0.3-10    lifecycle_1.0.3 
[33] vctrs_0.6.3      evaluate_0.22    glue_1.6.2       fansi_1.0.4     
[37] colorspace_2.1-0 tools_4.3.1      pkgconfig_2.0.3  htmltools_0.5.6 

Warning message:

Warning message in i.c == i.t:
“longer object length is not a multiple of shorter object length”

[eev38]

Unfortunately I'm still unable to use the costDistance function with the extent(-18, 142, -32, 44) that I require. The environment I created above ended up producing an error when it was originally just a warning:

Error in all.equal(..., check.attributes = FALSE) : all(i.c == i.t) is not TRUE Calls: costDistance ... isTRUE -> ae -> all.equal -> all.equal -> stopifnot In addition: Warning message: In i.c == i.t : longer object length is not a multiple of shorter object length

I created a new environment with the latest versions of r and g-distance (as well as a new conductance layer) and get the following error:

Error in h(simpleError(msg, call)) : error in evaluating the argument 'target' in selecting a method for function 'all.equal': sparseVector length cannot exceed 2^53 Calls: costDistance ... isTRUE -> ae -> all.equal -> .M2V -> .handleSimpleError -> h

Session info:
`R version 4.3.2 (2023-10-31)
Platform: x86_64-conda-linux-gnu (64-bit)
Running under: Red Hat Enterprise Linux 8.8 (Ootpa)

Matrix products: default
BLAS/LAPACK: /lib/libopenblasp-r0.3.25.so; LAPACK version 3.11.0

locale:
[1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
[3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
[5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
[7] LC_PAPER=en_US.UTF-8 LC_NAME=C
[9] LC_ADDRESS=C LC_TELEPHONE=C
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C

attached base packages:
[1] stats graphics grDevices utils datasets methods base

other attached packages:
[1] gdistance_1.6.4 Matrix_1.6-4 igraph_1.5.1 raster_3.6-26
[5] sp_2.1-2

loaded via a namespace (and not attached):
[1] compiler_4.3.2 magrittr_2.0.3 cli_3.6.2 Rcpp_1.0.11
[5] grid_4.3.2 rlang_1.1.2 pkgconfig_2.0.3 lattice_0.22-5
[9] terra_1.7-55 `

It appears to be an issue/limitation with the Matrix package and all.equal function. I haven't had any luck with this - any suggestions? Many thanks!

[andrewmarx]

Unfortunately, my current computer only has 32GB of ram, which isn't enough to actually run raster data with 100+ million cells to try and pinpoint exactly what's happening.

But here's a guess:

Given a 100 million cell raster, the transition matrix is basically going to be a 100M x 100M matrix. Rather than storing all the zeros in memory (which would make the matrix about 80 petabytes), we can use various schemes to just keep track of the location of non-zero values. A simple one for this explanation is just the row and column numbers. Those row and column numbers will never be larger than the original number of raster cells. In other words, a non-zero value at the very end of the matrix would be row 100M and column 100M.

The problem appears to be that all.equal() is converting the sparse matrix to a sparse vector for the comparison. There's nothing wrong with this approach, and it may even make the most sense practically, but it does mean that instead of representing the location of each non-zero with two index values, it's now going to do it with a single index value. In order to do that, a single index value will have to represent numbers up to the width x height of the matrix. So 100 million squared in this case, which is 1x10^16, or the index for a value in the very last spot of the matrix.

For some reason, the Matrix package is using the double data type (numeric in R on 64-bit systems) to represent this single index value for sparse vectors. The double/numeric data type only has 2^53 bits of precision (which is where this number in the error came from). 2^53 =9.007x10^15, which is the largest integer index a double/numeric data type can accurately represent before precision issues start to pop up. Unfortunately, that's slightly too small for trying to convert the 100Mx100M matrix to a sparse vector.

So basically, right around the 100M mark for number of raster cells, some things might start fall apart. The gdistance package wasn't really designed for such large rasters (the original author appeared to focus on readability and maintainability, which is totally reasonable). It might work for some things at this scale, but I can't guarantee anything.

Unfortunately, I don't know what alternatives might work. You could try looking at the costDist() function in the terra package. You could also look at GRASS GIS which has the r.cost tool, which I think is available in QGIS. I'm not sure if anyone has ever worked out or implemented a tiling approach for most of the metrics in gdistance (Omniscape sort of does for circuit flow)

[DemetriosG]

Thanks for the explanation Andrew. I had empirically determined that the largest working dimensionality was just shy of 95000 rows and columns. Given the 2^53 limitation, the exact number is 94906 x 94906. I completely understand that there was no provision for such large datasets at gdistance package inception time, but for many interesting applications where gdistance is rather useful, currently this is a considerable limitation. I am wondering how much work would it take to use Rcpp and package bit64 that offers a integer64 type to do the conversion from the sparce matrix to a vector of indices. If implemented, it will push the max integer that can be represented without loss of precision from 2^53 to 2^60-1. The latter translates to a square matrix of 1,7073,741,823 rows and columns. Many thanks for all your efforts. Demetrios From: Andrew Marx ***@***.***> Sent: Friday, December 15, 2023 7:08 PM To: AgrDataSci/gdistance ***@***.***> Cc: Gatziolis, Demetrios - FS, OR ***@***.***>; Manual ***@***.***> Subject: Re: [AgrDataSci/gdistance] raster size limit? (Issue #17) Unfortunately, my current computer only has 32GB of ram, which isn't enough to actually run raster data with 100+ million cells to try and pinpoint exactly what's happening. But here's a guess: Given a 100 million cell raster, the transition matrix is basically going to be a 100M x 100M matrix. Rather than storing all the zeros in memory (which would make the matrix about 80 petabytes), we can use various schemes to just keep track of the location of non-zero values. A simple one for this explanation is just the row and column numbers. Those row and column numbers will never be larger than the original number of raster cells. In other words, a non-zero value at the very end of the matrix would be row 100M and column 100M. The problem appears to be that all.equal() is converting the sparse matrix to a sparse vector for the comparison. There's nothing wrong with this approach, and it may even make the most sense practically, but it does mean that instead of representing the location of each non-zero with two index values, it's now going to do it with a single index value. In order to do that, a single index value will have to represent numbers up to the width x height of the matrix. So 100 million squared in this case, which is 1x10^16, or the index for a value in the very last spot of the matrix. For some reason, the Matrix package is using the double data type (numeric in R on 64-bit systems) to represent this single index value for sparse vectors. The double/numeric data type only has 2^53 bits of precision (which is where this number in the error came from). 2^53 =9.007x10^15, which is the largest integer index a double/numeric data type can accurately represent before precision issues start to pop up. Unfortunately, that's slightly too small for trying to convert the 100Mx100M matrix to a sparse vector. So basically, right around the 100M mark for number of raster cells, some things might start fall apart. The gdistance package wasn't really designed for such large rasters (the original author appeared to focus on readability and maintainability, which is totally reasonable). It might work for some things at this scale, but I can't guarantee anything. Unfortunately, I don't know what alternatives might work. You could try looking at the costDist() function in the terra package. You could also look at GRASS GIS which has the r.cost tool, which I think is available in QGIS. I'm not sure if anyone has ever worked out or implemented a tiling approach for most of the metrics in gdistance (Omniscape sort of does for circuit flow) - Reply to this email directly, view it on GitHub<#17 (comment)>, or unsubscribe<https://github.com/notifications/unsubscribe-auth/A6FSJWZTBIT3R2FSUWNB6QLYJUGCNAVCNFSM6AAAAAAVETE7T2VHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMYTQNJYGY4TKNBWGY>. You are receiving this because you are subscribed to this thread.Message ID: ***@***.******@***.***>> This electronic message contains information generated by the USDA solely for the intended recipients. Any unauthorized interception of this message or the use or disclosure of the information it contains may violate the law and subject the violator to civil or criminal penalties. If you believe you have received this message in error, please notify the sender and delete the email immediately.

[eev38]

Thank you for the detailed explanation and alternative suggestions Andrew. I understand that the package was not originally intended for such large raster layers. If there is any possibility of a working solution I have access to hpc resources and could attempt to help out on my end.

Currently, I am attempting to run the costDistance function with a smaller layer: extent (-20, 55, -40, 40) or (10, 35, -35, 30) but get 'Inf' as the output with no error message; I also have to set the crs for the raster layer as it is originally NA. I do know that the following extent (71, 135, -10, 41) has worked in the past for a colleague (gdistance version 1.2-2 or earlier), but when I replicate it with my current environment I get 'Inf' as the output. The function does work with the smallest I've tried (and has a set crs): (30, 55, 0, 15). I'm not quite sure if it's a bug or something wrong with my own environment. I would be grateful for any insight. Thank you again for your support.

R session Info:

R version 4.3.2 (2023-10-31)
Platform: x86_64-conda-linux-gnu (64-bit)
Running under: Red Hat Enterprise Linux 8.8 (Ootpa)
Matrix products: default

attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] gdistance_1.6.4 Matrix_1.6-4 igraph_1.5.1 raster_3.6-26
[5] sp_2.1-2
loaded via a namespace (and not attached):
[1] compiler_4.3.2 magrittr_2.0.3 cli_3.6.2 Rcpp_1.0.11
[5] codetools_0.2-19 grid_4.3.2 rlang_1.1.2 pkgconfig_2.0.3
[9] lattice_0.22-5 terra_1.7-55

[andrewmarx]

Currently, I am attempting to run the costDistance function with a smaller layer: extent (-20, 55, -40, 40) or (10, 35, -35, 30) but get 'Inf' as the output with no error message; I also have to set the crs for the raster layer as it is originally NA. I do know that the following extent (71, 135, -10, 41) has worked in the past for a colleague (gdistance version 1.2-2 or earlier), but when I replicate it with my current environment I get 'Inf' as the output. The function does work with the smallest I've tried (and has a set crs): (30, 55, 0, 15). I'm not quite sure if it's a bug or something wrong with my own environment. I would be grateful for any insight. Thank you again for your support.

If you're getting Inf's, it means that your to/from coords don't have a path connecting them in the transition matrix. Double check your input raster and transition function; if they are producing 0 for the connection between cells, then those cells are effectively not connected. Keep in mind NAs in the data.

If you can't find an issue with that, try looking at the values in your transition matrix before and after the geocorrection. If the geocorrection is changing reasonable values to unreasonable values (0 or effectively 0), then we need to determine if that is an issue related to things like parameters, crs, etc, or if it's an actual bug.

The costDistance() function is pretty simple, so I'm pretty sure there isn't a bug with it (unless something has changed with the igraph functions it uses).

Thank you for the detailed explanation and alternative suggestions Andrew. I understand that the package was not originally intended for such large raster layers. If there is any possibility of a working solution I have access to hpc resources and could attempt to help out on my end.

My main priority is just keeping the package functional and available on CRAN. In terms of enhancements, I might occasionally rewrite some calculations to improve performance if it isn't too involved. But if it requires significant changes, especially in package dependencies, then it probably won't happen. Partly because I want to keep the package easy to maintain and install, and partly because I'm not funded in any way to improve this package, so I can't spend a lot of time on it. I have another related connectivity package that I'm significantly more invested in, prioritizes performance, and I'm considering implementing similar equivalent functionality there. So If I am able to optimize these functions, they'll show up there.

[eev38]

Hi Andrew - thank you for getting back to me and apologies for the delayed response. I ended up using an older version of R (4.2.1) and gdistance (1.6.2) that seems to have worked. I created a raster layer with the following dimensions: 9600, 9000, 86400000 (nrow, ncol, ncell) and did not need to set the crs.

I couldn't resolve the issue with the Inf output. I didn't get the chance to test the geocorrection function but I did make two versions of the same raster layer (one with gdistance v. 1.6.4 and the other with 1.6.2). When I compared the transition matrix values, the older version had the expected values present while the new version had 0s and no crs set. I am happy to use the older version but can troubleshoot the geocorrection function if needed. Thank you for your support and I'm looking forward to seeing how the samc function progresses.

R session info below:

R version 4.2.1 (2022-06-23)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Red Hat Enterprise Linux 8.9 (Ootpa)

Matrix products: default

attached base packages:
[1] stats graphics grDevices utils datasets methods base

other attached packages:
[1] gdistance_1.6.2 Matrix_1.5-4 igraph_1.4.2 raster_3.6-20
[5] sp_1.6-0

loaded via a namespace (and not attached):
[1] compiler_4.2.1 magrittr_2.0.3 cli_3.6.0 Rcpp_1.0.10
[5] codetools_0.2-18 grid_4.2.1 rlang_1.1.1 pkgconfig_2.0.3
[9] lattice_0.20-45 terra_1.7-29

Edit: working version is 1.6.2 on R version 4.2.1

[andrewmarx]

I wouldn't spend your time troubleshooting the function, but if you could try a few different intermediate versions to see when the behavior changed, then I can take a look at what changed

[eev38]

I wouldn't spend your time troubleshooting the function, but if you could try a few different intermediate versions to see when the behavior changed, then I can take a look at what changed

I made a correction to my post above; I was using gdistance version 1.6.2 and R version 4.2.1 that was available on the server I'm working on. Not quite sure what it means for the package in that case.

[andrewmarx]

I wouldn't spend your time troubleshooting the function, but if you could try a few different intermediate versions to see when the behavior changed, then I can take a look at what changed

I made a correction to my post above; I was using gdistance version 1.6.2 and R version 4.2.1 that was available on the server I'm working on. Not quite sure what it means for the package in that case.

I added a significantly optimized version of transition() in 1.6.4 (for non-global rasters). It should produce the same results, but it may be something related to larger raster sizes that the unit tests aren't capturing. Instead of using transition() you could use gdistance:::.TfromR_new() and gdistance:::.TfromR_old() and see if they produce the same results

[eev38]

I added a significantly optimized version of transition() in 1.6.4 (for non-global rasters). It should produce the same results, but it may be something related to larger raster sizes that the unit tests aren't capturing. Instead of using transition() you could use gdistance:::.TfromR_new() and gdistance:::.TfromR_old() and see if they produce the same results

I tried making a new raster with both functions and still get the Inf output with both; the only difference was that I had to change spelling of symm to sym when using gdistance:::.TfromR_new().