Store time series data to big matrix

General
, ,
21

@Tazinho here are 3 reprex of the same code

22

Hi @Edward let's try to make your code a better Reprex first.

  1. As long as the code doesn't run in one run via the reprex function, let's workaround in the following way:
  • Close RStudio (don't save any data)
  • Restart R Studio (don't load any data from previous R session)
  • Run the code
  • Copy the code including the output here
23 
library(ncdf4)
library(bigmemory)
library(raster)
library(rgdal)
library(sp)

#Load the script needed to run do_pixel
source("~/Work/dopixel_2020.R")


#--load dates which are in R date format
ERA.dates <- readRDS("~/Work/ERA.dates.rds") #240 dates

#Load environmental variables

ts.evi <- readRDS("~/Work/ts.evi.rds")
ts.rel <- readRDS("~/Work/ts.rel.rds")
tile.xy <- readRDS("~/Work/tile.xy.rds")
ts.moist1 <- readRDS("~/Work/ts.moist1.rds")
ts.moist2 <- readRDS("~/Work/ts.moist2.rds")
ts.srad <- readRDS("~/Work/ts.srad.rds")
ts.tair <- readRDS("~/Work/ts.tair.rds")
ts.tsoil <- readRDS("~/Work/ts.tsoil.rds")


#------------------------------------------------------------------#
#--get date information from the MODIS time series-----------------#
#------------------------------------------------------------------#
# get the layer names provided by the call to extract evi:
layer.names<-labels(ts.evi)[[2]]   
# take the substring on these names, which give the dates:
date.text <-substr(layer.names,start=10,stop=16)
date.pos <- strptime(date.text,  format = "%Y%j", tz="UTC") #convert to R date format


#------------------------------------------------------------------#
#--create big matrix to store output of do_pixel-------------------#
#------------------------------------------------------------------#
pixel.date <- as.Date(seq.POSIXt(ISOdate(2000,6,1),by="week",length.out=1020,tz="UTC"))
#tail(pixel.date)
YEAR <- as.numeric(format(pixel.date, "%Y"))
YITS = length(unique(YEAR))
#--the column names for the big matrix
names<-c("lon","lat",
         "lq", "uq", "mean.evi", "sd.evi", "sum.evi", "amplitude",
         "peak.day", "trough.day",
         paste("td", 1:YITS, sep="."), paste("td.x", 1:YITS, sep="."), paste("td.evi", 1:YITS, sep="."),
         paste("pd", 1:YITS, sep="."), paste("pd.x", 1:YITS, sep="."), paste("pd.evi", 1:YITS, sep="."),
         paste("elon.m", 1:YITS, sep="."), paste("elon.m.x", 1:YITS, sep="."),
         paste("eoff.m", 1:YITS, sep="."), paste("eoff.m.x", 1:YITS, sep="."),
         paste("elon.f", 1:YITS, sep="."), paste("elon.f.x", 1:YITS, sep="."),
         paste("eoff.f", 1:YITS, sep="."), paste("eoff.f.x", 1:YITS, sep="."),
         paste("elon.l", 1:YITS, sep="."), paste("elon.l.x", 1:YITS, sep="."),
         paste("eoff.l", 1:YITS, sep="."), paste("eoff.l.x", 1:YITS, sep="."),
         paste("sum.evi.yr", 1:YITS, sep="."), paste("amp", 1:YITS, sep="."),
         paste("gsl", 1:YITS, sep="."), paste("gsl.peak", 1:YITS, sep="."), paste("gsl.long", 1:YITS, sep="."),
         paste("elon.f.evi", 1:YITS, sep="."), paste("elon.m.evi", 1:YITS, sep="."), paste("elon.l.evi", 1:YITS, sep="."),
         paste("eoff.f.evi", 1:YITS, sep="."), paste("eoff.m.evi", 1:YITS, sep="."), paste("eoff.l.evi", 1:YITS, sep="."),
         paste("tair.rank.td", 1:YITS, sep="."),paste("tair.rank.pd", 1:YITS, sep="."),
         paste("moist1.rank.td", 1:YITS, sep="."),paste("moist1.rank.pd", 1:YITS, sep=".")
) 
#--create the matrix using pixel.df to specifiy rows
#--and column names to speciify columns
statsMODIS<-big.matrix(nrow=length(tile.xy), ncol=length(names), type="double",
                       dimnames=list(row.names=NULL, col.names=names),
                       backingfile="statsMODIS.bin", descriptorfile="statsMODIS.desc")


#------------------------------------------------------------------#
#--start a for loop that will call do_pixel------------------------#
#------------------------------------------------------------------#
ERA.dates <- ERA.dates+14
pixel.df = data.frame(date = pixel.date)
for(p in 77910:77930)
  #for(p in 1:dim(ts.evi)[1]) #for all pixels in ts.evi
{
  print(p)
  #-0-only call dopixel if there is sufficient evi data in the pixel.
  if( length(which(is.na(ts.evi[p,])))  < 100 ) # there are  457 data points
  {
    #-1-create a data.frame for the pixel
    pixel.df = data.frame(date = pixel.date)
    #-2-create an approx function for each ERA variable
    af.tair = approxfun(ERA.dates,ts.tair[1,])
    af.tsoil = approxfun(ERA.dates,ts.tsoil[1,])
    af.srad = approxfun(ERA.dates,ts.srad[1,])
    af.moist1 = approxfun(ERA.dates,ts.moist1[1,])
    af.moist2 = approxfun(ERA.dates,ts.moist2[1,])
    #-3-create an approx function for the EVI data
    af.evi = approxfun(as.Date(date.pos),ts.evi[p,])
    #-4-use the approxfun's to fill the data frame data.frame
    pixel.df$tair <- af.tair(as.Date(pixel.date))
    pixel.df$tsoil <- af.tsoil(as.Date(pixel.date))
    pixel.df$srad <- af.srad(as.Date(pixel.date))
    pixel.df$moist1 <- af.moist1(as.Date(pixel.date))
    pixel.df$moist2 <- af.moist2(as.Date(pixel.date))
    pixel.df$evi <- af.evi(as.Date(pixel.date))
    pixel.df<-pixel.df[complete.cases(pixel.df),]
    #head(pixel.df)
    coords<-round(coordinates(tile.xy)[p,],2)
    statsMODIS[p,] <- dopixel(pixel.df,coords,PLOT=FALSE)
  } else 
  {
    coords<-round(coordinates(tile.xy)[p,],2)
    statsMODIS[p,] <- c(coords,rep(-9999,length=length(names)-2))
  }
}

console output:

[1] 77910
[1] 77911
[1] 77912
[1] 77913
[1] 77914
[1] 77915
[1] 77916
[1] 77917
[1] 77918
[1] 77919
[1] 77920
Error in SetRows.bm(x, i, value) :
number of items to replace is not a multiple of replacement length
In addition: There were 11 warnings (use warnings() to see them)
"#>"

24

Ok, now that we have a valid peace of code, we can carefully simplify it. We also need to find out where the error lies. If we know this, we can drop far more code to simplify our reprex again and so on.

  1. To start, you can change

    for(p in 77910:77930)

    into

    for(p in 77919:77920)

    and

    length(which(is.na(ts.evi[p,]))

    into

    sum(is.na(ts.evi[p,])

    Next, you should look in your code, which lines are not needed, or which data isn't needed. (Always check that the script produces the same error in the end).

  2. To see if the error lies in the if or the else part of your condition, please add

    print("true")

    in the first line of the if() part ...
    and

    print("false")

    in the first line of the else() part.

Afterwards, plesase post the simplified code together with the output (it would be good if you copy a "#>" in front of every line of the output; you can add the output to the markup. Pls also run warnings() in the end and copy it's output here.)

25 
for(p in 77919:77920)
  #for(p in 1:dim(ts.evi)[1]) #for all pixels in ts.evi, this will be slow!
{
  print(p)
  #-0-only call dopixel if there is sufficient evi data in the pixel.
  if( sum(is.na(ts.evi[p,]))  < 100 )  # there are circa 440 data points
  { print("true")
    #-1-create a data.frame for the pixel
    pixel.df = data.frame(date = pixel.date)
    #-2-create an approx function for each ERA variable
    af.tair = approxfun(ERA.dates,ts.tair[1,])
    af.tsoil = approxfun(ERA.dates,ts.tsoil[1,])
    af.srad = approxfun(ERA.dates,ts.srad[1,])
    af.moist1 = approxfun(ERA.dates,ts.moist1[1,])
    af.moist2 = approxfun(ERA.dates,ts.moist2[1,])
    #-3-create an approx function for the EVI data
    af.evi = approxfun(as.Date(date.pos),ts.evi[p,])
    #-4-use the approxfun's to fill the data frame data.frame
    pixel.df$tair <- af.tair(as.Date(pixel.date))
    pixel.df$tsoil <- af.tsoil(as.Date(pixel.date))
    pixel.df$srad <- af.srad(as.Date(pixel.date))
    pixel.df$moist1 <- af.moist1(as.Date(pixel.date))
    pixel.df$moist2 <- af.moist2(as.Date(pixel.date))
    pixel.df$evi <- af.evi(as.Date(pixel.date))
    pixel.df<-pixel.df[complete.cases(pixel.df),]
    #head(pixel.df)
    coords<-round(coordinates(tile.xy)[p,],2)
    statsMODIS[p,] <- dopixel(pixel.df,coords,PLOT=FALSE)
  }  else 
  { print("false")
    coords<-round(coordinates(tile.xy)[p,],2)
    statsMODIS[p,] <- c(coords,rep(-9999,length=length(names)-2))
  }
} 

warnings()

"#>" [1] 77919
"#>" [1] "true"
"#>" [1] 77920
"#>" [1] "true"
"#>" Error in SetRows.bm(x, i, value) :
"#>" number of items to replace is not a multiple of replacement length
"#>" In addition: Warning messages:
"#>" 1: In pd.evi - td.evi[-1] :
"#>" longer object length is not a multiple of shorter object length
"#>" 2: In pd.evi - td.evi[-1] :
"#>" longer object length is not a multiple of shorter object length
"#>" warnings()
"#>" Warning messages:
"#>" 1: In pd.evi - td.evi[-1] :
"#>" longer object length is not a multiple of shorter object length
"#>" 2: In pd.evi - td.evi[-1] :
"#>" longer object length is not a multiple of shorter object length

Is that what you mean?

26

Almost:

  1. You can add the output code within the markup (I mean within the ```. This will look nicer. The "#>" is basically just there, that someone else can copy from here directly to an editor. The ">" makes clear, that it's output and the "#" means that it's a comment. You don't need to add the quotes. Just #> is enough)
  2. From the print statements, we see that the code run within the else() statement is not evaluated in our for loop for our current reprex, so we can drop this code.
  3. By dropping code I don't mean that you should drop the code in the beginning. It is important that you include this code as well. Otherwise your minimal example is not self containing. Please readd this code.
  4. Instead I meant, that you should go carefully through your code and remove anything from it which is not necessary to reproduce the error (imagine you loaded a library which is not needed in the code or data that is not needed). You are the expert here and can make it much easier to let me (or others) help you here, if you see some unnecessary code (or data) which can be removed.
27 
library(bigmemory)


#Load the script needed to run do_pixel
source("~/Work/dopixel_2020.R")


#--load dates which are in R date format
ERA.dates <- readRDS("~/Work/ERA.dates.rds") #240 dates

#Load environmental variables

ts.evi <- readRDS("~/Work/ts.evi.rds")
tile.xy <- readRDS("~/Work/tile.xy.rds")
ts.moist1 <- readRDS("~/Work/ts.moist1.rds")
ts.moist2 <- readRDS("~/Work/ts.moist2.rds")
ts.srad <- readRDS("~/Work/ts.srad.rds")
ts.tair <- readRDS("~/Work/ts.tair.rds")
ts.tsoil <- readRDS("~/Work/ts.tsoil.rds")


#------------------------------------------------------------------#
#--get date information from the MODIS time series-----------------#
#------------------------------------------------------------------#
# get the layer names provided by the call to extract evi:
layer.names<-labels(ts.evi)[[2]]   
# take the substring on these names, which give the dates:
date.text <-substr(layer.names,start=10,stop=16)
date.pos <- strptime(date.text,  format = "%Y%j", tz="UTC") #convert to R date format


#------------------------------------------------------------------#
#--create big matrix to store output of do_pixel-------------------#
#------------------------------------------------------------------#
pixel.date <- as.Date(seq.POSIXt(ISOdate(2000,6,1),by="week",length.out=1020,tz="UTC"))
#tail(pixel.date)
YEAR <- as.numeric(format(pixel.date, "%Y"))
YITS = length(unique(YEAR))
#--the column names for the big matrix
names<-c("lon","lat",
         "lq", "uq", "mean.evi", "sd.evi", "sum.evi", "amplitude",
         "peak.day", "trough.day",
         paste("td", 1:YITS, sep="."), paste("td.x", 1:YITS, sep="."), paste("td.evi", 1:YITS, sep="."),
         paste("pd", 1:YITS, sep="."), paste("pd.x", 1:YITS, sep="."), paste("pd.evi", 1:YITS, sep="."),
         paste("elon.m", 1:YITS, sep="."), paste("elon.m.x", 1:YITS, sep="."),
         paste("eoff.m", 1:YITS, sep="."), paste("eoff.m.x", 1:YITS, sep="."),
         paste("elon.f", 1:YITS, sep="."), paste("elon.f.x", 1:YITS, sep="."),
         paste("eoff.f", 1:YITS, sep="."), paste("eoff.f.x", 1:YITS, sep="."),
         paste("elon.l", 1:YITS, sep="."), paste("elon.l.x", 1:YITS, sep="."),
         paste("eoff.l", 1:YITS, sep="."), paste("eoff.l.x", 1:YITS, sep="."),
         paste("sum.evi.yr", 1:YITS, sep="."), paste("amp", 1:YITS, sep="."),
         paste("gsl", 1:YITS, sep="."), paste("gsl.peak", 1:YITS, sep="."), paste("gsl.long", 1:YITS, sep="."),
         paste("elon.f.evi", 1:YITS, sep="."), paste("elon.m.evi", 1:YITS, sep="."), paste("elon.l.evi", 1:YITS, sep="."),
         paste("eoff.f.evi", 1:YITS, sep="."), paste("eoff.m.evi", 1:YITS, sep="."), paste("eoff.l.evi", 1:YITS, sep="."),
         paste("tair.rank.td", 1:YITS, sep="."),paste("tair.rank.pd", 1:YITS, sep="."),
         paste("moist1.rank.td", 1:YITS, sep="."),paste("moist1.rank.pd", 1:YITS, sep=".")
) 
#--create the matrix using pixel.df to specifiy rows
#--and column names to speciify columns
statsMODIS<-big.matrix(nrow=length(tile.xy), ncol=length(names), type="double",
                       dimnames=list(row.names=NULL, col.names=names),
                       backingfile="statsMODIS.bin", descriptorfile="statsMODIS.desc")


#------------------------------------------------------------------#
#--start a for loop that will call do_pixel------------------------#
#--this will be slow, but its useful for debugging-----------------#
#------------------------------------------------------------------#
ERA.dates <- ERA.dates+14
pixel.df = data.frame(date = pixel.date)
for(p in 77919:77920)
  #for(p in 1:dim(ts.evi)[1]) #for all pixels in ts.evi, this will be slow!
{
  print(p)
  #-0-only call dopixel if there is sufficient evi data in the pixel.
  if( sum(is.na(ts.evi[p,]))  < 100 )  # there are circa 440 data points
  { print("true")
    #-1-create a data.frame for the pixel
    pixel.df = data.frame(date = pixel.date)
    #-2-create an approx function for each ERA variable
    af.tair = approxfun(ERA.dates,ts.tair[1,])
    af.tsoil = approxfun(ERA.dates,ts.tsoil[1,])
    af.srad = approxfun(ERA.dates,ts.srad[1,])
    af.moist1 = approxfun(ERA.dates,ts.moist1[1,])
    af.moist2 = approxfun(ERA.dates,ts.moist2[1,])
    #-3-create an approx function for the EVI data
    af.evi = approxfun(as.Date(date.pos),ts.evi[p,])
    #-4-use the approxfun's to fill the data frame data.frame
    pixel.df$tair <- af.tair(as.Date(pixel.date))
    pixel.df$tsoil <- af.tsoil(as.Date(pixel.date))
    pixel.df$srad <- af.srad(as.Date(pixel.date))
    pixel.df$moist1 <- af.moist1(as.Date(pixel.date))
    pixel.df$moist2 <- af.moist2(as.Date(pixel.date))
    pixel.df$evi <- af.evi(as.Date(pixel.date))
    pixel.df<-pixel.df[complete.cases(pixel.df),]
    #head(pixel.df)
    coords<-round(coordinates(tile.xy)[p,],2)
    statsMODIS[p,] <- dopixel(pixel.df,coords,PLOT=FALSE)
  }  else 
  { print("false")
    coords<-round(coordinates(tile.xy)[p,],2)
    statsMODIS[p,] <- c(coords,rep(-9999,length=length(names)-2))
  }
} 

warnings()

I have thus removed unneeded libraries and objects.
Below is the output.

#> [1] 77919
#> [1] "true"
#> [1] 77920
#> [1] "true"
#> Error in SetRows.bm(x, i, value) : 
#> number of items to replace is not a multiple of replacement length
#> In addition: Warning messages:
#> 1: In pd.evi - td.evi[-1] :
#>  longer object length is not a multiple of shorter object length
#> 2: In pd.evi - td.evi[-1] :
#> longer object length is not a multiple of shorter object length
#> warnings()
#> Warning messages:
#> 1: In pd.evi - td.evi[-1] :
#>  longer object length is not a multiple of shorter object length
#> 2: In pd.evi - td.evi[-1] :
#>  longer object length is not a multiple of shorter object length
28

Ok.

  1. Please from now on ensure that you stick to this process. 1) Restart RStudio 2) Run the code 3) Post the code. (BTW: You don't need to split your input and output. Just post it in one chunk).

  2. You can remove the line

        else 
  { print("false")
    coords<-round(coordinates(tile.xy)[p,],2)
    statsMODIS[p,] <- c(coords,rep(-9999,length=length(names)-2))
  }

    because it's never evaluated.

  3. You can also drop the code for the if() condition, i.e.

    #-0-only call dopixel if there is sufficient evi data in the pixel.
  if( sum(is.na(ts.evi[p,]))  < 100 )  # there are circa 440 data points
  { print("true")

    as it's always true in the context of our two iterations.

  4. Next is a new step. We want to find the error. Before the line where you believe the error appears, please add a line print("before") and after the error line add a line print("after"). If you found the error, and added these lines, please post your code again.

29 
#Load the script needed to run do_pixel
source("~/Work/dopixel_2020.R")


#--load dates which are in R date format
ERA.dates <- readRDS("~/Work/ERA.dates.rds") #240 dates

#Load environmental variables

ts.evi <- readRDS("~/Work/ts.evi.rds")
tile.xy <- readRDS("~/Work/tile.xy.rds")
ts.moist1 <- readRDS("~/Work/ts.moist1.rds")
ts.moist2 <- readRDS("~/Work/ts.moist2.rds")
ts.srad <- readRDS("~/Work/ts.srad.rds")
ts.tair <- readRDS("~/Work/ts.tair.rds")
ts.tsoil <- readRDS("~/Work/ts.tsoil.rds")


#------------------------------------------------------------------#
#--get date information from the MODIS time series-----------------#
#------------------------------------------------------------------#
# get the layer names provided by the call to extract evi:
layer.names<-labels(ts.evi)[[2]]   
# take the substring on these names, which give the dates:
date.text <-substr(layer.names,start=10,stop=16)
date.pos <- strptime(date.text,  format = "%Y%j", tz="UTC") #convert to R date format


#------------------------------------------------------------------#
#--create big matrix to store output of do_pixel-------------------#
#------------------------------------------------------------------#
pixel.date <- as.Date(seq.POSIXt(ISOdate(2000,6,1),by="week",length.out=1020,tz="UTC"))
#tail(pixel.date)
YEAR <- as.numeric(format(pixel.date, "%Y"))
YITS = length(unique(YEAR))
#--the column names for the big matrix
names<-c("lon","lat",
         "lq", "uq", "mean.evi", "sd.evi", "sum.evi", "amplitude",
         "peak.day", "trough.day",
         paste("td", 1:YITS, sep="."), paste("td.x", 1:YITS, sep="."), paste("td.evi", 1:YITS, sep="."),
         paste("pd", 1:YITS, sep="."), paste("pd.x", 1:YITS, sep="."), paste("pd.evi", 1:YITS, sep="."),
         paste("elon.m", 1:YITS, sep="."), paste("elon.m.x", 1:YITS, sep="."),
         paste("eoff.m", 1:YITS, sep="."), paste("eoff.m.x", 1:YITS, sep="."),
         paste("elon.f", 1:YITS, sep="."), paste("elon.f.x", 1:YITS, sep="."),
         paste("eoff.f", 1:YITS, sep="."), paste("eoff.f.x", 1:YITS, sep="."),
         paste("elon.l", 1:YITS, sep="."), paste("elon.l.x", 1:YITS, sep="."),
         paste("eoff.l", 1:YITS, sep="."), paste("eoff.l.x", 1:YITS, sep="."),
         paste("sum.evi.yr", 1:YITS, sep="."), paste("amp", 1:YITS, sep="."),
         paste("gsl", 1:YITS, sep="."), paste("gsl.peak", 1:YITS, sep="."), paste("gsl.long", 1:YITS, sep="."),
         paste("elon.f.evi", 1:YITS, sep="."), paste("elon.m.evi", 1:YITS, sep="."), paste("elon.l.evi", 1:YITS, sep="."),
         paste("eoff.f.evi", 1:YITS, sep="."), paste("eoff.m.evi", 1:YITS, sep="."), paste("eoff.l.evi", 1:YITS, sep="."),
         paste("tair.rank.td", 1:YITS, sep="."),paste("tair.rank.pd", 1:YITS, sep="."),
         paste("moist1.rank.td", 1:YITS, sep="."),paste("moist1.rank.pd", 1:YITS, sep=".")
) 
#--create the matrix using pixel.df to specifiy rows
#--and column names to speciify columns
statsMODIS<-big.matrix(nrow=length(tile.xy), ncol=length(names), type="double",
                       dimnames=list(row.names=NULL, col.names=names),
                       backingfile="statsMODIS.bin", descriptorfile="statsMODIS.desc")


#------------------------------------------------------------------#
#--start a for loop that will call do_pixel------------------------#
#--this will be slow, but its useful for debugging-----------------#
#------------------------------------------------------------------#
ERA.dates <- ERA.dates+14


for(p in 77919:77920){
  
  
  print(p)
  
  #-1-create a data.frame for the pixel
  pixel.df = data.frame(date = pixel.date)
  #-2-create an approx function for each ERA variable
  af.tair = approxfun(ERA.dates,ts.tair[1,])
  af.tsoil = approxfun(ERA.dates,ts.tsoil[1,])
  af.srad = approxfun(ERA.dates,ts.srad[1,])
  af.moist1 = approxfun(ERA.dates,ts.moist1[1,])
  af.moist2 = approxfun(ERA.dates,ts.moist2[1,])
  #-3-create an approx function for the EVI data
  af.evi = approxfun(as.Date(date.pos),ts.evi[p,])
  #-4-use the approxfun's to fill the data frame data.frame
  pixel.df$tair <- af.tair(as.Date(pixel.date))
  pixel.df$tsoil <- af.tsoil(as.Date(pixel.date))
  pixel.df$srad <- af.srad(as.Date(pixel.date))
  pixel.df$moist1 <- af.moist1(as.Date(pixel.date))
  pixel.df$moist2 <- af.moist2(as.Date(pixel.date))
  pixel.df$evi <- af.evi(as.Date(pixel.date))
  pixel.df<-pixel.df[complete.cases(pixel.df),]
  
  coords<-round(coordinates(tile.xy)[p,],2)
  
  print("before")
  
  statsMODIS[p,] <- dopixel(pixel.df,coords,PLOT=FALSE)
  
  print("after")
  
} 

warnings()

#> [1] 77919
#> [1] "before"
#> [1] "after"
#> [1] 77920
#> [1] "before"
#> Error in SetRows.bm(x, i, value) : 
#>  number of items to replace is not a multiple of replacement length
#> In addition: Warning messages:
#> 1: In pd.evi - td.evi[-1] :
#>  longer object length is not a multiple of shorter object length
#> 2: In pd.evi - td.evi[-1] :
#>  longer object length is not a multiple of shorter object length

#> warnings()
#> Warning messages:
#> 1: In pd.evi - td.evi[-1] :
#>  longer object length is not a multiple of shorter object length
#> 2: In pd.evi - td.evi[-1] :
#>  longer object length is not a multiple of shorter object length
30

Ok. Can you add the following lines under the line print("before"):

print(dim(pixel.df))
print(dim(coords))
print(dim(statsMODIS[p,]))
31 
for(p in 77919:77920){
  
  
  print(p)
  
  #-1-create a data.frame for the pixel
  pixel.df = data.frame(date = pixel.date)
  #-2-create an approx function for each ERA variable
  af.tair = approxfun(ERA.dates,ts.tair[1,])
  af.tsoil = approxfun(ERA.dates,ts.tsoil[1,])
  af.srad = approxfun(ERA.dates,ts.srad[1,])
  af.moist1 = approxfun(ERA.dates,ts.moist1[1,])
  af.moist2 = approxfun(ERA.dates,ts.moist2[1,])
  #-3-create an approx function for the EVI data
  af.evi = approxfun(as.Date(date.pos),ts.evi[p,])
  #-4-use the approxfun's to fill the data frame data.frame
  pixel.df$tair <- af.tair(as.Date(pixel.date))
  pixel.df$tsoil <- af.tsoil(as.Date(pixel.date))
  pixel.df$srad <- af.srad(as.Date(pixel.date))
  pixel.df$moist1 <- af.moist1(as.Date(pixel.date))
  pixel.df$moist2 <- af.moist2(as.Date(pixel.date))
  pixel.df$evi <- af.evi(as.Date(pixel.date))
  pixel.df<-pixel.df[complete.cases(pixel.df),]
  
  coords<-round(coordinates(tile.xy)[p,],2)
  
  print("before")
  
  print(dim(pixel.df))
  print(dim(coords))
  print(dim(statsMODIS[p,]))
  
  statsMODIS[p,] <- dopixel(pixel.df,coords,PLOT=FALSE)
  
  print("after")
  
} 

warnings()

#> [1] 77919
#> [1] "before"
#> [1] 989   7
#> NULL
#> NULL
#> [1] "after"
#> [1] 77920
#> [1] "before"
#> [1] 970   7
#> NULL
#> NULL
#> Error in SetRows.bm(x, i, value) : 
#>  number of items to replace is not a multiple of replacement length
#> In addition: Warning messages:
#> 1: In pd.evi - td.evi[-1] :
#>  longer object length is not a multiple of shorter object length
#> 2: In pd.evi - td.evi[-1] :
#>  longer object length is not a multiple of shorter object length
32

Ok. Please always include the whole code. This makes it much easier to reason about it and I don't have to lookup things from the top of this thread.

It would be good if we can see what coords and statsMODIS[p, ] are. Can you replace dim with typeof() and class() and head() and post the results (if it doesn't get too long)?

Then let's look at it on another day (or hope someone else jumps in to the discussion).

If you see any way to reduce the amount of data further (maybe you don't need all the columns in the data), it would really help.

33 
library(bigmemory)


#Load the script needed to run do_pixel
source("~/Work/dopixel_2020.R")


#--load dates which are in R date format
ERA.dates <- readRDS("~/Work/ERA.dates.rds") #240 dates

#Load environmental variables

ts.evi <- readRDS("~/Work/ts.evi.rds")
tile.xy <- readRDS("~/Work/tile.xy.rds")
ts.moist1 <- readRDS("~/Work/ts.moist1.rds")
ts.moist2 <- readRDS("~/Work/ts.moist2.rds")
ts.srad <- readRDS("~/Work/ts.srad.rds")
ts.tair <- readRDS("~/Work/ts.tair.rds")
ts.tsoil <- readRDS("~/Work/ts.tsoil.rds")


#------------------------------------------------------------------#
#--get date information from the MODIS time series-----------------#
#------------------------------------------------------------------#
# get the layer names provided by the call to extract evi:
layer.names<-labels(ts.evi)[[2]]   
# take the substring on these names, which give the dates:
date.text <-substr(layer.names,start=10,stop=16)
date.pos <- strptime(date.text,  format = "%Y%j", tz="UTC") #convert to R date format


#------------------------------------------------------------------#
#--create big matrix to store output of do_pixel-------------------#
#------------------------------------------------------------------#
pixel.date <- as.Date(seq.POSIXt(ISOdate(2000,6,1),by="week",length.out=1020,tz="UTC"))
#tail(pixel.date)
YEAR <- as.numeric(format(pixel.date, "%Y"))
YITS = length(unique(YEAR))
#--the column names for the big matrix
names<-c("lon","lat",
         "lq", "uq", "mean.evi", "sd.evi", "sum.evi", "amplitude",
         "peak.day", "trough.day",
         paste("td", 1:YITS, sep="."), paste("td.x", 1:YITS, sep="."), paste("td.evi", 1:YITS, sep="."),
         paste("pd", 1:YITS, sep="."), paste("pd.x", 1:YITS, sep="."), paste("pd.evi", 1:YITS, sep="."),
         paste("elon.m", 1:YITS, sep="."), paste("elon.m.x", 1:YITS, sep="."),
         paste("eoff.m", 1:YITS, sep="."), paste("eoff.m.x", 1:YITS, sep="."),
         paste("elon.f", 1:YITS, sep="."), paste("elon.f.x", 1:YITS, sep="."),
         paste("eoff.f", 1:YITS, sep="."), paste("eoff.f.x", 1:YITS, sep="."),
         paste("elon.l", 1:YITS, sep="."), paste("elon.l.x", 1:YITS, sep="."),
         paste("eoff.l", 1:YITS, sep="."), paste("eoff.l.x", 1:YITS, sep="."),
         paste("sum.evi.yr", 1:YITS, sep="."), paste("amp", 1:YITS, sep="."),
         paste("gsl", 1:YITS, sep="."), paste("gsl.peak", 1:YITS, sep="."), paste("gsl.long", 1:YITS, sep="."),
         paste("elon.f.evi", 1:YITS, sep="."), paste("elon.m.evi", 1:YITS, sep="."), paste("elon.l.evi", 1:YITS, sep="."),
         paste("eoff.f.evi", 1:YITS, sep="."), paste("eoff.m.evi", 1:YITS, sep="."), paste("eoff.l.evi", 1:YITS, sep="."),
         paste("tair.rank.td", 1:YITS, sep="."),paste("tair.rank.pd", 1:YITS, sep="."),
         paste("moist1.rank.td", 1:YITS, sep="."),paste("moist1.rank.pd", 1:YITS, sep=".")
) 
#--create the matrix using pixel.df to specifiy rows
#--and column names to speciify columns
statsMODIS<-big.matrix(nrow=length(tile.xy), ncol=length(names), type="double",
                       dimnames=list(row.names=NULL, col.names=names),
                       backingfile="statsMODIS.bin", descriptorfile="statsMODIS.desc")


#------------------------------------------------------------------#
#--start a for loop that will call do_pixel------------------------#
#--this will be slow, but its useful for debugging-----------------#
#------------------------------------------------------------------#
ERA.dates <- ERA.dates+14


for(p in 77919:77920){
  
  
  print(p)
  
  #-1-create a data.frame for the pixel
  pixel.df = data.frame(date = pixel.date)
  #-2-create an approx function for each ERA variable
  af.tair = approxfun(ERA.dates,ts.tair[1,])
  af.tsoil = approxfun(ERA.dates,ts.tsoil[1,])
  af.srad = approxfun(ERA.dates,ts.srad[1,])
  af.moist1 = approxfun(ERA.dates,ts.moist1[1,])
  af.moist2 = approxfun(ERA.dates,ts.moist2[1,])
  #-3-create an approx function for the EVI data
  af.evi = approxfun(as.Date(date.pos),ts.evi[p,])
  #-4-use the approxfun's to fill the data frame data.frame
  pixel.df$tair <- af.tair(as.Date(pixel.date))
  pixel.df$tsoil <- af.tsoil(as.Date(pixel.date))
  pixel.df$srad <- af.srad(as.Date(pixel.date))
  pixel.df$moist1 <- af.moist1(as.Date(pixel.date))
  pixel.df$moist2 <- af.moist2(as.Date(pixel.date))
  pixel.df$evi <- af.evi(as.Date(pixel.date))
  pixel.df<-pixel.df[complete.cases(pixel.df),]
  
  coords<-round(coordinates(tile.xy)[p,],2)
  
  print("before")
  
  print(typeof(pixel.df))
  print(class(coords))
  print(head(statsMODIS[p,]))

  
  statsMODIS[p,] <- dopixel(pixel.df,coords,PLOT=FALSE)
  
  print("after")
  
} 

#> warnings()

#> [1] 77919
#> [1] "before"
#> [1] "list"
#> [1] "numeric"
#>         lon          lat           lq           uq     mean.evi       sd.evi 
#> 15.62000000 -18.16000000   0.11379423   0.25225500   0.16936279   0.05381537 
#> [1] "after"
#> [1] 77920
#> [1] "before"
#> [1] "list"
#> [1] "numeric"
#>   lon      lat       lq       uq mean.evi   sd.evi 
#>      0        0        0        0        0        0 
#> Error in SetRows.bm(x, i, value) : 
#>  number of items to replace is not a multiple of replacement length
#> In addition: Warning messages:
#> 1: In pd.evi - td.evi[-1] :
#>  longer object length is not a multiple of shorter object length
#> 2: In pd.evi - td.evi[-1] :
#>  longer object length is not a multiple of shorter object length

Thank you for making time to have a look at the code. Highly appreciated. Unfortunately, I need all the columns in the matrix.

34

Ok. Let's make a last try.

  1. I removed the for-loop and just wrote out the code for both iterations.
  2. In this way, we should get warnings (almost) immediately without running warnings() afterwards. (In case you want to debug warnings later you can always set options(warn = 2) which will convert them like errors; options(warn = 0) will get you back to the normal behaviour.)
  1. I added print statements to give us information about statsMODIS[p,] and dopixel(pixel.df,coords,PLOT=FALSE). It looks like in the dimensions of the object are correct in the first iteration but not in the second one. Maybe this gives us some good hint about what to change to make the second iteration work.

Please run the following code and post the results here:

library(bigmemory)


#Load the script needed to run do_pixel
source("~/Work/dopixel_2020.R")


#--load dates which are in R date format
ERA.dates <- readRDS("~/Work/ERA.dates.rds") #240 dates

#Load environmental variables

ts.evi <- readRDS("~/Work/ts.evi.rds")
tile.xy <- readRDS("~/Work/tile.xy.rds")
ts.moist1 <- readRDS("~/Work/ts.moist1.rds")
ts.moist2 <- readRDS("~/Work/ts.moist2.rds")
ts.srad <- readRDS("~/Work/ts.srad.rds")
ts.tair <- readRDS("~/Work/ts.tair.rds")
ts.tsoil <- readRDS("~/Work/ts.tsoil.rds")


#------------------------------------------------------------------#
#--get date information from the MODIS time series-----------------#
#------------------------------------------------------------------#
# get the layer names provided by the call to extract evi:
layer.names<-labels(ts.evi)[[2]]   
# take the substring on these names, which give the dates:
date.text <-substr(layer.names,start=10,stop=16)
date.pos <- strptime(date.text,  format = "%Y%j", tz="UTC") #convert to R date format


#------------------------------------------------------------------#
#--create big matrix to store output of do_pixel-------------------#
#------------------------------------------------------------------#
pixel.date <- as.Date(seq.POSIXt(ISOdate(2000,6,1),by="week",length.out=1020,tz="UTC"))
#tail(pixel.date)
YEAR <- as.numeric(format(pixel.date, "%Y"))
YITS = length(unique(YEAR))
#--the column names for the big matrix
names<-c("lon","lat",
         "lq", "uq", "mean.evi", "sd.evi", "sum.evi", "amplitude",
         "peak.day", "trough.day",
         paste("td", 1:YITS, sep="."), paste("td.x", 1:YITS, sep="."), paste("td.evi", 1:YITS, sep="."),
         paste("pd", 1:YITS, sep="."), paste("pd.x", 1:YITS, sep="."), paste("pd.evi", 1:YITS, sep="."),
         paste("elon.m", 1:YITS, sep="."), paste("elon.m.x", 1:YITS, sep="."),
         paste("eoff.m", 1:YITS, sep="."), paste("eoff.m.x", 1:YITS, sep="."),
         paste("elon.f", 1:YITS, sep="."), paste("elon.f.x", 1:YITS, sep="."),
         paste("eoff.f", 1:YITS, sep="."), paste("eoff.f.x", 1:YITS, sep="."),
         paste("elon.l", 1:YITS, sep="."), paste("elon.l.x", 1:YITS, sep="."),
         paste("eoff.l", 1:YITS, sep="."), paste("eoff.l.x", 1:YITS, sep="."),
         paste("sum.evi.yr", 1:YITS, sep="."), paste("amp", 1:YITS, sep="."),
         paste("gsl", 1:YITS, sep="."), paste("gsl.peak", 1:YITS, sep="."), paste("gsl.long", 1:YITS, sep="."),
         paste("elon.f.evi", 1:YITS, sep="."), paste("elon.m.evi", 1:YITS, sep="."), paste("elon.l.evi", 1:YITS, sep="."),
         paste("eoff.f.evi", 1:YITS, sep="."), paste("eoff.m.evi", 1:YITS, sep="."), paste("eoff.l.evi", 1:YITS, sep="."),
         paste("tair.rank.td", 1:YITS, sep="."),paste("tair.rank.pd", 1:YITS, sep="."),
         paste("moist1.rank.td", 1:YITS, sep="."),paste("moist1.rank.pd", 1:YITS, sep=".")
) 
#--create the matrix using pixel.df to specifiy rows
#--and column names to speciify columns
statsMODIS<-big.matrix(nrow=length(tile.xy), ncol=length(names), type="double",
                       dimnames=list(row.names=NULL, col.names=names),
                       backingfile="statsMODIS.bin", descriptorfile="statsMODIS.desc")


#------------------------------------------------------------------#
#--start a for loop that will call do_pixel------------------------#
#--this will be slow, but its useful for debugging-----------------#
#------------------------------------------------------------------#
ERA.dates <- ERA.dates+14

### 1st iteration

p <- 77919

print(p)

#-1-create a data.frame for the pixel
pixel.df = data.frame(date = pixel.date)
#-2-create an approx function for each ERA variable
af.tair = approxfun(ERA.dates,ts.tair[1,])
af.tsoil = approxfun(ERA.dates,ts.tsoil[1,])
af.srad = approxfun(ERA.dates,ts.srad[1,])
af.moist1 = approxfun(ERA.dates,ts.moist1[1,])
af.moist2 = approxfun(ERA.dates,ts.moist2[1,])
#-3-create an approx function for the EVI data
af.evi = approxfun(as.Date(date.pos),ts.evi[p,])
#-4-use the approxfun's to fill the data frame data.frame
pixel.df$tair <- af.tair(as.Date(pixel.date))
pixel.df$tsoil <- af.tsoil(as.Date(pixel.date))
pixel.df$srad <- af.srad(as.Date(pixel.date))
pixel.df$moist1 <- af.moist1(as.Date(pixel.date))
pixel.df$moist2 <- af.moist2(as.Date(pixel.date))
pixel.df$evi <- af.evi(as.Date(pixel.date))
pixel.df<-pixel.df[complete.cases(pixel.df),]

coords<-round(coordinates(tile.xy)[p,],2)

print("before")

print(paste("statsModis_t:"), typeof(statsMODIS[p,]))
print(paste("statsModis_c:"), class(statsMODIS[p,]))
print(paste("statsModis_d:"), dim(statsMODIS[p,]))
print(paste("statsModis_l:"), length(statsMODIS[p,]))
print(paste("statsModis_h:"), head(statsMODIS[p,]))

print(paste("dopixel_t:"), typeof(dopixel(pixel.df,coords,PLOT=FALSE)))
print(paste("dopixel_c:"), class(dopixel(pixel.df,coords,PLOT=FALSE)))
print(paste("dopixel_d:"), dim(dopixel(pixel.df,coords,PLOT=FALSE)))
print(paste("dopixel_l:"), length(dopixel(pixel.df,coords,PLOT=FALSE)))
print(paste("dopixel_h:"), head(dopixel(pixel.df,coords,PLOT=FALSE)))

statsMODIS[p,] <- dopixel(pixel.df,coords,PLOT=FALSE)

print("after")

### 2nd iteration

p <- 77919

print(p)

#-1-create a data.frame for the pixel
pixel.df = data.frame(date = pixel.date)
#-2-create an approx function for each ERA variable
af.tair = approxfun(ERA.dates,ts.tair[1,])
af.tsoil = approxfun(ERA.dates,ts.tsoil[1,])
af.srad = approxfun(ERA.dates,ts.srad[1,])
af.moist1 = approxfun(ERA.dates,ts.moist1[1,])
af.moist2 = approxfun(ERA.dates,ts.moist2[1,])
#-3-create an approx function for the EVI data
af.evi = approxfun(as.Date(date.pos),ts.evi[p,])
#-4-use the approxfun's to fill the data frame data.frame
pixel.df$tair <- af.tair(as.Date(pixel.date))
pixel.df$tsoil <- af.tsoil(as.Date(pixel.date))
pixel.df$srad <- af.srad(as.Date(pixel.date))
pixel.df$moist1 <- af.moist1(as.Date(pixel.date))
pixel.df$moist2 <- af.moist2(as.Date(pixel.date))
pixel.df$evi <- af.evi(as.Date(pixel.date))
pixel.df<-pixel.df[complete.cases(pixel.df),]

coords<-round(coordinates(tile.xy)[p,],2)

print("before")

print(paste("statsModis_t:"), typeof(statsMODIS[p,]))
print(paste("statsModis_c:"), class(statsMODIS[p,]))
print(paste("statsModis_d:"), dim(statsMODIS[p,]))
print(paste("statsModis_l:"), length(statsMODIS[p,]))
print(paste("statsModis_h:"), head(statsMODIS[p,]))

print(paste("dopixel_t:"), typeof(dopixel(pixel.df,coords,PLOT=FALSE)))
print(paste("dopixel_c:"), class(dopixel(pixel.df,coords,PLOT=FALSE)))
print(paste("dopixel_d:"), dim(dopixel(pixel.df,coords,PLOT=FALSE)))
print(paste("dopixel_l:"), length(dopixel(pixel.df,coords,PLOT=FALSE)))
print(paste("dopixel_h:"), head(dopixel(pixel.df,coords,PLOT=FALSE)))

statsMODIS[p,] <- dopixel(pixel.df,coords,PLOT=FALSE)

print("after")
35 
#> library(bigmemory)
#> 
#> 
#> #Load the script needed to run do_pixel
#> source("~/Work/dopixel_2020.R")
#> 
#> 
#> #--load dates which are in R date format
#> ERA.dates <- readRDS("~/Work/ERA.dates.rds") #240 dates
#> 
#> #Load environmental variables
#> 
#> ts.evi <- readRDS("~/Work/ts.evi.rds")
#> tile.xy <- readRDS("~/Work/tile.xy.rds")
#> ts.moist1 <- readRDS("~/Work/ts.moist1.rds")
#> ts.moist2 <- readRDS("~/Work/ts.moist2.rds")
#> ts.srad <- readRDS("~/Work/ts.srad.rds")
#> ts.tair <- readRDS("~/Work/ts.tair.rds")
#> ts.tsoil <- readRDS("~/Work/ts.tsoil.rds")
#> 
#> 
#> #------------------------------------------------------------------#
#> #--get date information from the MODIS time series-----------------#
#> #------------------------------------------------------------------#
#> # get the layer names provided by the call to extract evi:
#> layer.names<-labels(ts.evi)[[2]]   
#> # take the substring on these names, which give the dates:
#> date.text <-substr(layer.names,start=10,stop=16)
#> date.pos <- strptime(date.text,  format = "%Y%j", tz="UTC") #convert to R date format
#> 
#> 
#> #------------------------------------------------------------------#
#> #--create big matrix to store output of do_pixel-------------------#
#> #------------------------------------------------------------------#
#> pixel.date <- as.Date(seq.POSIXt(ISOdate(2000,6,1),by="week",length.out=1020,tz="UTC"))
#> #tail(pixel.date)
#> YEAR <- as.numeric(format(pixel.date, "%Y"))
#> YITS = length(unique(YEAR))
#> #--the column names for the big matrix
#> names<-c("lon","lat",
#>          "lq", "uq", "mean.evi", "sd.evi", "sum.evi", "amplitude",
#>          "peak.day", "trough.day",
#>          paste("td", 1:YITS, sep="."), paste("td.x", 1:YITS, sep="."), paste("td.evi", 1:YITS, sep="."),
#>          paste("pd", 1:YITS, sep="."), paste("pd.x", 1:YITS, sep="."), paste("pd.evi", 1:YITS, sep="."),
#>          paste("elon.m", 1:YITS, sep="."), paste("elon.m.x", 1:YITS, sep="."),
#>          paste("eoff.m", 1:YITS, sep="."), paste("eoff.m.x", 1:YITS, sep="."),
#>          paste("elon.f", 1:YITS, sep="."), paste("elon.f.x", 1:YITS, sep="."),
#>          paste("eoff.f", 1:YITS, sep="."), paste("eoff.f.x", 1:YITS, sep="."),
#>          paste("elon.l", 1:YITS, sep="."), paste("elon.l.x", 1:YITS, sep="."),
#>          paste("eoff.l", 1:YITS, sep="."), paste("eoff.l.x", 1:YITS, sep="."),
#>          paste("sum.evi.yr", 1:YITS, sep="."), paste("amp", 1:YITS, sep="."),
#>          paste("gsl", 1:YITS, sep="."), paste("gsl.peak", 1:YITS, sep="."), paste("gsl.long", 1:YITS, sep="."),
#>          paste("elon.f.evi", 1:YITS, sep="."), paste("elon.m.evi", 1:YITS, sep="."), paste("elon.l.evi", 1:YITS, sep="."),
#>          paste("eoff.f.evi", 1:YITS, sep="."), paste("eoff.m.evi", 1:YITS, sep="."), paste("eoff.l.evi", 1:YITS, sep="."),
#>          paste("tair.rank.td", 1:YITS, sep="."),paste("tair.rank.pd", 1:YITS, sep="."),
#>         paste("moist1.rank.td", 1:YITS, sep="."),paste("moist1.rank.pd", 1:YITS, sep=".")
#> ) 
#> #--create the matrix using pixel.df to specifiy rows
#> #--and column names to speciify columns
#> statsMODIS<-big.matrix(nrow=length(tile.xy), ncol=length(names), type="double",
#>                        dimnames=list(row.names=NULL, col.names=names),
#>                       backingfile="statsMODIS.bin", descriptorfile="statsMODIS.desc")
#> Loading required package: sp
#> 
#> 
#> #------------------------------------------------------------------#
#> #--start a for loop that will call do_pixel------------------------#
#> #--this will be slow, but its useful for debugging-----------------#
#> #------------------------------------------------------------------#
#> ERA.dates <- ERA.dates+14
#> 
#> ### 1st iteration
#> 
#> p <- 77919
#> 
#> print(p)
#> [1] 77919
#> 
#>  #-1-create a data.frame for the pixel
#>  pixel.df = data.frame(date = pixel.date)
#>  #-2-create an approx function for each ERA variable
#>  af.tair = approxfun(ERA.dates,ts.tair[1,])
#>  af.tsoil = approxfun(ERA.dates,ts.tsoil[1,])
#>  af.srad = approxfun(ERA.dates,ts.srad[1,])
#>  af.moist1 = approxfun(ERA.dates,ts.moist1[1,])
#>  af.moist2 = approxfun(ERA.dates,ts.moist2[1,])
#>  #-3-create an approx function for the EVI data
#>  af.evi = approxfun(as.Date(date.pos),ts.evi[p,])
#>  #-4-use the approxfun's to fill the data frame data.frame
#>  pixel.df$tair <- af.tair(as.Date(pixel.date))
#>  pixel.df$tsoil <- af.tsoil(as.Date(pixel.date))
#>  pixel.df$srad <- af.srad(as.Date(pixel.date))
#>  pixel.df$moist1 <- af.moist1(as.Date(pixel.date))
#>  pixel.df$moist2 <- af.moist2(as.Date(pixel.date))
#>  pixel.df$evi <- af.evi(as.Date(pixel.date))
#>  pixel.df<-pixel.df[complete.cases(pixel.df),]
#>  
#>  coords<-round(coordinates(tile.xy)[p,],2)
#>  
#>  print("before")
#> [1] "before"
#>  
#>  print(paste("statsModis_t:"), typeof(statsMODIS[p,]))
#> Error in print.default(paste("statsModis_t:"), typeof(statsMODIS[p, ])) : 
#>   invalid 'digits' argument
#> In addition: Warning message:
#> In print.default(paste("statsModis_t:"), typeof(statsMODIS[p, ])) :
#>   NAs introduced by coercion
#>  print(paste("statsModis_c:"), class(statsMODIS[p,]))
#> Error in print.default(paste("statsModis_c:"), class(statsMODIS[p, ])) : 
#>   invalid 'digits' argument
#> In addition: Warning message:
#> In print.default(paste("statsModis_c:"), class(statsMODIS[p, ])) :
#>   NAs introduced by coercion
#>  print(paste("statsModis_d:"), dim(statsMODIS[p,]))
#> [1] "statsModis_d:"
#>  print(paste("statsModis_l:"), length(statsMODIS[p,]))
#> Error in print.default(paste("statsModis_l:"), length(statsMODIS[p, ])) : 
#>   invalid 'digits' argument
#>  print(paste("statsModis_h:"), head(statsMODIS[p,]))
#> [1] "statsModis_h:"
#>  
#>  print(paste("dopixel_t:"), typeof(dopixel(pixel.df,coords,PLOT=FALSE)))
#> Error in print.default(paste("dopixel_t:"), typeof(dopixel(pixel.df, coords,  : 
#>   invalid 'digits' argument
#> In addition: Warning messages:
#> 1: In pd.evi - td.evi[-1] :
#>   longer object length is not a multiple of shorter object length
#> 2: In print.default(paste("dopixel_t:"), typeof(dopixel(pixel.df, coords,  :
#>   NAs introduced by coercion
#>  print(paste("dopixel_c:"), class(dopixel(pixel.df,coords,PLOT=FALSE)))
#> Error in print.default(paste("dopixel_c:"), class(dopixel(pixel.df, coords,  : 
#>   invalid 'digits' argument
#> In addition: Warning messages:
#> 1: In pd.evi - td.evi[-1] :
#>   longer object length is not a multiple of shorter object length
#> 2: In print.default(paste("dopixel_c:"), class(dopixel(pixel.df, coords,  :
#>   NAs introduced by coercion
#>  print(paste("dopixel_d:"), dim(dopixel(pixel.df,coords,PLOT=FALSE)))
#> [1] "dopixel_d:"
#> Warning message:
#> In pd.evi - td.evi[-1] :
#>   longer object length is not a multiple of shorter object length
#>  print(paste("dopixel_l:"), length(dopixel(pixel.df,coords,PLOT=FALSE)))
#> Error in print.default(paste("dopixel_l:"), length(dopixel(pixel.df, coords,  : 
#>   invalid 'digits' argument
#> In addition: Warning message:
#> In pd.evi - td.evi[-1] :
#>   longer object length is not a multiple of shorter object length
#>  print(paste("dopixel_h:"), head(dopixel(pixel.df,coords,PLOT=FALSE)))
#> [1] "dopixel_h:"
#> Warning message:
#> In pd.evi - td.evi[-1] :
#>   longer object length is not a multiple of shorter object length
#>  
#>  statsMODIS[p,] <- dopixel(pixel.df,coords,PLOT=FALSE)
#> Warning message:
#> In pd.evi - td.evi[-1] :
#>   longer object length is not a multiple of shorter object length
#>  
#>  print("after")
#> [1] "after"
#>  
#>  ### 2nd iteration
#>  
#>  p <- 77919
#> 
#>  print(p)
#> [1] 77919
#>  
#>  #-1-create a data.frame for the pixel
#>  pixel.df = data.frame(date = pixel.date)
#>  #-2-create an approx function for each ERA variable
#>  af.tair = approxfun(ERA.dates,ts.tair[1,])
#>  af.tsoil = approxfun(ERA.dates,ts.tsoil[1,])
#>  af.srad = approxfun(ERA.dates,ts.srad[1,])
#>  af.moist1 = approxfun(ERA.dates,ts.moist1[1,])
#>  af.moist2 = approxfun(ERA.dates,ts.moist2[1,])
#>  #-3-create an approx function for the EVI data
#>  af.evi = approxfun(as.Date(date.pos),ts.evi[p,])
#>  #-4-use the approxfun's to fill the data frame data.frame
#>  pixel.df$tair <- af.tair(as.Date(pixel.date))
#>  pixel.df$tsoil <- af.tsoil(as.Date(pixel.date))
#>  pixel.df$srad <- af.srad(as.Date(pixel.date))
#>  pixel.df$moist1 <- af.moist1(as.Date(pixel.date))
#>  pixel.df$moist2 <- af.moist2(as.Date(pixel.date))
#>  pixel.df$evi <- af.evi(as.Date(pixel.date))
#>  pixel.df<-pixel.df[complete.cases(pixel.df),]
#>  
#>  coords<-round(coordinates(tile.xy)[p,],2)
#>  
#>  print("before")
#> [1] "before"
#> 
#>  print(paste("statsModis_t:"), typeof(statsMODIS[p,]))
#> Error in print.default(paste("statsModis_t:"), typeof(statsMODIS[p, ])) : 
#>   invalid 'digits' argument
#> In addition: Warning message:
#> In print.default(paste("statsModis_t:"), typeof(statsMODIS[p, ])) :
#>   NAs introduced by coercion
#> print(paste("statsModis_c:"), class(statsMODIS[p,]))
#> Error in print.default(paste("statsModis_c:"), class(statsMODIS[p, ])) : 
#>   invalid 'digits' argument
#> In addition: Warning message:
#> In print.default(paste("statsModis_c:"), class(statsMODIS[p, ])) :
#>   NAs introduced by coercion
#>  print(paste("statsModis_d:"), dim(statsMODIS[p,]))
#> [1] "statsModis_d:"
#>  print(paste("statsModis_l:"), length(statsMODIS[p,]))
#> Error in print.default(paste("statsModis_l:"), length(statsMODIS[p, ])) : 
#>   invalid 'digits' argument
#>  print(paste("statsModis_h:"), head(statsMODIS[p,]))
#> [1] "statsModis_h:"
#>  
#>  print(paste("dopixel_t:"), typeof(dopixel(pixel.df,coords,PLOT=FALSE)))
#> Error in print.default(paste("dopixel_t:"), typeof(dopixel(pixel.df, coords,  : 
#>   invalid 'digits' argument
#> In addition: Warning messages:
#> 1: In pd.evi - td.evi[-1] :
#>   longer object length is not a multiple of shorter object length
#> 2: In print.default(paste("dopixel_t:"), typeof(dopixel(pixel.df, coords,  :
#>   NAs introduced by coercion
#>  print(paste("dopixel_c:"), class(dopixel(pixel.df,coords,PLOT=FALSE)))
#> Error in print.default(paste("dopixel_c:"), class(dopixel(pixel.df, coords,  : 
#>   invalid 'digits' argument
#> In addition: Warning messages:
#> 1: In pd.evi - td.evi[-1] :
#>   longer object length is not a multiple of shorter object length
#> 2: In print.default(paste("dopixel_c:"), class(dopixel(pixel.df, coords,  :
#>   NAs introduced by coercion
#>  print(paste("dopixel_d:"), dim(dopixel(pixel.df,coords,PLOT=FALSE)))
#> [1] "dopixel_d:"
#> Warning message:
#> In pd.evi - td.evi[-1] :
#>   longer object length is not a multiple of shorter object length
#>  print(paste("dopixel_l:"), length(dopixel(pixel.df,coords,PLOT=FALSE)))
#> Error in print.default(paste("dopixel_l:"), length(dopixel(pixel.df, coords,  : 
#>   invalid 'digits' argument
#> In addition: Warning message:
#> In pd.evi - td.evi[-1] :
#>   longer object length is not a multiple of shorter object length
#>  print(paste("dopixel_h:"), head(dopixel(pixel.df,coords,PLOT=FALSE)))
#> [1] "dopixel_h:"
#> Warning message:
#> In pd.evi - td.evi[-1] :
#>   longer object length is not a multiple of shorter object length
#>  
#>  statsMODIS[p,] <- dopixel(pixel.df,coords,PLOT=FALSE)
#> Warning message:
#> In pd.evi - td.evi[-1] :
#>   longer object length is not a multiple of shorter object length
#>  
#>  print("after")
#> [1] "after"
36

There were some issues with parentheses etc. I corrected it. Can you run it again?

library(bigmemory)


#Load the script needed to run do_pixel
source("~/Work/dopixel_2020.R")


#--load dates which are in R date format
ERA.dates <- readRDS("~/Work/ERA.dates.rds") #240 dates

#Load environmental variables

ts.evi <- readRDS("~/Work/ts.evi.rds")
tile.xy <- readRDS("~/Work/tile.xy.rds")
ts.moist1 <- readRDS("~/Work/ts.moist1.rds")
ts.moist2 <- readRDS("~/Work/ts.moist2.rds")
ts.srad <- readRDS("~/Work/ts.srad.rds")
ts.tair <- readRDS("~/Work/ts.tair.rds")
ts.tsoil <- readRDS("~/Work/ts.tsoil.rds")


#------------------------------------------------------------------#
#--get date information from the MODIS time series-----------------#
#------------------------------------------------------------------#
# get the layer names provided by the call to extract evi:
layer.names<-labels(ts.evi)[[2]]   
# take the substring on these names, which give the dates:
date.text <-substr(layer.names,start=10,stop=16)
date.pos <- strptime(date.text,  format = "%Y%j", tz="UTC") #convert to R date format


#------------------------------------------------------------------#
#--create big matrix to store output of do_pixel-------------------#
#------------------------------------------------------------------#
pixel.date <- as.Date(seq.POSIXt(ISOdate(2000,6,1),by="week",length.out=1020,tz="UTC"))
#tail(pixel.date)
YEAR <- as.numeric(format(pixel.date, "%Y"))
YITS = length(unique(YEAR))
#--the column names for the big matrix
names<-c("lon","lat",
         "lq", "uq", "mean.evi", "sd.evi", "sum.evi", "amplitude",
         "peak.day", "trough.day",
         paste("td", 1:YITS, sep="."), paste("td.x", 1:YITS, sep="."), paste("td.evi", 1:YITS, sep="."),
         paste("pd", 1:YITS, sep="."), paste("pd.x", 1:YITS, sep="."), paste("pd.evi", 1:YITS, sep="."),
         paste("elon.m", 1:YITS, sep="."), paste("elon.m.x", 1:YITS, sep="."),
         paste("eoff.m", 1:YITS, sep="."), paste("eoff.m.x", 1:YITS, sep="."),
         paste("elon.f", 1:YITS, sep="."), paste("elon.f.x", 1:YITS, sep="."),
         paste("eoff.f", 1:YITS, sep="."), paste("eoff.f.x", 1:YITS, sep="."),
         paste("elon.l", 1:YITS, sep="."), paste("elon.l.x", 1:YITS, sep="."),
         paste("eoff.l", 1:YITS, sep="."), paste("eoff.l.x", 1:YITS, sep="."),
         paste("sum.evi.yr", 1:YITS, sep="."), paste("amp", 1:YITS, sep="."),
         paste("gsl", 1:YITS, sep="."), paste("gsl.peak", 1:YITS, sep="."), paste("gsl.long", 1:YITS, sep="."),
         paste("elon.f.evi", 1:YITS, sep="."), paste("elon.m.evi", 1:YITS, sep="."), paste("elon.l.evi", 1:YITS, sep="."),
         paste("eoff.f.evi", 1:YITS, sep="."), paste("eoff.m.evi", 1:YITS, sep="."), paste("eoff.l.evi", 1:YITS, sep="."),
         paste("tair.rank.td", 1:YITS, sep="."),paste("tair.rank.pd", 1:YITS, sep="."),
         paste("moist1.rank.td", 1:YITS, sep="."),paste("moist1.rank.pd", 1:YITS, sep=".")
) 
#--create the matrix using pixel.df to specifiy rows
#--and column names to speciify columns
statsMODIS<-big.matrix(nrow=length(tile.xy), ncol=length(names), type="double",
                       dimnames=list(row.names=NULL, col.names=names),
                       backingfile="statsMODIS.bin", descriptorfile="statsMODIS.desc")


#------------------------------------------------------------------#
#--start a for loop that will call do_pixel------------------------#
#--this will be slow, but its useful for debugging-----------------#
#------------------------------------------------------------------#
ERA.dates <- ERA.dates+14

### 1st iteration

p <- 77919

print(p)

#-1-create a data.frame for the pixel
pixel.df = data.frame(date = pixel.date)
#-2-create an approx function for each ERA variable
af.tair = approxfun(ERA.dates,ts.tair[1,])
af.tsoil = approxfun(ERA.dates,ts.tsoil[1,])
af.srad = approxfun(ERA.dates,ts.srad[1,])
af.moist1 = approxfun(ERA.dates,ts.moist1[1,])
af.moist2 = approxfun(ERA.dates,ts.moist2[1,])
#-3-create an approx function for the EVI data
af.evi = approxfun(as.Date(date.pos),ts.evi[p,])
#-4-use the approxfun's to fill the data frame data.frame
pixel.df$tair <- af.tair(as.Date(pixel.date))
pixel.df$tsoil <- af.tsoil(as.Date(pixel.date))
pixel.df$srad <- af.srad(as.Date(pixel.date))
pixel.df$moist1 <- af.moist1(as.Date(pixel.date))
pixel.df$moist2 <- af.moist2(as.Date(pixel.date))
pixel.df$evi <- af.evi(as.Date(pixel.date))
pixel.df<-pixel.df[complete.cases(pixel.df),]

coords<-round(coordinates(tile.xy)[p,],2)

print("before")

print(paste(c("statsModis:", typeof(statsMODIS[p,]))))
print(paste(c("statsModis:", class(statsMODIS[p,]))))
print(paste(c("statsModis:", dim(statsMODIS[p,]))))
print(paste(c("statsModis:", length(statsMODIS[p,]))))
print(paste(c("statsModis:", head(statsMODIS[p,]))))

print(paste(c("dopixel:", typeof(dopixel(pixel.df,coords,PLOT=FALSE)))))
print(paste(c("dopixel:", class(dopixel(pixel.df,coords,PLOT=FALSE)))))
print(paste(c("dopixel:", dim(dopixel(pixel.df,coords,PLOT=FALSE)))))
print(paste(c("dopixel:", length(dopixel(pixel.df,coords,PLOT=FALSE)))))
print(paste(c("dopixel:", head(dopixel(pixel.df,coords,PLOT=FALSE)))))

statsMODIS[p,] <- dopixel(pixel.df,coords,PLOT=FALSE)

print("after")

### 2nd iteration

p <- 77919

print(p)

#-1-create a data.frame for the pixel
pixel.df = data.frame(date = pixel.date)
#-2-create an approx function for each ERA variable
af.tair = approxfun(ERA.dates,ts.tair[1,])
af.tsoil = approxfun(ERA.dates,ts.tsoil[1,])
af.srad = approxfun(ERA.dates,ts.srad[1,])
af.moist1 = approxfun(ERA.dates,ts.moist1[1,])
af.moist2 = approxfun(ERA.dates,ts.moist2[1,])
#-3-create an approx function for the EVI data
af.evi = approxfun(as.Date(date.pos),ts.evi[p,])
#-4-use the approxfun's to fill the data frame data.frame
pixel.df$tair <- af.tair(as.Date(pixel.date))
pixel.df$tsoil <- af.tsoil(as.Date(pixel.date))
pixel.df$srad <- af.srad(as.Date(pixel.date))
pixel.df$moist1 <- af.moist1(as.Date(pixel.date))
pixel.df$moist2 <- af.moist2(as.Date(pixel.date))
pixel.df$evi <- af.evi(as.Date(pixel.date))
pixel.df<-pixel.df[complete.cases(pixel.df),]

coords<-round(coordinates(tile.xy)[p,],2)

print("before")

print(paste(c("statsModis:", typeof(statsMODIS[p,]))))
print(paste(c("statsModis:", class(statsMODIS[p,]))))
print(paste(c("statsModis:", dim(statsMODIS[p,]))))
print(paste(c("statsModis:", length(statsMODIS[p,]))))
print(paste(c("statsModis:", head(statsMODIS[p,]))))

print(paste(c("dopixel:", typeof(dopixel(pixel.df,coords,PLOT=FALSE)))))
print(paste(c("dopixel:", class(dopixel(pixel.df,coords,PLOT=FALSE)))))
print(paste(c("dopixel:", dim(dopixel(pixel.df,coords,PLOT=FALSE)))))
print(paste(c("dopixel:", length(dopixel(pixel.df,coords,PLOT=FALSE)))))
print(paste(c("dopixel:", head(dopixel(pixel.df,coords,PLOT=FALSE)))))

statsMODIS[p,] <- dopixel(pixel.df,coords,PLOT=FALSE)

print("after")
37

BTW: Maybe you can reduce the size of statsMODIS by subsetting it for only the relevant rows. This SO-Post suggests to use the sub.big.matrix() function for that. Could you try to incorporate this into your code to make the reprex smaller.

38 
#> setwd("~/Work")
#> library(bigmemory)
#> 
#> 
#> #Load the script needed to run do_pixel
#> source("~/Work/dopixel_2020.R")
#> 
#> 
#> #--load dates which are in R date format
#> ERA.dates <- readRDS("~/Work/ERA.dates.rds") #240 dates
#> 
#> #Load environmental variables
#> 
#> ts.evi <- readRDS("~/Work/ts.evi.rds")
#> tile.xy <- readRDS("~/Work/tile.xy.rds")
#> ts.moist1 <- readRDS("~/Work/ts.moist1.rds")
#> ts.moist2 <- readRDS("~/Work/ts.moist2.rds")
#> ts.srad <- readRDS("~/Work/ts.srad.rds")
#> ts.tair <- readRDS("~/Work/ts.tair.rds")
#> ts.tsoil <- readRDS("~/Work/ts.tsoil.rds")
#> 
#> 
#> #------------------------------------------------------------------#
#> #--get date information from the MODIS time series-----------------#
#> #------------------------------------------------------------------#
#> # get the layer names provided by the call to extract evi:
#> layer.names<-labels(ts.evi)[[2]]   
#> # take the substring on these names, which give the dates:
#> date.text <-substr(layer.names,start=10,stop=16)
#> date.pos <- strptime(date.text,  format = "%Y%j", tz="UTC") #convert to R date format
#> 
#> 
#> #------------------------------------------------------------------#
#> #--create big matrix to store output of do_pixel-------------------#
#> #------------------------------------------------------------------#
#> pixel.date <- as.Date(seq.POSIXt(ISOdate(2000,6,1),by="week",length.out=1020,tz="UTC"))
#> #tail(pixel.date)
#> YEAR <- as.numeric(format(pixel.date, "%Y"))
#> YITS = length(unique(YEAR))
#> #--the column names for the big matrix
#> names<-c("lon","lat",
#>          "lq", "uq", "mean.evi", "sd.evi", "sum.evi", "amplitude",
#>          "peak.day", "trough.day",
#>          paste("td", 1:YITS, sep="."), paste("td.x", 1:YITS, sep="."), paste("td.evi", 1:YITS, sep="."),
#>          paste("pd", 1:YITS, sep="."), paste("pd.x", 1:YITS, sep="."), paste("pd.evi", 1:YITS, sep="."),
#>          paste("elon.m", 1:YITS, sep="."), paste("elon.m.x", 1:YITS, sep="."),
#>          paste("eoff.m", 1:YITS, sep="."), paste("eoff.m.x", 1:YITS, sep="."),
#>          paste("elon.f", 1:YITS, sep="."), paste("elon.f.x", 1:YITS, sep="."),
#>          paste("eoff.f", 1:YITS, sep="."), paste("eoff.f.x", 1:YITS, sep="."),
#>          paste("elon.l", 1:YITS, sep="."), paste("elon.l.x", 1:YITS, sep="."),
#>          paste("eoff.l", 1:YITS, sep="."), paste("eoff.l.x", 1:YITS, sep="."),
#>          paste("sum.evi.yr", 1:YITS, sep="."), paste("amp", 1:YITS, sep="."),
#>          paste("gsl", 1:YITS, sep="."), paste("gsl.peak", 1:YITS, sep="."), paste("gsl.long", 1:YITS, sep="."),
#>          paste("elon.f.evi", 1:YITS, sep="."), paste("elon.m.evi", 1:YITS, sep="."), paste("elon.l.evi", 1:YITS, sep="."),
#>          paste("eoff.f.evi", 1:YITS, sep="."), paste("eoff.m.evi", 1:YITS, sep="."), paste("eoff.l.evi", 1:YITS, sep="."),
#>          paste("tair.rank.td", 1:YITS, sep="."),paste("tair.rank.pd", 1:YITS, sep="."),
#>          paste("moist1.rank.td", 1:YITS, sep="."),paste("moist1.rank.pd", 1:YITS, sep=".")
#> ) 
#> #--create the matrix using pixel.df to specifiy rows
#> #--and column names to speciify columns
#> statsMODIS<-big.matrix(nrow=length(tile.xy), ncol=length(names), type="double",
#>                        dimnames=list(row.names=NULL, col.names=names),
#>                        backingfile="statsMODIS.bin", descriptorfile="statsMODIS.desc")
#Loading required package: sp
#> 
#> 
#> #------------------------------------------------------------------#
#> #--start a for loop that will call do_pixel------------------------#
#> #--this will be slow, but its useful for debugging-----------------#
#> #------------------------------------------------------------------#
#> ERA.dates <- ERA.dates+14
#> 
#> ### 1st iteration
#> 
#> p <- 77919
#> 
#> print(p)
#[1] 77919
#> 
#> #-1-create a data.frame for the pixel
#> pixel.df = data.frame(date = pixel.date)
#> #-2-create an approx function for each ERA variable
#> af.tair = approxfun(ERA.dates,ts.tair[1,])
#> af.tsoil = approxfun(ERA.dates,ts.tsoil[1,])
#> af.srad = approxfun(ERA.dates,ts.srad[1,])
#> af.moist1 = approxfun(ERA.dates,ts.moist1[1,])
#> af.moist2 = approxfun(ERA.dates,ts.moist2[1,])
#> #-3-create an approx function for the EVI data
#> af.evi = approxfun(as.Date(date.pos),ts.evi[p,])
#> #-4-use the approxfun's to fill the data frame data.frame
#> pixel.df$tair <- af.tair(as.Date(pixel.date))
#> pixel.df$tsoil <- af.tsoil(as.Date(pixel.date))
#> pixel.df$srad <- af.srad(as.Date(pixel.date))
#> pixel.df$moist1 <- af.moist1(as.Date(pixel.date))
#> pixel.df$moist2 <- af.moist2(as.Date(pixel.date))
#> pixel.df$evi <- af.evi(as.Date(pixel.date))
#> pixel.df<-pixel.df[complete.cases(pixel.df),]
#> 
#> coords<-round(coordinates(tile.xy)[p,],2)
#> 
#> print("before")
#> [1] "before"
#> 
#> print(paste(c("statsModis:", typeof(statsMODIS[p,]))))
#> [1] "statsModis:" "double"     
#> print(paste(c("statsModis:", class(statsMODIS[p,]))))
#> [1] "statsModis:" "numeric"    
#> print(paste(c("statsModis:", dim(statsMODIS[p,]))))
#> [1] "statsModis:"
#> print(paste(c("statsModis:", length(statsMODIS[p,]))))
#> [1] "statsModis:" "670"        
#> print(paste(c("statsModis:", head(statsMODIS[p,]))))
#>[1] "statsModis:" "0"           "0"           "0"           "0"           "0"           "0"          
#> 
#> print(paste(c("dopixel:", typeof(dopixel(pixel.df,coords,PLOT=FALSE)))))
#> [1] "dopixel:" "double"  
#> Warning message:
#> In pd.evi - td.evi[-1] :
#> longer object length is not a multiple of shorter object length
#> print(paste(c("dopixel:", class(dopixel(pixel.df,coords,PLOT=FALSE)))))
#> [1] "dopixel:" "numeric" 
#> Warning message:
#> In pd.evi - td.evi[-1] :
#>  longer object length is not a multiple of shorter object length
#> print(paste(c("dopixel:", dim(dopixel(pixel.df,coords,PLOT=FALSE)))))
#> [1] "dopixel:"
#> Warning message:
#> In pd.evi - td.evi[-1] :
#>  longer object length is not a multiple of shorter object length
#> print(paste(c("dopixel:", length(dopixel(pixel.df,coords,PLOT=FALSE)))))
#> [1] "dopixel:" "670"     
#> Warning message:
#> In pd.evi - td.evi[-1] :
#>  longer object length is not a multiple of shorter object length
#> print(paste(c("dopixel:", head(dopixel(pixel.df,coords,PLOT=FALSE)))))
#> [1] "dopixel:"           "15.62"              "-18.16"             "0.113794230769231"  "0.252255"             #>0.169362790269425"  "0.0538153702024827"
#> Warning message:
#> In pd.evi - td.evi[-1] :
#>  longer object length is not a multiple of shorter object length
#> 
#> statsMODIS[p,] <- dopixel(pixel.df,coords,PLOT=FALSE)
#> Warning message:
#> In pd.evi - td.evi[-1] :
#>  longer object length is not a multiple of shorter object length
#> 
#> print("after")
#> [1] "after"
#> 
#> ### 2nd iteration
#> 
#> p <- 77919
#> 
#> print(p)
#> [1] 77919
#> 
#> #-1-create a data.frame for the pixel
#> pixel.df = data.frame(date = pixel.date)
#> #-2-create an approx function for each ERA variable
#> af.tair = approxfun(ERA.dates,ts.tair[1,])
#> af.tsoil = approxfun(ERA.dates,ts.tsoil[1,])
#> af.srad = approxfun(ERA.dates,ts.srad[1,])
#> af.moist1 = approxfun(ERA.dates,ts.moist1[1,])
#> af.moist2 = approxfun(ERA.dates,ts.moist2[1,])
#> #-3-create an approx function for the EVI data
#> af.evi = approxfun(as.Date(date.pos),ts.evi[p,])
#> #-4-use the approxfun's to fill the data frame data.frame
#> pixel.df$tair <- af.tair(as.Date(pixel.date))
#> pixel.df$tsoil <- af.tsoil(as.Date(pixel.date))
#> pixel.df$srad <- af.srad(as.Date(pixel.date))
#> pixel.df$moist1 <- af.moist1(as.Date(pixel.date))
#> pixel.df$moist2 <- af.moist2(as.Date(pixel.date))
#> pixel.df$evi <- af.evi(as.Date(pixel.date))
#> pixel.df<-pixel.df[complete.cases(pixel.df),]
#> 
#> coords<-round(coordinates(tile.xy)[p,],2)
#> 
#> print("before")
#> [1] "before"
#> 
#> print(paste(c("statsModis:", typeof(statsMODIS[p,]))))
#> [1] "statsModis:" "double"     
#> print(paste(c("statsModis:", class(statsMODIS[p,]))))
#> [1] "statsModis:" "numeric"    
#> print(paste(c("statsModis:", dim(statsMODIS[p,]))))
#> [1] "statsModis:"
#> print(paste(c("statsModis:", length(statsMODIS[p,]))))
#> [1] "statsModis:" "670"        
#> print(paste(c("statsModis:", head(statsMODIS[p,]))))
#> [1] "statsModis:"        "15.62"              "-18.16"             "0.113794230769231"  "0.252255"           #>"0.169362790269425"  "0.0538153702024827"
#> 
#> print(paste(c("dopixel:", typeof(dopixel(pixel.df,coords,PLOT=FALSE)))))
#> [1] "dopixel:" "double"  
#> Warning message:
#> In pd.evi - td.evi[-1] :
#>  longer object length is not a multiple of shorter object length
#> print(paste(c("dopixel:", class(dopixel(pixel.df,coords,PLOT=FALSE)))))
#> [1] "dopixel:" "numeric" 
#> Warning message:
#> In pd.evi - td.evi[-1] :
#>  longer object length is not a multiple of shorter object length
#> print(paste(c("dopixel:", dim(dopixel(pixel.df,coords,PLOT=FALSE)))))
#> [1] "dopixel:"
#> Warning message:
#> In pd.evi - td.evi[-1] :
#>  longer object length is not a multiple of shorter object length
#> print(paste(c("dopixel:", length(dopixel(pixel.df,coords,PLOT=FALSE)))))
#> [1] "dopixel:" "670"     
#> Warning message:
#> In pd.evi - td.evi[-1] :
#>  longer object length is not a multiple of shorter object length
#> print(paste(c("dopixel:", head(dopixel(pixel.df,coords,PLOT=FALSE)))))
#> [1] "dopixel:"           "15.62"              "-18.16"             "0.113794230769231"  "0.252255"           #>"0.169362790269425"  "0.0538153702024827"
#>Warning message:
#> In pd.evi - td.evi[-1] :
#>  longer object length is not a multiple of shorter object length
#> 
#> statsMODIS[p,] <- dopixel(pixel.df,coords,PLOT=FALSE)
#> Warning message:
#> In pd.evi - td.evi[-1] :
#>  longer object length is not a multiple of shorter object length
#> 
#> print("after")
#> [1] "after"
39

Alright, I decided to check whether similar errors exist in other rows of the big matrix. statsMODIS[86016,] reports the same error suggesting that there are many rows of that nature. I am afraid all the rows are relevant as they contain data for a study area. I will send you results of a subset from row 77919 to row 86017.

40

I have created a subset of the big matrix and tile.xy then ran the code again.


#> library(bigmemory)
#> 
#> 
#> #Load the script needed to run do_pixel
#> source("~/Work/dopixel_2020.R")
#> 
#> 
#> #--load dates which are in R date format
#> ERA.dates <- readRDS("~/Work/ERA.dates.rds") #240 dates
#> 
#> #Load environmental variables
#> 
#> ts.evi <- readRDS("~/Work/ts.evi.rds")
#> tile.xy <- readRDS("~/Work/tile.xy.rds")
#> ts.moist1 <- readRDS("~/Work/ts.moist1.rds")
#> ts.moist2 <- readRDS("~/Work/ts.moist2.rds")
#> ts.srad <- readRDS("~/Work/ts.srad.rds")
#> ts.tair <- readRDS("~/Work/ts.tair.rds")
#> ts.tsoil <- readRDS("~/Work/ts.tsoil.rds")
#> 
#> 
#> #------------------------------------------------------------------#
#> #--get date information from the MODIS time series-----------------#
#> #------------------------------------------------------------------#
#> # get the layer names provided by the call to extract evi:
#> layer.names<-labels(ts.evi)[[2]]   
#> # take the substring on these names, which give the dates:
#> date.text <-substr(layer.names,start=10,stop=16)
#> date.pos <- strptime(date.text,  format = "%Y%j", tz="UTC") #convert to R date format
#> 
#> 
#> #------------------------------------------------------------------#
#> #--create big matrix to store output of do_pixel-------------------#
#> #------------------------------------------------------------------#
#> pixel.date <- as.Date(seq.POSIXt(ISOdate(2000,6,1),by="week",length.out=1020,tz="UTC"))
#> #tail(pixel.date)
#> YEAR <- as.numeric(format(pixel.date, "%Y"))
#> YITS = length(unique(YEAR))
#> #--the column names for the big matrix
#> names<-c("lon","lat",
#>          "lq", "uq", "mean.evi", "sd.evi", "sum.evi", "amplitude",
#>          "peak.day", "trough.day",
#>          paste("td", 1:YITS, sep="."), paste("td.x", 1:YITS, sep="."), paste("td.evi", 1:YITS, sep="."),
#>          paste("pd", 1:YITS, sep="."), paste("pd.x", 1:YITS, sep="."), paste("pd.evi", 1:YITS, sep="."),
#>          paste("elon.m", 1:YITS, sep="."), paste("elon.m.x", 1:YITS, sep="."),
#>          paste("eoff.m", 1:YITS, sep="."), paste("eoff.m.x", 1:YITS, sep="."),
#>          paste("elon.f", 1:YITS, sep="."), paste("elon.f.x", 1:YITS, sep="."),
#>          paste("eoff.f", 1:YITS, sep="."), paste("eoff.f.x", 1:YITS, sep="."),
#>          paste("elon.l", 1:YITS, sep="."), paste("elon.l.x", 1:YITS, sep="."),
#>          paste("eoff.l", 1:YITS, sep="."), paste("eoff.l.x", 1:YITS, sep="."),
#>          paste("sum.evi.yr", 1:YITS, sep="."), paste("amp", 1:YITS, sep="."),
#>          paste("gsl", 1:YITS, sep="."), paste("gsl.peak", 1:YITS, sep="."), paste("gsl.long", 1:YITS, sep="."),
#>          paste("elon.f.evi", 1:YITS, sep="."), paste("elon.m.evi", 1:YITS, sep="."), paste("elon.l.evi", 1:YITS, sep="."),
#>          paste("eoff.f.evi", 1:YITS, sep="."), paste("eoff.m.evi", 1:YITS, sep="."), paste("eoff.l.evi", 1:YITS, sep="."),
#>          paste("tair.rank.td", 1:YITS, sep="."),paste("tair.rank.pd", 1:YITS, sep="."),
#>          paste("moist1.rank.td", 1:YITS, sep="."),paste("moist1.rank.pd", 1:YITS, sep=".")
#> ) 
#> #--create the matrix using pixel.df to specifiy rows
#> #--and column names to speciify columns
#> statsMODIS<-big.matrix(nrow=length(tile.xy), ncol=length(names), type="double",
#>                        dimnames=list(row.names=NULL, col.names=names),
#>                        backingfile="statsMODIS.bin", descriptorfile="statsMODIS.desc")
#>Loading required package: sp
#> 
#> ###########Susbet the big matrix and tile.xy
#> 
#> statsMODIS <- sub.big.matrix(statsMODIS, firstRow = 77919, lastRow = 86017)
#> tile.xy <-tile.xy[c(77919:86017),]
#> dim(statsMODIS)
#>[1] 8099  670
#> #------------------------------------------------------------------#
#> #--start a for loop that will call do_pixel------------------------#
#> #--this will be slow, but its useful for debugging-----------------#
#> #------------------------------------------------------------------#
#> ERA.dates <- ERA.dates+14
#> 
#> ### 1st iteration
#> 
#> p <- 77919
#> 
#> print(p)
#> [1] 77919
#> 
#> #-1-create a data.frame for the pixel
#> pixel.df = data.frame(date = pixel.date)
#> #-2-create an approx function for each ERA variable
#> af.tair = approxfun(ERA.dates,ts.tair[1,])
#> af.tsoil = approxfun(ERA.dates,ts.tsoil[1,])
#> af.srad = approxfun(ERA.dates,ts.srad[1,])
#> af.moist1 = approxfun(ERA.dates,ts.moist1[1,])
#> af.moist2 = approxfun(ERA.dates,ts.moist2[1,])
#> #-3-create an approx function for the EVI data
#> af.evi = approxfun(as.Date(date.pos),ts.evi[p,])
#> #-4-use the approxfun's to fill the data frame data.frame
#> pixel.df$tair <- af.tair(as.Date(pixel.date))
#> pixel.df$tsoil <- af.tsoil(as.Date(pixel.date))
#> pixel.df$srad <- af.srad(as.Date(pixel.date))
#> pixel.df$moist1 <- af.moist1(as.Date(pixel.date))
#> pixel.df$moist2 <- af.moist2(as.Date(pixel.date))
#> pixel.df$evi <- af.evi(as.Date(pixel.date))
#> pixel.df<-pixel.df[complete.cases(pixel.df),]
#> 
#> coords<-round(coordinates(tile.xy)[p,],2)
#> Error in coordinates(tile.xy)[p, ] : subscript out of bounds
#> 
#> print("before")
#> [1] "before"
#> 
#> print(paste(c("statsModis:", typeof(statsMODIS[p,]))))
#>Error in GetRows.bm(x, i) : Illegal row index usage in extraction.
#> print(paste(c("statsModis:", class(statsMODIS[p,]))))
#>Error in GetRows.bm(x, i) : Illegal row index usage in extraction.
#> print(paste(c("statsModis:", dim(statsMODIS[p,]))))
#>Error in GetRows.bm(x, i) : Illegal row index usage in extraction.
#> print(paste(c("statsModis:", length(statsMODIS[p,]))))
#>Error in GetRows.bm(x, i) : Illegal row index usage in extraction.
#> print(paste(c("statsModis:", head(statsMODIS[p,]))))
#> Error in GetRows.bm(x, i) : Illegal row index usage in extraction.
#> 
#> print(paste(c("dopixel:", typeof(dopixel(pixel.df,coords,PLOT=FALSE)))))
#>Error in data.frame(coords[1], coords[2], lq, uq, mean.evi, sd.evi, sum.evi,  : 
#>  object 'coords' not found
#>In addition: Warning message:
#>In pd.evi - td.evi[-1] :
#> 
#> Error in data.frame(coords[1], coords[2], lq, uq, mean.evi, sd.evi, sum.evi,  : 
#>  object 'coords' not found > print(paste(c("dopixel:", class(dopixel(pixel.df,coords,PLOT=FALSE)))))
#>Error in data.frame(coords[1], coords[2], lq, uq, mean.evi, sd.evi, sum.evi,  : 
#>  object 'coords' not found
#>In addition: Warning message:
#>In pd.evi - td.evi[-1] :
#> 
#> Error in data.frame(coords[1], coords[2], lq, uq, mean.evi, sd.evi, sum.evi,  : 
#>  object 'coords' not found > print(paste(c("dopixel:", dim(dopixel(pixel.df,coords,PLOT=FALSE)))))
#>Error in data.frame(coords[1], coords[2], lq, uq, mean.evi, sd.evi, sum.evi,  : 
#>  object 'coords' not found
#>In addition: Warning message:
#>In pd.evi - td.evi[-1] :
 
#> Error in data.frame(coords[1], coords[2], lq, uq, mean.evi, sd.evi, sum.evi,  : 
#>  object 'coords' not found > print(paste(c("dopixel:", length(dopixel(pixel.df,coords,PLOT=FALSE)))))
#>Error in data.frame(coords[1], coords[2], lq, uq, mean.evi, sd.evi, sum.evi,  : 
#>  object 'coords' not found
#>In addition: Warning message:
#>In pd.evi - td.evi[-1] :
#> 
#> Error in data.frame(coords[1], coords[2], lq, uq, mean.evi, sd.evi, sum.evi,  : 
#>  object 'coords' not found > print(paste(c("dopixel:", head(dopixel(pixel.df,coords,PLOT=FALSE)))))
#>Error in data.frame(coords[1], coords[2], lq, uq, mean.evi, sd.evi, sum.evi,  : 
#>  object 'coords' not found
#>In addition: Warning message:
#>In pd.evi - td.evi[-1] :
#> 
#> Error in data.frame(coords[1], coords[2], lq, uq, mean.evi, sd.evi, sum.evi,  : 
#>  object 'coords' not found > 
#> statsMODIS[p,] <- dopixel(pixel.df,coords,PLOT=FALSE)
#>Error in data.frame(coords[1], coords[2], lq, uq, mean.evi, sd.evi, sum.evi,  : 
#>  object 'coords' not found
#>In addition: Warning message:
#>In pd.evi - td.evi[-1] :
#>  longer object length is not a multiple of shorter object length
#> 
#> print("after")
#>[1] "after"
#> 
#> ### 2nd iteration
#> 
#> p <- 77919
#> 
#> print(p)
#>[1] 77919
#> 
#> #-1-create a data.frame for the pixel
#> pixel.df = data.frame(date = pixel.date)
#> #-2-create an approx function for each ERA variable
#> af.tair = approxfun(ERA.dates,ts.tair[1,])
#> af.tsoil = approxfun(ERA.dates,ts.tsoil[1,])
#> af.srad = approxfun(ERA.dates,ts.srad[1,])
#> af.moist1 = approxfun(ERA.dates,ts.moist1[1,])
#> af.moist2 = approxfun(ERA.dates,ts.moist2[1,])
#> #-3-create an approx function for the EVI data
#> af.evi = approxfun(as.Date(date.pos),ts.evi[p,])
#> #-4-use the approxfun's to fill the data frame data.frame
#> pixel.df$tair <- af.tair(as.Date(pixel.date))
#> pixel.df$tsoil <- af.tsoil(as.Date(pixel.date))
#> pixel.df$srad <- af.srad(as.Date(pixel.date))
#> pixel.df$moist1 <- af.moist1(as.Date(pixel.date))
#> pixel.df$moist2 <- af.moist2(as.Date(pixel.date))
#> pixel.df$evi <- af.evi(as.Date(pixel.date))
#> pixel.df<-pixel.df[complete.cases(pixel.df),]
#> 
#> coords<-round(coordinates(tile.xy)[p,],2)
#>Error in coordinates(tile.xy)[p, ] : subscript out of bounds
#> 
#> print("before")
#>[1] "before"
#> 
#> print(paste(c("statsModis:", typeof(statsMODIS[p,]))))
#>Error in GetRows.bm(x, i) : Illegal row index usage in extraction.
#> print(paste(c("statsModis:", class(statsMODIS[p,]))))
#>Error in GetRows.bm(x, i) : Illegal row index usage in extraction.
#> print(paste(c("statsModis:", dim(statsMODIS[p,]))))
#>Error in GetRows.bm(x, i) : Illegal row index usage in extraction.
#> print(paste(c("statsModis:", length(statsMODIS[p,]))))
#>Error in GetRows.bm(x, i) : Illegal row index usage in extraction.
#> print(paste(c("statsModis:", head(statsMODIS[p,]))))
#>Error in GetRows.bm(x, i) : Illegal row index usage in extraction.
#> 
#> print(paste(c("dopixel:", typeof(dopixel(pixel.df,coords,PLOT=FALSE)))))
#>Error in data.frame(coords[1], coords[2], lq, uq, mean.evi, sd.evi, sum.evi,  : 
#>  object 'coords' not found
#>In addition: Warning message:
#>In pd.evi - td.evi[-1] :
#> 
#> Error in data.frame(coords[1], coords[2], lq, uq, mean.evi, sd.evi, sum.evi,  : 
#>  object 'coords' not found > print(paste(c("dopixel:", class(dopixel(pixel.df,coords,PLOT=FALSE)))))
#>Error in data.frame(coords[1], coords[2], lq, uq, mean.evi, sd.evi, sum.evi,  : 
#>  object 'coords' not found
#>In addition: Warning message:
#>In pd.evi - td.evi[-1] :
#> 
 #>Error in data.frame(coords[1], coords[2], lq, uq, mean.evi, sd.evi, sum.evi,  : 
#>  object 'coords' not found > print(paste(c("dopixel:", dim(dopixel(pixel.df,coords,PLOT=FALSE)))))
#>Error in data.frame(coords[1], coords[2], lq, uq, mean.evi, sd.evi, sum.evi,  : 
#>  object 'coords' not found
#>In addition: Warning message:
#>In pd.evi - td.evi[-1] :
#> 
#> Error in data.frame(coords[1], coords[2], lq, uq, mean.evi, sd.evi, sum.evi,  : 
#>  object 'coords' not found > print(paste(c("dopixel:", length(dopixel(pixel.df,coords,PLOT=FALSE)))))
#>Error in data.frame(coords[1], coords[2], lq, uq, mean.evi, sd.evi, sum.evi,  : 
#>  object 'coords' not found
#>In addition: Warning message:
#>In pd.evi - td.evi[-1] :
#> 
#> Error in data.frame(coords[1], coords[2], lq, uq, mean.evi, sd.evi, sum.evi,  : 
#>  object 'coords' not found > print(paste(c("dopixel:", head(dopixel(pixel.df,coords,PLOT=FALSE)))))
#>Error in data.frame(coords[1], coords[2], lq, uq, mean.evi, sd.evi, sum.evi,  : 
#>  object 'coords' not found
#>In addition: Warning message:
#>In pd.evi - td.evi[-1] :
#> 
#> Error in data.frame(coords[1], coords[2], lq, uq, mean.evi, sd.evi, sum.evi,  : 
#>  object 'coords' not found > 
#> statsMODIS[p,] <- dopixel(pixel.df,coords,PLOT=FALSE)
#>Error in data.frame(coords[1], coords[2], lq, uq, mean.evi, sd.evi, sum.evi,  : 
#>  object 'coords' not found
#>In addition: Warning message:
#>In pd.evi - td.evi[-1] :
#>  longer object length is not a multiple of shorter object length
#> 
#> print("after")
#>[1] "after"
#> dim(statsMODIS)
#>[1] 8099  670