Select to view content in your preferred language

# Exclude end points in polygon geometrie?

1172
8
10-31-2019 05:49 AM
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)
1 Solution

Accepted Solutions
MVP Emeritus

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

8 Replies
MVP Emeritus

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

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.

MVP Emeritus

Johannes

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

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]]‍``
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.

MVP Emeritus

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.

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.

Regular Contributor

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

Regular Contributor

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