Skip to content

Instantly share code, notes, and snippets.

@kalombos
Last active August 31, 2023 21:48
Show Gist options
  • Save kalombos/12ac7a25034ef9a3188af584bac22819 to your computer and use it in GitHub Desktop.
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
# -*- 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])
@kalombos
Copy link
Author

kalombos commented May 6, 2020

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