Using ArcGIS 10.1 with Spatial and 3D Analysts, and Python 2.7
I have a series of reclassified rasters in which every cell is either one or zero. I need to programmatically create a single raster which shows when each cell first had a value of one. The rasters contain the date in the filename, as well as in a date field in the raster. (More detailed description of the data below)
I am completely stymied. I've tried using the Plus function in Spatial Analyst, but I lose the date. I have very little experience working with rasters, so I'm hoping that the problem is my unfamiliarity with the terminology, and that someone can point me in the right direction.
If I can figure out how to do it 'by hand' in ArcMap, I'm confident I can create the scripts to automate the process. Of course, a ready to go script would be very nice :^)
The data is 35 years of sea-ice concentration shapefiles from the Canadian Ice Service. Clipped tifs were generated and reclassified so that each cell in the raster has a value of either one, if the ice concentration is below 50% (break-up), or zero if it is not (i.e. concentration greater than 50% or land). I then added the date, parsed from the filename, using Create and CalculateField_management. The objective is to create annual maps showing the progression of the ice break-up to use in biological analyses.
Because this analysis is local--the result in any cell depends only on the stack of values in that cell--we might look first to the Local Toolset for inspiration. It is best when such tools can be applied without any preliminary transformation of the data, because that avoids a potentially time-consuming loop over all the rasters.
Perhaps the most general solution available in this form uses Highest Position. This calculation applies to a sequence of rasters, which ought to be perfectly aligned and registered (so there is no complication entailing from resampling the values). It returns the index 1, 2, 3,… , etc., of the first raster at which the maximum value in the sequence is encountered. Thus, when the sequence is purely binary and there is a one, the position of its first occurrence is returned--and that's exactly what is needed.
When there is one raster per time unit (such as year), with no gaps, adding the year preceding the start year will convert this result into the raster date. Otherwise, reclassify the result to obtain the corresponding years.
Two cautions are in order.
First, if any of the input rasters has NoData in a cell, the output for that cell will be NoData. If that's a problem, the rasters will have to be pre-processed to replace the NoData values by suitable codes. Some constant negative value might work well. Use a conditional operator to make the replacements.
Second, if all of a cell's values are zeros, then
HighestPositionwill return the index of the first of those zeros--namely, 1. That is indistinguishable from the result of a cell where all values are ones. Although some conditional post-processing can be used to discriminate these two situations, a simple, elegant way to avoid this problem is to create a raster of zeros at the outset and list it first in the calculation. Now a 1 will be returned for any all-zero cell and a 2 will be returned for any all-one cell.