Closed micha-silver closed 2 years ago
@micha-silver, that is odd. So it produced a full surface on the 128 GB RAM machine but not on the 16 GB RAM machine? That does make it seem like some type of memory issue but I'm not entirely sure what is happening. What are the dimensions of your raster, how many grey levels are you using, what window size, and are you using na.rm? Once you send me that, I'll try to see if I can reproduce the issue on my end as my computer also has 16 GB of RAM though it may take a bit to run.
Thanks in advance for helping out.
Here's the code we are using:
class(stk_full)
[1] "SpatRaster"
attr(,"package")
[1] "terra"
grn <- stk_full$green
grn_quantized <- quantize_raster(r = grn, n_levels = 16,
method = "equal prob")
textures <- glcm_textures(r=grn_quantized,
w=c(5, 5),
shift=list(c(1,0), c(1,1), c(0, 1), c(-1,1)),
metrics = c("glcm_variance",
"glcm_contrast",
"glcm_dissimilarity",
"glcm_correlation",
"glcm_ASM"),
quantization = "none",
n_levels=16, na.rm = TRUE)
The raster dimensions:
> dim(stk_full)
[1] 16364 12509 13
and my session info:
> sessionInfo()
R version 4.1.2 (2021-11-01)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Debian GNU/Linux 11 (bullseye)
Matrix products: default
BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.13.so
locale:
[1] LC_CTYPE=en_GB.UTF-8 LC_NUMERIC=C
[3] LC_TIME=en_GB.UTF-8 LC_COLLATE=en_GB.UTF-8
[5] LC_MONETARY=en_GB.UTF-8 LC_MESSAGES=en_GB.UTF-8
[7] LC_PAPER=en_GB.UTF-8 LC_NAME=C
[9] LC_ADDRESS=C LC_TELEPHONE=C
[11] LC_MEASUREMENT=en_GB.UTF-8 LC_IDENTIFICATION=C
attached base packages:
[1] parallel stats graphics grDevices utils datasets methods
[8] base
other attached packages:
[1] terra_1.5-9 fusionImage_0.0.1 raster_3.5-12
[4] sp_1.4-6 GLCMTextures_0.3 supercells_0.6.3
[7] corrplot_0.92 doParallel_1.0.16 iterators_1.0.13
[10] foreach_1.5.1 ranger_0.13.1 caret_6.0-90
[13] lattice_0.20-45 remotes_2.4.2 exactextractr_0.7.2
[16] RStoolbox_0.2.6 sf_1.0-5 forcats_0.5.1
[19] stringr_1.4.0 dplyr_1.0.7 purrr_0.3.4
[22] readr_2.1.1 tidyr_1.1.4 tibble_3.1.6
[25] ggplot2_3.3.5 tidyverse_1.3.1
loaded via a namespace (and not attached):
[1] colorspace_2.0-2 ellipsis_0.3.2 class_7.3-19
[4] rprojroot_2.0.2 rgdal_1.5-28 fs_1.5.2
[7] rstudioapi_0.13 proxy_0.4-26 listenv_0.8.0
[10] prodlim_2019.11.13 fansi_0.5.0 lubridate_1.8.0
[13] xml2_1.3.3 codetools_0.2-18 splines_4.1.2
[16] R.methodsS3_1.8.1 jsonlite_1.7.2 pROC_1.18.0
[19] broom_0.7.10 dbplyr_2.1.1 R.oo_1.24.0
[22] rgeos_0.5-9 compiler_4.1.2 httr_1.4.2
[25] backports_1.4.1 assertthat_0.2.1 Matrix_1.3-4
[28] cli_3.1.0 prettyunits_1.1.1 tools_4.1.2
[31] gtable_0.3.0 glue_1.6.0 reshape2_1.4.4
[34] Rcpp_1.0.7 cellranger_1.1.0 vctrs_0.3.8
[37] nlme_3.1-153 timeDate_3043.102 gower_0.2.2
[40] ps_1.6.0 globals_0.14.0 rvest_1.0.2
[43] lifecycle_1.0.1 XML_3.99-0.8 future_1.23.0
[46] MASS_7.3-54 scales_1.1.1 ipred_0.9-12
[49] hms_1.1.1 curl_4.3.2 geosphere_1.5-14
[52] rpart_4.1-15 stringi_1.7.6 e1071_1.7-9
[55] pkgbuild_1.3.0 lava_1.6.10 rlang_0.4.12
[58] pkgconfig_2.0.3 recipes_0.1.17 processx_3.5.2
[61] tidyselect_1.1.1 parallelly_1.29.0 plyr_1.8.6
[64] magrittr_2.0.1 R6_2.5.1 generics_0.1.1
[67] DBI_1.1.2 pillar_1.6.4 haven_2.4.3
[70] withr_2.4.3 units_0.7-2 survival_3.2-13
[73] nnet_7.3-16 future.apply_1.8.1 modelr_0.1.8
[76] crayon_1.4.2 KernSmooth_2.23-20 utf8_1.2.2
[79] tzdb_0.2.0 grid_4.1.2 readxl_1.3.1
[82] data.table_1.14.2 callr_3.7.0 ModelMetrics_1.2.2.2
[85] reprex_2.0.1 digest_0.6.29 classInt_0.4-3
[88] R.utils_2.11.0 stats4_4.1.2 munsell_0.5.0
>
If you want to try the actual raster we are working with (be forewarned: 5 GB!), can be download from google drive
Regards, Micha
@micha-silver, awesome thanks. I'll run it through a couple different functions to see if it's something with the glcm_textures
function, terra::focalCpp
or some type of interaction between the 2.
@micha-silver, sorry I was not able to reproduce the issue on my computer which runs Windows 10 and has 16GB of RAM. Could you use debug
(as in the code I showed below) to see where the issue pops up. I would expect the issue to pop up in the loop defining out_list
where each element of the list is a multilayer SpatRaster
containing the texture measures for a single shift, or the loop where out_layer
is defined which averages the textures across shifts. So either this line out_list[[k]]<- terra::focalCpp(r, w=w, fun = C_glcm_textures_helper, w2=w, n_levels= n_levels, shift = shift[[k]], na_rm=na.rm)
or this line out_layer<- mean(do.call(c, lapply(out_list, terra::subset,j))) #Average across all shifts
. In debug mode you can run it line by line and then plot the intermediate surfaces to see where the clipping is happening.
library(terra) #‘1.5.9’
library(GLCMTextures) #‘0.3’
stk_full<- rast("Data/namibia_full.tif")
grn <- stk_full$green
grn_quantized <- quantize_raster(r = grn, n_levels = 16,
method = "equal prob")
debug(glcm_textures)
textures <- glcm_textures(r=grn_quantized,
w=c(5, 5),
shift=list(c(1,0), c(1,1), c(0, 1), c(-1,1)),
metrics = c("glcm_variance",
"glcm_contrast",
"glcm_dissimilarity",
"glcm_correlation",
"glcm_ASM"),
quantization = "none",
n_levels=16, na.rm = TRUE)
plot(textures)
sessionInfo()
# R version 4.0.4 (2021-02-15)
# Platform: x86_64-w64-mingw32/x64 (64-bit)
# Running under: Windows 10 x64 (build 19042)
#
# Matrix products: default
#
# locale:
# [1] LC_COLLATE=English_United States.1252 LC_CTYPE=English_United States.1252 LC_MONETARY=English_United States.1252 LC_NUMERIC=C
# [5] LC_TIME=English_United States.1252
#
# attached base packages:
# [1] stats graphics grDevices utils datasets methods base
#
# other attached packages:
# [1] GLCMTextures_0.3 terra_1.5-9
#
# loaded via a namespace (and not attached):
# [1] compiler_4.0.4 tools_4.0.4 sp_1.4-6 Rcpp_1.0.7 raster_3.5-12 codetools_0.2-18 grid_4.0.4 lattice_0.20-41
@ailich
I tried as you suggested, with debug(), and plotted the first two rasters from out_list[[k]]. Same problem ...
and here is the single band quantized raster that the glcm_textures() is based on:
I see that you are using R 4.0.4, whereas I have gone to 4.1.2. Can that have anything to do with this? (On the lab server, where this clipping did NOT happen, I have R 4.0.4) Let me know if there are additional tests I can try...
Cheers, Micha
Hi @micha-silver. After updating to R 4.1.2, I still am not running into this issue, even when lowering the memory fraction using terraOptions(memfrac=0.1)
to help ensure that the data can't fit fully in memory. The striping in your first set of plots and the separate chunk looks a lot like a tiling issue (see rspatial/terra#312 and rspatial/terra#452) in terra::focalCpp
though it's unclear to me why it's only occurring on one of the three machines. I may need to contact the author of the terra
package to see if he has any insight into what is happening. Maybe it's related to the operating system (Windows vs Linux)?
You can also try running a smaller data set that easily fits in memory, but use debug mode and in the line out_list[[k]]<- terra::focalCpp(r, w=w, fun = C_glcm_textures_helper, w2=w, n_levels= n_levels, shift = shift[[k]], na_rm=na.rm)
add wopt=list(steps=5)
so it looks like out_list[[k]]<- terra::focalCpp(r, w=w, fun = C_glcm_textures_helper, w2=w, n_levels= n_levels, shift = shift[[k]], na_rm=na.rm, wopt=list(steps=5))
, which should force it to tile even a small raster in 5 chunks so you can see if there's still a tiling issue. Additionally, you can try running Qfit
from my other R package on the large data as that function also relies on terra::focalCpp
and outputs many layers so it'd be helpful to see if the issue pops up there as well or if something specifically in C_glcm_textures_helper
causes memory issues with terra::focalCpp
.
Results from running just the first shift in debug mode:
sessionInfo()
R version 4.1.2 (2021-11-01)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 19042)
Matrix products: default
locale:
[1] LC_COLLATE=English_United States.1252 LC_CTYPE=English_United States.1252 LC_MONETARY=English_United States.1252
[4] LC_NUMERIC=C LC_TIME=English_United States.1252
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] GLCMTextures_0.3.1 terra_1.5-9
loaded via a namespace (and not attached):
[1] compiler_4.1.2 tools_4.1.2 sp_1.4-6 Rcpp_1.0.7 raster_3.5-11 codetools_0.2-18 grid_4.1.2
[8] lattice_0.20-45
@micha-silver, I have been able to reproduce this on an Ubuntu VM with ~10GB of RAM. My best guess as of now is that this is an issue in terra:focalCpp
when data cannot fit in memory on Linux systems as the same thing occurs in my Qfit
function in my other R package when using the Linux VM. I will report the issue there and hopefully it can be resolved soon.
library(terra)
library(GLCMTextures)
stk_full<- rast("Desktop/Test_data/namibia_full.tif")
grn <- stk_full$green
grn_quantized <- quantize_raster(r = grn, n_levels = 16,
method = "equal prob")
debug(glcm_textures)
textures <- glcm_textures(r=grn_quantized,
w=c(5, 5),
shift=list(c(1,0), c(1,1), c(0, 1), c(-1,1)),
metrics = c("glcm_variance",
"glcm_contrast",
"glcm_dissimilarity",
"glcm_correlation",
"glcm_ASM"),
quantization = "none",
n_levels=16, na.rm = TRUE)
# From Debug Mode
Browse[2]> plot(r)
Browse[2]> plot(out_list[[1]], col=grey.colors(100), add=TRUE)
rm(list=ls())
# Restarting R session...
library(terra)
library(MultiscaleDEM)
library(GLCMTextures)
stk_full<- rast("Desktop/Test_data/namibia_full.tif")
grn <- stk_full$green
grn_quantized <- quantize_raster(r = grn, n_levels = 16,
method = "equal prob")
out<- Qfit(grn_quantized, w=c(5,5), metrics = c(), na.rm=TRUE, return_params = TRUE)
plot(out)
plot(out$a, add=TRUE, col = rainbow(100))
sessionInfo()
# R version 4.1.2 (2021-11-01)
# Platform: x86_64-pc-linux-gnu (64-bit)
# Running under: Ubuntu 21.10
#
# Matrix products: default
# BLAS: /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.9.0
# LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.9.0
#
# locale:
# [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8 LC_MONETARY=en_US.UTF-8
# [6] LC_MESSAGES=en_US.UTF-8 LC_PAPER=en_US.UTF-8 LC_NAME=C 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
#
# loaded via a namespace (and not attached):
# [1] compiler_4.1.2 tools_4.1.2
This does not smell like a memory problem. It looks like a datatype/size error. I wonder if it can be reproduced with a small example like this:
library(GLCMTextures)
library(terra)
sessionInfo()
gdal()
r <- rast(volcano, extent= ext(2667400, 2667400 + ncol(volcano)*10, 6478700, 6478700 + nrow(volcano)*10), crs = "EPSG:27200")
rqt <- quantize_raster(r, n_levels = 16, method = "equal prob")
tx1 <- glcm_textures(rqt, w=c(5, 5), shift=list(c(1,0), c(1,1), c(0, 1), c(-1,1)),
metrics = c("glcm_variance", "glcm_contrast", "glcm_dissimilarity", "glcm_correlation", "glcm_ASM"), quantization="none", n_levels=16, na.rm=TRUE)
plot(tx1)
opt <- list(steps=10)
tx2 <- glcm_textures(rqt, w=c(5, 5), shift=list(c(1,0), c(1,1), c(0, 1), c(-1,1)),
metrics = c("glcm_variance", "glcm_contrast", "glcm_dissimilarity", "glcm_correlation", "glcm_ASM"), quantization = "none", n_levels=16, na.rm=TRUE, wopt=opt)
plot(tx2)
terraOptions(todisk=TRUE)
tx3 <- glcm_textures(rqt, w=c(5, 5), shift=list(c(1,0), c(1,1), c(0, 1), c(-1,1)),
metrics = c("glcm_variance", "glcm_contrast", "glcm_dissimilarity", "glcm_correlation", "glcm_ASM"), quantization = "none", n_levels=16, na.rm=TRUE, wopt=opt)
plot(tx3)
describe(sources(tx3))
@rhijmans, for the small example I'm not running into the issue. Only very small differences in values for the version written to disk.
library(GLCMTextures)
library(terra)
sessionInfo()
# R version 4.1.2 (2021-11-01)
# Platform: x86_64-pc-linux-gnu (64-bit)
# Running under: Ubuntu 21.10
#
# Matrix products: default
# BLAS: /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.9.0
# LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.9.0
#
# locale:
# [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8 LC_MONETARY=en_US.UTF-8
# [6] LC_MESSAGES=en_US.UTF-8 LC_PAPER=en_US.UTF-8 LC_NAME=C 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] terra_1.5-12 GLCMTextures_0.3.2
#
# loaded via a namespace (and not attached):
# [1] compiler_4.1.2 tools_4.1.2 sp_1.4-6 Rcpp_1.0.7 raster_3.5-12 codetools_0.2-18 grid_4.1.2 lattice_0.20-45
gdal()
# "3.2.2"
r <- rast(volcano, extent= ext(2667400, 2667400 + ncol(volcano)*10, 6478700, 6478700 + nrow(volcano)*10), crs = "EPSG:27200")
rqt <- quantize_raster(r, n_levels = 16, method = "equal prob")
tx1 <- glcm_textures(rqt, w=c(5, 5), shift=list(c(1,0), c(1,1), c(0, 1), c(-1,1)),
metrics = c("glcm_variance", "glcm_contrast", "glcm_dissimilarity", "glcm_correlation", "glcm_ASM"), quantization="none", n_levels=16, na.rm=TRUE)
plot(tx1)
opt <- list(steps=10)
tx2 <- glcm_textures(rqt, w=c(5, 5), shift=list(c(1,0), c(1,1), c(0, 1), c(-1,1)),
metrics = c("glcm_variance", "glcm_contrast", "glcm_dissimilarity", "glcm_correlation", "glcm_ASM"), quantization = "none", n_levels=16, na.rm=TRUE, wopt=opt)
plot(tx2)
terraOptions(todisk=TRUE)
tx3 <- glcm_textures(rqt, w=c(5, 5), shift=list(c(1,0), c(1,1), c(0, 1), c(-1,1)),
metrics = c("glcm_variance", "glcm_contrast", "glcm_dissimilarity", "glcm_correlation", "glcm_ASM"), quantization = "none", n_levels=16, na.rm=TRUE, wopt=opt)
plot(tx3)
describe(sources(tx3))
# [1] "Driver: GTiff/GeoTIFF"
# [2] "Files: /tmp/RtmpwrVqQZ/spat_W4sxz04CVkiAmHf_29719.tif"
# [3] "Size is 61, 87"
# [4] "Coordinate System is:"
# [5] "PROJCRS[\"NZGD49 / New Zealand Map Grid\","
# [6] " BASEGEOGCRS[\"NZGD49\","
# [7] " DATUM[\"New Zealand Geodetic Datum 1949\","
# [8] " ELLIPSOID[\"International 1924\",6378388,297,"
# [9] " LENGTHUNIT[\"metre\",1]]],"
# [10] " PRIMEM[\"Greenwich\",0,"
# [11] " ANGLEUNIT[\"degree\",0.0174532925199433]],"
# [12] " ID[\"EPSG\",4272]],"
# [13] " CONVERSION[\"New Zealand Map Grid\","
# [14] " METHOD[\"New Zealand Map Grid\","
# [15] " ID[\"EPSG\",9811]],"
# [16] " PARAMETER[\"Latitude of natural origin\",-41,"
# [17] " ANGLEUNIT[\"degree\",0.0174532925199433],"
# [18] " ID[\"EPSG\",8801]],"
# [19] " PARAMETER[\"Longitude of natural origin\",173,"
# [20] " ANGLEUNIT[\"degree\",0.0174532925199433],"
# [21] " ID[\"EPSG\",8802]],"
# [22] " PARAMETER[\"False easting\",2510000,"
# [23] " LENGTHUNIT[\"metre\",1],"
# [24] " ID[\"EPSG\",8806]],"
# [25] " PARAMETER[\"False northing\",6023150,"
# [26] " LENGTHUNIT[\"metre\",1],"
# [27] " ID[\"EPSG\",8807]]],"
# [28] " CS[Cartesian,2],"
# [29] " AXIS[\"(E)\",east,"
# [30] " ORDER[1],"
# [31] " LENGTHUNIT[\"metre\",1]],"
# [32] " AXIS[\"(N)\",north,"
# [33] " ORDER[2],"
# [34] " LENGTHUNIT[\"metre\",1]],"
# [35] " USAGE["
# [36] " SCOPE[\"Engineering survey, topographic mapping.\"],"
# [37] " AREA[\"New Zealand - North Island, South Island, Stewart Island - onshore.\"],"
# [38] " BBOX[-47.33,166.37,-34.1,178.63]],"
# [39] " ID[\"EPSG\",27200]]"
# [40] "Data axis to CRS axis mapping: 1,2"
# [41] "Origin = (2667400.000000000000000,6479570.000000000000000)"
# [42] "Pixel Size = (10.000000000000000,-10.000000000000000)"
# [43] "Metadata:"
# [44] " AREA_OR_POINT=Area"
# [45] "Image Structure Metadata:"
# [46] " INTERLEAVE=PIXEL"
# [47] "Corner Coordinates:"
# [48] "Upper Left ( 2667400.000, 6479570.000) (174d45'39.47\"E, 36d52'25.92\"S)"
# [49] "Lower Left ( 2667400.000, 6478700.000) (174d45'40.19\"E, 36d52'54.14\"S)"
# [50] "Upper Right ( 2668010.000, 6479570.000) (174d46' 4.09\"E, 36d52'25.52\"S)"
# [51] "Lower Right ( 2668010.000, 6478700.000) (174d46' 4.81\"E, 36d52'53.73\"S)"
# [52] "Center ( 2667705.000, 6479135.000) (174d45'52.14\"E, 36d52'39.83\"S)"
# [53] "Band 1 Block=61x6 Type=Float32, ColorInterp=Gray"
# [54] " Description = glcm_variance"
# [55] " Min=0.000 Max=3.136 "
# [56] " Minimum=0.000, Maximum=3.136, Mean=-9999.000, StdDev=-9999.000"
# [57] " NoData Value=nan"
# [58] " Metadata:"
# [59] " STATISTICS_MAXIMUM=3.1358886957169"
# [60] " STATISTICS_MEAN=-9999"
# [61] " STATISTICS_MINIMUM=0"
# [62] " STATISTICS_STDDEV=-9999"
# [63] "Band 2 Block=61x6 Type=Float32, ColorInterp=Undefined"
# [64] " Description = glcm_contrast"
# [65] " Min=0.000 Max=1.917 "
# [66] " Minimum=0.000, Maximum=1.917, Mean=-9999.000, StdDev=-9999.000"
# [67] " NoData Value=nan"
# [68] " Metadata:"
# [69] " STATISTICS_MAXIMUM=1.9166666679084"
# [70] " STATISTICS_MEAN=-9999"
# [71] " STATISTICS_MINIMUM=0"
# [72] " STATISTICS_STDDEV=-9999"
# [73] "Band 3 Block=61x6 Type=Float32, ColorInterp=Undefined"
# [74] " Description = glcm_dissimilarity"
# [75] " Min=0.000 Max=1.167 "
# [76] " Minimum=0.000, Maximum=1.167, Mean=-9999.000, StdDev=-9999.000"
# [77] " NoData Value=nan"
# [78] " Metadata:"
# [79] " STATISTICS_MAXIMUM=1.1666666679084"
# [80] " STATISTICS_MEAN=-9999"
# [81] " STATISTICS_MINIMUM=0"
# [82] " STATISTICS_STDDEV=-9999"
# [83] "Band 4 Block=61x6 Type=Float32, ColorInterp=Undefined"
# [84] " Description = glcm_correlation"
# [85] " Min=-0.058 Max=0.759 "
# [86] " Minimum=-0.058, Maximum=0.759, Mean=-9999.000, StdDev=-9999.000"
# [87] " NoData Value=nan"
# [88] " Metadata:"
# [89] " STATISTICS_MAXIMUM=0.75854831188917"
# [90] " STATISTICS_MEAN=-9999"
# [91] " STATISTICS_MINIMUM=-0.058159347623587"
# [92] " STATISTICS_STDDEV=-9999"
# [93] "Band 5 Block=61x6 Type=Float32, ColorInterp=Undefined"
# [94] " Description = glcm_ASM"
# [95] " Min=0.062 Max=1.000 "
# [96] " Minimum=0.062, Maximum=1.000, Mean=-9999.000, StdDev=-9999.000"
# [97] " NoData Value=nan"
# [98] " Metadata:"
# [99] " STATISTICS_MAXIMUM=1"
# [100] " STATISTICS_MEAN=-9999"
# [101] " STATISTICS_MINIMUM=0.062187500298023"
# [102] " STATISTICS_STDDEV=-9999"
all.equal(values(tx1),values(tx2))
#TRUE
all.equal(values(tx1),values(tx3))
# "Mean relative difference: 2.535531e-08"
Running the code on the nambia data with valgrind enabled I get:
==35459==
==35459== HEAP SUMMARY:
==35459== in use at exit: 288,656,068 bytes in 87,678 blocks
==35459== total heap usage: 9,294,801,906 allocs, 9,294,714,228 frees, 5,689,651,603,813 bytes allocated
==35459==
==35459== LEAK SUMMARY:
==35459== definitely lost: 2,048 bytes in 1 blocks
==35459== indirectly lost: 652 bytes in 255 blocks
==35459== possibly lost: 0 bytes in 0 blocks
==35459== still reachable: 288,653,368 bytes in 87,422 blocks
==35459== of which reachable via heuristic:
==35459== newarray : 4,264 bytes in 1 blocks
==35459== suppressed: 0 bytes in 0 blocks
==35459== Rerun with --leak-check=full to see details of leaked memory
==35459==
==35459== For lists of detected and suppressed errors, rerun with: -s
==35459== ERROR SUMMARY: 0 errors from 0 contexts (suppressed: 0 from 0)
Additionally, throughout execution I got some warnings all in this sort of form:
==35459== Warning: set address range perms: large range [0xbb643028, 0x11cffa638) (noaccess)
==35459== Warning: set address range perms: large range [0x1227fb040, 0x1841b2650) (undefined)
==35459== Warning: set address range perms: large range [0x59c8b028, 0xbb642638) (noaccess)
This does not show where things go wrong. Did you, or can you, run this with --leak-check=full
and --track-origins=yes
?
Thanks, @rhijmans. I'll rerun it with those flags checked (sorry, new to valgrind).
@rhijmans, the code finished running through valgrind (output below). This iteration looks a lot like what @micha-silver was getting (it's "clipped" rather than striped or offset)
R -d "valgrind --leak-check=full --track-origins=yes --log-file=valgrindout.txt"
> rm(list=ls())
> library(terra)
library(GLCMTextures)
stk_full<- rast("Desktop/Test_data/namibia_full.tif")
grn <- stk_full$green
grn_quantized <- quantize_raster(r = grn, n_levels = 16, method = "equal prob")
textures <- glcm_textures(r=grn_quantized,
w=c(5, 5),
shift=list(c(1,0), c(1,1)),
metrics = c("glcm_variance",
"glcm_contrast",
"glcm_dissimilarity",
"glcm_correlation",
"glcm_ASM"),
quantization = "none",
n_levels=16, na.rm = FALSE)
> writeRaster(textures, "valgrindtest2.grd")
> quit()
One layer of resultant surface (grey) overlain on original
Valgrind Log:
==322711== Memcheck, a memory error detector
==322711== Copyright (C) 2002-2017, and GNU GPL'd, by Julian Seward et al.
==322711== Using Valgrind-3.17.0 and LibVEX; rerun with -h for copyright info
==322711== Command: /usr/lib/R/bin/exec/R
==322711== Parent PID: 33891
==322711==
==322711== Warning: set address range perms: large range [0x59c8b040, 0xbb642620) (undefined)
==322711== Warning: set address range perms: large range [0xbb643040, 0x11cffa620) (undefined)
==322711== Warning: set address range perms: large range [0xbb643028, 0x11cffa638) (noaccess)
==322711== Warning: set address range perms: large range [0x1227fb040, 0x1841b2650) (undefined)
==322711== Warning: set address range perms: large range [0x59c8b028, 0xbb642638) (noaccess)
==322711== Warning: set address range perms: large range [0x59c8b040, 0xbb642650) (undefined)
==322711== Warning: set address range perms: large range [0x1227fb028, 0x1841b2668) (noaccess)
==322711== Warning: set address range perms: large range [0xbb643040, 0xec31eb60) (undefined)
==322711== Warning: set address range perms: large range [0xec31f040, 0x11cffab60) (undefined)
==322711== Warning: set address range perms: large range [0xbb643028, 0xec31eb78) (noaccess)
==322711== Warning: set address range perms: large range [0x1227fb040, 0x1534d6b68) (undefined)
==322711== Warning: set address range perms: large range [0x3be71040, 0x53375ff0) (undefined)
==322711== Warning: set address range perms: large range [0x1227fb028, 0x1534d6b80) (noaccess)
==322711== Warning: set address range perms: large range [0xbb643040, 0xea04cf70) (undefined)
==322711== Warning: set address range perms: large range [0x1227fb040, 0x139cffff0) (undefined)
==322711== Warning: set address range perms: large range [0x3be71028, 0x53376008) (noaccess)
==322711== Warning: set address range perms: large range [0x139d01040, 0x16870af70) (undefined)
==322711== Warning: set address range perms: large range [0xec31f028, 0x11cffab78) (noaccess)
==322711== Warning: set address range perms: large range [0x1227fb028, 0x139d00008) (noaccess)
==322711== Warning: set address range perms: large range [0xbb643028, 0xea04cf88) (noaccess)
==322711== Warning: set address range perms: large range [0x16870b040, 0x197114f90) (undefined)
==322711== Warning: set address range perms: large range [0x3be71040, 0x55ca2560) (undefined)
==322711== Warning: set address range perms: large range [0xbb643040, 0xd5474560) (undefined)
==322711== Warning: set address range perms: large range [0xbb643028, 0xd5474578) (noaccess)
==322711== Warning: set address range perms: large range [0x3be71028, 0x55ca2578) (noaccess)
==322711== Warning: set address range perms: large range [0x3be71040, 0x55ca2560) (undefined)
==322711== Warning: set address range perms: large range [0xbb643040, 0xd5474560) (undefined)
==322711== Warning: set address range perms: large range [0xbb643028, 0xd5474578) (noaccess)
==322711== Warning: set address range perms: large range [0x3be71028, 0x55ca2578) (noaccess)
==322711== Warning: set address range perms: large range [0x3be71040, 0x55ca2560) (undefined)
==322711== Warning: set address range perms: large range [0xd7c75040, 0xf1aa6560) (undefined)
==322711== Warning: set address range perms: large range [0xd7c75028, 0xf1aa6578) (noaccess)
==322711== Warning: set address range perms: large range [0x3be71028, 0x55ca2578) (noaccess)
==322711== Warning: set address range perms: large range [0x3be71040, 0x4fd946c0) (undefined)
==322711== Warning: set address range perms: large range [0xd7c75040, 0xebb986c0) (undefined)
==322711== Warning: set address range perms: large range [0xd7c75028, 0xebb986d8) (noaccess)
==322711== Warning: set address range perms: large range [0x3be71028, 0x4fd946d8) (noaccess)
==322711== Warning: set address range perms: large range [0x2b271040, 0x42bcbb90) (undefined)
==322711== Warning: set address range perms: large range [0xdb85a040, 0xf31b4b90) (undefined)
==322711== Warning: set address range perms: large range [0xdb85a028, 0xf31b4ba8) (noaccess)
==322711== Warning: set address range perms: large range [0x2b271028, 0x42bcbba8) (noaccess)
==322711== Warning: set address range perms: large range [0x2b271040, 0x42bcbb90) (undefined)
==322711== Warning: set address range perms: large range [0xdb85a040, 0xf31b4b90) (undefined)
==322711== Warning: set address range perms: large range [0xdb85a028, 0xf31b4ba8) (noaccess)
==322711== Warning: set address range perms: large range [0x2b271028, 0x42bcbba8) (noaccess)
==322711== Warning: set address range perms: large range [0x2b271040, 0x42bcbb90) (undefined)
==322711== Warning: set address range perms: large range [0xdb85a040, 0xf31b4b90) (undefined)
==322711== Warning: set address range perms: large range [0xdb85a028, 0xf31b4ba8) (noaccess)
==322711== Warning: set address range perms: large range [0x2b271028, 0x42bcbba8) (noaccess)
==322711== Warning: set address range perms: large range [0x2b271040, 0x42bcbb90) (undefined)
==322711== Warning: set address range perms: large range [0xdb85a040, 0xf31b4b90) (undefined)
==322711== Warning: set address range perms: large range [0xdb85a028, 0xf31b4ba8) (noaccess)
==322711== Warning: set address range perms: large range [0x2b271028, 0x42bcbba8) (noaccess)
==322711== Warning: set address range perms: large range [0x3230b040, 0x4432a238) (undefined)
==322711== Warning: set address range perms: large range [0x4432b040, 0x5634a238) (undefined)
==322711== Warning: set address range perms: large range [0x4432b028, 0x5634a250) (noaccess)
==322711== Warning: set address range perms: large range [0xdac5a040, 0xecc79238) (undefined)
==322711== Warning: set address range perms: large range [0xdac5a028, 0xecc79250) (noaccess)
==322711== Warning: set address range perms: large range [0x3230b028, 0x4432a250) (noaccess)
==322711== Warning: set address range perms: large range [0x3230b040, 0x4432a238) (undefined)
==322711== Warning: set address range perms: large range [0x4432b040, 0x5634a238) (undefined)
==322711== Warning: set address range perms: large range [0x4432b028, 0x5634a250) (noaccess)
==322711== Warning: set address range perms: large range [0xdac5a040, 0xecc79238) (undefined)
==322711== Warning: set address range perms: large range [0xdac5a028, 0xecc79250) (noaccess)
==322711== Warning: set address range perms: large range [0x3230b028, 0x4432a250) (noaccess)
==322711== Warning: set address range perms: large range [0x3230b040, 0x4432a238) (undefined)
==322711== Warning: set address range perms: large range [0x4432b040, 0x5634a238) (undefined)
==322711== Warning: set address range perms: large range [0x4432b028, 0x5634a250) (noaccess)
==322711== Warning: set address range perms: large range [0xdac5a040, 0xecc79238) (undefined)
==322711== Warning: set address range perms: large range [0xdac5a028, 0xecc79250) (noaccess)
==322711== Warning: set address range perms: large range [0x3230b028, 0x4432a250) (noaccess)
==322711== Warning: set address range perms: large range [0x3230b040, 0x4432a238) (undefined)
==322711== Warning: set address range perms: large range [0x4432b040, 0x5634a238) (undefined)
==322711== Warning: set address range perms: large range [0x4432b028, 0x5634a250) (noaccess)
==322711== Warning: set address range perms: large range [0xdac5a040, 0xecc79238) (undefined)
==322711== Warning: set address range perms: large range [0xdac5a028, 0xecc79250) (noaccess)
==322711== Warning: set address range perms: large range [0x3230b028, 0x4432a250) (noaccess)
==322711== Warning: set address range perms: large range [0x3230b040, 0x4432a238) (undefined)
==322711== Warning: set address range perms: large range [0x4432b040, 0x5634a238) (undefined)
==322711== Warning: set address range perms: large range [0x4432b028, 0x5634a250) (noaccess)
==322711== Warning: set address range perms: large range [0xdb85a040, 0xed879238) (undefined)
==322711== Warning: set address range perms: large range [0xdb85a028, 0xed879250) (noaccess)
==322711== Warning: set address range perms: large range [0x3230b028, 0x4432a250) (noaccess)
==322711== Warning: set address range perms: large range [0x2ae71040, 0x3e546de8) (undefined)
==322711== Warning: set address range perms: large range [0xbb643040, 0xced18de8) (undefined)
==322711== Warning: set address range perms: large range [0xbb643028, 0xced18e00) (noaccess)
==322711== Warning: set address range perms: large range [0xdac5a040, 0xee32fde8) (undefined)
==322711== Warning: set address range perms: large range [0xdac5a028, 0xee32fe00) (noaccess)
==322711== Warning: set address range perms: large range [0x2ae71028, 0x3e546e00) (noaccess)
==322711== Warning: set address range perms: large range [0x2ae71040, 0x3e546de8) (undefined)
==322711== Warning: set address range perms: large range [0xbb643040, 0xced18de8) (undefined)
==322711== Warning: set address range perms: large range [0xbb643028, 0xced18e00) (noaccess)
==322711== Warning: set address range perms: large range [0xdac5a040, 0xee32fde8) (undefined)
==322711== Warning: set address range perms: large range [0xdac5a028, 0xee32fe00) (noaccess)
==322711== Warning: set address range perms: large range [0x2ae71028, 0x3e546e00) (noaccess)
==322711== Warning: set address range perms: large range [0x2ae71040, 0x3e546de8) (undefined)
==322711== Warning: set address range perms: large range [0xbb643040, 0xced18de8) (undefined)
==322711== Warning: set address range perms: large range [0xbb643028, 0xced18e00) (noaccess)
==322711== Warning: set address range perms: large range [0xdac5a040, 0xee32fde8) (undefined)
==322711== Warning: set address range perms: large range [0xdac5a028, 0xee32fe00) (noaccess)
==322711== Warning: set address range perms: large range [0x2ae71028, 0x3e546e00) (noaccess)
==322711== Warning: set address range perms: large range [0x2ae71040, 0x3e546de8) (undefined)
==322711== Warning: set address range perms: large range [0xbb643040, 0xced18de8) (undefined)
==322711== Warning: set address range perms: large range [0xbb643028, 0xced18e00) (noaccess)
==322711== Warning: set address range perms: large range [0xdac5a040, 0xee32fde8) (undefined)
==322711== Warning: set address range perms: large range [0xdac5a028, 0xee32fe00) (noaccess)
==322711== Warning: set address range perms: large range [0x2ae71028, 0x3e546e00) (noaccess)
==322711== Warning: set address range perms: large range [0x2ae71040, 0x3e546de8) (undefined)
==322711== Warning: set address range perms: large range [0xbb643040, 0xced18de8) (undefined)
==322711== Warning: set address range perms: large range [0xbb643028, 0xced18e00) (noaccess)
==322711== Warning: set address range perms: large range [0xe1c5a040, 0xf532fde8) (undefined)
==322711== Warning: set address range perms: large range [0xe1c5a028, 0xf532fe00) (noaccess)
==322711== Warning: set address range perms: large range [0x2ae71028, 0x3e546e00) (noaccess)
==322711== Warning: set address range perms: large range [0x139d01028, 0x16870af88) (noaccess)
==322711== Warning: set address range perms: large range [0x59c8b028, 0xbb642668) (noaccess)
==322711== Warning: set address range perms: large range [0x16870b028, 0x197114fa8) (noaccess)
==322711==
==322711== HEAP SUMMARY:
==322711== in use at exit: 289,290,833 bytes in 73,362 blocks
==322711== total heap usage: 4,844,927,956 allocs, 4,844,854,594 frees, 3,052,199,909,815 bytes allocated
==322711==
==322711== 2,700 (2,048 direct, 652 indirect) bytes in 1 blocks are definitely lost in loss record 6,538 of 9,083
==322711== at 0x4848C73: realloc (in /usr/libexec/valgrind/vgpreload_memcheck-amd64-linux.so)
==322711== by 0xF57B20A: VSIReallocVerbose (in /usr/lib/libgdal.so.28.0.2)
==322711== by 0xF51C79D: CSLAddStringMayFail (in /usr/lib/libgdal.so.28.0.2)
==322711== by 0xF51C81C: CSLAddString (in /usr/lib/libgdal.so.28.0.2)
==322711== by 0xE210653: setCats(GDALRasterBand*, std::vector<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> >, std::allocator<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> > > >&) (write_gdal.cpp:90)
==322711== by 0xE2109F4: setBandCategories(GDALRasterBand*, SpatDataFrame&, std::vector<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> >, std::allocator<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> > > >) (write_gdal.cpp:163)
==322711== by 0xE21565A: SpatRaster::writeStartGDAL(SpatOptions&) (write_gdal.cpp:503)
==322711== by 0xE2160E8: SpatRaster::writeStart(SpatOptions&) (write.cpp:252)
==322711== by 0xE1ABF5C: SpatRaster::reclassify(std::vector<std::vector<double, std::allocator<double> >, std::allocator<std::vector<double, std::allocator<double> > > >, unsigned int, bool, bool, bool, bool, SpatOptions&) (raster_methods.cpp:3042)
==322711== by 0xE1AD1E6: SpatRaster::reclassify(std::vector<double, std::allocator<double> >, unsigned int, unsigned int, bool, bool, bool, bool, SpatOptions&) (raster_methods.cpp:3100)
==322711== by 0xE0E8568: Rcpp::CppMethod8<SpatRaster, SpatRaster, std::vector<double, std::allocator<double> >, unsigned int, unsigned int, bool, bool, bool, bool, SpatOptions&>::operator()(SpatRaster*, SEXPREC**) (Module_generated_CppMethod.h:783)
==322711== by 0xE11692B: Rcpp::class_<SpatRaster>::invoke_notvoid(SEXPREC*, SEXPREC*, SEXPREC**, int) (class.h:234)
==322711==
==322711== LEAK SUMMARY:
==322711== definitely lost: 2,048 bytes in 1 blocks
==322711== indirectly lost: 652 bytes in 255 blocks
==322711== possibly lost: 0 bytes in 0 blocks
==322711== still reachable: 289,288,133 bytes in 73,106 blocks
==322711== of which reachable via heuristic:
==322711== newarray : 4,264 bytes in 1 blocks
==322711== suppressed: 0 bytes in 0 blocks
==322711== Reachable blocks (those to which a pointer was found) are not shown.
==322711== To see them, rerun with: --leak-check=full --show-leak-kinds=all
==322711==
==322711== For lists of detected and suppressed errors, rerun with: -s
==322711== ERROR SUMMARY: 1 errors from 1 contexts (suppressed: 0 from 0)
Thanks. The valgrind error seems unrelated (and may be fixed by now). I still think this has to do with the size of integer values. Somehow, on each line, every other value is not read; or ignored after reading. Anyway, I am running some tests to see if I can find out.
I just cannot reproduce this. I wonder if this still gives the same bad results for you:
library(terra)
library(GLCMTextures)
r <- rast("namibia_full.tif")$green
grn_quantized <- quantize_raster(r, n_levels = 16, method = "equal prob", filename="quant.tif")
r <- rast("quant.tif")
e <- ext(c(320153, 320710, 7883350, 7883728))
r <- crop(r, e)
txt <- glcm_textures(r=r, w=c(5, 5), shift=list(c(1,0), c(1,1)),
metrics = c("glcm_variance", "glcm_contrast","glcm_dissimilarity", "glcm_correlation", "glcm_ASM"),
quantization = "none", n_levels=16, na.rm = FALSE, filename="text.tif")
# or with very low mem availability
terraOptions(todisk=T, steps=5, memmax=5)
That is, it would be nice to have much smaller example.
@rhijmans , I'll run some tests on smaller versions, with the parameters you suggested, and update terra. I was able to reproduce something similar on Windows though with a much smaller raster (rspatial/terra#519). It's still a fairly large raster but it runs in minutes rather than 10's of hours.
@rhijmans, the cropped example runs fine even with setting very low memory availability.
@micha-silver, I believe this is now fixed if you install the latest development version of terra
(rspatial/terra#519). I ran it on my Linux VM and it is no longer clipped.
Many thanks!
Meanwhile I had prepared a short script to create a synthetic raster and run the glcm_textures
. Now I can do some more testing...
# Packages for test
#-----------------------------------------------
pkg_list = c("microbenchmark", "glcm",
"raster", "terra",
"remotes" )
installed_packages <- pkg_list %in% rownames(installed.packages())
if (any(installed_packages == FALSE)) {
install.packages(pkg_list[!installed_packages])
}
# Load Packages
lapply(pkg_list, function(p) {
require(p, character.only = TRUE, quietly=TRUE)})
# Get recent GLCMTextures packages from Github
remotes::install_github("ailich/GLCMTextures", ref="arma_test")
library(GLCMTextures)
# Setup
#-----------------------------------------------
n <- 12000
m <- n*1.3
# Raster with NA in corners
fraction <- 0.25
to_na <- as.integer(fraction*m)
M <- matrix(runif(n*m), nrow=m)
for (i in 1:m) {
M[i, 1:(to_na-i*fraction)] <- NA
M[i, (n-(i*fraction)):n] <- NA
}
R <- rast(M)
# GLCM parameters
levels <- 16
wnd <- c(5, 5)
shft <- list(c(1,0))
stats <- c("glcm_variance")
cat("Window:", wnd, "\nShift:", as.character(shft),"\nGrey levels:", as.character(levels))
cat("\nRaster size:", as.character(dim(R)),
"\nGLCM statistics:", as.character(stats), "\n")
# Function
#----------------------------------------------
AddTextures <- function(x,
window = wnd,
grey_levels=levels,
shift=shft,
glcm_stats=stats) {
x_quantized <- quantize_raster(x, n_levels = grey_levels, method = "equal prob")
textures <- glcm_textures(r=x_quantized,
w=window,
n_levels=grey_levels,
shift=shift,
metrics = glcm_stats,
quantization = "none",
na.rm = TRUE)
return(textures)
} # End of AddTextures
# Micro benchmark
#-----------------------------------------------
print(microbenchmark(R.texture <- AddTextures(R),
# larger window
R.texture.15 <- AddTextures(R, window=c(15, 15), grey_levels=32),
# isotropic direction textures
R.texture.iso <- AddTextures(R,
shift=list(c(1,0), c(1,1), c(0,1), c(-1,1)),
glcm_stats = c("glcm_variance",
"glcm_homogeneity")),
times=1))
@micha-silver, thanks! I'm going to close this issue since it's resolved now, and we can continue the speed tests on the other issues thread.
@ailich Here are the results I am getting using the test code above
The original, synthetic raster:
A single, straightforward GLCM variable, 5x5 window
A larger window (15x15)
But when I request multiple variables, and full directionality (using the 4 shift directions) I get:
For completeness: (running on my laptop with 16GB RAM)
> sessionInfo()
R version 4.1.2 (2021-11-01)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Debian GNU/Linux 11 (bullseye)
Matrix products: default
BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.13.so
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] GLCMTextures_0.3.4 remotes_2.4.2 terra_1.5-17
[4] raster_3.5-15 sp_1.4-6 glcm_1.6.5
[7] microbenchmark_1.4.9
loaded via a namespace (and not attached):
[1] Rcpp_1.0.8 MASS_7.3-55 splines_4.1.2 lattice_0.20-45
[5] R6_2.5.1 multcomp_1.4-18 tools_4.1.2 pkgbuild_1.3.1
[9] grid_4.1.2 TH.data_1.1-0 cli_3.1.1 withr_2.4.3
[13] survival_3.2-13 rprojroot_2.0.2 crayon_1.4.2 Matrix_1.4-0
[17] processx_3.5.2 callr_3.7.0 ps_1.6.0 codetools_0.2-18
[21] curl_4.3.2 sandwich_3.0-1 compiler_4.1.2 prettyunits_1.1.1
[25] mvtnorm_1.1-3 zoo_1.8-9
Regards, Micha
@micha-silver, try updating terra
to the latest development version. I see you have 1.5-17 but it's up to 1.5-20 now. That should fix the striping issue.
We are using GLCMTextures to produce several texture layers for a large study ares - the extent is a bit more than 12,000 x 16,000. Running
quantize_raster()
in advance produces the correct 16 bit raster for our whole study area. But then, after runningglcm_textures()
the result is clipped somehow. See images below.I reran the function a few times on my home computer - 16GB RAM - with identical results. Then I tried on our dept server - with 128 GB RAM - and the result looked fine, the texture bands all covered the whole study area.
Is there some memory limit that would cause
glcm_textures
to do this clipping? (Just to be clear, on the low mem computer,quantize_raster
worked fine, resulting in a raster covering the whole region. Onlyglcm_textures
failed)Full study area:![Full study area](https://user-images.githubusercontent.com/3841768/147911547-bef22743-b4e9-4152-b5fc-6a76ebda53db.png)
Texture bands, after running glcm_textures:![GLCM textures result](https://user-images.githubusercontent.com/3841768/147911562-3f8af344-db84-4a18-a46a-eca30db50639.png)
Thanks, Micha