Exclude end points in polygon geometrie?

283
8
Jump to solution
10-31-2019 05:49 AM
JohannesBierer
Regular Contributor

Hello,

I try to find excactly identical points in a polygon feature class with the 10.4. basic license to search for polygon errors with python 2.7. The following script can't find end points in the row[1].lastPoint in line 18? I need to exclude start and end points. What's wrong or any ideas how to do better?

import arcpy

infc = r"R:\XXX\LSG_Fehler.shp"
outfc = r"C:\temp2\GleicheKoord.shp"
CoordList = list()

arcpy.env.overwriteOutput = True

# Enter for loop for each feature
#
for row in arcpy.da.SearchCursor(infc, ["OID@", "SHAPE@"]):
    # Print the current multipoint's ID
    #
    print("Feature {}:".format(row[0]))
    partnum = 0

    
    if row[1] == row[1].lastPoint:

        print row[1].lastPoint
        break

    # Step through each part of the feature
    
    else:
        for part in row[1]:
            # Print the part number
            #
            print("Part {}:".format(partnum))

            # Step through each vertex in the feature
            
                
            for pnt in part:

                if pnt:
                                   
                    # print("{}, {}".format(pnt.X, pnt.Y))
                    CoordList.append("{};{}".format(pnt.X, pnt.Y))
                else:
                    # If pnt is None, this represents an interior ring
                    #
                    print("Interior Ring:")
                    
            partnum += 1

# Find identical coordinates
output = set([x for x in CoordList if CoordList.count(x) > 1])
print(output)

spatial_reference = arcpy.SpatialReference(25832)

pt = arcpy.Point()
ptGeoms = []
for p in output:

    coord = p.split(";")
    
    pt.X = coord[0]
    print(coord[0])
    pt.Y = coord[1]
    print(coord[1])
    ptGeoms.append(arcpy.PointGeometry(pt, spatial_reference))

arcpy.CopyFeatures_management(ptGeoms, outfc)
arcpy.AddXY_management(outfc)‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍
0 Kudos
1 Solution

Accepted Solutions
DanPatterson_Retired
MVP Esteemed Contributor

Here is a clean multipart polygon shape.  Each part has an interior ring and the parts are nested.  The logic of line 18 seems to suggest that you might be looking to see where a segment of the polygon boundary crosses itself

# ---- s0, a polygon object

s0.__class__
<class 'arcpy.arcobjects.geometries.Polygon'

s0.pointCount
19

s0.partCount
2

