4

I am using ArcGIS 9.3 and python 2.5. My code is shown at the bottom of this post.

I have a raster containing fire intensity information at 200m cell resolution. I run the script and the first count cell part returns the size of the fire. This works fine.

I then clip the raster using a polygon representing a forested catchment. The new raster represents the area of the catchment that is burnt. When I try to use the same code to count the number of cells of the new raster, it always gives me an answer of zero, even if there is definitely a raster created by the clip tool.

I have tried every combination of every tool I can think of.

my code...

    fire = str(ignitionID) + "_" + str(weatherID) + "_Intensity.asc"
    print "fire: " + fire

    # ---------------------------------------------------------------------------
    # this section calculates the size of the fire for model calibration purposes
    # ---------------------------------------------------------------------------

    # Process: Build Raster Attribute Table...
    the_fire = input_fires + fire
    shutil.copy(the_fire, "C:\\FIRESTORM\\Workspace") # copy the fire to workspace
    the_fire = "C:\\FIRESTORM\\Workspace\\" + fire

    try:
      gp.BuildRasterAttributeTable_management(the_fire, "NONE")
      resultf = gp.GetCount_management(the_fire)  # count num of rows in raster
      countf = int(resultf.GetOutput(0))
      fire_area = countf * 4  # each cell has an area of 4 ha
      print "the fire area = " + str(fire_area) + " ha"

    except:  # the attribute table cannot be created coz the fire is too large

      fire_area = 300000 # default value for large fires 300000 ha
      print "the fire area EXCEEDS 262,140 ha"

    # -----------------------------------------------------------------------------
    # this section calculates the area of the catchments that are burnt by the fire
    # -----------------------------------------------------------------------------

    # Local variables...
    catch_extent = input_catch + "WS_Catchment_region.shp"  # catchment boundary
    fire_raster = "C:\\FIRESTORM\\Workspace\\fire_raster"
    fire_ascii = the_fire
    Output_raster = "C:\\FIRESTORM\\Workspace\\Extract_ASCI1"
    Output_ASCII = "C:\\FIRESTORM\\Workspace\\rastert_extract1.asc"
    catch_burn_area = 0  # reset catchment burn area

    # Process: ASCII to Raster...
    gp.ASCIIToRaster_conversion(fire_ascii, fire_raster, "INTEGER")

    # Process: Define Projection... based on projection of catch_extent
    desc = gp.Describe(catch_extent)
    spatialRef = desc.SpatialReference
    gp.DefineProjection_management(fire_raster, spatialRef)

    # Process: Extract by Mask...
    gp.ExtractByMask_sa(fire_raster, catch_extent, Output_raster)

    # Convert raster to ascii
    gp.RasterToASCII_conversion(Output_raster, Output_ASCII)

    # Process: Build Raster Attribute Table...
    gp.BuildRasterAttributeTable_management(Output_ASCII, "OVERWRITE")

    # count number of rows in raster
    resulti = gp.GetCount_management(Output_ASCII)
    counti = int(resulti.GetOutput(0))
    catch_burn_area = counti * 4  # each cell has an area of 4 ha

    print "catchment area burnt = " + str(catch_burn_area) + " ha"
  • 1
    why not use zonal statistics ? – radouxju Jan 20 '14 at 10:32
  • 1
    At a first glance, I would remove the conversion to ascii. – radouxju Jan 20 '14 at 10:34
  • 1
    Is the raster that is extracted by the mask turning into a floating point raster? If so the Build Raster Attribute Table tool will not create an attribute table and presumably the subsequent GetCount call is returning a zero? – Hornbydd Jan 20 '14 at 11:48
  • 1
    A zonal sum (using the polygon as a zone) would likely be a better solution in all respects: speed, flexibility, accuracy, data storage. – whuber Jan 20 '14 at 15:54
  • 1
    @Jens: You are correct, however, because the value in each of my cells is a fire intensity from 1 to 2,000,000 there are very few rows with the same value. GetCount gives me 99.9% accuracy and, to be honest, i didt know how else to count the number of cells in a raster. I will attempt to use a zonal sum for both parts though, thank you for your interest. – McVegas Jan 21 '14 at 12:15
1

couple of options that occur to me without getting too deep in your code. Both rely on a zone mask, which you could make by pushing your raster, which is presumably a mix of values, through a Con (doc here) so output has only values of 1 where there is data. Can take that mask/zone grid and use it w/:

  1. zonalgeometry (doc here)
  2. @whuber's suggestion of zonalsum the 1s and multiply by the area.

Alternatively, you could use a cursor to sniff the one and only VAT row and extract the COUNT field (not sure of the syntax for this back in 9.3). Example:


---EDIT---

csr = arcpy.SearchCursor(gridVAT)
row = csr.next()
min = row.getValue('COUNT')

Thinking COLUMNCOUNT + ROWCOUNT won't work on irregularly shaped (i.e., not square and oriented w/coord sys).

Agree w/@radouxju, don't really see the purpose of the ASCII representation. Seems wildly expensive, but can't see your whole workflow.

  • Thank you, this is really helpful. I have converted the raster into 1s and nodata using the Con function as you describe. At the expense of appearing really stupid, I use the zonalgeometry to calculate the raster area but the output of the ZonalGeometry function is another raster. I just want the Area as a number! Can you please tell me how to get a number out of the ZonalGeometry function? Thanks again for you help. – McVegas Feb 2 '14 at 2:06
  • Zonal stats will output a table if you ask it. BTW, McVegas, you now have enough reputation to upvote helpful answers like this one: doing so helps others in the future when they're searching for good solutions. – whuber Feb 3 '14 at 20:10

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.