6

I am currently struggling with a mathematical problem and can't seem to get it solved.

I have a bunch of rotated rectangular polygons (solution should also work for non-rectangular polygons) in different shapefiles. I want to (automatically) find the lower left corner of each of the polygons using python (arcpy).

First I tried to simply use the extent, but since they are rotated, the extent will bring false results. So how do I approach this? I googled for hours and the only thing I could come up with, is to somehow calculate the center of the polygon or somehow temporairly rotate it to 180 degrees and then use xmin and ymin. Here are a few pictures which may help to distract the attention from my bad english (sorry for that btw)

enter image description here enter image description here

I had no problems exporting the coordinates of the polygons by using the feature vertices to points tool. So far I extracted them once in a list of tupels and once in two different lists (one for x-coords and one for y-coords)

so here are a few coordinates to make it easier to test:

tupel-format:

tupelList=[[1792398.680577231, 4782539.85121522],
 [1792173.0363913027, 4780368.293228334],
 [1788935.7990357098, 4780713.732859781],
 [1789162.9530321995, 4782885.332685629],
 [1792398.680577231, 4782539.85121522]]

x-list/y-list format:

xList= [1792398.680577231, 1792173.0363913027, 1788935.7990357098, 1789162.9530321995, 1792398.680577231]

yList=[4782539.85121522, 4780368.293228334, 4780713.732859781, 4782885.332685629, 4782539.85121522]

I am currently using ArcGIS 10.3.1 width ally extensions and an advanced license. Version of Python is 2.7.something

  • 3
    Hmmm, take the extent, then find the vertex that's closest to the southwest corner of the extent? What vertex do you want if the polygon is V shaped? – mkennedy Nov 11 '15 at 23:12
  • Have a look at this gis.stackexchange.com/questions/121632/… It uses convex hull and rotation angle. I guess it's a good starting point for you. It's tuple anyway... – FelixIP Nov 11 '15 at 23:33
  • Even though you talk primarily about rectangles I think you should edit your question to show examples of what you mean by "non-rectangular" as this will dictate the solution – Hornbydd Nov 11 '15 at 23:41
  • 1
    @mkennedy this won't work for some polygons like parallelograms--> see a picture here (couldn't post more than two links, because I'm new here) oi67.tinypic.com/33wrytw.jpg (the brownish stuff is supposed to be the extent) – fry82 Nov 11 '15 at 23:47
  • @Hornbydd The polygons could as well be parallelograms or even non-symetrial polygons with four corners. I don't understand what difference it would make for the solution though – fry82 Nov 11 '15 at 23:50
7

Here's a very simple approach that offloads all the processing into the Sort GP tool. Since you have access to an Advanced license, sorting by shape and starting at the lower left corner gives quick results.

import os, arcpy

arcpy.env.overwriteOutput = True

inFC = r'<path>'
outFC = r'<path>'

# create output FC to hold points and field to link OID
spatref = arcpy.Describe(inFC).spatialReference
arcpy.CreateFeatureclass_management(*os.path.split(outFC), geometry_type="POINT",
                                    spatial_reference=spatref)
arcpy.AddField_management(outFC, "ID", "LONG")


with arcpy.da.SearchCursor(inFC, ["OID@", "SHAPE@"]) as sCursor:
    with arcpy.da.InsertCursor(outFC, ["ID", "SHAPE@"]) as iCursor:
        for oid, poly in sCursor:

            # using Geometry objects is very quick and also has the added
            # benefit of returning lists of geometries
            # Since we are sorting the vertices by LL, the first one is the answer
            verts = arcpy.FeatureVerticesToPoints_management(poly, arcpy.Geometry())
            sort = arcpy.Sort_management(verts, arcpy.Geometry(),
                                         [["SHAPE", "ASCENDING"]], "LL")[0]
            iCursor.insertRow([oid, sort])

enter image description here

From the explanation on spatial sorting, we see that N/S takes precedence over E/W:

Note that U gets priority over R. R is taken into considerations only when some features are at the same horizontal level.

  • Cool use of arcpy.Geometry() in arcpy.Sort_management. – Emil Brundage Nov 12 '15 at 18:35
4

Given the examples of the rectangles and parallelograms, and if I'm understanding your formulation of "lower-left" (i.e., "the most south-west") correctly, a naive solution:

You could sort the four vertices in order of ascending latitude(i.e., the south-most is first, north-most is last). If two are equal in latitude, their sequence doesn't matter. Then, among the two south-most points, select the one that is further west. This would give you the green-dotted corners in this illustration:

enter image description here

The following might be rare or non-existent cases, but the above falls apart with a quadrilateral where the two vertices between the south-most and north-most are of equal latitude, such as a "diamond" shape, or some deformation of it, or even when a perfect rectangle is rotated just so:

enter image description here

To deal robustly with these cases, you could use an approach that calculates the medial axis of the shape (https://en.wikipedia.org/wiki/Medial_axis), compares it's slope to 45 and 135 degrees, and determines the "bottom" and "top" of the shape compared to that. Then progress as above, with the west-most point of the two "bottom" points being the point to choose. That would work in two of the above illustrated cases. It would still fail in the case of a perfect square in "diamond" rotation - but then there's no "natural" bottom-left corner to that shape. Another way to deal with these cases where points 2 and 3 in the vertical sequence are at equal latitude is to actually take point 1, the south-most point; or the west-most point among 2 and 3.

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.