Login

compressing polygons for geodjango

Author:
showell
Posted:
June 9, 2009
Language:
Python
Version:
1.0
Tags:
geodjango polygons
Score:
0 (after 0 ratings)

The code shown allows you, in GeoDjango, to reduce the number of points in your polygons. It helps reduce storage needs and makes queries run faster, at the cost of some precision. It provides a variation on the simplify() method that comes with the GEOS API, allowing you to specify a number of points instead of a distance tolerance.

It is set up as a management command so that you can run it with python manage.py. See the django docs for how to set that up.

The example shown assumes a table called CountyBorders with fields "name" and "mpoly." It should be straightforward to adapt it for your needs. Look at the first three lines of the simplify() method in particular for customization. The rest of the code is pretty generic, and it should run fast enough in a one-time batch process for most needs.

The algorithm tries to keep the 75 points that provide the most definition for the shape. Each point in a polygon defines a triangle with its immediate neighbors. If that triangle has no area (the degenerate case), it is a midpoint on the segment between its neighbors and adds no value whatsoever. This principle is extended to say that the larger the triangle, the more value the point has in defining the shape. (You can find more refined algorithms, but this code seems to work fine by visual inspection.)

 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 6 years, 10 months ago
  2. Simple random file CAPTCHA by jeverling 3 years, 2 months ago
  3. Testing with unmanaged models by manfre 5 years, 11 months ago
  4. Variable._resolve_lookup monkeypatch by showell 5 years, 6 months ago
  5. Increase maximum number of changelist items for "Show all" link to appear by ramen 5 years, 6 months ago

Comments

Please login first before commenting.