References on simulating a raster with spatial dependence

I've simulated an N-by-N raster in the following way:

1. define a set $S$ containing a finite number $|S| = K$ of possible raster values (in my simulation, $K=3$ and the elements of $S$ are land uses / land covers)
2. Initialize an empty $N$-by-$N$ array of raster values; draw the $1, 1$ entry of the raster uniformly at random from $S$. I'll refer to the raster entries as $X_{i,j}$. At this point we have drawn a value for $X_{1, 1}$
3. Populate the first row of the raster using a conditional distribution $\Pr[X_{i,j} \,|\, X_{i,j-1}]$, i.e. populate $X_{1, j}$ for $j=2, 3, \ldots, N$
4. Use the same conditional distribution as above to populate the first column of the raster (i.e. $\Pr[X_{i,j} \,|\, X_{i-1,j}]$ is the same distribution as $\Pr[X_{i,j} \,|\, X_{i,j-1}]$); this allows us to populate $X_{i, 1}$ for $i=2, 3, \ldots, N$
5. Populate the remainder of the raster using a distribution that conditions on two neighbors: $\Pr[X_{i,j} \,|\, X_{i,j-1}, X_{i-1,j}]$. I assume that this distribution satisfies a symmetry property: the probability does not change if you switch the values of the two neighbors, i.e. $\Pr[X_{i,j} \,|\, X_{i,j-1} = r, X_{i-1,j} = s]$ = $\Pr[X_{i,j} \,|\, X_{i,j-1} = s, X_{i-1,j} = r]$, for any $s, r \in S$.

Those steps generate an $N$-by-$N$ raster whose cells take values in $S$.

I have several questions:

• The primitives of the model are the distributions $\Pr[X_{i,j} \,|\, X_{i,j-1}]$ and $\Pr[X_{i,j} \,|\, X_{i,j-1}, X_{i-1,j}]$. Do these objects imply any sort of stationary distribution over raster values $X_{i, j}$, and over the joint distribution of $X_{i, j}$ and its neighbors to the north, south, east, and west, i.e. the joint distribution of $\left(X_{i, j}, X_{i+1, j}, X_{i, j+1}, X_{i-1, j}, X_{i, j-1}\right)$?
• If the primitives do imply a stationary distribution, can the stationary distribution be calculated using the primitives (the same way we can calculate the stationary distribution of an one-dimensional, first-order Markov chain using an eigenvector decomposition of its transition matrix)?
• If there is a stationary joint distribution of $X_{i, j}$ and its neighbors to the north, south, east, and west, is that distribution "invariant to rotation", by which I mean that rotations of the neighbors around the central pixel do not change the probability? I suspect the answer is no: see the simulated raster below, which shows a sort of "streak" pattern from the 1, 1 pixel towards the top right (i.e., streaks moving in the direction in which I populated the raster in step 5).

Are there any references that answer these questions? I think what I'm simulating is a Markov random field on a grid, but I'm not sure what the proper terminology is.

R simulation:

library(ggplot2)

## Spatial dependence: distribution of ij raster pixel depends on its neighbors to north and west

S_sorted <- sort(c("crops", "forest", "pasture"))  # Set of land cover states
stopifnot(length(S_sorted) == 3)  # Assumed throughout

## Matrix with ij entry Pr[ S = s_j | S^{west} = s_i ]
## Symmetry assumption: this matrix is also used for Pr[ S = s_j | S^{north} = s_i ]
conditional_one_neighbor <- rbind(c(0.80, 0.10, 0.10),
c(0.05, 0.90, 0.05),
c(0.20, 0.10, 0.70))
stopifnot(isTRUE(all.equal(rowSums(conditional_one_neighbor),
rep(1.0, nrow(conditional_one_neighbor)))))  # Valid conditional probabilities

## Array with ijk entry Pr[ S = s_k | S^{north} = s_i, S^{west} = s_j ]
## Symmetry assumption: ijk entry also gives Pr[ S = s_k | S^{north} = s_j, S^{west} = s_i ],
## i.e. ijk entry and jik entry are always identical
conditional_two_neighbors <- array(0, dim=c(3, 3, 3))
conditional_two_neighbors[, , 1] <- rbind(c(0.95, 0.49, 0.29),
c(0.49, 0.01, 0.01),
c(0.29, 0.01, 0.02))  # Probability of crops
conditional_two_neighbors[, , 2] <- rbind(c(0.04, 0.50, 0.01),
c(0.50, 0.95, 0.40),
c(0.01, 0.40, 0.08))  # Probability of forest
conditional_two_neighbors[, , 3] <- rbind(c(0.01, 0.01, 0.70),
c(0.01, 0.04, 0.59),
c(0.70, 0.59, 0.90))  # Probability of pasture
for(k in seq_along(S_sorted)) {
stopifnot(all(conditional_two_neighbors[, , k] == t(conditional_two_neighbors[, , k])))  # Symmetry assumption
}

