1

I have a system that allows users to input data in various geometries, with which I do some processing using JTS and then store it in a MySQL database.

A few particular polygons that my users have submitted are registering as invalid by both MySQL and JTS. I have tried to open this polygon in various viewers and they display it without issue. I also glanced over the data myself and couldn't see anything wrong.

I am familiar with how a polygon can be invalid in general. What I can't figure out is what is invalid about my particular polygon. Any idea? Here is one of them:

POLYGON((-122.148755985 37.42723718100007,-122.148462038 37.42697222300006,-122.1476020649999 37.42622237100005,-122.14681093 37.42551726300007,-122.146799162 37.42550704500007,-122.1467579119999 37.42547357700005,-122.146714522 37.42544187700003,-122.146669111 37.42541203100006,-122.1466218039999 37.42538412200003,-122.146554117 37.42534915200002,-122.145898905 37.42500666200004,-122.146606488 37.42414940300006,-122.14718414 37.42346177100006,-122.147744212 37.42278424400007,-122.147784198 37.42273583200006,-122.1480067479999 37.42246395100005,-122.148371184 37.42202471400003,-122.148569196 37.42178912700007,-122.14798232 37.42249379400005,-122.1480067479999 37.42246395100005,-122.14949105 37.42068471600004,-122.149804428 37.42030734400004,-122.150144279 37.41989809300003,-122.150733848 37.41919052400004,-122.151078419 37.41877586800007,-122.151717718 37.41800663700008,-122.151912941 37.41777175300007,-122.152503005 37.41706315200003,-122.1531540729999 37.41628313600006,-122.153061958 37.41623509800007,-122.153290778 37.41595996300003,-122.153292466 37.41596085600003,-122.153309575 37.41594035100007,-122.157366324 37.41810165800007,-122.1496241789999 37.42741009900004,-122.149316868 37.42775509600006,-122.149258069 37.42770423800005,-122.149162514 37.42761216700006,-122.148755985 37.42723718100007))
Improve this question
2
  • 1
    Please convert it to a linestring and Edit the Question to include a picture, with the vertices numbered.
    – Vince
    Commented Nov 23 at 21:41
  • That's a good idea but I'm not quite sure how to do it. I tried loading the polygon into QGIS but I don't see an option to add numbering to the vertices. Any suggestion?
    – Dan
    Commented Nov 23 at 21:58

3 Answers 3

Reset to default
8

PostGIS 3.1.1, in PostgreSQL 13 (what I had handy), is able to report the issue --

WITH bad_data AS (
SELECT ST_GeomFromText('POLYGON((-122.148755985 37.42723718100007,-122.148462038 37.42697222300006,-122.1476020649999 37.42622237100005,-122.14681093 37.42551726300007,-122.146799162 37.42550704500007,-122.1467579119999 37.42547357700005,-122.146714522 37.42544187700003,-122.146669111 37.42541203100006,-122.1466218039999 37.42538412200003,-122.146554117 37.42534915200002,-122.145898905 37.42500666200004,-122.146606488 37.42414940300006,-122.14718414 37.42346177100006,-122.147744212 37.42278424400007,-122.147784198 37.42273583200006,-122.1480067479999 37.42246395100005,-122.148371184 37.42202471400003,-122.148569196 37.42178912700007,-122.14798232 37.42249379400005,-122.1480067479999 37.42246395100005,-122.14949105 37.42068471600004,-122.149804428 37.42030734400004,-122.150144279 37.41989809300003,-122.150733848 37.41919052400004,-122.151078419 37.41877586800007,-122.151717718 37.41800663700008,-122.151912941 37.41777175300007,-122.152503005 37.41706315200003,-122.1531540729999 37.41628313600006,-122.153061958 37.41623509800007,-122.153290778 37.41595996300003,-122.153292466 37.41596085600003,-122.153309575 37.41594035100007,-122.157366324 37.41810165800007,-122.1496241789999 37.42741009900004,-122.149316868 37.42775509600006,-122.149258069 37.42770423800005,-122.149162514 37.42761216700006,-122.148755985 37.42723718100007))',4326) as geom
)
SELECT ST_IsValidReason(geom)
FROM   bad_data;

Returns:

st_isvalidreason
Self-intersection[-122.148261898889 37.4221581012526]

Indeed, if you convert the feature to LINESTRING (by replacing the "POLYGON(" with "LINESTRING" and removing the trailing pattern), then to points (I used ArcGIS Pro and Feature Vertices to Points, but other solutions exist), label the points, and visually inspect them, the problem stands out:

sequenced vertices

Your coordinates double back upon themselves.

line overlap

The general solution for detecting this is calculating the angle at B across ABC for each triplet of sequential vertices, and looking for a very small (~0) value. Fixing it comes in two schools, dropping the later vertex, or the earlier one. You could also fold the removed vertices back into the line in an order that doesn't corrupt the shape. Evaluating which is most appropriate to your use case is left as an exercise.

