[R-sig-Geo] Problem with raster::calc

[Robert J. Hijmans]

Agus,

I think what you want to accomplish can be done like this:

x <- predict(SGCIRF56, sgcirRpred)

but if you want to use calc, I think it should be something like

sgcirRpred = calc(SGCIRF56, fun=function(x){coefficients(sgcirRlm)[1]
+ coefficients(sgcirRlm)[2:4] * x })

But in this case you need to make sure that the order of the layers in
SGCIRF56 is the same as in the coefficients.

subset(x,3)
is, it seems to me, an invalid statement
but following your approach, this might work:

sgcirRpred = calc(SGCIRF56,fun=function(x){coefficients(sgcirRlm)[1]+coefficients(sgcirRlm)[2]*x[1]
+ coefficients(sgcirRlm)[3]*x[2] + coefficients(sgcirRlm)[4]*x[3] } )

Best, Robert


On Fri, Apr 29, 2011 at 2:57 AM, Agustin Lobo <alobolistas at gmail.com> wrote:
> I've calculated a linear model using a subset of data from raster files
>> sgcirRlm = lm(SGRGBF40data[,4] ~ SGCIRF56data[,4]+SGCIRF56data[,5]+SGCIRF56data[,6]
>> coefficients(sgcirRlm)
>      (Intercept) SGCIRF56data[, 4] SGCIRF56data[, 5] SGCIRF56data[, 6]
>    -2075.0432230         0.9229117         0.9963710        -0.8310794
>
> and want to apply the coefficients to the multiband raster
>> show(SGCIRF56)
> class       : RasterBrick
> dimensions  : 1760, 2640, 3  (nrow, ncol, nlayers)
> resolution  : 1, 1  (x, y)
> extent      : 0, 2640, -1759, 1  (xmin, xmax, ymin, ymax)
> projection  : +proj=utm +zone=31 +ellps=intl +units=m +no_defs
> values      : /media/Iomega_HDD/UAVetal/CALIBRACIONRADIOM/TESTCASA/CALSEL/SGCIR/SGCIRWBPS125F56v2.tif
> min values  : 0 0 0
> max values  : 65535 65535 65535
>
> to calculate the predicted raster, for which I'm using subset() within calc():
>> sgcirRpred = calc(SGCIRF56,fun=function(x){coefficients(sgcirRlm)[1]+coefficients(sgcirRlm)[2]*subset(x,1) + coefficients(sgcirRlm)[3]*subset(x,2) + coefficients(sgcirRlm)[4]*subset(x,3)})
>
> but this results into an error:
> Error in .local(x, fun, ...) : cannot use this function
> Calls: calc -> calc -> .local
>
> is this because subset() cannot be used within the function to be
> provided to calc()?
>
> Thanks
>> sessionInfo()
> R version 2.13.0 (2011-04-13)
> Platform: x86_64-pc-linux-gnu (64-bit)
>
> 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=en_US.UTF-8
>  [9] LC_ADDRESS=en_US.UTF-8        LC_TELEPHONE=en_US.UTF-8
> [11] LC_MEASUREMENT=en_US.UTF-8    LC_IDENTIFICATION=en_US.UTF-8
>
> attached base packages:
> [1] grid      stats     graphics  grDevices utils     datasets  methods
> [8] base
>
> other attached packages:
> [1] gridExtra_0.7 ggplot2_0.8.9 rgdal_0.6-33  raster_1.8-16 sp_0.9-76
> [6] reshape_0.8.4 plyr_1.4      proto_0.3-8   rkward_0.5.6
>
> loaded via a namespace (and not attached):
> [1] lattice_0.19-26 tools_2.13.0
>
> Agus
>