Reputation: 1
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