Zonal statistics for raster using a matching mask file

5304
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
14 Replies
MorganHarris
New Contributor III

Sorry if I'm being a hassle! I converted some of my rasters to TIF just to mess with before I got your reply and have been just trying to get this to run, but I am getting the error shown below that says my input value raster isn't the right format or doesn't exist. The script is a little different since I have tweaked some things, but if you could take a look and see if you can tell what the problem is that would be great!

>>> import arcpy, os 

... from arcpy import env 

... from arcpy.sa import * 

...  

... rasterWorkspace = r"C:\Users\AmyShipley\Desktop\Offline_2014_GIS\Exp_Height_test_TIFF" 

... polygonWorkspace = r"C:\Users\AmyShipley\Desktop\Offline_2014_GIS\shp_files\Distance_Buffers\Span_IDx10m_test_group" 

...  

... 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, "FID", rasterName, "MINIMUM", "DATA") 

...             outZonalStats.save("zs_" + shapefileName)

...            

Runtime error  Traceback (most recent call last):   File "<string>", line 19, in <module>   File "c:\program files (x86)\arcgis\desktop10.2\arcpy\arcpy\sa\Functions.py", line 6145, in ZonalStatistics     ignore_nodata)   File "c:\program files (x86)\arcgis\desktop10.2\arcpy\arcpy\sa\Utils.py", line 47, in swapper     result = wrapper(*args, **kwargs)   File "c:\program files (x86)\arcgis\desktop10.2\arcpy\arcpy\sa\Functions.py", line 6138, in Wrapper     ignore_nodata)   File "c:\program files (x86)\arcgis\desktop10.2\arcpy\arcpy\geoprocessing\_base.py", line 498, in <lambda>     return lambda *args: val(*gp_fixargs(args, True)) ExecuteError: Failed to execute. Parameters are not valid. ERROR 000860: Input value raster: is not the type of Composite Geodataset, or does not exist. Failed to execute (ZonalStatistics). 

0 Kudos
JakeSkinner
Esri Esteemed Contributor

Add some print statements as shown previously.  It looks like it is failing when trying to run the ZonalStatistics for the intended raster.  We will need to see what raster is being passed to this function.  Adding "print rasterName" under "rasterName = raster.split(".")[0]" should help with this.

MorganHarris
New Contributor III

I added the print statements and it printed 10_100l for both shapefileName and rasterName, but then I still get the same error

>>> import arcpy, os 

... from arcpy import env 

... from arcpy.sa import * 

...  

... rasterWorkspace = r"C:\Users\AmyShipley\Desktop\Offline_2014_GIS\Exp_Height_test_TIFF" 

... polygonWorkspace = r"C:\Users\AmyShipley\Desktop\Offline_2014_GIS\shp_files\Distance_Buffers\Span_IDx10m_test_group" 

...  

... env.workspace = polygonWorkspace 

... shapefiles = arcpy.ListFeatureClasses("*", "Polygon") 

...  

... env.workspace = rasterWorkspace 

... rasters = arcpy.ListRasters("*", "TIF") 

...  

... 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, "FID", rasterName, "MINIMUM", "DATA") 

...             print "Saving raster"

...             outZonalStats.save("zs_" + shapefileName)

...            

10_100l

10_100l

Runtime error  Traceback (most recent call last):   File "<string>", line 21, in <module>   File "c:\program files (x86)\arcgis\desktop10.2\arcpy\arcpy\sa\Functions.py", line 6145, in ZonalStatistics     ignore_nodata)   File "c:\program files (x86)\arcgis\desktop10.2\arcpy\arcpy\sa\Utils.py", line 47, in swapper     result = wrapper(*args, **kwargs)   File "c:\program files (x86)\arcgis\desktop10.2\arcpy\arcpy\sa\Functions.py", line 6138, in Wrapper     ignore_nodata)   File "c:\program files (x86)\arcgis\desktop10.2\arcpy\arcpy\geoprocessing\_base.py", line 498, in <lambda>     return lambda *args: val(*gp_fixargs(args, True)) ExecuteError: Failed to execute. Parameters are not valid. ERROR 000860: Input value raster: is not the type of Composite Geodataset, or does not exist. Failed to execute (ZonalStatistics). 

0 Kudos
JakeSkinner
Esri Esteemed Contributor

That was my mistake.  We will need to add the full path to the shapefile.  Replace the outZonalStats with the following:

outZonalStats = ZonalStatistics(polygonWorkspace + os.sep + shapefile, "FID", rasterName, "MINIMUM", "DATA")

The os.sep is simply python for a backslash.

MorganHarris
New Contributor III

Alright! Got it to work! I had to change the outZonalStats line to have the raster represented by the full path name too. I also ended up adding a line to get zonal stats as table and it worked out fine.

Thanks for your help!

0 Kudos