Hi again GeoNet community!
I have a question and the answer may involve script. Just a warning, I have very rarely dealt with scripting and so if the answer is a script code then I may need some extra guidance. Anyway, my problem is this:
I have about 1300 raster layers that I need to iterate through. For each raster layer, I need to run zonal statistics using a polygon layer. Each polygon layer corresponds to each raster in space and they also have the same file name (i.e. the raster name is 110_10R and the polygon file name is 110_10R as well, they are just saved in different folders. How can I make it so that 110_10R (raster) uses 110_10R (polygon) for its zonal statistics? Remember this needs to iterate 1300 times (actually 1300 x 12) so it needs to be 110_10R (raster) with 110_10R (polygon) then it moves on a runs 111_10R (raster) with 111_10R (polygon)...etc...
Remember, scripting novice, but I realize that it may be the only way to solve this issue
Thanks!
Solved! Go to Solution.
Something like below should get you started:
import arcpy, os from arcpy import env from arcpy.sa import * rasterWorkspace = r"C:\Users\AmyShipley\Desktop\Offline_2014_GIS\Exposure_Height" polygonWorkspace = r"C:\Users\AmyShipley\Desktop\Offline_2014_GIS\shp_files\Distance_Buffers" env.workspace = polygonWorkspace shapefiles = arcpy.ListFeatureClasses("*", "Polygon") env.workspace = rasterWorkspace rasters = arcpy.ListRasters("*", "TIF") for shapefile in shapefiles: shapefileName = shapefile.split(".")[0] for raster in rasters: rasterName = raster.split(".")[0] if rasterName == shapefileName: outZonalStats = ZonalStatistics(shapefileName, "value", rasterName, "SUM", "DATA") outZonalStats.save("zoneStats_" + shapefileName)
The above code will iterate through the raster workspace and find all TIFF files. If you are working with a different raster type, take a look at the ListRasters to change to the correct type, or simply remove "TIF" to use all rasters.
Hi Morgan,
Where are the polygon feature classes/shapefiles located? Are they in the same directory as the rasters with the same name?
They are not in geodatabases or anything. They are just in different folders the pathnames may help you visualize how it is arranged in the directory?
Rasters:
C:\Users\AmyShipley\Desktop\Offline_2014_GIS\Exposure_Height
Polygons:
C:\Users\AmyShipley\Desktop\Offline_2014_GIS\shp_files\Distance_Buffers
Something like below should get you started:
import arcpy, os from arcpy import env from arcpy.sa import * rasterWorkspace = r"C:\Users\AmyShipley\Desktop\Offline_2014_GIS\Exposure_Height" polygonWorkspace = r"C:\Users\AmyShipley\Desktop\Offline_2014_GIS\shp_files\Distance_Buffers" env.workspace = polygonWorkspace shapefiles = arcpy.ListFeatureClasses("*", "Polygon") env.workspace = rasterWorkspace rasters = arcpy.ListRasters("*", "TIF") for shapefile in shapefiles: shapefileName = shapefile.split(".")[0] for raster in rasters: rasterName = raster.split(".")[0] if rasterName == shapefileName: outZonalStats = ZonalStatistics(shapefileName, "value", rasterName, "SUM", "DATA") outZonalStats.save("zoneStats_" + shapefileName)
The above code will iterate through the raster workspace and find all TIFF files. If you are working with a different raster type, take a look at the ListRasters to change to the correct type, or simply remove "TIF" to use all rasters.
Great! Thanks, I will work with this and see what I can do with it!
Actually, one question before I get started working with this. After shapefile.split, what does (".") mean? Also, what should fill in the (" * ") space when defining the workspace? Does that stay as an asterisk so that it iterates through?
Since you are working with Shapefiles, the name will be <name>.shp. The split method splits the shapefile name at the period and returns two values "<name>" and "shp" in a list. The [0] indicates to return the first value in the list. The same is true when dealing with a raster.
The "*" is a wildcard and will return all results when searching for shapefiles or rasters. If you just wanted to find a shapefile that began with "Pa", you could use:
shapefiles = arcpy.ListFeatureClasses("Pa*", "Polygon")
Ok, great, thanks for describing those for me
Ok, I've been trying to work with this code, but I must be doing something wrong. I have filled in the spaces with my own information and everything seems good to go, but when I hit enter nothing really happens. It doesn't give me errors or anything, but it also seems like it isn't doing anything. I have some specific questions below:
1) When I look at my raster files outside of arcmap they are each their own folder and have many files inside each folder that are .adf. Could this be a problem?
2) Is there a way to define where the data saves, and where is it saving right now so I can see if it did anything?
Thanks!
1) It appears you are working with ESRI GRIDs. You will want to change:
rasters = arcpy.ListRasters("*", "TIF")
to:
rasters = arcpy.ListRasters("*", "GRID")
2) The rasters are saving to your env.workspace. You set this to your rasterWorkspace variable. So, the rasters should save here once once you make the change mentioned above.
Also, it's helpful to add print statements throughout your code. For example:
for shapefile in shapefiles: shapefileName = shapefile.split(".")[0] print shapefileName for raster in rasters: rasterName = raster.split(".")[0] print rasterName if rasterName == shapefileName: outZonalStats = ZonalStatistics(shapefileName, "value", rasterName, "SUM", "DATA") print "Saving raster" outZonalStats.save("zoneStats_" + shapefileName)