s0.firstPoint
<Point (300010.0, 5000010.0, #, #)>

s0.lastPoint
<Point (300006.0, 5000007.0, #, #)>

s0.JSON
'{"rings":[[[300010,5000010],[300010,5000000],[300000,5000000],[300000,5000010],[300010,5000010]],[[300003,5000009],[300003,5000003],[300009,5000003],[300009,5000009],[300003,5000009]],[[300008,5000008],[300008,5000004],[300004,5000004],[300004,5000008],[300008,5000008]],[[300006,5000007],[300005,5000005],[300007,5000005],[300006,5000007]]],"spatialReference":{"wkid":2146,"latestWkid":2951}}'

s0.__geo_interface__['coordinates']
[
 [[(300010.0, 5000010.0), (300010.0, 5000000.0), (300000.0, 5000000.0),
   (300000.0, 5000010.0), (300010.0, 5000010.0)],
  [(300003.0, 5000009.0), (300003.0, 5000003.0), (300009.0, 5000003.0),
   (300009.0, 5000009.0), (300003.0, 5000009.0)]],

 [[(300008.0, 5000008.0), (300008.0, 5000004.0), (300004.0, 5000004.0),
   (300004.0, 5000008.0), (300008.0, 5000008.0)], [(300006.0, 5000007.0),
   (300005.0, 5000005.0), (300007.0, 5000005.0), (300006.0, 5000007.0)]]
 ]

You might want to check line 18 since comparing the whole shape to the last point will always fail

Inspecting the parts and reassembling might help

arr = s0.getPart(0)

arr
<Array [<Point (300010.0, 5000010.0, #, #)>, <Point (300010.0, 5000000.0, #, #)>,
        <Point (300000.0, 5000000.0, #, #)>, <Point (300000.0, 5000010.0, #, #)>,
        <Point (300010.0, 5000010.0, #, #)>,
        None,
        <Point (300003.0, 5000009.0, #, #)>, <Point (300003.0, 5000003.0, #, #)>,
        <Point (300009.0, 5000003.0, #, #)>, <Point (300009.0, 5000009.0, #, #)>,
        <Point (300003.0, 5000009.0, #, #)>]>

[(p.X, p.Y) for p in arr if p]
[(300010.0, 5000010.0), (300010.0, 5000000.0), (300000.0, 5000000.0),
 (300000.0, 5000010.0), (300010.0, 5000010.0), (300003.0, 5000009.0),
 (300003.0, 5000003.0), (300009.0, 5000003.0), (300009.0, 5000009.0),
 (300003.0, 5000009.0)]

Notice that None has been removed and you could check for a crossing point if this were a 'bad' shape

View solution in original post

8 Replies
DanPatterson_Retired
MVP Esteemed Contributor

Here is a clean multipart polygon shape.  Each part has an interior ring and the parts are nested.  The logic of line 18 seems to suggest that you might be looking to see where a segment of the polygon boundary crosses itself

# ---- s0, a polygon object

s0.__class__
<class 'arcpy.arcobjects.geometries.Polygon'

s0.pointCount
19

s0.partCount
2

s0.firstPoint
<Point (300010.0, 5000010.0, #, #)>

s0.lastPoint
<Point (300006.0, 5000007.0, #, #)>

s0.JSON
'{"rings":[[[300010,5000010],[300010,5000000],[300000,5000000],[300000,5000010],[300010,5000010]],[[300003,5000009],[300003,5000003],[300009,5000003],[300009,5000009],[300003,5000009]],[[300008,5000008],[300008,5000004],[300004,5000004],[300004,5000008],[300008,5000008]],[[300006,5000007],[300005,5000005],[300007,5000005],[300006,5000007]]],"spatialReference":{"wkid":2146,"latestWkid":2951}}'

s0.__geo_interface__['coordinates']
[
 [[(300010.0, 5000010.0), (300010.0, 5000000.0), (300000.0, 5000000.0),
   (300000.0, 5000010.0), (300010.0, 5000010.0)],
  [(300003.0, 5000009.0), (300003.0, 5000003.0), (300009.0, 5000003.0),
   (300009.0, 5000009.0), (300003.0, 5000009.0)]],

 [[(300008.0, 5000008.0), (300008.0, 5000004.0), (300004.0, 5000004.0),
   (300004.0, 5000008.0), (300008.0, 5000008.0)], [(300006.0, 5000007.0),
   (300005.0, 5000005.0), (300007.0, 5000005.0), (300006.0, 5000007.0)]]
 ]

You might want to check line 18 since comparing the whole shape to the last point will always fail

Inspecting the parts and reassembling might help

arr = s0.getPart(0)

arr
<Array [<Point (300010.0, 5000010.0, #, #)>, <Point (300010.0, 5000000.0, #, #)>,
        <Point (300000.0, 5000000.0, #, #)>, <Point (300000.0, 5000010.0, #, #)>,
        <Point (300010.0, 5000010.0, #, #)>,
        None,
        <Point (300003.0, 5000009.0, #, #)>, <Point (300003.0, 5000003.0, #, #)>,
        <Point (300009.0, 5000003.0, #, #)>, <Point (300009.0, 5000009.0, #, #)>,
        <Point (300003.0, 5000009.0, #, #)>]>

[(p.X, p.Y) for p in arr if p]
[(300010.0, 5000010.0), (300010.0, 5000000.0), (300000.0, 5000000.0),
 (300000.0, 5000010.0), (300010.0, 5000010.0), (300003.0, 5000009.0),
 (300003.0, 5000003.0), (300009.0, 5000003.0), (300009.0, 5000009.0),
 (300003.0, 5000009.0)]

Notice that None has been removed and you could check for a crossing point if this were a 'bad' shape

JohannesBierer
Regular Contributor

How do you get from a shapefile (string) to a polygon object?

Yes, I'm looking for the points where a polygon boundary crosses itself. The idea was to compare the coordinates because in this points the coordinates of two points must be identical.

0 Kudos
DanPatterson_Retired
MVP Esteemed Contributor

Johannes

Checking and repairing geometries—Data Management toolbox | ArcGIS Desktop 

As part of the process, why not try a FeatureClassToFeatureClass first, apparently it repairs geometry issues during the load process.

Check Geometry—Data Management toolbox | ArcGIS Desktop 

Are you sure that the self intersections are where they are?

Why not provide the coordinates of a 'bad' shape so people can have a look at them. The self-intersection is not going to be a point, but two lines that cross

self_crossing = [[0,0], [1,10], [0,10], [10,10], [10,10], [0,0]]
0 Kudos
JohannesBierer
Regular Contributor

I have tried checking and repairing geometries - and topology must not overlap - with no result. For ArcGIS it's not a mistake? The script I postet found the point because of double coordinates but as well every starting and end point of the shape? It's a very small area you nearly can't see if you have a look at the shape.

0 Kudos
DanPatterson_Retired
MVP Esteemed Contributor

Don't know about the topology, but you keep saying shapefile.  FeatureClassToFeatureclass.  Use it to get the shapefile into a gdb first, that will remove geometry errors for an individual shape.  Overlapping polygons isn't a geometry error, it is a topology error, a subtle but important distinction.

JohannesBierer
Regular Contributor

Doesn't help to load to a fgdb. I have attached two shapes (PolygonError.zip) in the first question with the polygon shape and the error location.

0 Kudos
JohannesBierer
Regular Contributor

There's a second mistake in the shape of the same type?

0 Kudos
JohannesBierer
Regular Contributor

arr = s0.getPart(0) was the right hint to a solution, thank you

0 Kudos