Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

I have problems while using the ZonalHistogram function in Python. It seems to give strange results, while, with the same code run in arcgis' python window everything works fine. I am trying to study the snow cover per elevation ranges, thus including in the function a MODIS image (showing snow cover) as the input raster, and a reclissified DEM as the input value raster. I tried to modify the reclassification many times, but still not solving the problem... Moreover, inverting the input raster with the input value raster (meaning, inverting the function parameters), everything works fine, but it's a mess since I'm working on thousands of images and having the MODIS not always the same number of columns (while being the DEM always the same). Hope someone could help me. Here is my script, so that you can have a look at it in order to better understand my request:

import arcpy, os, sys, shutil

newpath= r'G:\\Chile_2012\\MODIS10A2.V005\\provvisoria\\Zonal_Histogram\\temporary'
if not os.path.exists (newpath):
    os.makedirs (newpath)

arcpy.env.workspace= 'G:\\Chile_2012\\MODIS10A2.V005\\Cobertura_nival\\Cuencas_Copia\\campione'
arcpy.env.overwriteOutput= True
from arcpy.sa import *

listRaster= arcpy.ListRasters ()
elev_ranges= 'G:\\Chile_2012\\MODIS10A2.V005\\DEM\\elev_ranges'
out_dir= 'G:\\Chile_2012\\MODIS10A2.V005\\provvisoria\\Zonal_Histogram\\'

try:
    for raster in listRaster:
        arcpy.CheckOutExtension ('Spatial')
        ZonalHistogram (raster, 'Value', elev_ranges, newpath+ '\\'+ raster)

    arcpy.ClearEnvironment ('workspace')
    arcpy.env.workspace= newpath
    listTable = arcpy.ListTables ()
    print listTable

    arcpy.Merge_management (listTable, out_dir+ 'Rio_Aconcagua_Merge.dbf')

    shutil.rmtree (newpath)

except:
    print arcpy.GetMessages ()

I would add that my idea would be to work with layers (by means of functions that could transform rasters in layers) since the table returned in arcgis has a "LABEL" field (coming from the classification in symbology) which is the thing I would like to include in the final table and which is not included while working out of arcgis (in pywin), in order to know on what I am working on.

share|improve this question
add comment

1 Answer

up vote 1 down vote accepted

I solved the problem. Strangely, it seems the Zonal_Histogram function works with the symbology of a file, and not only with the attribute table. I found it out by comparing my results by assigning both a stretched symbology, and a classified one to the DEM. Arcgis everytime considered the representation of the file, giving a field called "LABEL" when there was a classified symbology, while giving the result I wanted to avoid (the same I had with pywin) with the stretched representation. By right-clicking on the file in the Table View and choosing the "SaveAsLayerFile" botton, the symbology is made permanent, so everytime I open the file, the layer is always represented with the classified symbology I chose (and which I don't have to reselect everytime), thus giving as output the attribute table I need. Python, working on the layer file that was assigned the permanent classified symbology, gives me the right result now! Guess why it's not written in the ESRI help...

share|improve this answer
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.