8

This geometry is invalid:

POLYGON ((234107.4999285005 581107.0000356995, 234250.1399285011 581055.220035702, 234347.470028501 581007.8300357014, 234427.899928499 580959.1000357009, 234483.1350285001 580921.6430357024, 234483.9310285002 580921.1030356996, 234680.0683284998 580765.4714356996, 234603.1971285008 580681.0272357017, 234581.3791285008 580657.0597357005, 234555.6352284998 580633.0777356997, 234534.345028501 580615.3550356999, 234506.4820284992 580596.5460356995, 234484.0913285017 580583.7058357, 234458.9995285012 580571.5055356994, 234357.9406284988 580530.2466357015, 234307.3619285002 580510.5507357009, 234294.4100284986 580545.3700357005, 234280.2700285017 580597.4400357008, 234188.8799284995 580868.2600357011, 234163.6399285011 580943.5600357018, 234163.2110285014 580944.8316357024, 234144.8799284995 580999.1800356992, 234110.3599284999 581098.8800357021, 234107.5000285022 581107.0000356995, 234107.4999285005 581107.0000356995))

from shapely.validation import explain_validity
explain_validity(geometry) # 'Self-intersection[234107.500043162 581106.999994076]'

I tried the following but cannot find the solution:

First: Create a multipoint from the polygon and find the intersection. The multipoint is not valid due to self-intersection but except for the first and last point there are no intersections according to this check: 

 multipoint = shapely.geometry.MultiPoint(list(geometrie.exterior.coords))
    s = 0
    for s in range(len(multipoint)):
        for n, point in enumerate(multipoint):
             if n != s:
                   if point.intersects(multipoint[s]):
                        print('{}: {}'.format(n, point.intersects(multipoint[s])))
             s += 1

Although there seem to be 2 points that intersect visually: enter image description here

Second: I applied a buffer, this makes the polygon or multipoint valid. 

    geometry = geometry.buffer(0.0)
    geometry.is_valid #True

However, the resulting polygon is almost None:

geometry.buffer(0.01).area # 0.0004149480225530316

Unfortunately I do not understand why. I also do not understand the x,y coordinates that my geometries look like, it is hard to find the right search terms for this. Can you please help me?

PS the same neighborhood with slightly different coordinates from an alternative shape file (second in the picture, green is valid):

POLYGON ((234107.5 581107, 234250.1400146484 581055.2200317383, 234347.4700317383 581007.8300170898, 234427.9000244141 580959.1000366211, 234483.1350097656 580921.6430053711, 234483.9310302734 580921.1030273438, 234680.0684204102 580765.4716186523, 234603.1972045898 580681.0272216797, 234581.3792114258 580657.0598144531, 234555.6351928711 580633.0778198242, 234534.3450317383 580615.3550415039, 234506.4819946289 580596.5462036133, 234484.0914306641 580583.7059936523, 234458.9996337891 580571.5056152344, 234357.9407958984 580530.2465820313, 234307.3619995117 580510.5508422852, 234294.4100341797 580545.3699951172, 234280.2700195313 580597.4400024414, 234188.8800048828 580868.2600097656, 234163.6400146484 580943.5599975586, 234163.2109985352 580944.8317871094, 234144.8800048828 580999.1799926758, 234110.3599853516 581098.8800048828, 234107.5 581107))

enter image description here

PS2 Datasource: shapefile with neighbourhoods from the Netherlands crs:

{'proj': 'sterea',
 'lat_0': 52.15616055555555,
 'lon_0': 5.38763888888889,
 'k': 0.9999079,
 'x_0': 155000,
 'y_0': 463000,
 'ellps': 'bessel',
 'units': 'm',
 'no_defs': True}
Improve this question
4
  • 1
    I believe the area is given in degrees, not meters. If you want this value in square meters you should reproject this feature to any metric reference system.
    – dmh126
    May 6, 2019 at 8:11
  • Thanks! This pointed me to "if you work with a crs.unit = degree all calculations are wrong" gis.stackexchange.com/questions/166675/… -> I will change the crs to an appropriate format in meters, hopefully that solves my issue!
    – Wouter
    May 6, 2019 at 9:13
  • 1
    No, that is wrong and totally unrelated. Shapely operates on coordinates in any format given, it only looks at the numbers. In this case your perfectly fine local, apparently metric coordinates. Areas in shapely are always in whatever unit your coordinates are in.
    – bugmenot123
    May 6, 2019 at 9:15
  • Thanks, indeed, I have now changed the crs to meters and the issues remain. And indeed I cannot find anything in the Shapely docs that says that the unit of measure has to be in meters... Any other ideas please?
    – Wouter
    May 6, 2019 at 9:28

1 Answer 1

Reset to default
11

Unfortunately your issue is not just a duplicate coordinate. If you zoom in to the self-intersection, you can see a tiny extra area (this is probably the result of your buffer(0) btw):

qgis validity check

Now this is ambitious. Is there one point too many (the one on the left)? Or are the lines connected wrong here?

Look at surrounding data

I would recommend looking at the surrounding data to see if you can infer the correct topology and then manually fix this.

Manual guessing

On the other hand, this is a really really tiny corner of the polygon so maybe just remove the very last defining coordinate (234107.5000285022 581107.0000356995) and be happy.

Exploit the apparent convexity of the polygon

Or, if you feel safe and want an automated approach, you can calculate the convex hull as this polygon looks like it should be convex.

To get a MultiPoint geometry you need to exclude either the first or last coordinate. They are repeated in the LinearRings that defines the rings of a Polygon.

points = MultiPoint(polygon.exterior.coords[1:])
convex_hull = points.convex_hull
convex_hull.is_valid  # -> True
convex_hull.is_simple  # -> True

pistachio is the "fixed" polygon, mortadella is the broken one:

enter image description here

Improve this answer
4
  • Many thanks! And sorry, but how did you zoom in on the self-intersection? Also, the geography (polygon) that I shared was original from the source, without buffer. And: I found a second geography in a different shape file for the same neighborhood which is valid. This solves my issue and I really want to learn more about this so have added it to my question hoping it will inspire to better understand the issue with the first polygon. Notice the different color in my Jupyter notebook.
    – Wouter
    May 6, 2019 at 10:16
  • 1
    I used QGIS to look at it.
    – bugmenot123
    May 6, 2019 at 10:55
  • 7
    +1 for the detailed explanation with pictures and for showing color mortadella to us.
    – jjmontes
    Mar 24, 2020 at 2:11
  • Further to this, is there a way of moving the lines (let's say 5 meters) left and right of the intersection point to remove this intersection point? Thanks
    – user1655130
    Apr 15, 2021 at 13:30

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.