Tuesday, 23 February 2016

geodjango - Certain MultiPolygons cause “BOOM! Could not generate outside point!” in PostGIS?


I've just posted this at StackOverflow and someone suggested I try here. Original: https://stackoverflow.com/q/5863859/246265



I'm trying to represent a rectangular area which crosses 180 degrees longitude. For more background see https://stackoverflow.com/q/5737506/246265


Here's my test case:


from django.contrib.gis.geos import Polygon, MultiPolygon
from my_project.my_app.models import Photo

a = Polygon.from_bbox((30, -80, 180, 80)) # the part to the east of 180
b = Polygon.from_bbox((-180, -80, -160, 80)) # a part to the west of 180
c = Polygon.from_bbox((-180, -80, -100, 80)) # a larger part to the west of 180

ok = MultiPolygon(a,b)

ok2 = MultiPolygon(c)
boom = MultiPolygon(a,c)

# This works
Photo.objects.filter(location__coveredby=ok)[:1]
# This also works so c is ok
Photo.objects.filter(location__coveredby=ok2)[:1]
# This gives "BOOM! Could not generate outside point!"
Photo.objects.filter(location__coveredby=boom)[:1]


# splitting c doesn't help
c1 = Polygon.from_bbox((-180, -80, -140, 80))
c2 = Polygon.from_bbox((-140, -80, -100, 80))
test = MultiPolygon(a,c1,c2)
Photo.objects.filter(location__coveredby=test)[:1]
# BOOM! Could not generate outside point!

By changing the numbers I can make this error come and go. (-180, -80, x, 80) works where x <= -140 for example. For every number there is a threshold like this but I can't find a pattern. For boxes with the same area, some will work and others not. For boxes with the same width some will work and others not.


I can look at the SQL being generated but the areas are represented in binary (EWKB) and I'm not sure how to read it.


Can anyone explain this?




Answer



I think I've figured this out, please correct me if my assumptions are wrong.


The MultiPolygon(a,c) covers more than half the earth's surface so postGIS thinks I'm referring to the area outside my MultiPolygon, because it's smaller. This area includes the poles, as Mike Toews pointed out, postGIS doesn't like these extreme values.


My solution is to use an SQL OR to join each Polygon, rather than putting them in a MultiPolygon. My area is split up in the first place to allow for searching areas bigger than half the earth's surface.


No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...