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])
@hirowatari
Copy link

hirowatari commented May 5, 2020

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:

    def crosses_prime_meridian?(min, max)
      min = min + 360 if min < 0
      max = max + 360 if max < 0
      min < 180 && max >180 ||
        min < 180 && max < min ||
        min > 180 && max > 180 && min > max
    end

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.

@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