Offset polygonZ in arcpy

ToddHawley · ‎10-17-2012

Arcinfo 10.0 SP5
Python 2.6
Win 7 Ultimate SP1
Hello,

I have been trying to find a way to offset a polygonZ in arcmap. My application involves earthwork where a single bounding polygonZ is to be offset inward at 200ft horizontally and vertically at 100ft after which there will be a flat 20 ft bench and the process repeats. My ultimate goal is to create a TIN that I can Contour.

I created the code below which creates something close to what I need however there is a problem. I have included three files "ArcScene, PlanView, CrossSection'. The first two are 3D and 2D viewpoints of the output from the included code. The last is a cross section of what the output should look like. Upon inspection of the contours in the first two files you will see that the benches are sloped across their width which needs to be flat. The benches however need to be sloped along their length. This is due to the fact that I was forced to interpolate each buffer to the DEM and then raise it to the appropriate elevation using the Plus3D() tool. The bounding polygon mentioned above was interpolated by the same dem which is represented in the PlanView drawing by the labeled contour lines which shows how the DEM is a plane on a slight slope. So in short I need to start with a 3D line or polygon and offset while preserving the 3D values so the benches can be flat across their width. I am also open to any raster based methods anyone may have to offer.

If anyone has any recommendations of how this could be accomplished in ArcMap I would greatly appreciate your ideas and suggestions.

Thanks
Todd Hawley

import arcpy, sys, string
from arcpy import env
from arcpy import AddMessage
from arcpy import SearchCursor
from arcpy import UpdateCursor
from arcpy import AddError
env.workspace = "c:/1/Program.gdb"

env.overwriteOutput = True

#Input Polygon
P_CropLine = sys.argv[1]
#DEM
P_DEM = sys.argv[2]
# Distance Values for Berm, Ditch, and Benches
P_BermHeight = float(sys.argv[3])
P_BermWidth = float(sys.argv[4])
P_DitchDepth = float(sys.argv[5])
P_DitchWidth = float(sys.argv[6])
P_VertInter = float(sys.argv[7])
P_BenchWidth = float(sys.argv[8])
P_TotalHeight = float(sys.argv[9])
#Variable for Generalize Tool
P_Gen = float(sys.argv[10])

P_Output = sys.argv[11]

# Begin Function section
#Used to Create new Fields in a table
def FieldNameExists(tablename, fieldname, fieldtype):

    try:
        fieldList = arcpy.ListFields(tablename)
        fieldPresent = 0
        for field in fieldList:
            if field.baseName == fieldname:
                fieldPresent = 1

        if fieldPresent ==0:
            arcpy.AddField_management(tablename, fieldname, fieldtype,"","","","","NULLABLE", "")
        elif fieldPresent ==1:
            arcpy.DeleteField_management(tablename, fieldname)
            arcpy.AddField_management(tablename, fieldname, fieldtype,"","","","", "NULLABLE" , "")
        del fieldList
        del field

    except Exception as e:
        fException('Error in FieldNameExists Function',e)

#Called during when exception is raised
def fException (string, ex):
    try:
        
        AddError("\n!!!" + string + "!!!\n")
        import traceback, sys
        tb = sys.exc_info()[2]
        AddError("Line " + str(tb.tb_lineno))
        AddMessage('\n' + str(ex.message))
        sys.exit(0)
    except Exception as e:
        AddMessage("!!!Error in fException Function!!!\n")
        import traceback, sys
        tb = sys.exc_info()[2]
        AddMessage("Line \n" + str(tb.tb_lineno))
        AddMessage("\n" + str(e.message))
        sys.exit(0)
#add various ID fields
def fAddPointID(tbl):
    try:
        FieldNameExists(tbl,'PT_ID','short')
        Ucur_PT_ID = UpdateCursor(tbl)
        i = 0
        for item in Ucur_PT_ID:
            i += 1
            item.setValue('PT_ID',i)
            Ucur_PT_ID.updateRow(item)
        del item
        del Ucur_PT_ID
    except Exception as e:
        fException('Error while creating PT_ID',e)
        
def fAddPartID(tbl):
    try:
        FieldNameExists(tbl,'Prt_ID','short')
        Ucur_Prt_ID = UpdateCursor(tbl)
        i = 0
        for item in Ucur_Prt_ID:
            i += 1
            item.setValue('Prt_ID',i)
            Ucur_Prt_ID.updateRow(item)
        del item
        del Ucur_Prt_ID
    except Exception as e:
        fException('Error while creating PT_ID',e)

def fOutput(infc):
    
    try:
        arcpy.CopyFeatures_management(infc, P_Output)

    except Exception as e:
        fException('Error in fOutput Function', e)        
# Clear features from feature class
def fClearFeat(tbl):
    try:
        arcpy.DeleteFeatures_management(tbl)
        
       
    
    except Exception as e:
        fException('Error spliting crop line',e)    
# uses a feature class that new buffer outputs are appended to.
def fCollect(tbl,outfeat,IND):
    try:
        
        Ucur_Part = UpdateCursor(tbl)
        for item in Ucur_Part:
            item.setValue('Part_ID',IND)
            Ucur_Part.updateRow(item)
      
        del Ucur_Part
                
       
        Scur_Col = SearchCursor(tbl)
        Icur_Col = arcpy.InsertCursor(outfeat)
        
        for item in Scur_Col:
            
            Icur_Col.insertRow(item)
        del Icur_Col
 
        
        del Scur_Col
       
    except Exception as e:
        fException('Error adding to collection', e)
    
#End Function Section

