Monday 23 May 2016

python - How to calculate total number of data cells in raster?


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"


No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...