Loading libraries
rm(list = ls())
library("gstlearn")
##
## Attaching package: 'gstlearn'
## The following objects are masked from 'package:base':
##
## message, toString
library("tidyr") # essentially, for drop_na function
Setting global variables
verbose = TRUE
flagGraphic = TRUE
# space dimension
ndim = 2
defineDefaultSpace(ESpaceType_RN(),ndim)
## NULL
Setting seed for the random number generator
law_set_random_seed(5584)
## NULL
set.seed(13155)
A Poisson data set will be generated and various neighborhoods will be explored around a specific node of a regular grid
dxref = 0.1
grid = DbGrid_create(nx = c(10,10),
dx=c(dxref, dxref))
xlim = grid$getExtrema(0)
ylim = grid$getExtrema(1)
Setting parameters for the dbgrid hosting the data
coormin = grid$getCoorMinimum()
coormax = grid$getCoorMaximum()
nech = 100
data = Db_createFromBox(nech, coormin, coormax)
Plotting the data (if needed for more control over the graph parameters, use ggplot instead)
ggDefaultGeographic()+
plot.point(data, flagCst=True)+
plot.decoration(xlab="Easting", ylab="Northing")
Defining a standard moving neighborhood
nmini = 1
nmaxi = 15
radius = 0.3
nsect = 8
nsmax = 3
neigh = NeighMoving_create(flag_xvalid = FALSE,
nmaxi = nmaxi,
radius = radius,
nmini = nmini,
nsect = nsect,
nsmax = nsmax)
neigh
##
## Moving Neighborhood
## ===================
## Minimum number of samples = 1
## Maximum number of samples = 15
## Number of angular sectors = 8
## Maximum number of points per sector = 3
## Maximum horizontal distance = 0.3
Attaching the neighborhood to the grid and selecting the neighbors retained
node = 55
neigh$attach(data, grid)
## [1] 0
ranks = VectorInt()
neigh$select(node, ranks)
## NULL
dataSel = data$clone()
dum = dataSel$addSelectionByRanks(ranks)
Plotting the result
ggDefaultGeographic()+
plot.point(data,
color = "red")+
plot.point(dataSel,
color = "blue") +
plot.neigh(neigh = neigh,
grid = grid,
node = node,
flagCell = FALSE)+
plot.decoration(title="Standard Neighborhood",
xlab="Easting",
ylab="Northing")
## Warning in names(guides$guides)[to_change] <- paste0(names(guides$guides), :
## number of items to replace is not a multiple of replacement length
In the following chunk, variables are generated in the grid containing the cell extension. This cell extension replaces the constant value of the mesh. The new values define blocks around each sample point that will be represented below
# defining the parameters
nech = grid$getSampleNumber()
mini = 0.5
maxi = 2.5
blx = numeric()
bly = numeric()
# adding the random values
for (i in 1:nech){
blx = append(blx, runif(n = 1, min = mini, max = maxi))
bly = append(bly, runif(n = 1, min = mini, max = maxi))
}
blx = dxref*blx
bly = dxref*bly
# incorporating the columns into the data base :
dum = grid$addColumns(blx, "X-ext", ELoc_BLEX(), 0)
dum = grid$addColumns(bly, "Y-ext", ELoc_BLEX(), 1)
Displaying each block and its cell extension Until the next gstlearn version release, this needs to be done manually thanks to ggplot
df = data[]
colnames(df) = c("rank", "x1", "x2")
extensions = data.frame(x = NA, y = NA, node = NA)
for (node in 1:nech){
edges = data.frame(grid$getCellEdges(node-1))
edges = cbind(edges,NA)
colnames(edges) = c("x","y", "node")
extensions = rbind(extensions, edges[1:4,])
extensions$node[(4*node-2):(4*node+1)] = node
}
extensions = drop_na(extensions)
p = ggplot()+
geom_point(data=df,
aes(x = x1, y = x2),
color = "black",
size = 1)+
geom_polygon(data=extensions,
aes(x=x, y=y, group=node),
color = "black",
linewidth = 0.2,
fill = NA)
p
Choosing a specific cell again and determining the standard block neighborhood
node = 56
neigh$attach(data, grid)
## [1] 0
ranks = VectorInt()
neigh$select(node, ranks)
## NULL
dataSel = data$clone()
dum = dataSel$addSelectionByRanks(ranks)
Plotting (TODO: this currently generates a warning message)
ggDefaultGeographic()+
plot.point(data,
color = "red")+
plot.point(dataSel,
color = 'blue')+
plot.neigh(neigh = neigh,
grid = grid,
node = node,
flagCell=TRUE)+
plot.decoration(title="Standard Neighborhood")
## Warning in names(guides$guides)[to_change] <- paste0(names(guides$guides), :
## number of items to replace is not a multiple of replacement length
Using the cell neighborhood to force the selection of all samples belonging to the block
nmini = 1
neigh = NeighCell_create(FALSE, nmini)
neigh
##
## Cell Neighborhood
## =================
## Reject samples which do not belong to target Block
Attaching the neighborhood to the grid and selecting the neighbors retained
node = 56
neigh$attach(data, grid)
## [1] 0
ranks = VectorInt()
neigh$select(node, ranks)
## NULL
dataSel = data$clone()
dum = dataSel$addSelectionByRanks(ranks)
Plotting the result
ggDefaultGeographic()+
plot.point(data,
color = "red")+
plot.point(dataSel,
color='blue')+
plot.neigh(neigh,
grid,
node,
flagCell=TRUE)+
plot.decoration(title="Standard Neighborhood")
## Warning in names(guides$guides)[to_change] <- paste0(names(guides$guides), :
## number of items to replace is not a multiple of replacement length
This part is dedicated to data organized along lines, under the assumption that the typical distance between two consecutive flight lines is much larger than the typical distance between two observation along a line. To illustrate this case, data will be generated at the nodes of a regular grid with dx>>dy. The file is converted into a point file to make the rest of the procedure more flexible
Creating the grid and the point file containing the data
grid = DbGrid_create(nx=c(20,100), dx=c(10,1))
iuid = grid$addColumns(VectorHelper_simulateUniform(grid$getSampleNumber()),"z", ELoc_Z())
point = Db_createReduce(grid)
ggDefaultGeographic()+
plot.point(grid,size=0.3)