compressing polygons for geodjango

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
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

More like this

  1. Enumeration field by nail.xx 5 years, 9 months ago
  2. Simple random file CAPTCHA by jeverling 2 years, 1 month ago
  3. Testing with unmanaged models by manfre 4 years, 10 months ago
  4. Variable._resolve_lookup monkeypatch by showell 4 years, 5 months ago
  5. Increase maximum number of changelist items for "Show all" link to appear by ramen 4 years, 5 months ago

Comments

(Forgotten your password?)