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
- Template tag - list punctuation for a list of items by shapiromatron 11 months, 2 weeks ago
- JSONRequestMiddleware adds a .json() method to your HttpRequests by cdcarter 11 months, 2 weeks ago
- Serializer factory with Django Rest Framework by julio 1 year, 6 months ago
- Image compression before saving the new model / work with JPG, PNG by Schleidens 1 year, 7 months ago
- Help text hyperlinks by sa2812 1 year, 7 months ago
Comments
Please login first before commenting.