3

In ArcMap 10, I have a few dozen tiffs with cells that represent habitat area that has been burned. Each burned cell has a Julian date of the fire associated with it.

I want to end up with polygons that represent a single fire event each, defined as (8-way) contiguous groups of cells where no bordering burned cells differ by >1 Julian day; i.e., where burned cells differ by 2 or more days, these would be classed as separate fires. If I convert the tiffs to shp files based on the date field, it groups all the contiguous cells for a single date.

How can I merge only the polygons that are both adjacent and differing by one day?

Improve this question
3
  • You could reclassify the original rasters to "event" rather than "date" -- is this something you want to automate or semi-automate rather than do manually?
    – Erica
    Aug 13, 2014 at 13:19
  • In general this cannot be done, at least not in a unique manner. Consider a fire that progressed for weeks, moving about as it did: there may be no pixels that satisfy these criteria. A correct solution to this question would thereby completely miss such a major event or else incorrectly classify it as a huge number of separate small events. What do you want to do in such cases?
    – whuber
    Aug 13, 2014 at 15:46
  • Whuber is correct that the approach I'm going for is imperfect when a fire takes a few days to jump cells. But, we know from looking at the data that this is a rare situation. Usually fires have at least a few contiguous cells. Erica, reclassifying to "event" is, I think, just what we're after. And, yes, ideally there's a tool that leads to this solution. We could do it manually, but it's a lot of work for 12 months X 14 years.
    – ecologist1234
    Aug 13, 2014 at 17:21

2 Answers 2

Reset to default
1

It seems that solution here is using raster itself. I've tried to reproduce raster with 'similar' cell values in different zones using set of polygons (142 total), distance to their boundaries and maximum distance inside individual polygon. Max distance minus actuall distance produced areas with breaks in values at the boundary. Raster cellsize is 1. enter image description here

It is just a matter of iteration through raster cells and comparing values with 8 neighbours, skipping already assigned. This picture compares original polygons with cells grouping derived: enter image description here

Original boundaries shown in white, derived cell groups shown by unique colour. Output raster contains 99 groups only. Obviously algorithm I've used in attempt to reproduce your data fails to produce true break lines at the boundary of 2 or more polygons with similar or close distance pattern. This is why 'fire' managed to 'sneak' into other catchment through cell with 'similar' value in it. Raster is 682*542 in size and it took 2 minutes to process on my out of date home PC. I've used

arcpy.RasterToNumPyArray(raster,"","","",-9999)

Not sure how it will perform on larger set due to memory limitations with this function. Let me know if you are interested in script itself and I'll try to post it. Frankly I found it very frustrating experience to post code sample on this site...

Improve this answer
2
  • Thanks. Promising to see it's doable, but alas I'm a python neophyte at best.
    – ecologist1234
    Aug 17, 2014 at 15:07
  • I'll post script as soon as possible. The only thing you have to learn is how to add new script to custom toolbox. You'll have to change 2 lines in my code because they are hard coded
    – FelixIP
    Aug 18, 2014 at 2:16
1 
# GROUPS adjacent slightly different cells
# extension - Spatial Analyst

import arcpy, traceback, sys
from arcpy import env
env.overwriteOutput = True
import numpy
try:
  def showPyMessage():
    arcpy.AddMessage(str(time.ctime()) + " - " + message)

#       load input raster to array and define it's size, origin point, cell size
#       replace line below. Consider converting tif to integer grid
    raster=arcpy.Raster('C:/FELIX_DATA/theDates')
    myArray = arcpy.RasterToNumPyArray(raster,"","","",-9999)
    nRows,nCols=myArray.shape
    cellsTotal=nCols*nRows
    d=arcpy.Describe(raster)
    origin=d.extent.lowerLeft
    cSize=raster.meanCellHeight
#       raster to store assigned group number. Set to -1
    blankR=raster-raster-1
    blank = arcpy.RasterToNumPyArray(blankR,"","","",-9999)
#       directions to find cell neighbour
    dRow=(0,-1,-1,-1,0,1,1,1)
    dCol=(1,1,0,-1,-1,-1,0,1)
#       group counter
    groupNo=0
    arcpy.SetProgressor("step", "", 0, cellsTotal)
    for nRow in range(nRows):
            for nCol in range (nCols):
#                       unassigned cell value = -1
                    inGroup=blank[nRow,nCol]
                    if inGroup==-1:
                            groupNo+=1
                            blank[nRow,nCol]=groupNo
                            lstInit=[(nRow,nCol)]
                            lst2grow=[]
                            while len(lstInit)>0:
                                    for cell in lstInit:
                                            nR,nC=cell
                                            v=myArray[nR,nC]
                                            for i in xrange(8):
                                                    curRow=nR+dRow[i]
                                                    curCol=nC+dCol[i]
                                                    if curRow in range(nRows) and curCol in range(nCols):
                                                            inGroup=blank[curRow,curCol]
                                                            if inGroup==-1:
                                                                    vCur=myArray[curRow,curCol]
                                                                    if vCur in xrange(v-1,v+2):
                                                                            blank[curRow,curCol]=groupNo
                                                                            lst2grow.append((curRow,curCol))
                                                                            arcpy.SetProgressorPosition()
                                    lstInit=lst2grow[:]
                                    lst2grow=[]
                            arcpy.SetProgressorLabel("Groups found %i" %groupNo)

    del myArray
    myRaster = arcpy.NumPyArrayToRaster(blank,origin,cSize,cSize)
    del blank
#       replace line below
    myRaster.save('C:/FELIX_DATA/groups')
except:
  message = "\n*** PYTHON ERRORS *** "; showPyMessage()
  message = "Python Traceback Info: " + traceback.format_tb(sys.exc_info()[2])[0]; showPyMessage()
  message = "Python Error Info: " +  str(sys.exc_type)+ ": " + str(sys.exc_value) + "\n"; showPyMessage()
Improve this answer
2

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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