4

I'm using ArcMap 10.2.2 and I have a raster that has Nodata around the edges. How can I extract the tightest rectangle from this raster that does not include Nodata? The screenshot may help illustrate my question.

This is my raster and Nodata is indicated as yellow. What tool/script can I run/write that will extract the raster based on the red bounding rectangle? The purpose of this is to NOT use manual/visual inspection to get the red bounding rectangle extent values. I intend to use whatever method that works to include in a larger Python script.

clip

3
  • Perhaps you could cut into very small chunks, test each chunk by horizontal/vertical for existence of NoData and then merge. – Michael Stimson Mar 25 '15 at 3:09
  • 1
    Sorry about the screenshot... moved it to another place to host and lost it for a while. Should be working – nokalake Apr 2 '15 at 20:57
  • Screenshot is broken again. Perhaps re-edit and upload to the site rather than hosting it somewhere else. – Fezter Jan 17 '16 at 23:06
5

This is potentially a very difficult problem when the borders are ragged. A brute-force search of the optimum could require computational time that is proportional to the square of the number of cells in the image (a value that often will be in the billions, trillions, or greater).

One promising approach is to relax the conditions a little bit and actually allow a few NoData values to sneak in: simply penalize their inclusion. This affords some welcome flexibility, too: a small string of NoData values penetrating from the outside (such as a skinny gap) could be included rather than cutting off most of the image.

This suggests maximizing the area of the remaining image, subject to a penalty that depends on how many NoData cells are left in the remaining image. To solve this problem it helps to interpolate positions within the image, so that it changes from a large discrete problem (where only integral row and column indexes are allowed) to a continuous problem. Bilinear interpolation should be fine.

Only four variables are needed to represent the coordinates of the lower left and upper right of the final included image. This makes the problem quite tractable.

To expedite the calculations of the numbers of included NoData cells, first compute a cumulative sum array Y by summing along columns and then accumulating those sums along rows. To obtain the number of NoData cells included in the rectangle spanning the coordinates (i,k) through (j,l), just calculate

Y[k,l] - Y[k,j] - Y[i,l] + Y[i,j]

That's much faster than re-summing them each time.

A canned optimizer (in R) had some difficulties finding solutions. I was able to tweak it by starting with a small penalty factor and gradually increasing it, restarting each new search at the value of the preceding one. Occasionally this fails and needs to be restarted with a different value, but most of the time it works beautifully.

