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:


## 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))
                           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)[1], dim(conditional_two_neighbors)[2]))))

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,
                               prob=conditional_two_neighbors[neighbor_north_index, neighbor_west_index, ])

raster_df <- expand.grid(x=seq_len(nrow_raster),

## 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",
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)

simulated raster

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?


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.