R Under development (unstable) (2023-10-29 r85433 ucrt) -- "Unsuffered Consequences"
Copyright (C) 2023 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> library(intamap)
Loading required package: sp
> library(automap)
> library(gstat)
> set.seed(13531)
> 
> plotFigs = FALSE
> npts = 1000
> pts=SpatialPoints(cbind(runif(npts),runif(npts)))
> d = SpatialPointsDataFrame(cbind(0,0),data.frame(z=1))
> observations = krige(z~1,d,pts,vgm(1, "Sph", .5,anis=c(90,0.5)),nsim=1,nmax=50,beta=0)
[using conditional Gaussian simulation]
> 
> #spplot(sim)
> 
> predictionLocations = spsample(observations, 1000, "regular")
> 
> # We dont know the projection of the data at this stage, assume it is
> # somehow metric
> 
> 
> 
> formulaString = as.formula(sim1~1)
> krigingObject = createIntamapObject(
+ 	observations = observations,
+ 	predictionLocations = predictionLocations,
+   formulaString = formulaString
+ )
> class(krigingObject) = c("automap")
> 
> checkSetup(krigingObject)
Checking object ... OK
> object = preProcess(krigingObject)
> objTemp = estimateAnisotropy(object)
> #rotate Data
> objTemp$observations = rotateAnisotropicData(objTemp$observations, objTemp$anisPar)
> #Estimate Variogram Model
> vario = autofitVariogram(objTemp$formulaString, objTemp$observations, model = "Sph")$var_model
> objTemp$anisPar
$ratio
[1] 1.473469

$direction
[1] 10.79148

$Q
         Q11      Q22       Q12
[1,] 52.7867 108.0044 -10.92161

$doRotation
[1] TRUE

> if (plotFigs) {
+   spplot(object$observations, "sim1", col.regions=bpy.colors())
+   spplot(objTemp$observations, "sim1", col.regions=bpy.colors())
+   plot(variogram(sim1~1, object$observations, alpha=c(0, 90)), vario)
+ }
> print(vario, digits = 3)
  model   psill range
1   Nug 0.00423 0.000
2   Sph 0.76877 0.341
> 
> 
> vmod = vgm(1, "Sph", 1, anis = c(90, 0.5))
> krigingObject$observations = krige(z~1, d, pts, vmod, nsim = 1, nmax = 50, beta = 0)
[using conditional Gaussian simulation]
> object = preProcess(krigingObject)
> objTemp = estimateAnisotropy(object)
> objTemp$anisPar
$ratio
[1] 1.463776

$direction
[1] 0.4183025

$Q
          Q11      Q22       Q12
[1,] 30.88621 66.17211 -0.257632

$doRotation
[1] TRUE

> objTemp$observations = rotateAnisotropicData(objTemp$observations, objTemp$anisPar)
> vario = autofitVariogram(objTemp$formulaString, objTemp$observations, model="Sph")$var_model
> print(vario, digits = 3)
  model psill range
1   Nug 0.000 0.000
2   Sph 0.943 0.671
> vmod
  model psill range ang1 anis1
1   Sph     1     1   90   0.5
> 
> vmod = vgm(1, "Sph", 2, anis=c(45, 0.2))
> krigingObject$observations = krige(z~1, d, pts, vmod, nsim = 1, nmax = 50, beta = 0)
[using conditional Gaussian simulation]
> object = preProcess(krigingObject)
> objTemp = estimateAnisotropy(object)
> objTemp$anisPar
$ratio
[1] 2.932813

$direction
[1] 42.98085

$Q
          Q11      Q22       Q12
[1,] 54.28156 60.69178 -45.39902

$doRotation
[1] TRUE

> objTemp$observations = rotateAnisotropicData(objTemp$observations, objTemp$anisPar)
> vario = autofitVariogram(objTemp$formulaString, objTemp$observations, model = "Sph")$var_model
> print(vario, digits = 3)
  model   psill range
1   Nug 0.00102  0.00
2   Sph 0.89779  1.24
> vmod
  model psill range ang1 anis1
1   Sph     1     2   45   0.2
> 
> 
> vmod = vgm(1,"Sph",1,anis=c(135,0.5))
> krigingObject$observations = krige(z~1,d,pts,vmod,nsim=1,nmax=50,beta=0)
[using conditional Gaussian simulation]
> object = preProcess(krigingObject)
> objTemp=estimateAnisotropy(object)
> objTemp$anisPar
$ratio
[1] 1.370417

$direction
[1] -49.22412

$Q
          Q11      Q22      Q12
[1,] 34.70327 31.72582 10.02324

