6

I've made a model in ModelBuilder in ArcMap 10.3, what it basically does is creating polygons over open areas such as farm-lands, grass-lands, marshes etc that has a low slope.

My challenge is that I don't want polygons that are smaller than a specific area. lets say in this example that the area can not be smaller than 160m2. This is easy by creating an expression removing polygons with area less than that. But I'm interested in the actual size, it needs to be at least 40x40m. (diameter of 40m if it is doable with circle vs square).

Examples

Figure 1 and 2 are OK because they both have an area inside which is 40x40m, but figure 3 does not and should thus be removed/deleted. It also visualizes the trouble with square meters, as the figure 3 has the same area as figure 1.

As long as the polygon have one area of at least 40x40 the remainder "noise"/bulbs in the polygon does not matter.


@FelixIP gave one solution in the first answer on this post that almost worked. But after some testing I found that there were some flaws in the output after running the two scripts. It solved my problem since, luckily enough, no errors occurred in my analysis extent. Here are two pictures that depicts the two errors i found:

error1

error2

As you can see here, if we had moved the point from which the buffer was made, it would be possible to arrange it so that the buffer would fit inside the polygon boundaries.

  • To clarify, the x and y dimensions must both be >= 40m? – Paul Oct 12 '15 at 18:21
  • Does 20m-radius circle suffice? If yes see gis.stackexchange.com/questions/147790/… or you really talking of fitting 40*40 square into polygon? – FelixIP Oct 12 '15 at 19:02
9

Use my script from Checking if polygon fits inside another polygon using ArcGIS or QGIS?

to create points inside your polygons. Each point is a centre of maximum inscribed circle of the parcel. Please note, script tested on shapefiles only and it does not handle polygons with holes, i.e. donut like polygons.

Table of centre point contains field theDist, that is the radius of above circle. enter image description here

Select all points with theDist=>20, and create 20m buffers around them: enter image description here Use Minimum Bounding Geometry to create rectangles from buffers enter image description here Attach this script to tool with no parameters. It is a tool to run from mxd

import arcpy, traceback, os, sys,math
from math import radians,sin,cos

mxd = arcpy.mapping.MapDocument("CURRENT")
parent=arcpy.mapping.ListLayers(mxd,"parent")[0]
child=arcpy.mapping.ListLayers(mxd,"child")[0]

# NUMERIC ! common field
linkField="PAR_ID"
d=arcpy.Describe(parent)
SR=d.spatialReference
g=arcpy.Geometry()
gr=(math.sqrt(5)-1)/2

try:
    def showPyMessage():
        arcpy.AddMessage(str(time.ctime()) + " - " + message)

    #golden section to find minimum
    def gss(a,b,tol):
        c=b-gr*(b-a)
        d=a+gr*(b-a)
        while abs(c-d)>tol:       
            fc=f(c);fd=f(d)
            if fc<fd:
                b=d
                d=c
                c=b-gr*(b-a)
            else:
                a=c
                c=d
                d=a+gr*(b-a)
        return (b+a)/2

    # rotate polygon
    def ShapeMake(pGon,angle):
        ar=arcpy.Array()
        a=radians(angle)
        part=pGon.getPart(0)
        for p in part:
            x,y=p.X-pGon.centroid.X,p.Y-pGon.centroid.Y
            xN=cos(a)*x+sin(a)*y
            yN=-sin(a)*x+cos(a)*y
            pN=arcpy.Point(xN+pGon.centroid.X,yN+pGon.centroid.Y)
            ar.add(pN)
        pgonRotated=arcpy.Polygon(ar,SR)
        return pgonRotated

    #function to minimise
    def f(a):
        pgonRot=ShapeMake(square,a)
        intR=pgonRot.difference(bigPolygon)
        return intR.area

    result=arcpy.GetCount_management(child)
    nF=int(result.getOutput(0))
    initFidList=[]
    arcpy.SetProgressor("step", "", 0, nF,1)
    with arcpy.da.UpdateCursor(child, ("SHAPE@",linkField)) as rows:
        for row in rows:
            square=row[0]
            PID=row[1]
            quer='%s%s%s=%i'%('"',linkField,'"',PID)
            parent.definitionQuery=quer
            bigPolygon=arcpy.CopyFeatures_management(parent,g)[0]
            intR=square.difference(bigPolygon)

            # square fits inside parent without rotation
            if intR.area==0:initFidList.append(PID)
            # find minimum area outside parent at differenr rotations
            else:
                angle=gss(0.0,90,1e-3)
                # square fits if rotated
                if f(angle)==0:
                    row[0]=ShapeMake(square,angle)
                    rows.updateRow(row)
                    initFidList.append(PID)
            arcpy.SetProgressorPosition()

    quer='"FID" IN '+str(initFidList)
    quer='%s%s%s in %s'%('"',linkField,'"',str(tuple(initFidList)))
    parent.definitionQuery=""
    arcpy.SelectLayerByAttribute_management(child, "NEW_SELECTION", quer)


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

Script works on shapefiles and assumes:

  1. that your original polygons named PARENT
  2. that your squares polygons named CHILD
  3. They have common numeric field "PAR_ID"

It child fits completely into original, script adds this square to selection. If not it attempts to fit it by rotating around centre points, and adds to selection if successful

enter image description here

  • Ok I have to ask... what's the deal with the golden ratio? – nickves Oct 13 '15 at 20:04
  • 2
    It is one of many techniques to find extremum of unimodal function. en.wikipedia.org/wiki/Golden_section_search In this case we are trying to find rotation at which geometric difference between square and original polygon is minimum. – FelixIP Oct 13 '15 at 20:18
  • Mr FelixIP, always one step ahead! :-) – Hornbydd Oct 13 '15 at 23:50
  • I think this answer is great, but to offer a simpler, but far less robust alternative that might work in this case: if the polygons are oriented vertically in the projection (as in the example image) where sides either run parallel or perpendicular to the horizontal (if that makes sense), you could use cursor and access the extent properties of the Shape attribute, and then just do some basic subtraction to determine if the side lengths are long enough. It doesn't handle rotation like this answer does though. – nicksan Oct 14 '15 at 1:07
  • @Hornbydd this was easy to assemble, I've made it out of 3 ready elements (centre, golden section, polygon rotation). This one gis.stackexchange.com/questions/166326/… IS hard case though. Obviously your area - river cross-sections, although I am hydrologist by trade also – FelixIP Oct 14 '15 at 1:19
3

Interesting question!

OK thinking out aloud here and untested so may not capture all possible shapes. It looks like you are mapping to a cell size of 10m and that your rectangular polygons are raster based. The logic I am suggesting only works if the original raster based polygons do not touch each other. In @FelixIP answer he has polygons sharing boundaries. If this is truly the case how about this approach?

  1. Convert your polygons to a binary grid where a cell covered by your polygon becomes 1, everything else a 0.
  2. Run the Focal Statistics tool with a circular window of 2 pixels with SUM as the option.
  3. Any pixel surrounded by 4x4 1's will sum up to 16.
  4. You could then convert only pixels that are 16 into points and use those to select the polygons as the required minimum dimension. You just need to invert the selection to find the ones you don't want.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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