How to stratify a raster object in R

Question

My goal in the current work is to sumulate a database containning n variables : (1)latitude of the location, (2)longitude of the location, (3...n-1) covariates, (n)stratum. In stratum variable, I'm looking for the number of the compartment from which each points comes from. But i'm struggling to get the strata. I've created two rasters R1 and R2.

#'Creating the window of domaine D
D <- c(300, 300)  # Square Domaine D of side 300

Win <- owin(xrange =c(0, D[1]), yrange =c(0,D[2])) 

spatstat.options(npixel=c(D[1],D[2]))

# Creating Rasters
library(raster)
# First raster
ext <- extent(Win$xrange, Win$yrange)
v1 <- rnorm(n=seq(from = -10, to=50,length.out = 100), mean = 0, sd=1)
R1 <- raster(matrix(v1,byrow=T,nrow=10, ncol=10))
extent(R1) <- ext
crs(R1) <- "+proj=lcc +lat_1=48 +lat_2=33 +lon_0=-100 +datum=WGS84"
names(R1) <- 'temp sim'
R1
# Second Raster
v2 <- rnorm(n=seq(from = 100, to=2000,length.out = 100), mean = 0, sd=1)
R2 <- raster(matrix(v1,byrow=T,nrow=10, ncol=10))
extent(R2) <- ext
crs(R2) <- "+proj=lcc +lat_1=48 +lat_2=33 +lon_0=-100 +datum=WGS84"
names(R2) <- 'precip sim'
R2
# Raster List
Rasters.group <- stack(R1, R2)

After doing that I've convert the rasters into .im objects

#' Creating .im objects based on the previous rasters. 
#' Each .im object will represent an environmental condition
library(maptools)
IM1 <- as.im.RasterLayer(R1)
IM1 <- as.im.RasterLayer(R2)

#We can plot them as well
library(fields)
par(mfrow=c(2,2))
image.plot(list(x=IM1$xcol, y=IM1$yrow, z=IM1$v), main= "temp sim", asp=1)
image.plot(list(x=IM2$xcol, y=IM2$yrow, z=IM2$v), main= "precip sim", asp=1)

Now I've combined the .im objects in order to generate spacial points from them.

# Combination  of invironments

Com <- IM1 + IM1
image.plot(list(x=Com$xcol, y=Com$yrow, z=Com$v), main= "Global env", asp=1, col=c("green", "red", "blue", "yellow", "White", "Black"))

# Generating points Log - Gaussian Cox process

lgc.process <- rLGCP('matern', mu=Com, var=0.5, scale=0.05, nu=1)
lgc.process

xycoor <- cbind(lgc.process$x, lgc.process$y)

lgc.data <- cbind(xycoor, extract(Rasters.group, xycoor))

The following is the result:

> head(lgc.data)
                           temp.sim precip.sim
[1,] 194.13086 169.76961 1.82615774 1.82615774
[2,]  31.07026 271.56473 0.09453324 0.09453324
[3,] 109.21386 232.53282 2.30704591 2.30704591
[4,]  27.16373 121.88652 0.83505001 0.83505001
[5,]  97.95982 105.71822 1.85318918 1.85318918
[6,] 167.63728  61.16835 1.14231913 1.14231913

Now it remains the last column in which i'm expecting to see from which compatment(stratum) each of the points comes from. I need at least 10 strata. All what I've tried donnot lead me to what i'm looking for. I've used split function, but I'm not getting result.

How should i do?

Answer 1

There are several things in your code that seem strange. E.g. both when you create v1 and v2 you have a special vector for the argument n, but only the length is used, so it is equivalent to:

v1 <- rnorm(100)
v2 <- rnorm(100)

In the beginning you say something about longitude and latitude, but you should be aware that spatstat only works with planar projected coordinates. It is great to go through sf/maptools/raster/stars/terra if you need to work with different projections and convert between them. In your code in the question however you never do any of that and following code is somewhat similar to what you do:

library(spatstat)
Win <- square(300)
IM1 <- as.im(matrix(v1, 10, 10), W = Win)
IM2 <- as.im(matrix(v2, 10, 10), W = Win)
Com <- IM1 + IM1
lgc.process <- rLGCP('matern', mu=Com, var=0.5, scale=0.05, nu=1)
lgc.process # Note this has a lot of points
#> Planar point pattern: 792607 points
#> window: binary image mask
#> 10 x 10 pixel array (ny, nx)
#> enclosing rectangle: [0, 300] x [0, 300] units

Divide domain into nx*ny quadrats (a tessellation) and give them labels (marks)

nx <- 4
ny <- 3
qtess <- quadrats(Win, nx, ny)
marks(qtess) <- 1:(nx*ny)
qtess
#> Tessellation
#> Tiles are equal rectangles, of dimension 75 x 100 units 
#> 4 by 3 grid of tiles
#> Tessellation is marked
#> window: rectangle = [0, 300] x [0, 300] units
plot(qtess) #' marks/labels not visible on this type of plot

Convert into a factor valued image with the mark labels as the value of each region:

qim <- as.im(qtess, dimyx = c(ny, nx))
qim
#> factor-valued pixel image
#> factor levels:
#>  [1] "1"  "2"  "3"  "4"  "5"  "6"  "7"  "8"  "9"  "10" "11" "12"
#> 3 x 4 pixel array (ny, nx)
#> enclosing rectangle: [0, 300] x [0, 300] units
plot(qim, valuesAreColours = FALSE)

Finally extract the values of the image at the locations of the point pattern:

val <- qim[lgc.process]
xycor <- coords(lgc.process)
lgc.data <- cbind(xycor, temp = IM1[lgc.process], precip = IM2[lgc.process], val = val)
head(lgc.data)
#>           x         y      temp     precip val
#> 1 260.93114 167.25626 1.6887990 -0.1590026   8
#> 2 253.33075 242.38444 2.1093809 -0.9992986   4
#> 3  23.57170 244.86029 2.3924167  0.3739248   1
#> 4  14.27467 263.48256 2.3924167  0.3739248   1
#> 5  91.47319  82.87186 1.6035174  1.8282477  10
#> 6 157.93998  26.90952 0.6159271  0.7130826  11