Improve this answer
1
  • 1
    Thanks I've never actually used PostGIS before and didn't know it could do that. As mentioned my production system is MySQL but maybe I should spin up my own personal little PostGIS for situations like this.
    – Dan
    Commented Nov 23 at 23:12
6

If you would prefer an all Java solution then you should look at the JTS IsValidOp

So a program like:

String polyStr = "POLYGON((-122.148755985 37.42723718100007,-122.148462038 37.42697222300006,-122.1476020649999 37.42622237100005,-122.14681093 37.42551726300007,-122.146799162 37.42550704500007,-122.1467579119999 37.42547357700005,-122.146714522 37.42544187700003,-122.146669111 37.42541203100006,-122.1466218039999 37.42538412200003,-122.146554117 37.42534915200002,-122.145898905 37.42500666200004,-122.146606488 37.42414940300006,-122.14718414 37.42346177100006,-122.147744212 37.42278424400007,-122.147784198 37.42273583200006,-122.1480067479999 37.42246395100005,-122.148371184 37.42202471400003,-122.148569196 37.42178912700007,-122.14798232 37.42249379400005,-122.1480067479999 37.42246395100005,-122.14949105 37.42068471600004,-122.149804428 37.42030734400004,-122.150144279 37.41989809300003,-122.150733848 37.41919052400004,-122.151078419 37.41877586800007,-122.151717718 37.41800663700008,-122.151912941 37.41777175300007,-122.152503005 37.41706315200003,-122.1531540729999 37.41628313600006,-122.153061958 37.41623509800007,-122.153290778 37.41595996300003,-122.153292466 37.41596085600003,-122.153309575 37.41594035100007,-122.157366324 37.41810165800007,-122.1496241789999 37.42741009900004,-122.149316868 37.42775509600006,-122.149258069 37.42770423800005,-122.149162514 37.42761216700006,-122.148755985 37.42723718100007))";
WKTReader reader = new WKTReader();
Polygon poly = (Polygon) reader.read(polyStr);
IsValidOp valid = new IsValidOp(poly);
if(!valid.isValid()){
    System.out.println(valid.getValidationError());
}

which gives the result:

Ring Self-intersection at or near point (-122.1480067479999, 37.42246395100005, NaN)
Improve this answer
1
  • Ah right I forgot about getValidationError, in fact that gives me an idea about how I can better communicate these kinds of issues to my users.
    – Dan
    Commented Nov 24 at 17:29
2

Another option using Python:

from shapely.wkt import loads
from shapely.validation import explain_validity
from shapely.geometry import Point

wkt = r"POLYGON((-122.148755985 37.42723718100007,-122.148462038 37.42697222300006,-122.1476020649999 37.42622237100005,-122.14681093 37.42551726300007,-122.146799162 37.42550704500007,-122.1467579119999 37.42547357700005,-122.146714522 37.42544187700003,-122.146669111 37.42541203100006,-122.1466218039999 37.42538412200003,-122.146554117 37.42534915200002,-122.145898905 37.42500666200004,-122.146606488 37.42414940300006,-122.14718414 37.42346177100006,-122.147744212 37.42278424400007,-122.147784198 37.42273583200006,-122.1480067479999 37.42246395100005,-122.148371184 37.42202471400003,-122.148569196 37.42178912700007,-122.14798232 37.42249379400005,-122.1480067479999 37.42246395100005,-122.14949105 37.42068471600004,-122.149804428 37.42030734400004,-122.150144279 37.41989809300003,-122.150733848 37.41919052400004,-122.151078419 37.41877586800007,-122.151717718 37.41800663700008,-122.151912941 37.41777175300007,-122.152503005 37.41706315200003,-122.1531540729999 37.41628313600006,-122.153061958 37.41623509800007,-122.153290778 37.41595996300003,-122.153292466 37.41596085600003,-122.153309575 37.41594035100007,-122.157366324 37.41810165800007,-122.1496241789999 37.42741009900004,-122.149316868 37.42775509600006,-122.149258069 37.42770423800005,-122.149162514 37.42761216700006,-122.148755985 37.42723718100007))"

polygon = loads(wkt)

print(polygon.is_valid)
#False
print(explain_validity(polygon))
#'Self-intersection[-122.148261898889 37.4221581012526]'

#Plot the polygon and the point of the self intersection
import geopandas as gpd
ax = gpd.GeoSeries([polygon], crs=4326).plot(figsize=(10,10), color="yellow", alpha=0.1, zorder=0)
gpd.GeoSeries([polygon], crs=4326).boundary.plot(ax=ax, color="black", linewidth=0.5, zorder=1)
gpd.GeoSeries([Point(-122.148262, 37.422127)], crs=4326).plot(ax=ax, marker="+", color="red", 
                                                              markersize=400, zorder=2)

enter image description here

Improve this answer

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.