from django.core.management.base import BaseCommand
from geo.models import CountyBorders
from django.contrib.gis.geos.collections import MultiPolygon
from django.contrib.gis.geos.polygon import Polygon
from django.contrib.gis.geos.linestring import LinearRing
import time
class Command(BaseCommand):
option_list = BaseCommand.option_list
help = 'Simplify polygons'
args = ''
def handle(self, *app_labels, **options):
simplify()
print 'DONE!'
def simplify():
precision = 75
klasses = [CountyBorders] # change this or parameterize it
max = 2000 # for whole batch make this None
for klass in klasses:
objs = klass.objects.all()
if max:
objs = objs[:max]
for obj in objs:
fix_region(obj, precision)
def fix_region(obj, precision):
new_mpoly = fix_poly(obj.mpoly, precision, obj.id, obj.name)
obj.mpoly = new_mpoly
obj.save()
def fix_poly(mpoly, precision, id, name):
new_polies = []
for poly in mpoly:
new_rings = []
for ring in poly:
points = simplify_points(list(ring), precision, id, name)
new_rings.append(LinearRing(points))
new_polies.append(Polygon(*new_rings))
# new_polies.append(poly) # so we can see diff
new_mpoly = MultiPolygon(*new_polies)
return new_mpoly
def benefit_function(midpt, pt1, pt2):
# This corresponds to the area of the triangle
# that a point adds to a polygon. It favors
# either midpoints that are between two neighbors
# that are far away from each other or midpoints
# that are orthogonally far the neighbors' connecting
# segment.
xm, ym = midpt
x1, y1 = pt1
x2, y2 = pt2
return abs((x2 -x1)* (y1 - ym) - (x1-xm) * (y2 -y1))
def simplify_points(points, max, id, name):
# 338 points takes 0.012 secs to compress
tups = []
t = time.time()
for i, point in enumerate(points):
if i == 0 or i == len(points) - 1:
benefit = None
else:
benefit = benefit_function(point, points[i-1], points[i+1])
tups.append((point, benefit))
removals = 0
while len(tups) > max + 1:
min_benefit = tups[1][1]
min_i = 1
for i, (point, benefit) in enumerate(tups):
# print i, benefit
if benefit and benefit < min_benefit:
min_benefit = benefit
min_i = i
del tups[min_i]
removals += 1
for i in [min_i-1, min_i]:
point = tups[i][0]
if i == 0 or i == len(tups) - 1:
benefit = None
else:
benefit = benefit_function(point, tups[i-1][0], tups[i+1][0])
tups[i] = (point, benefit)
points = [point for point, benefit in tups]
if removals > 600:
print id, name, removals, 'removals took', time.time() - t
return points
Comments