[R-sig-Geo] Comparing 2 rasterbricks using functions from TSdist package

Loïc Dutrieux |o|c@dutr|eux @end|ng |rom c|r@d@|r
Thu May 28 12:49:25 CEST 2020


Hi Jackson,

calc is meant to apply a function pixel wise on a single raster object. 
Its equivalent for multiple raster objects is overlay.

Kind regards,
Loïc

On 05/28/2020 09:05 AM, Bede-Fazekas Ákos wrote:
> Dear Jackson,
> 
> I think this is what you are searching for:
> new_s3 <- setValues(r, sapply(1:ncell(r), function(i) 
> fun(as.vector(s[i]), as.vector(s2[i]))))
> plot(new_s3, colNA = "red")
> 
> HTH,
> Ákos Bede-Fazekas
> Hungarian Academy of Sciences
> 
> 
> 
> 2020.05.28. 0:18 keltezéssel, Jackson Rodrigues írta:
>> Dear all,
>>
>> Few days ago I got a important help from this list and I am really 
>> thankful
>> for that.
>> Today I have a similar challenge with 2 rasterbrick files.
>>
>> I would like to compare 2 rasterbricks using any similarity function 
>> (let's
>> take CorDistance) from package TSdist.
>> In my real data, both rasterbricks came from the same region but using
>> different measurement techniques.
>>
>> I am trying to adapt the suggestion I got from here rather than 
>> looping on
>> each cell, but no matter what I try I always get an error related to NA.
>> It seems to be very simples but I am stuck since last Monday.
>>
>> I have played with several conditional ("ifelse", "if", "if else" etc) 
>> and
>> logical tests ("|", "||", "&" etc) but I got nothing so far.
>>
>> Could someone help me and show me a path to go through?
>> Below there is the cleanest code I have tried.
>>
>> Best regards
>>
>> Jackson
>>
>> ##############
>> library(raster)
>> library(TSdist)
>>
>> set.seed(12)
>> r <- raster(nrow=10, ncol=10)
>> s <- lapply(1:200, function(i) setValues(r, rnorm(ncell(r),
>>                                                    sample.int(5,1), 
>> 0.5)))
>> s <- stack(s)
>> s[s < 0] <- NA
>> s2<-s^2
>>
>> # Visualize that some pixels have NAs and other don't
>> hasna <- stackApply(s, indices = 1, fun = function(x, na.rm){anyNA(x)})
>> hasna2 <- stackApply(s2, indices = 1, fun = function(x, na.rm){anyNA(x)})
>>
>> par(mfrow=c(1,2))
>> plot(hasna);plot(hasna2)
>>
>> # CorDistance function does not handle NAs
>> CorDistance(as.vector(s[1]),as.vector(s2[1]))
>> [1] 0.2192382
>>
>> CorDistance(as.vector(s[2]),as.vector(s2[2]))
>> <simpleError in .common.ts.sanity.check(x): NA in the series>
>>    [1] NA
>>
>> fun <- function(x,y){
>>    out <- ifelse(anyNA(x) | anyNA(y),
>>                  yes = NA,
>>                  no = unname(CorDistance(x,y)))
>>    return(out)
>> }
>>
>> # Check that it works
>> fun(as.vector(s[1]),as.vector(s2[1]))
>> [1] 0.2192382
>>
>> # Check that it does not work
>> fun(as.vector(s[2]),as.vector(s2[2]))
>> [1] NA
>>
>> new_s1 <- fun(s,s2)
>> <simpleError in if (any(is.na(x))) { stop("NA in the series")}:
>>    argument is not interpretable as logical>
>>
>> # Check that it does not work with calc
>> new_s2 <- calc(s,s2, fun = fun)
>> Error in .calcTest(x[1:5], fun, na.rm, forcefun, forceapply) :
>>    cannot use this function
>>
>> plot(s_out)
>> #############
>>
> 
> 
> 
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo using r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> 



More information about the R-sig-Geo mailing list