Inside Buffer

[LornaMurison]

Hi everyone,
I'm using 9.2 SP6.
I have a a fairly complex lakes feature class (shown in attached picture).  I am trying to create a buffer 20 meters inside each feature and none of the methods seems to be working.
The buffer tool with a negative distance produces some kind of error before the process can be completed.
The buffer wizard tool produces the result shown in brown in the other attached image.  First of all it's the blue area that I want to keep, so I would have to add the extra erase step, and second, it completely ignored the more complicated and larger polygon on top.
I used ArcMap advanced settings utility to change the buffer wizard method to "feature optimized coordinate system" which produced the same results.

Does anyone know how, using a script, tool, or anything, how to create a buffer inside a polygon, keeping in mind that this process will have to be run MANY times.

Thank-You!

[ChrisSnyder]

Try:

1. Breaking the input FC into singlepart geometry using the MultipartToSinglepart tool
2. Running the RepairGeometry tool. Are there any "errors"?
3. Then run the negative buffer.

Does it work then?

The issue of course is the extreemly complex geometry of the lake polygons. If the buffer tools fails on the original FC, maybe try writing a script that will loop through the features and make a buffer for each feature in a try/except block. That way you will be able to figure out the offending features (maybe it's only a few of them?).

[Anonymous User]

Original User: Lady_Jane

Thank-you for your reply, Chris.

I have found 2 polygons that the buffer tool is ignoring.  If I select all polygons except for those 2, then the process works.  Of course, they are the two most complicated polygons in the shapefile.
Does anyone know of a way to simplify the polygons before buffering, or of allowing complex polygons to be buffered?

[Anonymous User]

Original User: clm42

Well there is a very tedious method of the trace tool. Just set it to off set the buffer distance and trace around the lake. Than finish only that part of the polygon and trace around again with an off set of zero. When you finish the polygon it will be "buffered". I just had to do this very thing to the entire coastline of my county because it wouldnt buffer the crazy geometry.

[DanPatterson_Retired]

run a Check Geometry and if needed a Repair Geometry on the input file to rule out geometry errors.
If the shapes are multipart shapes, try using the multiparttosinglepart tool on the file prior to interior buffering

[Anonymous User]

Original User: Lady_Jane

Unfortunately it is not a geometry issue.  I have checked geometry, found no errors, ran repair geometry anyway, and the large complex polygon is still causing an error.
I am going to try converting to a coverage first.
Thanks everyone, still looking for ideas 🙂

[ChrisSnyder]

If workstation fails (try the regionbuffer tool), you might try this in ArcGIS:

1. Convert your lake polys to a line (FeatureToLine tool)
2. Buffer the lines 20m (on both sides, round ends). If the shoreline buffer fails try:
   a. Build a fishnet polygon (CreateFishnet tool and then FeatureToPolygon) of sufficient # of cells to break the large contiguous line features into smaller segments. This might take a couple tries to get right.
   b Buffer the "chopped up" lines.
3. Erase the buffer from the lake poly. If the erase fails try:
   a. Intersect the lake polys with the fishnet you created for the buffers
   b. Try the intersect again with the "chopped up" lake polys

Generally, the issue is that ArcGIS chokes when you try to do some sort of geoprocessing operation with a feature that has a kajillion vertices. The reason being that all the feature's x,y, pairs need to get loaded into memory in order to do the operation, and often the number of vertrices in a complex feature exceed the limits of the ArcGIS algorithm. The solution is to break your input features into small enough pieces to allow the algorithms to run without maxing out the memory limitations.

[Anonymous User]

Original User: Lady_Jane

I don't know what you mean by fishnet, I've never used that tool before, I will have to look into it.
Right now I have run the simplify polygon tool and the output still has sufficient detail.  I am now attempting to buffer the simplified features which has been going on for almost 3 hours now.  I take this as a good sign though, since previously it would crash after a minute or so.
Provided this buffer works, my plan will be to calculate the number of vertices on each feature, select those over a given threshold and simplify them before performing the buffer.
If this doesn't work I will look into your suggestions.
Thank-You!

[ChrisSnyder]

The fishnet tool is basically a way to make a polygon layer that looks like, well, a fishnet that has adjacent rectangualr (or square) polygons (like a checkerboard). I do a lot of large-dataset processing in Python, and the Fishnet tool is always my step #1 in being able to break up large datasets to get them to work with the ESRI tools.

For example:

# Description
# -----------
# This creates an area of interest feature class, a few buffers of the area of interest, and a fishnet polygon to use for tile-by-tile psudo-parallel processing. 
# Written for Python version: 2.5.1 (PythonWin)
# Written for ArcGIS version: 9.3.1 (should be forward compatible).
# Author: Chris Snyder, WA Department of Natural Resources, chris.snyder(at)wadnr.gov
# Date created:
# Last Updated: 20090923, csny490

try:
    #Import system modules and creates the Geoprocessor object (the v9.2 way)
    import sys, string, os, glob, shutil, time, traceback, arcgisscripting
    gp = arcgisscripting.create(9.3)
        
    #Defines some functions
    def showGpMessage():
        try:
            print >> open(logFile, 'a'), "\n" + gp.GetMessages() + "\n"
            print "\n" + gp.GetMessages() + "\n"
        except:
            pass
    def showPyMessage():
        try:
            print >> open(logFile, 'a'), str(time.ctime()) + " - " + str(message)
            print str(time.ctime()) + " - " + str(message)
        except:
            pass
        
    #Sets a date/time stamp variable in the format yyyymmddhhmmss (example: 20060330132345)
    dateTimeStamp = time.strftime('%Y%m%d%H%M%S')

    #Specifies the root variable, makes the logFile variable, and does some error checking...
    root = sys.argv[1]
    if os.path.exists(root)== False:
        print "Specified root directory: " + root + " does not exist... Bailing out!"
        sys.exit()
    scriptName = sys.argv[0].split("\\")[len(sys.argv[0].split("\\")) - 1][0:-3] #Gets the name of the script without the .py extension  
    logFile = root + "\\log_files\\" + scriptName + "_" + dateTimeStamp[:8] + ".log" #Creates the logFile variable
    if os.path.exists(root + "\\log_files") == False:
        os.mkdir(root + "\\log_files")
    if os.path.exists(logFile)== True:
        os.remove(logFile)
        message = "Deleting log file with the same name and datestamp... Recreating " + logFile; showPyMessage()
    workspaceDir = "overlay_index"
    if os.path.exists(root + "\\" + workspaceDir)== False:
        message = "Creating layer directory: " + root + "\\" + workspaceDir; showPyMessage()
        os.mkdir(root + "\\" + workspaceDir)

    #Attempts to check out the highest grade license available...
    if gp.CheckProduct("ArcInfo") == "Available":
        gp.SetProduct("ArcInfo")
    elif gp.CheckProduct("ArcEditor") == "Available":
        gp.SetProduct("ArcEditor")
    elif gp.CheckProduct("ArcView") == "Available":
        gp.SetProduct("ArcView")
    else:
        messsage =  "No ArcGis licesnses are available... Exiting script!"; showPyMessage(); sys.exit()
    message =  "Selected an " + gp.ProductInfo() + " license"; showPyMessage()

    #Sets some environment variables    
    gp.overwriteoutput = True
    gp.loghistory = False
    sr = 'PROJCS["NAD_1983_HARN_StatePlane_Washington_South_FIPS_4602_Feet",GEOGCS["GCS_North_American_1983_HARN",DATUM["D_North_American_1983_HARN",SPHEROID["GRS_1980",6378137.0,298.257222101]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Lambert_Conformal_Conic"],PARAMETER["False_Easting",1640416.666666667],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",-120.5],PARAMETER["Standard_Parallel_1",45.83333333333334],PARAMETER["Standard_Parallel_2",47.33333333333334],PARAMETER["Latitude_Of_Origin",45.33333333333334],UNIT["Foot_US",0.3048006096012192]]'
    gp.OutputCoordinateSystem = sr
    gp.XYTolerance = "0.001 METERS"
    gp.XYResolution = "0.0005 METERS"
    gp.workspace = root + "\\" + workspaceDir

    #Process: Creates a FGDB to hold all the spatial layers
    fgdbName = "overlay_index"
    fgdbPath = gp.workspace + "\\" + fgdbName + ".gdb"
    gp.CreateFileGDB_management(gp.workspace, fgdbName); showGpMessage()

    #Process: Makes a FC of the PLS polys we want
    plsFC = root + "\\gis_layers\\pls.gdb\\pls"
    plsFL = "pls_feature_layer"
    gp.MakeFeatureLayer_management(plsFC, plsFL, "SUR_OWN_CD > 0 OR TIM_OWN_CD > 0", "", ""); showGpMessage()
    dnrMgmtAllFC = fgdbPath + "\\dnr_mgmt_all"
    gp.Dissolve_management(plsFL, dnrMgmtAllFC, "", "", "SINGLE_PART"); showGpMessage()

    #Process: Creates a fishnet layer (the basis of the index polygons)
    fishnetLine1FC = fgdbPath + "\\fishnet_line1"
    dsc = gp.describe(dnrMgmtAllFC)
    xMin = dsc.extent.xmin
    yMin = dsc.extent.ymin
    xMax = dsc.extent.xmax
    yMax = dsc.extent.ymax
    bufferWidth = 1 #Buffer width expands the the extent of the fishnet by that many map units (to make sure it extends past DNR lands a bit)
    numberOfRows = 6
    numberOfColumns = 6
    gp.CreateFishnet_management(fishnetLine1FC, str(xMin - bufferWidth) + " " + str(yMin - bufferWidth), str(xMin - bufferWidth) + " " + str(yMin), "0", "0", numberOfRows, numberOfColumns, str(xMax + bufferWidth) + " " + str(yMax + bufferWidth), "NO_LABELS", ""); showGpMessage()
    fishnet1FC = fgdbPath + "\\fishnet_poly1"
    gp.FeatureToPolygon_management(fishnetLine1FC, fishnet1FC, "", "", ""); showGpMessage()
    gp.Delete_management(fishnetLine1FC, ""); showGpMessage()

    #Process: Creates a 2nd fishnet layer (that is more dense than the original)
    fishnetLine2FC = fgdbPath + "\\fishnet_line2"
    densificationFactor = 12 #this densification factor ABSOLUTELY NEEDS to be a multiple of 2 (e.g. 2,4,6,8,etc)
    gp.CreateFishnet_management(fishnetLine2FC, str(xMin - bufferWidth) + " " + str(yMin - bufferWidth), str(xMin - bufferWidth) + " " + str(yMin), "0", "0", str(numberOfRows * densificationFactor), str(numberOfColumns * densificationFactor), str(xMax + bufferWidth) + " " + str(yMax + bufferWidth), "NO_LABELS", ""); showGpMessage()
    fishnet2FC = fgdbPath + "\\fishnet_poly2"
    gp.FeatureToPolygon_management(fishnetLine2FC, fishnet2FC, "", "", ""); showGpMessage()
    gp.Delete_management(fishnetLine2FC, ""); showGpMessage()

    #Process: Updates fishnet1FC with fishnet2FC for specified areas
    hcpUnitFC = root + "\\gis_layers\\hcpunit.gdb\\hcpunit"
    oesfFC = fgdbPath + "\\oesf"
    gp.Select_analysis(hcpUnitFC, oesfFC, "HCPUNIT_ID = 1", "", ""); showGpMessage()
    fishnet1FL = "fishnet_1_feature_layer"
    gp.MakeFeatureLayer_management(fishnet1FC, fishnet1FL, "", "", ""); showGpMessage()
    gp.SelectLayerByLocation_management(fishnet1FL, "INTERSECT", oesfFC, "", "NEW_SELECTION"); showGpMessage()
    fishnet2FL = "fishnet_2_feature_layer"
    gp.MakeFeatureLayer_management(fishnet2FC, fishnet2FL, "", "", ""); showGpMessage()
    gp.SelectLayerByLocation_management(fishnet2FL, "HAVE_THEIR_CENTER_IN", fishnet1FL, "", "NEW_SELECTION"); showGpMessage()
    gp.Delete_management(oesfFC, ""); showGpMessage()
    fishnetFC = fgdbPath + "\\fishnet_poly"
    gp.Update_analysis(fishnet1FC, fishnet2FL, fishnetFC, "BORDERS", ""); showGpMessage()

    #Process: Make a line FC version of fishnetFC
    fishnetLineFC = fgdbPath + "\\fishnet_lines"
    gp.FeatureToLine_management(fishnetFC, fishnetLineFC, "", "NO_ATTRIBUTES"); showGpMessage()
    
    #Process: Unions everything together
    union1FC = fgdbPath + "\\union1"
    gp.Union_analysis(fishnetFC + ";" + dnrMgmtAllFC, union1FC, "ONLY_FID", "", "GAPS"); showGpMessage()

    #BLAH BLAH BLAH - I had to trim this part out to post it to the forum (stupid 10,000 charater limit!!!)


    #Process: Compacts the PGDB
    gp.Compact_management(fgdbPath); showGpMessage()

    #Process: Creates a .txt file that indicates to other scripts how many tiles there are (this is used so that the run_union_master script doesn't have to use the gp and waste memory)    
    print >> open(root + "\\log_files\\" + "therearethismanytiles_" + str(numberOfTiles) + ".txt", 'a'), "therearethismanytiles_" + str(numberOfTiles)
    
    message = "ALL DONE!"; showPyMessage()
    print >> open(root + "\\log_files\\" + scriptName + "_isalldone.txt", 'a'), scriptName + "_isalldone.txt"
except:
    print >> open(root + "\\log_files\\" + scriptName + "_bombed.txt", 'a'), scriptName + "_bombed.txt"
    message = "\n*** LAST GEOPROCESSOR MESSAGE (may not be source of the error)***"; showPyMessage()
    showGpMessage()
    message = "\n*** PYTHON ERRORS *** "; showPyMessage()
    message = "Python Traceback Info: " + traceback.format_tb(sys.exc_info()[2])[0]; showPyMessage()
    message = "Python Error Info: " +  str(sys.exc_type)+ ": " + str(sys.exc_value) + "\n"; showPyMessage()
    message = "\n*** PYTHON LOCAL VARIABLE LIST ***"; showPyMessage()
    variableCounter = 0                      
    while variableCounter < len(locals()):
        message =  str(list(locals())[variableCounter]) + " = " + str(locals()[list(locals())[variableCounter]]); showPyMessage()
        variableCounter = variableCounter + 1
    showPyMessage()

[Anonymous User]

Original User: whuber

Thanks everyone, still looking for ideas 🙂

Part of the complexity of the large lakes lies in their many islands, Lorna, so here's a thought: instead of negatively buffering the lakes, why not positively buffer their complement?  Normally the complement is easily obtained: create a feature surrounding the lakes (a rectangle will do, but a feature drawn to closely surround the problem lakes would be better), form its union with the lakes, and remove the problem lake features from the union.  What's left is their complement.  In this case, each island will be a separate feature in the complement and will therefore get separately buffered, which might be easier for ArcGIS to do.  The big question is whether ArcGIS will be able to perform the union at all, but there is a hope: my experience is that buffering tends to have more bugs than other geoprocessing operations (it uses algorithms that differ in some important ways from the more basic set-theoretic operations of union and intersection).  Of course if you succeed in externally buffering the complement, you will then proceed to form its complement to produce the internal lake buffer you're looking for.

If this doesn't succeed, another approach is to convert the lakes to their polyline boundaries and attempt to buffer the polylines.  A subsequent union of that buffer with the original lake polygons--if it succeeds--will produce a mess, but in that mess will be the negative buffers you want.

Raster analysis offers another route that is likely to succeed, but its solutions will have some grid discretization error.  Whether that's acceptable or not depends on what you'll be doing with these buffers.