I have a land cover classification that consists of three files. One file is the LCC for the western half of my study area, another covers the eastern half (these two overlap in the middle), and the third represents a singe vegetation type which is present in both of the other two images (overlaps both). I would like to stack them western, eastern, single veg type, from bottom to top respectively, and end up with a raster that keeps only the top most cell values. Basically, I want a single image with the same appearance as the three when stacked. I'm working with ArcMap 10.2 and Erdas Imagine 2013 – so solutions in either software would be helpful.
[GIS] Stacking overlapping rasters, keeping cell values on top
arcgis-10.2arcgis-desktoperdas-imaginerasterstack
Related Solutions
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.
I had the same issue some months ago.
Use gdal_merge to generate a new file from the 3 independent rasters.
In OSGeo4W command line you can do this:
gdal_merge.bat -separate -of GTiff -o output.tif input1.tif input2.tif input3.tif
In QGIS you can do the same with a GUI in the raster plugin "merge" tool.
Best Answer
There are several methods to combine rasters in ArcGis; generally ERDAS is faster at raster functions than ArcGis but I'm so out-of-touch with ERDAS that I can't even remember how to mosaic.
The first thing to consider is order of precedence. Lets assume that west is more important than east. Using the Mosaic to New Raster tool you can combine the three rasters into one. The only concerns are to get the rasters in the right order and use the appropriate overlap method, assuming you are using this as a tool and not in a python script either:
Use overlap first and select the vegetation type raster, west then east.
Use overlap last and select east, west then vegetation type raster.
NoData never overwrites data... so to ensure that top rasters don't mask lower rasters you must set the NoData value using the SetNull tool; for example if your raster contains classified data but 0 for nothing use SetNull to make the 0's into NoData.
Another option is to create a mosaic dataset, which doesn't combine the rasters but rather draws them in turn.. you can adjust the drawing order after creation until you're happy with the result then use Copy Raster to make it into a single raster dataset, this tool can also be launched by right clicking on the mosaic dataset and selecting raster to different format (or something like that).