ailich / GLCMTextures

This R package calculates the most common gray-level co-occurrence matrix (GLCM) texture metrics used for spatial analysis on raster data.
https://ailich.github.io/GLCMTextures/
GNU General Public License v3.0
12 stars 4 forks source link

add ... or wopt #15

Closed rhijmans closed 2 years ago

rhijmans commented 2 years ago

With glcm_textures and quantize_raster one cannot pass options to terra. That would be useful for debugging; to use less common filetypes; and to manage memory use without setting global options.

As these functions do not have ellipses for other purposes, you could perhaps use them for this (and pass them on as wopt=list(...)

ailich commented 2 years ago

Thanks @rhijmans for the suggestion. I added wopt to those two functions with a default of wopt=list(), and passed it to what I believe are the relevant terra functions: classify, as.int, writeRaster, focalCpp, subset, and app. Should that be sufficient?

rhijmans commented 2 years ago

That sounds good. I think that fun(, filename="", ...) is cleaner and easier to use, but that is up to you, and you need to consider if you want to use the ellipses in the future for other purposes.

ailich commented 2 years ago

If wopt is specified in the function is it still ok to perform standard mathematical and trigonometric operations, or should I attempt to wrap them into a call to app or lapp so that I can specify wopt in the calculations? For example in my other R package I have code such as:

 slp<- atan(sqrt(params$d^2 + params$e^2))
 slp<- slp*(180/pi)
 asp<- atan2(params$e,params$d)
 eastness<- sin(asp)
 rdmv<- (r - localmean)/(localsd) #r, localmean, and localsd are single layer SpatRasters
rhijmans commented 2 years ago

Sorry, I missed this one. It is a bit tricky. If you do one or two of these "build-in" functions it might be faster to do them separately. But if there are many app or lapp might be faster if the files are large such that temp files are written (and that is when speed matters most).

slpaslp <- lapp(params[[c("d", "e")]], function(x, y) {
        slp <- atan(sqrt(x^2 + y^2))
        slp <- slp*(180/pi)
        asp <- atan2(y, x)
        cbind(slp, sin(asp))
    }   
)

wopt can be useful for different things, including for debugging, but the most important place is the final terra::method were the output is written to disk.

ailich commented 2 years ago

Thanks @rhijmans but I have actually found app to be substantially slower than the generic functions for large data. For example calculating the product of a 6 layer raster with 16,000 rows and 12,500 columns is over 32 times faster with prod than with app. There's some functions such as prod are listed as in the Raster Algebra online documentation but not included as a method in terra::math (also atan2 works but is listed in neither). Using the function in a call to terra::math would allow specification of wopt. Additionally, given the differences in speed maybe it would be useful to have terra::app use terra::math if it recognizes the function name?


library(terra)
r<- rast(nrow=16000, ncol=12500, nlyr=6)
set.seed(5)
values(r)<- sample(c(0,1),size = ncell(r)*nlyr(r), replace = TRUE, prob = c(0.1,0.9))

st1<- Sys.time()
v1<- prod(r)
end1<- Sys.time()
print(difftime(end1, st1, units="secs"))
#Time difference of 15.01893 secs

st2<- Sys.time()
v2<- app(r, fun=prod)
end2<- Sys.time()
print(difftime(end2, st2, units="secs"))
#Time difference of 489.1865 secs
rhijmans commented 2 years ago

It is much faster when you use "prod". Here using a smaller variation of your example:

library(terra)
r<- rast(nrow=1600, ncol=12500, nlyr=6)
set.seed(5)
values(r)<- sample(c(0,1),size = ncell(r)*nlyr(r), replace = TRUE, prob = c(0.1,0.9))

system.time(v1 <- prod(r))
#   user  system elapsed 
#   0.41    0.22    0.62 

system.time(v2 <- app(r, fun="prod"))
#   user  system elapsed 
#   0.49    0.11    0.59 

So app did not recognize prod as being equivalent to "prod", as it did for, e.g. mean, min max. That has now been corrected and I get

system.time(v2 <- app(r, fun=prod))
#   user  system elapsed 
#   0.43    0.17    0.61 
rhijmans commented 2 years ago

prod is not "Math" generic, it is a "Summary" generic and you can pass wopt, according to the documentation (that is a bit hard to find); see ?terra::`Summary-methods`

ailich commented 2 years ago

@rhijmans, thanks for explanation, clarification, and fix! I guess just atan2 currently doesn't have the wopt argument.

rhijmans commented 2 years ago

It does not. See ?atan2. I could perhaps implement a fast version with wopt etc via lapp

rhijmans commented 2 years ago

I have added atan.2 for this. That is simpler and easier to discover.

ailich commented 2 years ago

Awesome, thank you. I'll make changes based on your suggestions and the new additions.

rhijmans commented 2 years ago

R CMD check does not like atan.2 (it thinks it is the atan S3 method for objects of class "2") so I renamed it to atan_2. Also, ?terra::prod now works.

ailich commented 2 years ago

@rhijmans, not really essential for me as I just import terra so I can just call prod without specifying it's from terra, but terra::prod appears to have been removed from the namespace as I now get Error: 'prod' is not an exported object from 'namespace:terra'

rhijmans commented 2 years ago

prod is defined as a member of the Summary-Generics. See ?S4groupGeneric
Therefore, terra::prod does not exist, but you can use prod , same for sum, min.

Specifying the package:: for generic methods is not that meaningful anyway, because the method used depends on the argument. For example, there is no difference between terra::crop(x, e) and raster::crop(x, e). The method used depends on the class of x.

ailich commented 2 years ago

Gotcha, thanks!