Hi,
I am trying to run a script that stores time series vegetation and reanalysis data based on geographic points in a big matrix. But I have two problems: (1) the code fails at some point (when running the for loop at row 77920) and (2) a reprex of the code displays an error message.
The code successfully works over a smaller study area but I get the following error message when applied over a larger area: Error in SetRows.bm(x, i, value) :
number of items to replace is not a multiple of replacement length
A reprex of the code displays the following error message: Rendering reprex...
Error: This reprex appears to crash R
Call reprex() again with std_out_err = TRUE to get more info
Herewith is the code, it's a bit long...
###This script stores time series vegetation and reanalysis data based on geographic points in a big matrix
#ts.evi is matrix that contains 19-year vegetation data
#tile.xy are geographic/spatial points of length 705184 (dimension 705184 by 457)
#ts.moist1 is a matrix that contains 19-year moisture data at each geographic point (dimensions 705184 by 240)
#ts.moist2 is a matrix that contains 19-year moisture data at each geographic point, but different to ts.moist2 (dimensions 705184 by 240)
#ts.srad is a matrix that contains 19-year solar radiation data at each geographic point (dimensions 705184 by 240)
#ts.tair is a matrix that contains 19-year air temperature data at each geographic point (dimensions 705184 by 240)
#ts.tsoil is a matrix that contains 19-year soil temperature data at each geographic point (dimensions 705184 by 240)
#The vegetation data is available every 16 days of the time series while the reanalysis data is available on a monthly basis
library(reprex)
library(raster)
library(ncdf4)
library(bigmemory)
library(foreach)
library(doMC)
library(doParallel)
library(raster)
library(ncdf4)
library(bigmemory)
library(rgdal)
library(sp)
cores<-12
registerDoMC(cores=cores)
source("/home/muhoko/Biomes_of_Southern_Africa/Tests/Kruger/dopixel_2020_EVI_kruger_negative1_removed.R")#Code used to run dopixel in the for loop
#--open the ERA file extract time
setwd("/home/fbiome/Attribution")
ncERA = nc_open("TTRforceTime.nc")
lon <- ncvar_get(ncERA, "longitude")
lat <- ncvar_get(ncERA, "latitude", verbose = F)
t <- ncvar_get(ncERA, "time")
#--put dates in R date format
ERA.dates <- as.Date(t/24, origin = '1900-01-01', tz="UTC") #240 dates
#load saved objects
ts.evi <- readRDS("/home/muhoko/Biomes_of_Southern_Africa/Protected_areas/Protected_areas_NDVI_EVI/cropped/export/ts.evi.rds")
tile.xy <- readRDS("~/Biomes_of_Southern_Africa/Protected_areas/Protected_areas_NDVI_EVI/Phenological metrics of Southern Africa/Saved_objects/Extracted_data/tile.xy.rds")
ts.moist1 <- readRDS("~/Biomes_of_Southern_Africa/Protected_areas/Protected_areas_NDVI_EVI/Phenological metrics of Southern Africa/Saved_objects/Extracted_data/ts.moist1.rds")
ts.moist2 <- readRDS("~/Biomes_of_Southern_Africa/Protected_areas/Protected_areas_NDVI_EVI/Phenological metrics of Southern Africa/Saved_objects/Extracted_data/ts.moist2.rds")
ts.srad <- readRDS("~/Biomes_of_Southern_Africa/Protected_areas/Protected_areas_NDVI_EVI/Phenological metrics of Southern Africa/Saved_objects/Extracted_data/ts.srad.rds")
ts.tair <- readRDS("~/Biomes_of_Southern_Africa/Protected_areas/Protected_areas_NDVI_EVI/Phenological metrics of Southern Africa/Saved_objects/Extracted_data/ts.tair.rds")
ts.tsoil <- readRDS("~/Biomes_of_Southern_Africa/Protected_areas/Protected_areas_NDVI_EVI/Phenological metrics of Southern Africa/Saved_objects/Extracted_data/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"))
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 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( length(which(is.na(ts.evi[p,]))) < 100 ) # there are circa 440 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))
}
}
'''