Permalink
Browse files

added initial quadtree implementation

  • Loading branch information...
1 parent 4458895 commit b99656454d738d1cb5ffa1fd507949b3dfab98b8 @trey0 trey0 committed May 5, 2012
View
@@ -4,8 +4,12 @@
*~
/example/dev.db
/example/example.db
+/example/database.db
/docs/_build/
/pip-log.txt
/static/js/*.r*.js
/static/css/*.r*.css
/geocamLayer/build
+/geocamLayer/static/geocamLayer/tiles
+/example/sourceme.sh
+/data
View
Binary file not shown.
View
@@ -5,13 +5,12 @@
# All Rights Reserved.
# __END_LICENSE__
-from django.core.management import execute_manager
-try:
- import settings # Assumed to be in the same directory.
-except ImportError:
- import sys
- sys.stderr.write("Error: Can't find the file 'settings.py' in the directory containing %r. It appears you've customized things.\nYou'll have to run django-admin.py, passing it your settings module.\n(If the file settings.py does indeed exist, it's causing an ImportError somehow.)\n" % __file__)
- sys.exit(1)
+import os
+import sys
if __name__ == "__main__":
- execute_manager(settings)
+ os.environ.setdefault("DJANGO_SETTINGS_MODULE", "example.settings")
+
+ from django.core.management import execute_from_command_line
+
+ execute_from_command_line(sys.argv)
View
@@ -9,8 +9,8 @@
DEBUG = True
TEMPLATE_DEBUG = DEBUG
import os, sys
-APP = os.path.abspath(os.path.dirname(os.path.dirname(__file__)))
-PROJ_ROOT = os.path.abspath(os.path.dirname(__file__))
+APP = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))
+PROJ_ROOT = os.path.dirname(os.path.abspath(__file__))
sys.path.append(APP)
ADMINS = (
@@ -22,7 +22,7 @@
DATABASES = {
'default': {
'ENGINE': 'django.db.backends.sqlite3',
- 'NAME': 'database.db',
+ 'NAME': os.path.join(PROJ_ROOT, 'database.db'),
'USER': '',
'PASSWORD': '',
'HOST': '',
View
@@ -4,4 +4,25 @@
# All Rights Reserved.
# __END_LICENSE__
-from django.contrib.gis import admin
+from django.contrib import admin
+
+from geocamLayer.models import Feature, QuadTreeCell
+
+class FeatureAdmin(admin.ModelAdmin):
+ list_display = ('lat',
+ 'lng',
+ 'name',
+ 'description',
+ 'cell')
+
+class QuadTreeCellAdmin(admin.ModelAdmin):
+ list_display = ('zoom',
+ 'x',
+ 'y',
+ 'count',
+ 'isLeaf',
+ 'lng',
+ 'lat')
+
+admin.site.register(Feature, FeatureAdmin)
+admin.site.register(QuadTreeCell, QuadTreeCellAdmin)
@@ -0,0 +1,62 @@
+#!/usr/bin/env python
+# __BEGIN_LICENSE__
+# Copyright (C) 2008-2010 United States Government as represented by
+# the Administrator of the National Aeronautics and Space Administration.
+# All Rights Reserved.
+# __END_LICENSE__
+
+import sys
+import itertools
+
+from geocamLayer.noaaExampleParser import readNoaaWeatherStations
+from geocamLayer.quadTree import QuadTree
+from geocamLayer.models import Feature, QuadTreeCell
+
+
+def importFeatures(tree, features, maxNumFeatures):
+ if maxNumFeatures:
+ features = itertools.islice(features, maxNumFeatures)
+ for feature in features:
+ tree.addFeature(feature)
+ sys.stdout.write('.')
+ sys.stdout.flush()
+
+
+def importData(opts):
+ if opts.clean:
+ print 'deleting old data'
+ # for loops are a work around for bulk delete problem in sqlite
+ # with objects.all().delete()
+ for f in Feature.objects.all():
+ f.delete()
+ for q in QuadTreeCell.objects.all():
+ q.delete()
+
+ tree = QuadTree()
+ if opts.noaa:
+ sys.stdout.write('noaa import')
+ importFeatures(tree, readNoaaWeatherStations(opts.noaa), opts.maxNumFeatures)
+ print
+ print
+ tree.debugStats()
+
+
+def main():
+ import optparse
+ parser = optparse.OptionParser('usage: %prog')
+ parser.add_option('-n', '--maxNumFeatures',
+ type='int',
+ help='Cap on number of features to import')
+ parser.add_option('-c', '--clean',
+ action='store_true', default=False,
+ help='[USE CAUTION!] Delete old data before import')
+ parser.add_option('--noaa',
+ help='Path to NOAA weather stations example file')
+ opts, args = parser.parse_args()
+ if args:
+ parser.error('expected no arguments')
+ importData(opts)
+
+
+if __name__ == '__main__':
+ main()
View
@@ -4,10 +4,13 @@
# All Rights Reserved.
# __END_LICENSE__
-#from django.contrib.gis.db import models
-#from django.contrib.gis.geos import Point
+import datetime
+import time
+import random
+import math
+
from django.db import models
-import datetime, time, random
+
class BaseFeature(models.Model):
# demo feature model demonstrating
@@ -44,6 +47,10 @@ class Feature(BaseFeature):
timespan = models.FloatField(null=True, blank=True)
name = models.CharField(max_length=80)
description = models.TextField()
+ cell = models.ForeignKey('QuadTreeCell', null=True, blank=True)
+
+ def __unicode__(self):
+ return u'Feature "%s" (%.6f, %.6f)' % (self.name, self.lng, self.lat)
@staticmethod
def randomFeature():
@@ -62,3 +69,85 @@ def getTimeSpan(self): return self.timespan
def getName(self): return self.name
def getDescriptionHTML(self): return self.description
def getProperties(self): return {} # self.properties
+
+class QuadTreeCell(models.Model):
+ zoom = models.PositiveIntegerField()
+ x = models.PositiveIntegerField()
+ y = models.PositiveIntegerField()
+ count = models.PositiveIntegerField(default=0)
+ isLeaf = models.BooleanField(default=True)
+ # centroid: lng, lat
+ lng = models.FloatField(default=0.0)
+ lat = models.FloatField(default=0.0)
+ # bbox: west, south, east, north
+ west = models.FloatField(null=True, blank=True)
+ south = models.FloatField(null=True, blank=True)
+ east = models.FloatField(null=True, blank=True)
+ north = models.FloatField(null=True, blank=True)
+
+ class Meta:
+ ordering = ('zoom', 'x', 'y')
+
+ def __unicode__(self):
+ return u'QuadTreeCell (%d, %d, %d) count=%d' % (self.zoom, self.x, self.y, self.count)
+
+ @staticmethod
+ def getSizeForZoom(zoom):
+ return 360.0 / (2.0 ** zoom)
+
+ @staticmethod
+ def getIndexAtLonLat(zoom, lonLat):
+ size = QuadTreeCell.getSizeForZoom(zoom)
+ lng, lat = lonLat
+ x = int((lng - (-180)) / size)
+ y = int((lat - (-90)) / size)
+ return (zoom, x, y)
+
+ @staticmethod
+ def getCellAtIndex(coords):
+ zoom, x, y = coords
+ cell, _created = QuadTreeCell.objects.get_or_create(zoom=zoom, x=x, y=y)
+ return cell
+
+ @staticmethod
+ def getCellAtLonLat(zoom, lonLat):
+ return (QuadTreeCell.getCellAtIndex
+ (QuadTreeCell.getIndexAtLonLat(zoom, lonLat)))
+
+ def getSize(self):
+ return QuadTreeCell.getSizeForZoom(self.zoom)
+
+ def getMinCorner(self):
+ size = self.getSize()
+ return (-180 + self.x * size,
+ -90 + self.y * size)
+
+ def getBounds(self):
+ size = self.getSize()
+ west, south = self.getMinCorner()
+ return (west, south, west + size, south + size)
+
+ def updateStats(self, feature):
+ self.lng = (self.count * self.lng + feature.lng) / (self.count + 1)
+ self.lat = (self.count * self.lat + feature.lat) / (self.count + 1)
+
+ if self.count:
+ self.west = min(self.west, feature.lng)
+ self.south = min(self.south, feature.lat)
+ self.east = max(self.east, feature.lng)
+ self.north = max(self.north, feature.lat)
+ else:
+ self.west = feature.lng
+ self.south = feature.lat
+ self.east = feature.lng
+ self.north = feature.lat
+
+ self.count += 1
+
+ def getDiameter(self):
+ if self.count:
+ dx = self.east - self.west
+ dy = self.north - self.south
+ return math.sqrt(dx**2 + dy**2)
+ else:
+ return 0
@@ -0,0 +1,44 @@
+# __BEGIN_LICENSE__
+# Copyright (C) 2008-2010 United States Government as represented by
+# the Administrator of the National Aeronautics and Space Administration.
+# All Rights Reserved.
+# __END_LICENSE__
+
+from geocamLayer.models import Feature
+
+
+def parseDegrees(data, positiveDirection):
+ direction = data[-1]
+ direction = 1 if direction == positiveDirection else -1
+ data = data[:-1]
+ fields = [float(x) for x in data.split('-')]
+ degrees, minutes = fields[:2]
+ val = degrees + minutes / 60.0
+ if len(fields) == 3:
+ seconds = fields[2]
+ val += seconds / 3600.0
+ return direction * val
+
+
+def parseLatitude(latStr):
+ return parseDegrees(latStr, 'N')
+
+
+def parseLongitude(latStr):
+ return parseDegrees(latStr, 'E')
+
+
+def readNoaaWeatherStations(fname):
+ for line in open(fname):
+ data = line.strip().split(';')
+
+ blockNumber, stationNumber, icaoLocation, placeName,\
+ state, countryName, wmoRegion, stationLatitude,\
+ stationLongitude, upperAirLatitude, upperAirLongitude,\
+ stationElevation, upperAirElevation, rbsnIndicator\
+ = data
+
+ yield Feature(lat=parseLatitude(stationLatitude),
+ lng=parseLongitude(stationLongitude),
+ name=placeName,
+ description=countryName)
View
@@ -0,0 +1,70 @@
+# __BEGIN_LICENSE__
+# Copyright (C) 2008-2010 United States Government as represented by
+# the Administrator of the National Aeronautics and Space Administration.
+# All Rights Reserved.
+# __END_LICENSE__
+
+import math
+
+from geocamLayer.models import Feature, QuadTreeCell
+
+MAX_FEATURES_PER_CELL = 10
+
+# beyond zoom MAX_ZOOM - 1, corresponding to cells approximately 1 meter
+# in size, all points are considered to be "in the same place" and we
+# don't split any further.
+MAX_ZOOM = 26
+
+
+class QuadTree(object):
+ def __init__(self):
+ self.root = QuadTreeCell.getCellAtIndex((0, 0, 0))
+
+ def addFeature(self, feature):
+ self.addFeatureToCell(feature, self.root)
+
+ def addFeatureToCell(self, feature, cell):
+ cell.updateStats(feature)
+
+ if cell.isLeaf:
+ feature.cell = cell
+ feature.save()
+
+ if cell.count >= MAX_FEATURES_PER_CELL and cell.zoom < MAX_ZOOM - 1:
+ self.splitCell(cell)
+ else:
+ self.addFeatureToZoom(feature, cell.zoom + 1)
+
+ cell.save()
+
+ def addFeatureToZoom(self, feature, zoom):
+ cell = (QuadTreeCell.getCellAtLonLat
+ (zoom, (feature.lng, feature.lat)))
+ self.addFeatureToCell(feature, cell)
+
+ def createCell(self, zoom, x, y):
+ return QuadTreeCell.objects.create(zoom=zoom, x=x, y=y)
+
+ def splitCell(self, cell):
+ cell.isLeaf = False
+ for feature in Feature.objects.filter(cell=cell):
+ self.addFeatureToZoom(feature, cell.zoom + 1)
+
+ def debugStats(self):
+ print '=== DEBUG STATS ==='
+ zoom = 0
+ while 1:
+ cellsAtZoom = QuadTreeCell.objects.filter(zoom=zoom)
+ if not cellsAtZoom:
+ break
+ numCellsAtZoom = cellsAtZoom.count()
+ numFeatures = 0
+ diameterSum = 0.0
+ for cell in cellsAtZoom:
+ numFeatures += cell.count
+ diameterSum += cell.getDiameter()
+ meanDiameter = diameterSum / numCellsAtZoom
+ maxNumCellsAtZoom = math.ceil(0.5 * 4 ** zoom)
+ print ('zoom=%d numCells=%d maxNumCells=%d numFeatures=%d meanDiameter=%.1f'
+ % (zoom, numCellsAtZoom, maxNumCellsAtZoom, numFeatures, meanDiameter))
+ zoom += 1

0 comments on commit b996564

Please sign in to comment.