[R-sig-Geo] trouble with the 'focal' function in 'raster'

[Robert J. Hijmans]

Marcia,

This is a bug that is caused by
> mean(NA, na.rm=T)
[1] NaN

Because of this, "newr" will have NaN values. These are equivalent to
NA in the R "is.na" function but not in the C macro R_ISNA that is
used by focal.  Anyway, here is a work-around for the current version
of raster:

library(raster)
r <- raster(ncols=60, nrows=60, xmn=0)
r[] <- 1:(60*60)
values(r)[5]<-NA
values(r)[(30*30):(60*60)]<-NA
values(r)[round(runif(3000,1,3600))]<-NA
w = matrix(1,3,3)

for (i in 1:40){
newr <- focal(r, w=w, fun=mean, na.rm=T, pad=TRUE, padValue=NA, NAonly=TRUE)
r <- calc(newr, function(x){ x[is.nan(x)] = NA; x})
}

On Mon, Jun 24, 2013 at 7:14 AM, Koen Hufkens <koen.hufkens at gmail.com> wrote:
> This works on a system on which your code doesn't (ubuntu x64 / R 3.0.1 /
> raster 2.1-37).
>
> library(raster)
> r <- raster(ncols=60, nrows=60, xmn=0)
> r[] <- 1:(60*60)
> values(r)[5]<-NA
> values(r)[(30*30):(60*60)]<-NA
> values(r)[round(runif(3000,1,3600))]<-NA
>
> for (i in 1:40){
>   if (i==1){
>     writeRaster(focal(r,w=3, fun=mean, na.rm=F, pad=TRUE, padValue=NA,
> NAonly=TRUE),'tmp.tif',overwrite=TRUE)
>   }else{
>   t <- raster('tmp.tif')
>   writeRaster(focal(t,w=3, fun=mean, na.rm=F, pad=TRUE, padValue=NA,
> NAonly=TRUE),'tmp.tif',overwrite=TRUE)
>   }
>
> } # end loop i
>
> plot(raster('tmp.tif'))
>
>
> Somehow the history of the file gets in the way of your iterative approach,
> writing to file decouples the original from the new file. This is either a
> safeguard, or the file structure/content isn't properly updated, Robert
> Hijmans should know the details on which is which.
>
> Cheers,
> Koen
>
>
>
> On Mon, Jun 24, 2013 at 3:00 PM, Marcia Macedo <marcia.n.macedo at gmail.com>wrote:
>
>> Hi All,
>>
>> I am trying to use the focal function to iteratively fill NA's in a raster,
>> but am having trouble running it on my system. A colleague who was kindly
>> helping me troubleshoot my code noticed that it ran as expected with an
>> earlier version of R and the 'raster' package.
>>
>> The following example runs well on his setup (i.e., NA's are gradually
>> filled in with each iteration), but not on mine (i.e., the raster only
>> changes in the first iteration).
>>
>> #################################################################
>> # A test raster with plenty of NAs
>>
>> r <- raster(ncols=60, nrows=60, xmn=0)
>> r[] <- 1:(60*60)
>> values(r)[5]<-NA
>> values(r)[(30*30):(60*60)]<-NA
>> values(r)[round(runif(3000,1,3600))]<-NA
>>
>> for (i in 1:40){
>>
>> newr <- focal(r,w=3, fun=mean, na.rm=T, pad=TRUE, padValue=NA, NAonly=TRUE)
>>
>>  print(table(is.na(r[]),is.na(newr[])))
>>
>> r<-newr
>>
>> } # end loop i
>>
>> #################################################################
>> *The above code works on:*
>>
>> (Windows 7)
>> R version 2.15.0
>> raster 1.9-99 (1-June-2012)
>>
>> *It doesn't work on:*
>>
>> (Windows 7)
>> R version 3.0.1
>> raster 2.1-38
>>
>> The latter setup also gives the following warning message with each
>>  iteration:
>>
>> *In .getW(w): the computation of the weights matrix has changed in version
>> 2.1-35. The sum of weights is now 1.*
>>
>> Is it possible that the arguments for the function have been changed in the
>> most recent update, but not in the documentation? Any thoughts about what
>> might be going on would much appreciated!
>>
>> Best,
>> Marcia
>>
>> ---------------------------------------
>> *Marcia Macedo*
>> Postdoctoral Fellow
>> Woods Hole Research Center
>> mmacedo at whrc.org
>> 508.444.1538
>>
>>         [[alternative HTML version deleted]]
>>
>> _______________________________________________
>> R-sig-Geo mailing list
>> R-sig-Geo at r-project.org
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>
>
>         [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo