topic ZonalStatisticsAsTable in Geoprocessing Questions

ZonalStatisticsAsTable

IrèneDaubard — Tue, 28 Jul 2020 19:19:13 GMT

Hello everyone! I'm starting with Python and I hope someone will help me see through the mistakes I might have made.

I'm making a zona lstatistics as table loop with python and i would like to make ane average of all TIF files pixels according to a grid (.shp). When I run my script I've the error message bellow: ERROR 000876: Output raster: E:\climat_data\tb_51N_Precip_2014-01-15.dbf's extension is invalid for the output raster format.
Failed to execute (ZonalStatistics).

I'm leaving my script here... maybe someone can help me figure it out:

arcpy.CheckOutExtension("Spatial")
#Set local variables
in_zone_data = 'grid51N.shp'
zone_field = 'id'

raster = arcpy.ListRasters("51N*","ALL")
for stat in raster:
print (stat)

outDir = r"E:\climat_data"

for stat in raster:
        raster_name = os.path.basename(stat).rstrip(os.path.splitext(stat)[1])
        outTable = os.path.join(outDir,"tb_"+raster_name+".dbf")
        arcpy.gp.ZonalStatistics_sa(in_zone_data, zone_field, stat, outTable, "MEAN", "NODATA")
        print 'Done'+raster_name

Thank you for any help.

Re: ZonalStatisticsAsTable

DanPatterson — Tue, 28 Jul 2020 19:39:42 GMT

could you throw code formatting in

/blogs/dan_patterson/2016/08/14/script-formatting

And I think you will need to put some print statements in since *.dbf files can only go into a folder and it isn't clear where you are trying to send it given your code above

Re: ZonalStatisticsAsTable

IrèneDaubard — Sat, 11 Dec 2021 15:36:58 GMT

#The main data folder contains climat data (precipitation, temperature and vegetation index) from 2014.
#The script has the objective to reproject all the climat data to the right projection for Philippines witch is "WGS 84 UTM zone 51N" and at the same time, resampling the cell size to 4km*4km.
#Then we create another list with other attributes witch are "51N_..." to make zonal statistics as table process with a grid having IDs for each cell. We want to have the mean of all the climatic data for each celle of the grid (GRID51N.shp).
#Finally, we have to join by ID, all the climat variable in order to get one layer with all the climat data.
#Different steps are commented bellow:
#________________________________________________________________________________________


#IMPORT MODULES
import arcpy, os, sys
from arcpy import env
from arcpy.sa import *
arcpy.env.workspace=r"E:\climat_data"

arcpy.env.overwriteOutput = True 
##arcpy.AlterField_management(grid51N_shp, "id", "loc", "TEXT")

#Adding the right spatial reference for the reproject loop
sr = arcpy.SpatialReference(32651)#"WGS 84 UTM zone 51N" spatial reference

#Getting only the TIFF files of the workspace
##lstFiles = []
##for root, dirs, files in os.walk(path):
##    files = arcpy.ListRasters("*", "TIF")  
##    for fichero in files:  
##        (nombreFichero, extension) = os.path.splitext(fichero)  
##        if(extension == ".tif"):  
##            lstFiles.append(os.path.join(path, fichero))
##print (lstreproj)

#####################Reprojection and resampling

reprojList = arcpy.ListRasters("*","TIF")#Another method to get a list in order to get only the TIFF files
for TIFFtoReproj in reprojList:
    print (TIFFtoReproj)
    
for data in reprojList:
	print "Processing reprojection and resampling: "+data
	outRaster ="51N_"+data
	arcpy.ProjectRaster_management(data,outRaster,sr, cell_size="4000 4000",
        resampling_type="NEAREST")
	print "Rejoction and resampling complete"
	
####################Processing zonal statistics on a new list

arcpy.CheckOutExtension("Spatial")
#Set local variables
in_zone_data = 'grid51N.shp'
zone_field = 'id'

raster = arcpy.ListRasters("51N*","ALL")
for stat in raster:
    print (stat)

outDir = r"E:\climat_data"

for stat in raster:
        raster_name = os.path.basename(stat).rstrip(os.path.splitext(stat)[1])
        outTable = os.path.join(outDir,"tb_"+raster_name+".dbf")
        arcpy.gp.ZonalStatistics_sa(in_zone_data, zone_field, stat, outTable, "MEAN", "NODATA")
        print 'Done'+raster_name

##Traceback (most recent call last):
##  File "E:\script.py", line 60, in <module>
##    arcpy.gp.ZonalStatistics_sa(in_zone_data, zone_field, stat, outTable, "MEAN", "NODATA")
##  File "C:\Program Files (x86)\ArcGIS\Desktop10.7\ArcPy\arcpy\geoprocessing\_base.py", line 510, in <lambda>
##    return lambda *args: val(*gp_fixargs(args, True))
##ExecuteError: Failed to execute. Parameters are not valid.
##ERROR 000876: Output raster: E:\climat_data\tb_51N_Precip_2014-01-15.dbf's extension is invalid for the output raster format.
##Failed to execute (ZonalStatistics).


for table in arcpy.ListTables():
    print (table)

#####Joining the ".dbf" files to the attribute table of the GRID
##joinList = arcpy.ListeTables("*",dBASE)#The new list contains the previous outputs witch are database (.dbf) 
##for data in joinList:         
##    print (data)                         
##
##for data in joinList:
##    outField = "IDjoined_"+data            

##    print "Processing join field:"+data
##    arcpy.JoinField_management(in_data=grid51N, join_table=data, join_flied="id", fields="")
##    print "Join to the grid complete"
‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

Actually I've indicated my workspace as the place where I get all my input witch is a folder. I would have liked my .dbf documents to go in that same folder but it doesn't work.
The code worked when I removed the +".dbf" on line 59 but I only got .tif in the end. Which makes sense.

Thank you for your help.

Re: ZonalStatisticsAsTable

DanPatterson — Tue, 28 Jul 2020 21:49:35 GMT

Are you using arcmap or pro? your print statements aren't consistent with python 3.

Justd as an FYI..., esri grids can't begin with a number

"51N_"+data

so that would have been an issue without the *.tif file extension.

Re: ZonalStatisticsAsTable

IrèneDaubard — Sat, 11 Dec 2021 15:37:01 GMT

I'm using arcgis pro...

I made a model builder and export the script just to know the different tools I've to use.

# -*- coding: utf-8 -*-
"""Generated by ArcGIS ModelBuilder on: 2020-07-24 11:41:30
All ModelBuilder functionality may not be exported. Edits may be required for equivalency with the original model.
"""

import arcpy

# To allow overwriting the outputs change the overwrite option to true.
arcpy.env.overwriteOutput = False

# Local variables:
grid51N_shp__4_ = r"H:\CG_mosq\grid\grid51N.shp"
grid51N_shp__10_ = grid51N_shp__4_
grid51N_shp = r"H:\CG_mosq\grid\grid51N.shp"
hourlyPrecipRate_2014_01_01_2014_01_14_tif__3_ = r"H:\CG_mosq\climat_data\precipitation\hourlyPrecipRate_2014-01-01_2014-01-14\hourlyPrecipRate_2014-01-01_2014-01-14.tif"
hourlyPrecipRate_20140101_20 = r"H:\CG_mosq\MyProject_arcgispro_modelbuilder\MyProject_arcgispro_modelbuilder.gdb\hourlyPrecipRate_20140101_202"
hourlyPrecipRate_20140101_201 = r"H:\CG_mosq\MyProject_arcgispro_modelbuilder\MyProject_arcgispro_modelbuilder.gdb\hourlyPrecipRate_20140101_203"
precip_01_01_reproj_resampled_grid_dbf = r"H:\CG_mosq\climat_data\precipitation\precip-01-01_reproj_resampled_grid.dbf"
grid51N_shp__5_ = r"H:\CG_mosq\grid\grid51N.shp"
grid51N_shp__7_ = grid51N_shp__5_
grid51N_shp__2_ = r"H:\CG_mosq\grid\grid51N.shp"
Tair_f_inst_2014_01_01_2014_01_14_tif = r"H:\CG_mosq\climat_data\air_temperature\Tair_f_inst_2014-01-01_2014-01-14\Tair_f_inst_2014-01-01_2014-01-14.tif"
Tair_f_inst_20140101_2014011 = r"H:\CG_mosq\MyProject_arcgispro_modelbuilder\MyProject_arcgispro_modelbuilder.gdb\Tair_f_inst_20140101_2014011"
Tair_f_inst_20140101_20140111 = r"H:\CG_mosq\MyProject_arcgispro_modelbuilder\MyProject_arcgispro_modelbuilder.gdb\Tair_f_inst_20140101_20140111"
tair_01_01_reproj_resampled_grid_dbf = r"H:\CG_mosq\climat_data\air_temperature\tair-01-01_reproj_resampled_grid.dbf"
grid51N_shp__6_ = r"H:\CG_mosq\grid\grid51N.shp"
grid51N_shp__9_ = grid51N_shp__6_
grid51N_shp__3_ = r"H:\CG_mosq\grid\grid51N.shp"
v3f0acee6d99dbb08f5a84d3bdddf49c4_EVI_median_tif__2_ = r"H:\CG_mosq\climat_data\vegetation_index\EVI_2014-01-01_2014-01-14\3f0acee6d99dbb08f5a84d3bdddf49c4.EVI_median.tif"
c3f0acee6d99dbb08f5a84d3bddd = r"H:\CG_mosq\MyProject_arcgispro_modelbuilder\MyProject_arcgispro_modelbuilder.gdb\c3f0acee6d99dbb08f5a84d3bddd"
c3f0acee6d99dbb08f5a84d3bddd1 = r"H:\CG_mosq\MyProject_arcgispro_modelbuilder\MyProject_arcgispro_modelbuilder.gdb\c3f0acee6d99dbb08f5a84d3bddd1"
evi_01_01_reproj_resampled_grid_dbf = r"H:\CG_mosq\climat_data\vegetation_index\evi-01-01_reproj_resampled_grid.dbf"
sr = arcpy.SpatialReference(32651)#"WGS 84 UTM zone 51N"

# Process: Project Raster (2)
arcpy.ProjectRaster_management(in_raster=hourlyPrecipRate_2014_01_01_2014_01_14_tif__3_, out_raster=hourlyPrecipRate_20140101_20, out_coor_system="PROJCS['WGS_1984_UTM_Zone_51N',GEOGCS['GCS_WGS_1984',DATUM['D_WGS_1984',SPHEROID['WGS_1984',6378137.0,298.257223563]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]],PROJECTION['Transverse_Mercator'],PARAMETER['False_Easting',500000.0],PARAMETER['False_Northing',0.0],PARAMETER['Central_Meridian',123.0],PARAMETER['Scale_Factor',0.9996],PARAMETER['Latitude_Of_Origin',0.0],UNIT['Meter',1.0]]", resampling_type="NEAREST", cell_size="11171.1925206363 11171.1925206363", geographic_transform="", Registration_Point="", in_coor_system="GEOGCS['GCS_WGS_1984',DATUM['D_WGS_1984',SPHEROID['WGS_1984',6378137.0,298.257223563]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]]", vertical="NO_VERTICAL")

# Process: Resample (2)
arcpy.Resample_management(in_raster=hourlyPrecipRate_20140101_20, out_raster=hourlyPrecipRate_20140101_201, cell_size="4000 4000", resampling_type="NEAREST")

# Process: Zonal Statistics as Table (2)
arcpy.gp.ZonalStatisticsAsTable_sa(in_zone_data=grid51N_shp, zone_field="id", in_value_raster=hourlyPrecipRate_20140101_201, out_table=precip_01_01_reproj_resampled_grid_dbf, ignore_nodata="DATA", statistics_type="MEAN")

# Process: Join Field
arcpy.JoinField_management(in_data=grid51N_shp__4_, in_field="id", join_table=precip_01_01_reproj_resampled_grid_dbf, join_field="MEAN", fields="")

# Process: Project Raster
arcpy.ProjectRaster_management(in_raster=Tair_f_inst_2014_01_01_2014_01_14_tif, out_raster=Tair_f_inst_20140101_2014011, out_coor_system=sr)

# Process: Resample
arcpy.Resample_management(in_raster=Tair_f_inst_20140101_2014011, out_raster=Tair_f_inst_20140101_20140111, cell_size="4000 4000", resampling_type="NEAREST")

# Process: Zonal Statistics as Table
arcpy.gp.ZonalStatisticsAsTable_sa(in_zone_data=grid51N_shp__2_, zone_field="id", in_value_raster=Tair_f_inst_20140101_20140111, out_table=tair_01_01_reproj_resampled_grid_dbf, ignore_nodata="DATA", statistics_type="MEAN")

 Process: Join Field (2)
arcpy.JoinField_management(in_data=grid51N_shp__5_, in_field="id", join_table=tair_01_01_reproj_resampled_grid_dbf, join_field="id", fields="")

# Process: Project Raster (3)
arcpy.ProjectRaster_management(in_raster=v3f0acee6d99dbb08f5a84d3bdddf49c4_EVI_median_tif__2_, out_raster=c3f0acee6d99dbb08f5a84d3bddd, out_coor_system="PROJCS['WGS_1984_UTM_Zone_51N',GEOGCS['GCS_WGS_1984',DATUM['D_WGS_1984',SPHEROID['WGS_1984',6378137.0,298.257223563]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]],PROJECTION['Transverse_Mercator'],PARAMETER['False_Easting',500000.0],PARAMETER['False_Northing',0.0],PARAMETER['Central_Meridian',123.0],PARAMETER['Scale_Factor',0.9996],PARAMETER['Latitude_Of_Origin',0.0],UNIT['Meter',1.0]]", resampling_type="NEAREST", cell_size="464.671253451976 464.671253451976", geographic_transform="", Registration_Point="", in_coor_system="GEOGCS['GCS_WGS_1984',DATUM['D_WGS_1984',SPHEROID['WGS_1984',6378137.0,298.257223563]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]]", vertical="NO_VERTICAL")

# Process: Resample (3)
arcpy.Resample_management(in_raster=c3f0acee6d99dbb08f5a84d3bddd, out_raster=c3f0acee6d99dbb08f5a84d3bddd1, cell_size="4000 4000", resampling_type="NEAREST")

# Process: Zonal Statistics as Table (3)
arcpy.gp.ZonalStatisticsAsTable_sa(in_zone_data=grid51N_shp__3_, zone_field="id", in_value_raster=c3f0acee6d99dbb08f5a84d3bddd1, out_table=evi_01_01_reproj_resampled_grid_dbf, ignore_nodata="DATA", statistics_type="MEAN")

# Process: Join Field (3)
arcpy.JoinField_management(in_data=grid51N_shp__6_, in_field="id", join_table=evi_01_01_reproj_resampled_grid_dbf, join_field="id", fields="")
‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

I'll switch the "51N_"+data and get back to you just in case... I currently no longer have access to my office.

Thank you

Re: ZonalStatisticsAsTable

Luke_Pinner — Fri, 31 Jul 2020 21:47:04 GMT

arcpy.sa.ZonalStatistics outputs a raster not a table. You need arcpy.sa.ZonalStatisticsAsTable