[GIS] Creating raster by randomly choosing cell value from multiple overlapping rasters

arcgis-10.0arcgis-desktopmap-algebrarasterspatial-analyst

I'm using ArcGIS Desktop 10 with its Spatial Analyst extension.

How do I combine multiple rasters into one, always choosing randomly from the values of overlapping cells?

I have an image that may explain this better:

example

Best Answer

Pick was created for problems like this. Think of it as the "switch" (or "case") version of "con", which is the map algebra implementation of "if...else."

If there are 3 overlapping rasters, for instance, the (Python) syntax would look like

inPositionRaster = 1 + int(3 * CreateRandomRaster())
Pick(inPositionRaster, [inRas01, inRas02, inRas03])

Note that pick starts indexing at 1, not 0.

Edit

(see the comment thread)

To cope with NoData values, first you need to turn off ArcGIS's NoData handling. Do this by creating grids that have a special (but valid) value in place of NoData, such as 99999 (or whatever: but make sure to choose a value that is larger than any valid number that can appear: this will be handy later). This requires use of the IsNull request, as in

p01 = Con(IsNull(inRas01), 99999, inRas01)
p02 = Con(IsNull(inRas02), 99999, inRas01)
p03 = Con(IsNull(inRas03), 99999, inRas01)

For example, consider the case of these one-row grids (NoData is shown as "*"):

inRas01:  1  2 19  4  *  *  *  *
inRas02:  9  2  *  * 13 14  *  *
inRas03: 17  *  3  * 21  * 23  *

The result is to put a 99999 in place of each "*".

Next, imagine all these rasters as flat arrays of wooden blocks with NoData corresponding to missing blocks (holes). When you vertically stack these rasters, blocks will fall into any holes beneath them. We need that behavior to avoid picking NoData values: we don't want any vertical gaps in the stacks of blocks. The order of the blocks in each tower doesn't really matter. To this end, we can obtain each tower by ranking the data:

q01 = Rank(1, [p01, p02, p03])
q02 = Rank(2, [p01, p02, p03])
q03 = Rank(3, [p01, p02, p03])

In the example, we obtain

q01:      1     2     3     4    13    14    23 99999
q02:      9     2    19 99999    21 99999 99999 99999
q03:     17 99999 99999 99999 99999 99999 99999 99999

Note that the ranks are from lowest to highest, so that q01 contains the lowest values at each location, q02 contains the second-lowest, etc. The NoData codes don't start showing up until all the valid numbers are collected, because those codes are larger than any valid numbers.

In order to avoid picking these NoData codes during the random selection, you need to know how many blocks are stacked at each location: this tells us how many valid values occur. One way to handle this is to count the number of NoData codes and subtract that from the total number of selection grids:

n0 = 3 - EqualToFrequency(99999, [q01, q02, q03])

This yields

n0:       3 2 2 1 2 1 1 0

To handle the cases where n=0 (so there's nothing available to select), set them to NoData:

n = SetNull(n0 == 0, n0)

Now

n:        3 2 2 1 2 1 1 *

This will also guarantee that your (temporary) NoData codes disappear in the final calculation. Generate random values between 1 and n:

inPositionRaster = 1 + int(n * CreateRandomRaster())

For instance, this raster might look like

inPositionRaster: 3 2 1 1 2 1 1 *

All its values lie between 1 and the corresponding value in [n].

Select values exactly as before:

selection = Pick(inPositionRaster, [q01, q02, q03])

This would result in

selection:       17  2  3  4 21 14 23  *

To check that all is ok, try selecting all output cells that have the NoData code (99999 in this example): there shouldn't be any.

Although this running example uses only three grids to select from, I have written it in a way that generalizes easily to any number of grids. With lots of grids, writing a script (to loop over the repeated operations) will be invaluable.

Related Question