try:
    #Variable to provide ID for all polygons at the same vertical stage
    partID = 1
    #Temp feature classes
    IM_Input = arcpy.CreateFeatureclass_management('in_memory','input','polygon')
    IM_OS = arcpy.CreateFeatureclass_management('in_memory', 'OS', 'polygon')
    IM_Pri =  arcpy.CreateFeatureclass_management('in_memory', 'Pri', 'polygon','','','disabled')
    IM_3D = arcpy.CreateFeatureclass_management('in_memory', 'F3D','polygon')
    
    IM_Collect = arcpy.CreateFeatureclass_management('in_memory','Collect', 'polygon','','','enabled')
    IM_Out = arcpy.CreateFeatureclass_management('in_memory', 'out','polyline')

    
    
   
    FieldNameExists(IM_Collect,'Part_ID', 'short')
    FieldNameExists(IM_Input,'Part_ID', 'short')
    FieldNameExists(IM_Pri,'Part_ID,','short')

    
    #Copy Input, InterpolateShap and add to collection
    
    arcpy.CopyFeatures_management(P_CropLine, IM_Pri)
    
    arcpy.InterpolateShape_3d(P_DEM,IM_Pri,IM_3D,0.01,'','linear','TRUE')

    fClearFeat(IM_Pri)
    FieldNameExists(IM_3D,'Part_ID','short')
    
    fCollect(IM_3D,IM_Collect,partID)
    
    arcpy.CopyFeatures_management(IM_3D,IM_Pri)
  
    fClearFeat(IM_3D)
    
    # Create outside top of berm
    
    partID += 1
    arcpy.Buffer_analysis(IM_Pri,IM_OS,P_BermHeight *-2)
    fClearFeat(IM_Pri)
  
    arcpy.Plus_3d(P_DEM,P_BermHeight,'rast1')
    arcpy.InterpolateShape_3d('rast1',IM_OS,IM_3D,0.01,'','linear','TRUE')

    fClearFeat(IM_OS)
    fCollect(IM_3D,IM_Collect,partID)
    
    arcpy.CopyFeatures_management(IM_3D,IM_Pri)
    fClearFeat(IM_3D)
    
    #Create inside top of berm
    
    partID += 1
    arcpy.Buffer_analysis(IM_Pri,IM_OS,P_BermWidth*-1)
    fClearFeat(IM_Pri)
    
   
    arcpy.InterpolateShape_3d('rast1',IM_OS,IM_3D,0.01,'','linear','TRUE')
    
    fClearFeat(IM_OS)
    fCollect(IM_3D,IM_Collect,partID)
    arcpy.management.Delete('rast1')
    
    arcpy.CopyFeatures_management(IM_3D, IM_Pri)
    fClearFeat(IM_3D)
    
    #Create outside bottom of ditch
    
    partID +=1
    arcpy.Buffer_analysis(IM_Pri,IM_OS,P_BermHeight * -2)
    fClearFeat(IM_Pri)
  
    
    arcpy.InterpolateShape_3d(P_DEM,IM_OS,IM_3D,0.01,'','linear','TRUE')
    
    fClearFeat(IM_OS)
    fCollect(IM_3D,IM_Collect,partID)
    
    arcpy.CopyFeatures_management(IM_3D, IM_Pri)
    fClearFeat(IM_3D)
    
    #Create inside bottom of ditch
    
    partID +=1
    arcpy.Buffer_analysis(IM_Pri,IM_OS,P_DitchWidth * -1)
    fClearFeat(IM_Pri)
  
    arcpy.InterpolateShape_3d(P_DEM,IM_OS,IM_3D,0.01,'','linear','TRUE')
    
    fClearFeat(IM_OS)
    fCollect(IM_3D,IM_Collect,partID)
    
    arcpy.CopyFeatures_management(IM_3D, IM_Pri)
    fClearFeat(IM_3D)   
 
    # Begin slope section
    CumHeight = 0
    Lift = 1
    level = 0
    #Buffer polygon and interpolate to appropriate height until elevation reaches P_Total_Height Var
    while CumHeight < P_TotalHeight:
        # Buffer until 100 lift is created then create flat bench
        while level < 100:
            partID+= 1
            level += P_VertInter
            CumHeight += P_VertInter
            Lift += 1
            arcpy.Buffer_analysis(IM_Pri,IM_OS,P_VertInter * -2)
            fClearFeat(IM_Pri)
          
            arcpy.Plus_3d(P_DEM,Lift * P_VertInter,'rast1')
            arcpy.InterpolateShape_3d('rast1',IM_OS,IM_3D,0.01,'','linear','TRUE')
            arcpy.management.Delete('rast1')
            fClearFeat(IM_OS)
            fCollect(IM_3D,IM_Collect,partID)
          
            arcpy.CopyFeatures_management(IM_3D, IM_Pri)
            fClearFeat(IM_3D)
            
        #Create flat bench
        partID +=1
        level = 0
        arcpy.Buffer_analysis(IM_Pri,IM_OS,P_BenchWidth * -1)
        fClearFeat(IM_Pri)
     
        arcpy.Plus_3d(P_DEM,Lift * P_VertInter,'rast1')
        arcpy.InterpolateShape_3d('rast1', IM_OS, IM_3D,0.01,'','linear','TRUE')
        arcpy.management.Delete('rast1')
        fClearFeat(IM_OS)
        fCollect(IM_3D, IM_Collect,partID)
        arcpy.CopyFeatures_management(IM_3D,IM_Pri)
        fClearFeat(IM_3D)

    arcpy.FeatureToLine_management(IM_Collect, IM_Out)
    
    fOutput(IM_Out)

except Exception as e:
    fException('Error spliting crop line',e)