Iterate many raster and many shapefiles in python

4912
9
Jump to solution
10-18-2014 09:36 AM
NathanMietkiewicz
New Contributor

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

0 Kudos
1 Solution

Accepted Solutions
JoshuaBixby
MVP Esteemed Contributor

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)

View solution in original post

9 Replies
DanPatterson_Retired
MVP Emeritus

Are there errors?  Are you looking for style comments? Anything specific?

0 Kudos
NathanMietkiewicz
New Contributor

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

0 Kudos
curtvprice
MVP Esteemed Contributor

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")

0 Kudos
NathanMietkiewicz
New Contributor

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

0 Kudos
curtvprice
MVP Esteemed Contributor

What does "not print i" mean?

0 Kudos
NathanMietkiewicz
New Contributor

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

0 Kudos
curtvprice
MVP Esteemed Contributor

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

0 Kudos
JoshuaBixby
MVP Esteemed Contributor

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)

NathanMietkiewicz
New Contributor

Joshua,

That did the trick!  Thank you both for your help with this, it is greatly appreciated!

Best,

Nate

0 Kudos