R Under development (unstable) (2024-02-19 r85946 ucrt) -- "Unsuffered Consequences"
Copyright (C) 2024 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("geostatsp")
Loading required package: Matrix
Loading required package: terra
terra 1.7.71
> 
> param = c(range=1, shape=1.5,	anisoRatio=2, anisoAngleDegrees=-25)
> 
> matern(c(0, 0.001, 100000), param=param)
[1] 1.000000 0.999994 0.000000
attr(,"param")
            range             shape          variance            nugget 
        1.0000000         1.5000000         1.0000000         0.0000000 
       anisoRatio anisoAngleRadians anisoAngleDegrees 
        2.0000000        -0.4363323       -25.0000000 
> 
> #x=c(0, 0.001, 100000);param=c(param, variance=1)
> 
> #resultFull = .C("matern", as.double(x), as.integer(length(x)),
> #		as.double(param["range"]), as.double(param["shape"]),
> #		as.double(param["variance"]))
> 
> 
> # example with raster
> myraster = rast(nrows=40,ncols=60,xmin=-3,xmax=3,ymin=-2,ymax=2)
> 
> # plot correlation of each cell with the origin
> myMatern = matern(myraster, y=c(0,0), param=param)
> myMatern[1:3,1:3]
        lyr.1
1 0.002385719
2 0.003084067
3 0.003973211
4 0.002349589
5 0.003047536
6 0.003940191
7 0.002302637
8 0.002996124
9 0.003886871
>  
> 
> 
> bob = function(x) {
+ thepar = attributes(x)$param
+ pdf(tempfile("matern", tmpdir=".", fileext=".pdf"))
+ plot(x, main=
+ 				paste(
+ 				paste(names(thepar), thepar, sep="="),
+ 				collapse=", "),cex.main=0.5
+ )
+ dev.off()
+ }
> 
> bob(myMatern)
null device 
          1 
> 
> 
> bob(matern(myraster, y=c(1,-0.5), 
+ 				param =	c(range=1, shape=1.5,	anisoRatio=2, anisoAngleDegrees=25)			
+ 						)
+ )
null device 
          1 
> 
> bob(matern(myraster,y= c(0,0), 
+ 				param =	c(range=1, shape=25.1,	anisoRatio=2, anisoAngleDegrees=-25)			
+ 		)
+ )
null device 
          1 
> 
> bob(matern(myraster,y= c(0,0), 
+ 				param =	c(range=0, shape=1.5,	anisoRatio=2, anisoAngleDegrees=-25)			
+ 		)
+ )
null device 
          1 
> 
> bob(matern(myraster, y=c(0,0), 
+ 				param =	c(range=100000, shape=1.5,	anisoRatio=2, anisoAngleDegrees=-25)			
+ 		)
+ )
null device 
          1 
> 
> 
> 				
> 				
> 
> # correlation matrix for all cells with each other
> myraster = rast(nrows=4,ncols=6,xmin=-3,xmax=3,ymin=-2,ymax=2)
> myMatern = matern(myraster, param=c(range=0, shape=2))
> 
> dim(myMatern)
[1] 24 24
> myMatern[1:3,1:3]
3 x 3 Matrix of class "dsyMatrix"
     [,1] [,2] [,3]
[1,]    1    0    0
[2,]    0    1    0
[3,]    0    0    1
> 
> param = c(range=0.2, shape=1.5)
> set.seed(0)
> 
> mypoints = vect(cbind(runif(10), runif(10)), 
+ 		crs = "epsg:2000",
+ 		atts=data.frame(id=1:10))
> 
> myDist = forceSymmetric(as.matrix(distance(mypoints)))
> 
> myDist3 = myDist[lower.tri(myDist)]
> 
> myMatern1 = matern(mypoints, param)
> myMatern2 = matern(myDist, param)
> myMatern3 = matrix(NA, length(mypoints),length(mypoints))
> myMatern3[lower.tri(myMatern3)] = matern(myDist3, param)
> 
> 
> 
> class(myMatern1)
[1] "dsyMatrix"
attr(,"package")
[1] "Matrix"
> class(myMatern2)
[1] "dsyMatrix"
attr(,"package")
[1] "Matrix"
> class(myMatern3)
[1] "matrix" "array" 
> 
> 
> (myMatern1 - myMatern2)[1:8, 1:4]
8 x 4 Matrix of class "dgeMatrix"
              [,1]         [,2]         [,3]         [,4]
[1,]  0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
[2,]  0.000000e+00 1.355887e-08 1.491364e-07 1.096928e-08
[3,]  0.000000e+00 1.491364e-07 9.684880e-02 1.787965e-04
[4,]  0.000000e+00 1.096928e-08 1.787965e-04 3.393112e-07
[5,] -4.163336e-17 4.066149e-06 2.550393e-05 2.008934e-06
[6,]  0.000000e+00 9.514759e-11 2.488279e-04 3.295310e-07
[7,]  0.000000e+00 7.405549e-07 1.106377e-05 4.935520e-07
[8,]  0.000000e+00 2.307813e-08 1.273513e-06 1.863947e-08
> sum((myMatern1 - myMatern2)^2)
[1] 0.01802854
> (myMatern1 - myMatern3)[1:8, 1:4]
8 x 4 Matrix of class "dgeMatrix"
              [,1]         [,2] [,3]          [,4]
[1,]            NA           NA   NA            NA
[2,]  0.000000e+00           NA   NA            NA
[3,]  0.000000e+00 1.665335e-16   NA            NA
[4,]  0.000000e+00 0.000000e+00    0            NA
[5,] -4.163336e-17 0.000000e+00    0  0.000000e+00
[6,]  0.000000e+00 0.000000e+00    0 -6.938894e-18
[7,]  0.000000e+00 0.000000e+00    0  0.000000e+00
[8,]  0.000000e+00 0.000000e+00    0  0.000000e+00
> sum((myMatern1 - myMatern3)^2, na.rm=TRUE)
[1] 2.951492e-32
> 
> 
> 
> proc.time()
   user  system elapsed 
   4.20    0.20    4.35