[R-sig-Geo] Error in zapply?

Loïc loic.dutrieux at wur.nl
Sat Sep 20 16:02:09 CEST 2014


Hi Srivinas,

The as.POSIXct() is probably the reason why it falls back in the 
previous year for January 1st. If dates is sufficient, then use as.Date().

See reproducible example:

library(raster)
library(zoo)

time <- as.Date("2004-01-01") + seq(0,365,8)
r <- raster(nc=10, nr=10)
s <- brick(lapply(1:length(time), function(x) setValues(r, 
G0dm[x]+100*rnorm(ncell(r)) )))
s <- setZ(s, time)

zApply(s, by = as.yearmon, fun = max)

Hope this helps,
Loïc


On 19/09/2014 10:44, Srinivas V wrote:
> Hi,
>
> I have a raster brick of MODIS 8 day LST data and I'm trying to generate
> a brick with monthly maximum LST values for the time period 2000-2013.
> I'm using the zoo and raster packages to get monthly averages, however
> it appears that the zapply function does not handle 1st January very
> well.To illustrate the problem I have provided outputs for a single
> year, instead of returning a brick with 12 layers it returns a brick
> with 13 layers and the first layer is assigned to December of the
> previous year.
>
> I have also manually changed the time stamp from 1st January to 3rd
> January and it works fine. However to do this across multiple tiles and
> across years is going to problematic and I cannot automate the
> processing.  Is this a bug and what is the alternative?
>
> I would appreciate any advice on dealing with this issue.
>
> Thanks!
>
> library(raster)
> library(zoo)
>
>
>   > dy #dates extracted from file names
>    [1] "2010-01-01 IST" "2010-01-09 IST" "2010-01-17 IST" "2010-01-25
> IST" "2010-02-02 IST" "2010-02-10 IST" "2010-02-18 IST" "2010-02-26 IST"
> "2010-03-06 IST" "2010-03-14 IST" "2010-03-22 IST" "2010-03-30 IST"
> [13] "2010-04-07 IST" "2010-04-15 IST" "2010-04-23 IST" "2010-05-01 IST"
> "2010-05-09 IST" "2010-05-17 IST" "2010-05-25 IST" "2010-06-02 IST"
> "2010-06-10 IST" "2010-06-18 IST" "2010-06-26 IST" "2010-07-04 IST"
> [25] "2010-07-12 IST" "2010-07-20 IST" "2010-07-28 IST" "2010-08-05 IST"
> "2010-08-13 IST" "2010-08-21 IST" "2010-08-29 IST" "2010-09-06 IST"
> "2010-09-14 IST" "2010-09-22 IST" "2010-09-30 IST" "2010-10-08 IST"
> [37] "2010-10-16 IST" "2010-10-24 IST" "2010-11-01 IST" "2010-11-09 IST"
> "2010-11-17 IST" "2010-11-25 IST" "2010-12-03 IST" "2010-12-11 IST"
> "2010-12-19 IST" "2010-12-27 IST"
>
>   > as.yearmon(dy)# no problem here, as.yearmon() is a function from
> package zoo
>
>    [1] "Jan 2010" "Jan 2010" "Jan 2010" "Jan 2010" "Feb 2010" "Feb 2010"
> "Feb 2010" "Feb 2010" "Mar 2010" "Mar 2010" "Mar 2010" "Mar 2010" "Apr
> 2010" "Apr 2010" "Apr 2010" "May 2010" "May 2010" "May 2010" "May 2010"
> [20] "Jun 2010" "Jun 2010" "Jun 2010" "Jun 2010" "Jul 2010" "Jul 2010"
> "Jul 2010" "Jul 2010" "Aug 2010" "Aug 2010" "Aug 2010" "Aug 2010" "Sep
> 2010" "Sep 2010" "Sep 2010" "Sep 2010" "Oct 2010" "Oct 2010" "Oct 2010"
> [39] "Nov 2010" "Nov 2010" "Nov 2010" "Nov 2010" "Dec 2010" "Dec 2010"
> "Dec 2010" "Dec 2010"
>
>   > LSTWG<-setZ(LSTWG,as.POSIXct(dy))#setZ used to set time stamp
>
>   > LSTWG #raster brick for year 2010
> class       : RasterBrick
> dimensions  : 3334, 2294, 7648196, 46  (nrow, ncol, ncell, nlayers)
> resolution  : 1000.16, 999.8986  (x, y)
> extent      : -728812.9, 1565553, -1033.987, 3332628  (xmin, xmax, ymin,
> ymax)
> coord. ref. : +proj=utm +no_defs +zone=43 +a=6378137 +rf=298.257223563
> +towgs84=0.000,0.000,0.000 +to_meter=1
> data source : in memory
> names       : X2010.01.01, X2010.01.09, X2010.01.17, X2010.01.25,
> X2010.02.02, X2010.02.10, X2010.02.18, X2010.02.26, X2010.03.06,
> X2010.03.14, X2010.03.22, X2010.03.30, X2010.04.07, X2010.04.15,
> X2010.04.23, ...
> min values  :       13279,       13136,       13159, 13201,
> 13116,       13179,       13379,       13522, 13622,       13758,
> 13901,       13865,       14132, 14007,       13915, ...
> max values  :       15612,       15722,       15766, 15734,
> 15808,       16005,       16145,       16223, 16409,       16493,
> 16333,       16496,       16546, 16694,       16716, ...
> time        : 2010-01-01, 2010-12-27 (min, max)
>
>   > LSTWGm<-zApply(LSTWG,by=as.yearmon,fun=max)
>
>   > LSTWGm #instead of returning a brick with 12 layers it returns a
> brick with 13 layers and the first layer is assigned to December of the
> previous year
> class       : RasterBrick
> dimensions  : 3334, 2294, 7648196, 13  (nrow, ncol, ncell, nlayers)
> resolution  : 1000.16, 999.8986  (x, y)
> extent      : -728812.9, 1565553, -1033.987, 3332628  (xmin, xmax, ymin,
> ymax)
> coord. ref. : +proj=utm +no_defs +zone=43 +a=6378137 +rf=298.257223563
> +towgs84=0.000,0.000,0.000 +to_meter=1
> data source :
> /tmp/R_raster_karthik/raster_tmp_2014-09-09_114458_2361_81754.grd
> names       : Dec.2009, Jan.2010, Feb.2010, Mar.2010, Apr.2010,
> May.2010, Jun.2010, Jul.2010, Aug.2010, Sep.2010, Oct.2010, Nov.2010,
> Dec.2010
> min values  :    13279,    13201,    13522,    13901,    14132,
> 14360,    13865,    12649,    12048,    14142,    14031, 13868,    13667
> max values  :    15612,    15766,    16223,    16496,    16716,
> 16647,    16671,    16357,    16439,    16467,    16382, 15886,    15671
>               : Dec 2009, Jan 2010, Feb 2010, Mar 2010, Apr 2010, May
> 2010, Jun 2010, Jul 2010, Aug 2010, Sep 2010, Oct 2010, Nov 2010, Dec 2010
>
>   >fix(dy)#changed time stamp of first tile to 3rd January
>
>   > dy
>    [1] "2010-01-03 IST" "2010-01-09 IST" "2010-01-17 IST" "2010-01-25
> IST" "2010-02-02 IST" "2010-02-10 IST" "2010-02-18 IST" "2010-02-26 IST"
> "2010-03-06 IST" "2010-03-14 IST" "2010-03-22 IST" "2010-03-30 IST"
> [13] "2010-04-07 IST" "2010-04-15 IST" "2010-04-23 IST" "2010-05-01 IST"
> "2010-05-09 IST" "2010-05-17 IST" "2010-05-25 IST" "2010-06-02 IST"
> "2010-06-10 IST" "2010-06-18 IST" "2010-06-26 IST" "2010-07-04 IST"
> [25] "2010-07-12 IST" "2010-07-20 IST" "2010-07-28 IST" "2010-08-05 IST"
> "2010-08-13 IST" "2010-08-21 IST" "2010-08-29 IST" "2010-09-06 IST"
> "2010-09-14 IST" "2010-09-22 IST" "2010-09-30 IST" "2010-10-08 IST"
> [37] "2010-10-16 IST" "2010-10-24 IST" "2010-11-01 IST" "2010-11-09 IST"
> "2010-11-17 IST" "2010-11-25 IST" "2010-12-03 IST" "2010-12-11 IST"
> "2010-12-19 IST" "2010-12-27 IST"
>
>   > LSTWG<-setZ(LSTWG,as.POSIXct(dy))
>
>   > LSTWG #raster brick for year 2010, time stamp manually changed.
>
> class       : RasterBrick
> dimensions  : 3334, 2294, 7648196, 46  (nrow, ncol, ncell, nlayers)
> resolution  : 1000.16, 999.8986  (x, y)
> extent      : -728812.9, 1565553, -1033.987, 3332628  (xmin, xmax, ymin,
> ymax)
> coord. ref. : +proj=utm +no_defs +zone=43 +a=6378137 +rf=298.257223563
> +towgs84=0.000,0.000,0.000 +to_meter=1
> data source : in memory
> names       : X2010.01.01, X2010.01.09, X2010.01.17, X2010.01.25,
> X2010.02.02, X2010.02.10, X2010.02.18, X2010.02.26, X2010.03.06,
> X2010.03.14, X2010.03.22, X2010.03.30, X2010.04.07, X2010.04.15,
> X2010.04.23, ...
> min values  :       13279,       13136,       13159, 13201,
> 13116,       13179,       13379,       13522, 13622,       13758,
> 13901,       13865,       14132, 14007,       13915, ...
> max values  :       15612,       15722,       15766, 15734,
> 15808,       16005,       16145,       16223, 16409,       16493,
> 16333,       16496,       16546, 16694,       16716, ...
> time        : 2010-01-03, 2010-12-27 (min, max)
>
>   > LSTWGm<-zApply(LSTWG,by=as.yearmon,fun=max)
>
>   > LSTWGm
> class       : RasterBrick
> dimensions  : 3334, 2294, 7648196, 12  (nrow, ncol, ncell, nlayers)
> resolution  : 1000.16, 999.8986  (x, y)
> extent      : -728812.9, 1565553, -1033.987, 3332628  (xmin, xmax, ymin,
> ymax)
> coord. ref. : +proj=utm +no_defs +zone=43 +a=6378137 +rf=298.257223563
> +towgs84=0.000,0.000,0.000 +to_meter=1
> data source :
> /tmp/R_raster_karthik/raster_tmp_2014-09-09_114723_2361_46029.grd
> names       : Jan.2010, Feb.2010, Mar.2010, Apr.2010, May.2010,
> Jun.2010, Jul.2010, Aug.2010, Sep.2010, Oct.2010, Nov.2010, Dec.2010
> min values  :    13279,    13522,    13901,    14132,    14360,
> 13865,    12649,    12048,    14142,    14031,    13868,    13667
> max values  :    15766,    16223,    16496,    16716,    16647,
> 16671,    16357,    16439,    16467,    16382,    15886,    15671
>               : Jan 2010, Feb 2010, Mar 2010, Apr 2010, May 2010, Jun
> 2010, Jul 2010, Aug 2010, Sep 2010, Oct 2010, Nov 2010, Dec 2010
>
>
> Regards
>



More information about the R-sig-Geo mailing list