## Sum over ij should be 1, i.e. valid probability distribution
stopifnot(isTRUE(all.equal(apply(conditional_two_neighbors, MARGIN=c(1, 2), FUN=sum),
matrix(1, dim(conditional_two_neighbors), dim(conditional_two_neighbors)))))

nrow_raster <- 500
raster <- matrix(NA, nrow_raster, nrow_raster)  # Simulate a square land cover raster
raster[1, 1] <- sample(S_sorted, size=1)  # First pixel chosen uniformly at random

## Populate (rest of) first row using conditional_one_neighbor
for(i in seq(2, nrow_raster)) {
neighbor_index <- which(S_sorted == raster[1, i-1])  # Index into S_sorted
raster[1, i] <- sample(S_sorted, size=1, prob=conditional_one_neighbor[neighbor_index, ])
}

## Populate (rest of) first column using conditional_one_neighbor
for(i in seq(2, nrow_raster)) {
neighbor_index <- which(S_sorted == raster[i-1, 1])  # Index into S_sorted
raster[i, 1] <- sample(S_sorted, size=1, prob=conditional_one_neighbor[neighbor_index, ])
}

## Populate rest of raster using conditional_two_neighbors (conditional on north and west neighbors)
for(i in seq(2, nrow_raster)) {
for(j in seq(2, nrow_raster)) {
neighbor_north_index <- which(S_sorted == raster[i-1, j])
neighbor_west_index <- which(S_sorted == raster[i, j-1])
raster[i, j] <- sample(S_sorted,
size=1,
prob=conditional_two_neighbors[neighbor_north_index, neighbor_west_index, ])
}
}

raster_df <- expand.grid(x=seq_len(nrow_raster),
y=seq_len(nrow_raster))
raster_df$land_cover=as.vector(raster) dim(raster_df) ## Raster origin (1, 1) is at bottom-left p <- (ggplot(raster_df, aes(x=x, y=y, fill=land_cover)) + geom_raster() + scale_fill_manual("land cover", values=c("crops"="#ff7f00", "forest"="#4daf4a", "pasture"="#a6cee3"))) p ggsave("simulated_raster.png", width=12, height=10) ## Convergence to stationary distribution as we move away from the origin, to the top-right? table(raster_df$land_cover) / nrow(raster_df)  # Empirical marginal distribution

## Distribution conditional on one neighbor (empirical as above)
tables_list <- lapply(seq_len(nrow_raster), function(i) {
return(table(raster[i, seq(1, nrow_raster-1)], raster[i, seq(2, nrow_raster)]))
})
tables_sum <- Reduce("+", tables_list)
joint_distrib <- tables_sum / sum(tables_sum)
empirical_conditional_one_neighbor <- joint_distrib / matrix(rep(rowSums(joint_distrib), 3), 3, 3)
rowSums(empirical_conditional_one_neighbor)  # Sanity check, all 1

## Compare (empirical) empirical_conditional_one_neighbor to conditional_one_neighbor: they're clearly different,
## one is an initial distribution and the other is something like a stationary distribution
round(empirical_conditional_one_neighbor, 3)
round(conditional_one_neighbor, 3) When I say that I see "streaks" in the figure above, what I mean is that there are many clumps of forest that are long in the [down left]-to-[up right] direction, and not many clumps of forest that are long in the [up left]-to-[down right] direction. My intuition is that this is a consequence of the way that I simulated the raster, and that it means that, if there is a stationary joint distribution of $\left(X_{i, j}, X_{i+1, j}, X_{i, j+1}, X_{i-1, j}, X_{i, j-1}\right)$, i.e. of a pixel and its four neighbors, that distribution is not invariant to rotation. Is that correct?

Edit: has anyone reading this question read Introduction to Markov Random Fields by Andrew Blake and Pushmeet Kohli (chapter 1 available at https://mitpress.mit.edu/sites/default/files/titles/content/9780262015776_sch_0001.pdf)? If so, do you know whether I can find an answer to my question(s) in that book?