Because of the interpolation, an extra row or column of pixels might be stripped away unnecessarily. A little post-processing could detect this circumstance and fix it up. (I haven't implemented that.)

Here is the progress of one calculation involving an 85 by 119 array with randomly generated NoData cells within 15 pixels of its border all around. The upper row plots the masked image (starting with the original at the upper left). The bottom row, from left to right, depicts the cumulative sum array Y and the successive masks that were found. The final one includes 4717 pixels, which is slightly less than the correct value of 4895 pixels (because of the aforementioned interpolation issue). This calculation took 0.3 seconds.

It scales beautifully: on a similar configuration with 100 times as many pixels the optimization step took under one second, only three times longer. Another factor of 100 (an image of over 100,000,000 pixels) brought the time up to 80 seconds. (You have to turn off the display of the intermediate images--R is particularly slow at that.)

Figure


The R code that produced this example is a bit of a hack--it's not bulletproof nor is it optimized for speed. It is offered to illustrate the details and perhaps to serve as a point of departure for anyone who would like to improve or port it.

#
# Create data.
#
border <- 15
m <- 55+2*border
n <- 89+2*border
set.seed(17)
x <- rep(0, m*n)
e <- runif(m*n) < 1/3
i <- as.vector(outer(1:m, 1:n, 
                     function(a,b) (a < 1+border | a > m-border) | 
                       (b < 1+border | b > n-border)))
x[i] <- e[i]
x <- matrix(x, m, n)

par(mfcol=c(2,4))
#
# Display the original.
#
n.bad <- sum(x)
image(x, col=c("Gray", "Red"), main=paste("Bad pixels =", n.bad))
#
# Compute the cumulative sums of NA data.
#
y <- apply(apply(x, 1, cumsum), 1, cumsum)
image(y, col=rainbow(max(m,n)), main="Cumulative bad pixels")
#
# Perform bilinear interpolation into an array.
#
access <- function(y, ij) {
  bounds <- function(i, n) {
    lower <- pmin(n-1, pmax(1, floor(i)))
    upper <- lower+1
    delta <- pmax(0, pmin(1, c(upper-i, i-lower)))
    return (list(lower=lower, upper=upper, delta=delta))
  }
  b.1 <- bounds(ij[1], dim(y)[1])
  b.2 <- bounds(ij[2], dim(y)[2])
  values <- c(y[b.1$lower,b.2$lower], y[b.1$upper,b.2$lower], 
              y[b.1$lower,b.2$upper], y[b.1$upper,b.2$upper])
  weights <- as.vector(outer(b.1$delta, b.2$delta, '*'))
  return (sum(values*weights))
}

# z <- matrix(NA, floor(pi*m), floor(pi*n))
# for (i in 1:dim(z)[1]) for (j in 1:dim(z)[2]) z[i,j] <- access(y, c(i,j)/pi)
# image(z)

#
# Define the objective function (to be minimized).
#
f <- function(indexes, y, rho=1, rho.boundary=1) {
  clamp <- function(x, d) pmin(d, pmax(1, x))
  m <- dim(y)[1]; n <- dim(y)[2]
  i0 <- min(indexes[1:2])
  k0 <- max(indexes[1:2])
  j0 <- min(indexes[3:4])
  l0 <- max(indexes[3:4])
  i <- clamp(i0, m)
  k <- clamp(k0, m)
  j <- clamp(j0, n)
  l <- clamp(l0, n)
  area <- abs((k-i+1) * (l-j+1))
  penalty <- access(y, c(k,l)) - access(y, c(k,j)) - access(y, c(i,l)) + access(y, c(i,j))
  penalty.b <- sum((c(i,j,k,l) - c(i0,j0,k0,l0))^2)
  return (-area + rho*penalty + rho.boundary*penalty.b)
}

# f.r <- matrix(NA, floor(pi*m), floor(pi*n))
# for (i in 1:dim(z)[1]) for (j in 1:dim(z)[2]) f.r[i,j] <- f(c(15,i,7,j), z)
# image(f.r)

#
# Sneak up on a solution by starting with small penalties and 
# increasing them gradually.
#
theta.0 <- c(1,m, 1,n)
rho <- 1/2
n.max <- 16
n.bad.old <- n.bad
repeat {
  if (n.max <= 0) break
  sol <- nlm(function(theta) f(theta, y, rho, 1), theta.0, steptol=0.1)
  mask <- matrix(FALSE, m, n)
  ij <- sol$estimate
  i0 <- ceiling(min(ij[1:2]))
  k0 <- floor(max(ij[1:2]))
  j0 <- ceiling(min(ij[3:4]))
  l0 <- floor(max(ij[3:4]))
  mask[i0:k0, j0:l0] <- TRUE
  x.mask <- x
  x.mask[!mask] <- NA
  #image(x)

  n.bad <- sum(x.mask, na.rm=TRUE)
  if (n.bad < n.bad.old) {
    image(x.mask, col=c("Gray", "Red"), 
          main=paste("Bad pixels =", n.bad), sub=paste("Penalty =", round(rho, 1)))
    image(mask, col=c("Orange", "White"), main="Mask")
    n.bad.old <- n.bad
  }
  if (n.bad == 0) break
  n.max <- n.max - 1
  rho <- rho*1.5
  theta.0 <- sol$estimate
}
#
# Compare the number of displayed pixels to the best value.
#
sum(mask)
(m-2*border)*(n-2*border)
1

This is untested, but my thinking would be:

  1. Convert No Data to Polygons.
  2. Get the bounds of each polygon and use them to create a rectangle. That is, for the left polygon, the rightmost X value will be the Left Bounding Coordinate of the red rectangle. For the right polygon, the leftmost X value will be the Right Bounding Coordinate of the red rectangle.
  3. Repeat for top and bottom No-Data polygons.
  4. Create rectangle using the collected coordinates.

Note, this might not work if you get L-shaped polygons if the top and sides touch. Or if you get O-shaped polygons if all sides touch. But you might be able to figure out how to handle those cases.

Your Answer

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

Not the answer you're looking for? Browse other questions tagged or ask your own question.