Simulations on Sphere
D. Renard
17 août 2018
Introduction
rm(list=ls())
library(RGeostats)
library(rgl)
library(knitr)
knit_hooks$set(webgl = hook_webgl)
variety.define(flag.sphere=TRUE)##
## Parameters for Variety Definition
## ---------------------------------
## The Spherical Variety is defined
## - Radius of the Sphere = 6371.000000constant.define("nproc",NA)This file is meant to demonstrate the procedure used to perform Simulations on the Sphere. This function requires the definition of the change of space; the resulting file corresponds to a regular grid file, where nodes are regularly located along longitude and latitude regular angles. Finally for commodity, the visualization is performed using RGL library.
Preliminary tasks
First we need to create the output grid where the nodes vary along the longitude from 0 to 360 by mesh of 1 degree (hence 361 nodes); and from -90 to +90 degrees in latitude by mesh of 1 degree (hence 181 meshes).
grid = db.create(nx=c(361,181),x0=c(0,90),dx=c(1,1))For the simulations, we use a Spatial isotropic Model which is defined below:
model = model.create(vartype="cubic",range=0.5)Some parameters are defined and will be kept constant for all tests (unless overwritten). This is the case in particular for:
- The number of basic functions nbfdef
nbfdef = 10000In general, this function enables simulating any covariance model (defined on the Sphere). It requires the preliminary calculation of the covariance spectrum. Two specific constructions have been coded: the spectrum of Chentsov Covariance and the one for Exponential Covariance. In addition any standard 3-D model (defined in the Euclidean space) can be converted into the corresponding Model on the Sphere using the following Yadrenko formula:
\[ C(\alpha) = C_3 \left(2 \sin \frac{\alpha}{2} \right) for \, 0<= \alpha <= \pi \]
Checking Spectrum
In a first step, we will pay attention on the generated spectrum of the covariance. In particular, we make sure that the coefficients of the spectrum are all defined, positive, monotonously decreasing. Moreover the series are truncated if the coefficients are constantly smaller than a given tolerance and as soon as there sum is reasonably close to one. However, their number is limited to a maximum number.
Unless we are using a specific construction mode, the coefficients of the spectrum are calculated as follows. We start with the Schoenberg formula:
\[ C(\alpha) = \sum_{n=0}^{\infty} a_n P_n \left( cos \alpha \right) \]
Using the orthogonality of the \(P_n\) polynomials, we derive:
\[ a_n = \int_{0}^{\pi} C(\alpha) \tilde{P}_n(cos \alpha) sin \alpha \sqrt{\frac{2n+1}{2}} d \alpha \]
- Particular case of Chentsov Covariance:
set.keypair("Simsph_Shunt",1)
dummy = simsph(grid,model,special=1,verbose=TRUE)##
## Generation of Covariance spectrum in Spherical Coordinates
## ----------------------------------------------------------
## Using Chentsov construction
## Spectrum Tolerance = 1e-05
## Random generation seed = 232131
## Number of frequencies = 10000
## Number of basic functions = 100
## Cumulated Spectrum = 0.999936
## Sum of negative weights = 0.000000freqs = get.keypair("Simsph_Spectrum_Frequencies")
plot(freqs, main="Chentsov Covariance Spectrum",
type="l",ylab="Frequency")
plot(freqs, type="b", xlim=c(1,20),
main="Chentsov Covariance Spectrum (zoom)",ylab="Frequency")
- Particular case of Exponential covariance. The spectrum depends upon the range of the structure: this range is inherited from the range of the first basic structure contained in the Model (passed as argument).
set.keypair("Simsph_Shunt",1)
dummy = simsph(grid,model,special=2,verbose=TRUE)##
## Generation of Covariance spectrum in Spherical Coordinates
## ----------------------------------------------------------
## Using Exponential construction
## Spectrum Tolerance = 1e-05
## Random generation seed = 232131
## Number of frequencies = 10000
## Number of basic functions = 100
## Cumulated Spectrum = 0.999800
## Sum of negative weights = 0.000000freqs = get.keypair("Simsph_Spectrum_Frequencies")
plot(freqs, main="Exponential Covariance Spectrum",
type="l",ylab="Frequency")
plot(freqs, type="b", xlim=c(1,20),
main="Exponential Covariance Spectrum (zoom)",ylab="Frequency")
- Standard case with the current Model definition:
set.keypair("Simsph_Shunt",1)
dummy = simsph(grid,model,special=0,verbose=TRUE)##
## Generation of Covariance spectrum in Spherical Coordinates
## ----------------------------------------------------------
## Number of discretization = 360
## Spectrum Tolerance = 1e-05
## Random generation seed = 232131
## Number of frequencies = 41
## Number of basic functions = 100
## Cumulated Spectrum = 0.999567
## Sum of negative weights = 0.000000freqs = get.keypair("Simsph_Spectrum_Frequencies")
plot(freqs, main="Standard Model Spectrum",type="l",ylab="Frequency")
plot(freqs, type="b", xlim=c(1,50),
main="Standard Model Spectrum (zoom)",ylab="Frequency")
Checking the Random draw of Basic Functions
Then a number of basic functions is generated (their number is equal to the nbf argument). Each basic structure is characterized by:
- its degree \(N\) according to the probability law induced by the spectrum
- its order \(\tilde{P}\) which must lie in the interval \([-N;+N]\)
- its phase
Given the spectrum, we can check the distribution of all these parameters. We plot the experimental histogram of the generated degree values (in black) and overlay the theoretical discretized spectrum (in red).
For the degree
This is perform ed using the internal function display.spectrum.degrees.
display.spectrum.degrees <- function(special=0,model,nbf=nbfdef,
dslice=20,title="Undefined")
{
set.keypair("Simsph_Shunt",2)
dummy = simsph(grid,model,special=special,verbose=FALSE,nbf=nbf)
degrees = get.keypair("Simsph_Array_Degree")
freqs = get.keypair("Simsph_Spectrum_Frequencies")
fmax = range(freqs)[2]
dmax = range(degrees)[2]
h.result = hist(degrees,breaks=seq(-0.5,0.5+dmax),plot=FALSE)
nval = length(h.result$mids)
plot(h.result$mids,h.result$density,xlim=c(0,dslice),ylim=c(0,fmax),
type="b",col="black",pch=19,
xlab="Degree N",ylab="Frequency",main=title)
lines(x=h.result$mids,y=freqs[1:nval],add=T,
type="b",col="red",pch=19)
set.keypair("Simsph_Shunt",0)
invisible()
}The following tests are performed for different models:
- for the Chentsov Covariance
display.spectrum.degrees(special=1,model,
title="Chentsov Covariance Spectrum")
- for the Exponential covariance
display.spectrum.degrees(special=2,model,
title="Exponential Covariance Spectrum")
- for the standard Model
display.spectrum.degrees(special=0,model,dslice=40,
title="Standard Model Spectrum")
For the order
This is perform ed using the internal function display.spectrum.order.
display.spectrum.order <- function(special=0,model,nbf=nbfdef,
title="Undefined")
{
set.keypair("Simsph_Shunt",2)
dummy = simsph(grid,model,special=special,verbose=TRUE,nbf=nbf)
degrees = get.keypair("Simsph_Array_Degree")
orders = get.keypair("Simsph_Array_Order")
freqs = get.keypair("Simsph_Spectrum_Frequencies")
dmax = range(degrees)[2]
fmax = range(freqs)[2]
h.result = hist(orders,breaks=seq(-0.5-dmax,0.5+dmax),plot=FALSE)
nval = floor(length(h.result$mids) / 2) + 1
plot(h.result$mids,h.result$density,xlim=c(-20,20),
type="b",col="black",pch=19,
xlab="Order K",ylab="Frequency",main=title)
set.keypair("Simsph_Shunt",0)
invisible()
}The following tests are performed for different models:
- for the Chentsov Covariance
display.spectrum.order(special=1,model,
title="Chentsov Covariance Spectrum")##
## Generation of Covariance spectrum in Spherical Coordinates
## ----------------------------------------------------------
## Using Chentsov construction
## Spectrum Tolerance = 1e-05
## Random generation seed = 232131
## Number of frequencies = 10000
## Number of basic functions = 10000
## Cumulated Spectrum = 0.999936
## Sum of negative weights = 0.000000
## Maximum degree = 5747
## Minimum order = -2218
## Maximum order = 2742
- for the Exponential Covariance
display.spectrum.order(special=2,model,
title="Exponential Covariance Spectrum")##
## Generation of Covariance spectrum in Spherical Coordinates
## ----------------------------------------------------------
## Using Exponential construction
## Spectrum Tolerance = 1e-05
## Random generation seed = 232131
## Number of frequencies = 10000
## Number of basic functions = 10000
## Cumulated Spectrum = 0.999800
## Sum of negative weights = 0.000000
## Maximum degree = 8094
## Minimum order = -3804
## Maximum order = 3862
- for the standard Model
display.spectrum.order(special=0,model,
title="Standard Model Spectrum")##
## Generation of Covariance spectrum in Spherical Coordinates
## ----------------------------------------------------------
## Number of discretization = 360
## Spectrum Tolerance = 1e-05
## Random generation seed = 232131
## Number of frequencies = 41
## Number of basic functions = 10000
## Cumulated Spectrum = 0.999567
## Sum of negative weights = 0.000000
## Maximum degree = 40
## Minimum order = -34
## Maximum order = 39
For the phase
This is perform ed using the internal function display.spectrum.phase.
display.spectrum.phase <- function(special=0,model,nbf=nbfdef,
title="Undefined")
{
set.keypair("Simsph_Shunt",2)
dummy = simsph(grid,model,special=special,verbose=FALSE,nbf=nbf)
phases = get.keypair("Simsph_Array_Phase")
h.result = hist(phases,breaks=100,
main="Histogram of phases",xlab="")
set.keypair("Simsph_Shunt",0)
invisible()
}The simulation of phase is independent from the type of covariance; therefore only one type is used (i.e. Chentsov covariance). The histogram is build with 100 classes.
display.spectrum.phase(special=1,model,
title="Chentsov Covariance Spectrum")
Checking the Legendre Functions
The simulation of the basic function (rank j) over the sphere is obtained through the following formula:
\[ Z_{j}(\theta,\phi) = 2 \tilde{P}_{N_j,K_j} \left(cos \theta \right) cos \left( K \phi + \psi_j \right) \]
where:
- \({\theta, \phi}\) stand for the longitude and latitude angles
- \(N_j\) is an integer number drawn at random according to the covariance spectrum
- \(K_j\) is an integer number drawn according to a uniform law within \([-N,+N]\)
- \(\psi_j\) is a random phase drawn within \([0,2\pi]\)
- \(P_{N,K}(x)\) is a normalized associated Legendre function of degree \(N\) and order \(K\)
Finally, the outcome of the GRF at any point on the sphere (characterized by its longitude and latitude) is obtained as follows:
\[ Z(\theta, \phi) = \frac {\sum_{j}^{n} Z_j(\theta,\phi)}{\sqrt{n}}\]
In this paragraph, we will check the Legendre Function of various degrees and orders. The result should be represented graphically. Unfortunately, the output according to RGL interface does not allow including the results in a standard PDF or HTML file. Therefore the next chunk should be operated several time, in a trial-and-error manner, with different values for the testing parameters.
a = simsph(grid,model,special=1,flag.test=TRUE,
test.degree=1,test.order=0,test.phase=0)
#plot(a)Simulations on the Sphere
In this paragraph, all the plotting activities are commented out at it does not seem to be feasible to hook Rgl library within knitr procedure of Rstudio.
Now that all the elements have been tested, it is time to perform the simulation on sphere. We start with a Cubic Model with a range of 0.5.
model = model.create(vartype="cubic",range=0.5,sill=1)
a = simsph(grid,model,nbf=1000)
#plot(a)We also use a periodic Model such as the Matern one with a range of 0.03 and a third parameter set to 0.01.
model = model.create(vartype="J-Bessel",range=0.03,param=0.01)
a = simsph(grid,model,nbf=1000)
#plot(a)