Login

compressing polygons for geodjango

Author:
showell
Posted:
June 9, 2009
Language:
Python
Version:
1.0
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. Template tag - list punctuation for a list of items by shapiromatron 10 months, 1 week ago
  2. JSONRequestMiddleware adds a .json() method to your HttpRequests by cdcarter 10 months, 2 weeks ago
  3. Serializer factory with Django Rest Framework by julio 1 year, 5 months ago
  4. Image compression before saving the new model / work with JPG, PNG by Schleidens 1 year, 6 months ago
  5. Help text hyperlinks by sa2812 1 year, 6 months ago

Comments

Please login first before commenting.