topic Re: Inside Buffer in Geoprocessing Questions

Inside Buffer

LornaMurison — Wed, 07 Jul 2010 13:29:19 GMT

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!

Re: Inside Buffer

ChrisSnyder — Wed, 07 Jul 2010 16:43:43 GMT

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?).

Re: Inside Buffer

Anonymous User — Wed, 07 Jul 2010 18:43:19 GMT

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?

Re: Inside Buffer

Anonymous User — Wed, 07 Jul 2010 19:30:30 GMT

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.

Re: Inside Buffer

DanPatterson_Retired — Wed, 07 Jul 2010 19:49:14 GMT

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

Re: Inside Buffer

Anonymous User — Thu, 08 Jul 2010 11:42:55 GMT

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 🙂

Re: Inside Buffer

ChrisSnyder — Thu, 08 Jul 2010 15:17:22 GMT

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.

Re: Inside Buffer

Anonymous User — Thu, 08 Jul 2010 16:18:38 GMT

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!

Re: Inside Buffer

ChrisSnyder — Sat, 11 Dec 2021 16:16:41 GMT

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()

Re: Inside Buffer

Anonymous User — Fri, 09 Jul 2010 13:15:36 GMT

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.

Re: Inside Buffer

LornaMurison — Fri, 09 Jul 2010 13:41:32 GMT

FYI
This problem is part of the larger problem I am trying to solve in this thread:
http://forums.arcgis.com/threads/7695-Challenge-Lake-vs.-River

I have had success creating a negative buffer from a simplified polygon by converting the polygon to a line, buffering the line, and erasing the line buffer from the simplified polygon.
Creating the complement of the lakes is a really interesting idea, but I have such a huge number of lakes to do this for that I don't think it is feasible.

As part of a process further down the line I need to re-buffer my buffers to get a more accurate reading of lake area and to help with creating pour points for catchment delineation. This second buffer process is giving me errors as well.

Re: Inside Buffer

Anonymous User — Fri, 09 Jul 2010 20:08:51 GMT

Original User: whuber

Creating the complement of the lakes is a really interesting idea, but I have such a huge number of lakes to do this for that I don't think it is feasible.

You can do it all in one step, Lorna. I wrote about focusing on just the problem lakes because they are problems! But if buffering them succeeds, then you should be able to buffer all the lakes by computing their mutual complement (one operation) and buffering it (one operation).

As part of a process further down the line I need to re-buffer my buffers to get a more accurate reading of lake area and to help with creating pour points for catchment delineation. This second buffer process is giving me errors as well.

Yes, it will give you problems, because all the "round" parts of buffers contain many, many vertices. That creates a huge complex feature to buffer, potentially causing difficulties (or at least taking a lot of time). Usually it's better not to re-buffer a buffer but to re-buffer the original feature by the difference in buffer radii.

Re: Inside Buffer

LornaMurison — Mon, 12 Jul 2010 17:38:55 GMT

You can do it all in one step, Lorna. I wrote about focusing on just the problem lakes because they are problems! But if buffering them succeeds, then you should be able to buffer all the lakes by computing their mutual complement (one operation) and buffering it (one operation).

You're right, I could do it in one step more or less. I'm satisfied, however, with converting polygons to lines, buffering the lines, and erasing the buffers.
This does a good job of getting rid of narrow areas, or splitting polygons where they narrow. As you can see in the attached picture. The red line represents the original lake outline, and the purple is the buffered lake (26 meter buffer).

The reason why I was wanting to re-buffer the negative buffer is to get a more accurate reading of the lake area. By re-buffering the lakes, the rivers would still be gone, except the area calculations would be more accurate.

I am now hoping that I can just add a field to my negatively buffered lakes with an estimation of lake area. But that is another question entirely. Is there a way to estimate the area of a buffer polygon if you know the buffer distance and the area of the original polygon?

Re: Inside Buffer

Anonymous User — Tue, 13 Jul 2010 12:06:10 GMT

Original User: whuber

Is there a way to estimate the area of a buffer polygon if you know the buffer distance and the area of the original polygon?

No, but if you're willing to do the (easy) calculation of perimeter, then there is a way: the additional buffer area is approximately equal to the buffer radius times the perimeter. The approximation is good when the radius is small compared to the polygon's size (as measured by the square root of its area, say) and the polygon's boundary is not too tortuous. In particular, it should work well for situations like the one you illustrated. (If you add to the approximation the area of a circle whose radius equals the buffer radius, this approximation often gives exact answers.)

For example: suppose the polygon area is 10,000 square meters, its perimeter is 500 meters, and the buffer radius is 1 meter. Then the (externally) buffered polygon area is estimated as

[INDENT]10,000 + 500 * 1 [+ Pi * 1^2]
= 10,500 [or 10,503.14] square meters.

[/INDENT]We expect this to work well in the example because (a) sqrt(10,000) : 1 = 100:1 shows the polygon is much larger than the buffer radius and (b) 500 : sqrt(10,000) = 5:1, a number that measures the tortuosity, is small.

Re: Inside Buffer

LornaMurison — Tue, 13 Jul 2010 12:44:12 GMT

For example: suppose the polygon area is 10,000 square meters, its perimeter is 500 meters, and the buffer radius is 1 meter. Then the (externally) buffered polygon area is estimated as

[INDENT]10,000 + 500 * 1 [+ Pi * 1^2]
= 10,500 [or 10,503.14] square meters.
[/INDENT]

Thank you so much, this is exactly what I am looking for, it's going to be a huge help!
Could you please clarify the formula though?
I am not sure what the buffer radius is being multiplied by and what pi is being added to. Is it 1(10,000 + 500)+(pi*1^2) ?
Thanks

*edit
I have tested out the formula and now understand it to be BufferRadius(Perimeter)+Area+(BufferRadius^2*pi)
It works extremely well, thank you.