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)
Solved! Go to Solution.
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
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
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.
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]]
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.
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.
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.
There's a second mistake in the shape of the same type?
arr = s0.getPart(0) was the right hint to a solution, thank you