$doRotation
[1] TRUE

> objTemp$observations=intamap:::rotateAnisotropicData(objTemp$observations,objTemp$anisPar)
> vario = autofitVariogram(objTemp$formulaString,objTemp$observations,model="Sph")$var_model
> print(vario, digits = 3)
  model   psill range
1   Nug 0.00589  0.00
2   Sph 0.63573  0.67
> vmod
  model psill range ang1 anis1
1   Sph     1     1  135   0.5
> 
> 
> 
> vmod = vgm(1,"Sph",.3,anis=c(90,0.5))
> krigingObject$observations = krige(z~1,d,pts,vmod,nsim=1,nmax=100,beta=0)
[using conditional Gaussian simulation]
> object = preProcess(krigingObject)
> objTemp=estimateAnisotropy(object)
> objTemp$anisPar
$ratio
[1] 1.244549

$direction
[1] -10.91575

$Q
          Q11      Q22      Q12
[1,] 92.11681 138.1476 9.220176

$doRotation
[1] TRUE

> objTemp$observations=intamap:::rotateAnisotropicData(objTemp$observations,objTemp$anisPar)
> vario = autofitVariogram(objTemp$formulaString,objTemp$observations,model="Sph")$var_model
> print(vario, digits = 3)
  model psill range
1   Nug 0.000 0.000
2   Sph 0.805 0.176
> vmod
  model psill range ang1 anis1
1   Sph     1   0.3   90   0.5
> if (plotFigs) {
+   p1 = plot(variogram(sim1~1,object$observations,alpha=c(0,90)),vmod,ylim = c(0,1.2),xlim=c(0,0.6),main="orig,orig")
+   p2 = plot(variogram(sim1~1,object$observations,alpha=c(0,90)),vario,ylim = c(0,1.2),xlim=c(0,0.6),main="orig,fitted")
+   p3 = plot(variogram(sim1~1,objTemp$observations,alpha=c(0,90)),vmod,ylim = c(0,1.2),xlim=c(0,0.6),main="rot,orig")
+   p4 = plot(variogram(sim1~1,objTemp$observations,alpha=c(0,90)),vario,ylim = c(0,1.2),xlim=c(0,0.6),main = "rot,fitted")
+ 
+   print(p1,position = c(0,0.5,0.5,1),more = TRUE)
+   print(p2,position = c(0.5,0.5,1,1),more = TRUE)
+   print(p3,position = c(0,0,0.5,0.5),more = TRUE)
+   print(p4,position = c(0.5,0,1,0.5))
+ 
+   plot(variogram(sim1~1,objTemp$observations,alpha=c(seq(0,150,30))),vmod,ylim = c(0,1.2),xlim=c(0,0.6),main="rot,orig")
+ 
+   spplot(object$observations,"sim1",col.regions=bpy.colors())
+   spplot(objTemp$observations,"sim1",col.regions=bpy.colors())
+ }
> 
> 
> 
> 
> data(sic2004)
> coordinates(sic.val)=~x+y
> sic.val$value=sic.val$dayx
> x = sic.test$x
> y = sic.test$y
> 
> coordinates(sic.test)=~x+y
> 
> stest = sic.test[(x > -10000 & x < 140000 & y > 100000 & y < 240000),]
>  
> obj<-createIntamapObject(formulaString = "joker~1",
+    observations=sic.val,
+    predictionLocations = stest,
+    params = list(doAnisotropy = TRUE),
+    class = "automap" )
> obj = preProcess(obj)
> obj = estimateParameters(obj)
> obj$anisPar
$ratio
[1] 1.350924

$direction
[1] 25.93149

$Q
              Q11         Q22           Q12
[1,] 9.948346e-06 1.43261e-05 -2.787868e-06

$doRotation
[1] TRUE

> obj = spatialPredict(obj)
> obj = postProcess(obj)
> summary(as.data.frame(obj$outputTable), digits = 3)
       x                y               mean        variance    
 Min.   : -8351   Min.   :103106   Min.   :121   Min.   : 4535  
 1st Qu.: 20626   1st Qu.:136147   1st Qu.:146   1st Qu.:19001  
 Median : 63885   Median :167801   Median :153   Median :19935  
 Mean   : 59128   Mean   :170171   Mean   :173   Mean   :18875  
 3rd Qu.: 94458   3rd Qu.:210897   3rd Qu.:157   3rd Qu.:20282  
 Max.   :131632   Max.   :239246   Max.   :544   Max.   :20320  
> 
> proc.time()
   user  system elapsed 
  10.84    0.93   11.73