Zonal statistics for raster using a matching mask file

4184
14
Jump to solution
01-27-2015 05:19 AM
MorganHarris
New Contributor III

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!

0 Kudos
1 Solution

Accepted Solutions
JakeSkinner
Esri Esteemed Contributor

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.

View solution in original post

0 Kudos
14 Replies
JakeSkinner
Esri Esteemed Contributor

Hi Morgan,

Where are the polygon feature classes/shapefiles located?  Are they in the same directory as the rasters with the same name?

0 Kudos
MorganHarris
New Contributor III

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

0 Kudos
JakeSkinner
Esri Esteemed Contributor

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.

View solution in original post

0 Kudos
MorganHarris
New Contributor III

Great! Thanks, I will work with this and see what I can do with it!

0 Kudos
MorganHarris
New Contributor III

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?

0 Kudos
JakeSkinner
Esri Esteemed Contributor

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

MorganHarris
New Contributor III

Ok, great, thanks for describing those for me

0 Kudos
MorganHarris
New Contributor III

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!

0 Kudos
JakeSkinner
Esri Esteemed Contributor

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)