Hi all,
I am trying to iterate zonal statistics on many rasters (+200) over many shapefiles (+30), and step through the shapefiles such that each output table will be calculated across all raster for one shapefile. For instance, shapefile #1 processes all zonal statistics for all rasters and merges it to one table, then shapefile #2 does the same, and so on until all shapefiles have been processed. Any help/advice would be greatly appreciated! Below is the code I am currently working with.
# Import system modules
import arcpy, sys, os, string, glob
# Check out any necessary licenses
arcpy.CheckOutExtension("spatial")
# Overwrite pre-existing files
arcpy.env.overwriteOutput = True
#Get the raster datasets in the input workspace and loop through them from the start
def zonal_climate_fire():
for ws in listWS:
arcpy.env.workspace = ws
rasterlist = arcpy.ListRasters()
shplist = arcpy.ListFeatureClasses()
for k in rasterlist:
for i in shplist:
out_tbl = tbl + "\\" + k + "_ZStats"
print '==============================================='
print 'Zonal Statistics based on fire occurence : ' + k[6:-3]
print '......'
# Preform zonal statistics on the polygon file within the variable 'znlyr'
z = arcpy.gp.ZonalStatisticsAsTable_sa(i, InZnValFld, k, out_tbl, "DATA", "MIN_MAX_MEAN")
# Create a new field to add the filename information
arcpy.AddField_management(z, fieldname1, "DOUBLE")
# Create a new field to add the filename information
arcpy.AddField_management(z, fieldname2, "DOUBLE")
# Add the Year value from the filename to the new field
arcpy.CalculateField_management(z, fieldname1,k[6:-3])
# Add the Year value from the filename to the new field
arcpy.CalculateField_management(z, fieldname2,k[11:])
arcpy.Append_management(z, Template, "","","")
#--------------------Main to set variables--------------------------------------------------
znshp = "C:\\tmp\\znshp\\"
grd = "C:\\tmp\\grd\\"
tbl = "C:\\tmp\\tbl\\"
InZnValFld = "ID"
fieldname1= "Year"
fieldname2= "Month"
Template = "C:\\tmp\\tbl\\Template.dbf"
#Creates a list of multiple workspaces
listWS = [grd, znshp]
zonal_climate_fire ()
Thank you all in advance,
Nate
Solved! Go to Solution.
Does your shapefile workspace only have shapefiles and your raster workspace only have rasters? That is how it looks with znshp and grd. If that is the case, your listing of shapefiles will return nothing in the first pass of the ws loop, which means the second loop won't execute. What about changing zonal_climate_fire to accept the shapefile and raster workspaces as arguments:
def zonal_climate_fire(rasWS, shpWS):
arcpy.env.workspace = rasWS
rasterlist = arcpy.ListRasters()
arcpy.env.workspace = shpWS
shplist = arcpy.ListFeatureClasses()
for k in rasterlist:
for i in shplist:
out_tbl = tbl + "\\" + k + "_ZStats"
print '==============================================='
print 'Zonal Statistics based on fire occurence : ' + k[6:-3]
print '......'
# Preform zonal statistics on the polygon file within the variable 'znlyr'
z = arcpy.gp.ZonalStatisticsAsTable_sa(i, InZnValFld, k, out_tbl, "DATA", "MIN_MAX_MEAN")
# Create a new field to add the filename information
arcpy.AddField_management(z, fieldname1, "DOUBLE")
# Create a new field to add the filename information
arcpy.AddField_management(z, fieldname2, "DOUBLE")
# Add the Year value from the filename to the new field
arcpy.CalculateField_management(z, fieldname1,k[6:-3])
# Add the Year value from the filename to the new field
arcpy.CalculateField_management(z, fieldname2,k[11:])
arcpy.Append_management(z, Template, "","","")
zonal_climate_fire(grd, znshp)
Are there errors? Are you looking for style comments? Anything specific?
Hi Dan,
Im looking for some advice on stepping through the shapefiles. I am able to list the rasters and the shapefiles in their respective workspaces, but I cannot seem to call each shapefile, process all the rasters, then step onto the next shapefile and repeat raster processing. I hope this makes sense.
Thanks in advance.
Nate
Looks to me like you are trying to do too many things at once. If you don't have the main zonal stats code working, make a separate script to solve that. Then, try to get a loop printing your component names and variables to work right, like this:
def zonal_climate_fire():
for ws in listWS:
arcpy.env.workspace = ws
rasterlist = arcpy.ListRasters()
shplist = arcpy.ListFeatureClasses()
for k in rasterlist:
for i in shplist:
out_table = os.path.join(tbl, "{}_ZStats".format(k))
code1 = ras[6:-3]
code2 = ras[11:]
print("Raster: {} Shape: {} Code1: {} Code2: {}".format(
k,i, code1, code2))
print(out_table)
Then you can put the two pieces together.
BTW this is the usual format for zonal statistics As Table, you don't need to capture the result object as already know the results are written to out_table.
arcpy.ZonalStatisticsAsTable_sa(zone_dataset, zone_field, raster,
out_table, "DATA", "MIN_MAX_MEAN")
Hi Curtis,
Thank you for the quick response. I was able to get the zonal statistics loop work fine if I did not call the "for i in shplist" but instead fed one shapefile at a time manually - which is what I am trying to overcome.
Unfortunately I am not able to get your naming script to work. It stops during the "for i in shplist" command, which was happening to me originally as well. I am able to print the rasterlist, shplist, and k, but not i.
Any idea what may be happening?
Thanks,
Nate
What does "not print i" mean?
Hi Curtis,
Sorry for the lack of description there. What I was referring to was the script crashes when it encounters the "for i in shplist:" loop. A small test of trying to print the variable 'i' confirmed that. I have confirmed that "shplist = arcpy.ListFeatureClasses() is properly working and storing the appropriate shapefiles. Additionally, the first for loop (for k in rasterlist:) is working fine in that I am able to call the individual raster files (k) without a problem. I have a remedial knowledge of python and do not understand why this nested loop:
for k in rasterlist:
for i in shplist:
would not be working - logically it seems as if it should be an appropriate command.
I hope that clears things up. Thank you again for your help with this.
Best,
Nate
Beats me too. This works fine for me (ArcGIS 10.1 SP 1, ArcMap python prompt):
>>> fcs = arcpy.ListFeatureClasses()
>>> grids = arcpy.ListRasters()
>>> for k in grids:
... for i in fcs:
... print (k)
... print (i)
...
Proposal.jpg
centroid_points.shp
Proposal.jpg
County_Data_polygons1.shp
Does your shapefile workspace only have shapefiles and your raster workspace only have rasters? That is how it looks with znshp and grd. If that is the case, your listing of shapefiles will return nothing in the first pass of the ws loop, which means the second loop won't execute. What about changing zonal_climate_fire to accept the shapefile and raster workspaces as arguments:
def zonal_climate_fire(rasWS, shpWS):
arcpy.env.workspace = rasWS
rasterlist = arcpy.ListRasters()
arcpy.env.workspace = shpWS
shplist = arcpy.ListFeatureClasses()
for k in rasterlist:
for i in shplist:
out_tbl = tbl + "\\" + k + "_ZStats"
print '==============================================='
print 'Zonal Statistics based on fire occurence : ' + k[6:-3]
print '......'
# Preform zonal statistics on the polygon file within the variable 'znlyr'
z = arcpy.gp.ZonalStatisticsAsTable_sa(i, InZnValFld, k, out_tbl, "DATA", "MIN_MAX_MEAN")
# Create a new field to add the filename information
arcpy.AddField_management(z, fieldname1, "DOUBLE")
# Create a new field to add the filename information
arcpy.AddField_management(z, fieldname2, "DOUBLE")
# Add the Year value from the filename to the new field
arcpy.CalculateField_management(z, fieldname1,k[6:-3])
# Add the Year value from the filename to the new field
arcpy.CalculateField_management(z, fieldname2,k[11:])
arcpy.Append_management(z, Template, "","","")
zonal_climate_fire(grd, znshp)
Joshua,
That did the trick! Thank you both for your help with this, it is greatly appreciated!
Best,
Nate