Last active
August 31, 2023 21:48
-
-
Save kalombos/12ac7a25034ef9a3188af584bac22819 to your computer and use it in GitHub Desktop.
This is a solution for polygon selection in google maps. Some polygons can across antimeridian which will be getting wrong by postgis when you make ST_Within query, for example. So, this function allows to split polygon to multipolygon, parts of wich are contained in both hemispheres if needs
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
# -*- coding: utf-8 -*- | |
from django.contrib.gis.geos import MultiPolygon, GEOSGeometry | |
from django.db import connection | |
def is_crossed_antimeridian(polygon): | |
crossed = False | |
coords = [] | |
for polygon_set in polygon: | |
for (lng, lat) in polygon_set: | |
coords.append(lng) | |
current_lng = coords.pop(0) | |
for coord in coords: | |
if abs(current_lng-coord) > 180: | |
return True | |
current_lng = coord | |
return crossed | |
def split_polygon_by_antimeridian(polygon): | |
""" | |
Function split polygon to Multipolygon by antimeridian | |
:param polygon: Polygon | |
:return: Multipolygon | |
""" | |
if not is_crossed_antimeridian(polygon): | |
return MultiPolygon(polygon, srid=polygon.srid) | |
srid = polygon.srid, | |
wkt = polygon.wkt | |
with connection.cursor() as cursor: | |
cursor.execute(""" | |
SELECT ST_Union( | |
ST_Multi( | |
ST_Intersection( | |
ST_MakeEnvelope(0, -90, 180, 90, %s), | |
ST_MakeValid(ST_ShiftLongitude(ST_GeomFromText(%s, %s))) | |
) | |
) | |
, | |
ST_Translate( | |
ST_Multi( | |
ST_Intersection( | |
ST_MakeEnvelope(180, -90, 360, 90, %s), | |
ST_MakeValid(ST_ShiftLongitude(ST_GeomFromText(%s, %s))) | |
) | |
), -360, 0 | |
) | |
) | |
""", (srid, wkt, srid, srid, wkt, srid)) | |
row = cursor.fetchone() | |
return GEOSGeometry(row[0]) |
Thanks for you feedback it is nice to know that your work has helped somebody.
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Thanks for this. It mostly worked for me except for the whether it crosses the prime meridian. I think your function will work if it crosses the prime meridian on a short segment, but not if it wraps around much of the world.
Here's what I ended up with:
My problem space was slightly different since I knew I was always moving left to right and only dealing with a line segment not a polygon.
Anyway, thanks for posting this. Saved me a bunch of time.