vasudha chaturvedi
vasudha chaturvedi

Reputation: 1

krigeST error with STIDF and STFDF object

I want to produce daily, weekly and monthly kriged mapsvar1.pred is the output I am looking for.I created STIDF object with temporal Units in the variogram model following link for temporal Units. I get an error when running krigeST. Error is: Error in chol.default(A) : the leading minor of order 3 is not positive I know this error is because of duplicates as per issues associated to error. But I need duplicates because one of the parameter varies which is shown with light condition.

For weekly and monthly averages I created STFDF object which has no duplicates but shows another error. Error is:Error in na.fail.default(list(value = c(8.1925753837777, 22.0429913536433, : missing values in object which is because of NA values, but omitting NA values removes all the values from STFDF object.

timelag for STIDF object with empirical variogram shows values such as 15,45 empirical variogram output of STIDF objectwhereas for STFDF object it shows days empirical variogram output of STIDF aggregated to weeks an STFDF object.Why does it mean?

Data Structure

dput(head(time,50))
structure(list(ID = c(783, 784, 845, 846, 907, 908, 966, 1024, 
1081, 1, 2, 3, 61, 62, 63, 121, 122, 123, 181, 182, 183, 237, 
238, 239, 293, 294, 295, 349, 350, 351, 352, 353, 426, 427, 428, 
429, 430, 495, 496, 497, 498, 499, 571, 572, 573, 574, 575, 645, 
646, 647), time = c("11:00:00", "15:00:00", "11:00:00", "15:00:00", 
"11:00:00", "15:00:00", "11:00:00", "11:00:00", "11:00:00", "11:00:00", 
"15:00:00", "19:00:00", "11:00:00", "15:00:00", "19:00:00", "11:00:00", 
"15:00:00", "19:00:00", "11:00:00", "15:00:00", "19:00:00", "11:00:00", 
"15:00:00", "19:00:00", "11:00:00", "15:00:00", "19:00:00", "07:00:00", 
"11:00:00", "15:00:00", "19:00:00", "24:00:00", "07:00:00", "11:00:00", 
"15:00:00", "19:00:00", "24:00:00", "07:00:00", "11:00:00", "15:00:00", 
"19:00:00", "24:00:00", "07:00:00", "11:00:00", "15:00:00", "19:00:00", 
"24:00:00", "07:00:00", "11:00:00", "15:00:00"), date. = structure(c(19544, 
19544, 19544, 19544, 19544, 19544, 19544, 19544, 19544, 19545, 
19545, 19545, 19545, 19545, 19545, 19545, 19545, 19545, 19545, 
19545, 19545, 19545, 19545, 19545, 19545, 19545, 19545, 19547, 
19547, 19547, 19547, 19547, 19547, 19547, 19547, 19547, 19547, 
19547, 19547, 19547, 19547, 19547, 19547, 19547, 19547, 19547, 
19547, 19547, 19547, 19547), class = "Date"), light.condition = c("light", 
"light", "light", "light", "light", "light", "dark", "dark", 
"dark", "light", "light", "light", "light", "light", "light", 
"light", "light", "light", "dark", "dark", "dark", "dark", "dark", 
"dark", "dark", "dark", "dark", "light", "light", "light", "light", 
"light", "light", "light", "light", "light", "light", "light", 
"light", "light", "light", "light", "dark", "dark", "dark", "dark", 
"dark", "dark", "dark", "dark"), value = c(9.71401740645347, 
6.67113336110194, 22.2683202886655, 21.8176624186212, 1.10941911224629, 
1.06500494060907, 8.14032847129926, 21.4166772666492, 1.18754247027769, 
0.0550587583781332, 0.0908023190127035, -0.0287235863761574, 
-0.669267461955569, -0.0999384615406973, -0.467014583537735, 
-0.753432013505277, -0.0399606378578826, -0.365409650446191, 
-0.279498941786479, -0.0769385885223229, -0.252105706894821, 
-0.650782191920685, -0.906986174981546, -1.05863160747964, -0.485452126906809, 
-0.925896081018024, -1.04278402652926, 0.415677741263201, 0.277695366178094, 
0.266598920484192, 0.296489601928946, 0.450302004553194, 0.18961979977428, 
0.186103097004661, 0.181460098328402, 0.219629416931167, 0.201007497364092, 
0.188124540529136, 0.14345259055834, 0.153335645962729, 0.134674763692412, 
0.167406242060629, 0.275140881978198, 0.319478289598817, 0.288869076687038, 
0.342902199465362, 0.347413584337766, 0.158633196058699, 0.181909102924666, 
0.20224766775325), lat = c(68.623823, 68.623823, 68.623877, 68.623877, 
68.623781, 68.623781, 68.623823, 68.623877, 68.623781, 68.623899, 
68.623899, 68.623899, 68.623843, 68.623843, 68.623843, 68.623944, 
68.623944, 68.623944, 68.623899, 68.623899, 68.623899, 68.623843, 
68.623843, 68.623843, 68.623944, 68.623944, 68.623944, 68.624034, 
68.624034, 68.624034, 68.624034, 68.624034, 68.624093, 68.624093, 
68.624093, 68.624093, 68.624093, 68.62398, 68.62398, 68.62398, 
68.62398, 68.62398, 68.624034, 68.624034, 68.624034, 68.624034, 
68.624034, 68.624093, 68.624093, 68.624093), long = c(-149.584395, 
-149.584395, -149.584574, -149.584574, -149.584725, -149.584725, 
-149.584395, -149.584574, -149.584725, -149.583522, -149.583522, 
-149.583522, -149.583303, -149.583303, -149.583303, -149.583483, 
-149.583483, -149.583483, -149.583522, -149.583522, -149.583522, 
-149.583303, -149.583303, -149.583303, -149.583483, -149.583483, 
-149.583483, -149.58195, -149.58195, -149.58195, -149.58195, 
-149.58195, -149.581698, -149.581698, -149.581698, -149.581698, 
-149.581698, -149.581512, -149.581512, -149.581512, -149.581512, 
-149.581512, -149.58195, -149.58195, -149.58195, -149.58195, 
-149.58195, -149.581698, -149.581698, -149.581698), class = "data.frame")

Time object

datetime_str <- paste(time$date., time$time)
starttime <- as.POSIXct(datetime_str, tz = '', tryFormats = "%Y-%m-%d %H:%M:%OS")
time$startime = starttime

STIDF object

grd = data.frame(time)
coordinates(grd) = ~long + lat
proj4string(grd) <- CRS("+init=EPSG:6335")
valueSP = SpatialPoints(grd@coords,CRS("+init=epsg:6335"))
valueDF = data.frame(value=grd$value)
valueTM = grd$startime
stidf = STIDF(valueSP,valueTM,valueDF)
stplot(stidf,number=5,main="random days spatio-temporal values",as.table=T,colorkey=T )

Duplicates zerodist(stidf@sp) duplicates in STIDF object

Variogram Model

metric <- vgmST("metric", joint = vgm(psill=182,model="Per",nugget=30,range=1), stAni=3)
attr(metric, "temporal unit")="days"

Kriging

pred.grd = SpatialPixels(SpatialPoints(makegrid(valueSP,n=69105)),proj4string = proj4string(valueSP))
n=10
tm.grd=xts(1:n,seq(min(index(stidf)),max(index(stidf)),length=n))
prediction.grd = STF(pred.grd,tm.grd)
pred <- krigeST(value~1, data=stidf, newdata=prediction.grd,metric,nmax=20,computeVar=TRUE) 

Aggregate data to weeks/months STFDF object x=aggregate(stidf,"weeks",mean,na.rm=TRUE)

Duplicates for STFDF object zerodist(x@sp) no duplicates for STFDF object

Empirical variogram STFDF object

var_stidf_weeks= variogramST(value ~ 1,data=x,tunit="weeks",tlags=0:5,assumeRegular=F)
sum(is.na(var_stidf_weeks))
var_stidf_weeks1=na.omit(var_stidf_weeks)
var_stidf_weeks1$dist=var_stidf_weeks1$dist*20
var_stidf_weeks1$gamma=var_stidf_weeks1$gamma/3
var_stidf_weeks1$timelag=var_stidf_weeks1$timelag*4
var_stidf_weeks1$avgDist=var_stidf_weeks1$avgDist*20

Fitting variogram

sumMetric <- vgmST("sumMetric", space = vgm(psill=46,model="Per",nugget=1,range=4),
                   time = vgm(psill=0.09,model="Exc",nugget=0.01,range=2),
                   joint = vgm(psill=1,model="Per",nugget=1,range=1), stAni=15)
sumMetric_Vgm <- fit.StVariogram(var_stidf_weeks1,sumMetric,stAni=3)
attr(sumMetric_Vgm, "temporal unit")="weeks"

Kriging prediction

tm.grd=xts(1:n,seq(min(index(x)),max(index(x)),length=n))
pred <- krigeST(value~1, data=x, STF(pred.grd,tm.grd),sumMetric_Vgm) 

Just for comparison empirical variogram of STIDF

var_stidf_daily = variogramST(value~1,data=stidf,tunit="days",tlags=0:5,assumeRegular=F,na.omit=T, cores=4)
var_stidf_daily
sum(is.na(var_stidf_daily))

Upvotes: 0

Views: 45

Answers (0)

Related Questions