AnsweredAssumed Answered

Creating multipart polygons with ArcPy? Dissolve Failure?

Question asked by Chad84 on Aug 24, 2015
Latest reply on Jan 7, 2016 by Chad84

I need the functionality of the editor > merge command for use in a python script tool. I realize this function is not directly accessible through ArcPy, but I am hoping someone can help with a work around?

 

The Issue: I have a series of single part polygon features which do not touch or overlap, but have common field values. I need all these features to exist as multipart features for the random point generation to operate properly and consider all these regions equally. However, the dissolve function inconsistently fails leaving single part polygons.

 

Essentially, I am trying to go from single part geometries to multipart by creating feature layer selections and dissolving on that feature layer selection.

 

I have included a code snippet below. The code works and does exactly as required, but it is inconsistent. The dissolve fails when preceded by other scripts or when significant edits with debugging have been made to other parts of this script. However, when called from a new .mxd document it runs consistently well.

 

Any advice and support would be gratefully accepted!
Thanks!

 

#Custom exceptions.
class LicenseException(Exception):
    pass


#Import system modules.
import arcpy
import sys
import os
import traceback
from arcpy import env
from arcpy.sa import *


#Enable geoprocessor logging.
arcpy.SetLogHistory(True)


try:
    #---Check the installed ArcGIS version.---
    iv = arcpy.gp.GetInstallInfo()


    version = " ArcGIS %s : %s" % ('Version', iv.get('Version','Version Unknown'))


    if version.find("10.") > 0:
        ArcGIS10 = True
    else:
        ArcGIS10 = False
        arcpy.AddError("ERROR: This tool requires ArcGIS version 10.1 or greater...exiting script.")
        sys.exit("")


    del iv, version, ArcGIS10


    #---Check the ArcGIS administrative license.---
    if arcpy.CheckProduct("ArcInfo") == "AlreadyInitialized" or arcpy.CheckProduct("ArcInfo") == "Available":  
        pass
    else:
        raise LicenseException()


    #---Check out extension licenses.---
    if arcpy.gp.CheckExtension("spatial") == "Available":
        arcpy.CheckOutExtension("spatial")
    else:
        arcpy.AddError("ERROR: Spatial Analyst extension is unavailable. Please enable Spatial analyst and try again...exiting script.")
        sys.exit("")


    #---Get tool share folder file path.---
    mxd = arcpy.mapping.MapDocument("Current")
    mxdpath = mxd.filePath
    (filepath, filename) = os.path.split(mxdpath)
    MCP_Path = str(filepath)


    #---Set the geoprocessing environment.---
    workspace = MCP_Path + "\OutputData\ToolOut.gdb"
    scratchspace = MCP_Path + "\OutputData\ScratchFolder\Scratch.gdb"


    env.scratchWorkspace = scratchspace
    env.workspace = workspace
    env.overwriteOutput = True


    #---Read user input.---
    BI = arcpy.GetParameterAsText(0)


    BIdesc = arcpy.Describe(BI) 


    #---Execute Reclass.---
    DS_RC_filename = BIdesc.name + "_DS_RC"
    DS_RC_out = os.path.join(workspace, DS_RC_filename)


    DS_Reclass = Reclassify(BI, "Value", 
    RemapRange([[-1000,-75,8],[-75,-50,7],[-50,-35,6],[-35,-20,5],[-20,-12,4],[-12,-6,3],[-6,-3,2],[-3,-1,1],[-1,0,"NODATA"]]), "NODATA")
    DS_Reclass.save(DS_RC_out)


    #---Execute raster to polygon conversion.---
    DS_Polygon_filename = BIdesc.name + "_DS_Polygon"
    DS_Polygon = os.path.join(workspace, DS_Polygon_filename)


    arcpy.RasterToPolygon_conversion(DS_RC_out, DS_Polygon, "SIMPLIFY", "VALUE")


    if arcpy.Exists(DS_RC_out):
        arcpy.Delete_management(DS_RC_out)
    del DS_RC_filename, DS_Polygon_filename 


    #---Select polygons of areas >= 1000m2.---
    DS_PS_out_filename = BIdesc.name + "_DS_PSelect"
    DS_PS_out = os.path.join(workspace, DS_PS_out_filename)


    sql = '{0} >= {1}'.format("Shape_Area",'1000')
    if arcpy.Exists("DS_PS"):
        arcpy.Delete_management("DS_PS")
    DS_PS = arcpy.management.MakeFeatureLayer(DS_Polygon,"DS_PS",sql)
    arcpy.CopyFeatures_management("DS_PS", DS_PS_out) #Done to hopefully stabilize the following selection


    arcpy.Delete_management(DS_PS)
    del DS_PS_out_filename, sql


    #---------------------Perform Dissolve----------------------------------    
    DS_F_out_filename = BIdesc.name + "_DS"
    DS_F_out = os.path.join(workspace, DS_F_out_filename)


    arcpy.Dissolve_management(DS_PS_out,DS_F_out,"grid_code","#","MULTI_PART","DISSOLVE_LINES") 


    del BsMmaxS, sql


    #---Check to ensure dissolve function worked properly.---
    sql = '{0} = {1}'.format("grid_code",1)
    arcpy.MakeTableView_management(DS_F_out, "DS_F_out_TV",sql)
    count = int(arcpy.GetCount_management("DS_F_out_TV").getOutput(0))
    if count > 1:
        del sql
        if arcpy.Exists(DS_F_out):
            arcpy.Delete_management(DS_F_out)
        arcpy.Delete_management("DS_F_out_TV")
        arcpy.AddError("Geoprocessing Error due to memory limitations with 32bit processing. Assuming 64bit processing is unavaiable, please close your ArcMap instance, open it again and reprocess...exiting script.")
        #Above error prompt is my current best guess since the polygons can be quite complex and can contain many verticies. 
        sys.exit("")
    else:
        arcpy.Delete_management("DS_F_out_TV")
        del sql


    #Continue code here..

Outcomes