- June 9, 2009
- geodjango polygons
- 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 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] if i == 0 or i == len(tups) - 1: benefit = None else: benefit = benefit_function(point, tups[i-1], tups[i+1]) 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