Merge pull request #56 from openaddresses/preview-actual

Preview Actual [Experimental]

openaddr/preview.py

@@ -28,16 +28,9 @@
 def render(src_filename, png_filename, width, resolution, mapbox_key):
     '''
     '''
-    _, points_filename = mkstemp(prefix='points-', suffix='.bin')
-
     try:
-        _L.info('Writing from {} to {}...'.format(src_filename, points_filename))
-        points = project_lonlats(iterate_file_lonlats(src_filename))
-        write_points(points, points_filename)
-
-        xmin, ymin, xmax, ymax = calculate_bounds(points_filename)
+        xmin, ymin, xmax, ymax = calculate_bounds(src_filename)
     except:
-        os.remove(points_filename)
         raise
 
     surface, context, scale = make_context(xmin, ymin, xmax, ymax, width, resolution)
@@ -47,9 +40,10 @@ def render(src_filename, png_filename, width, resolution, mapbox_key):
     # Map units per reference pixel (http://www.w3.org/TR/css3-values/#reference-pixel)
     muppx = resolution / scale
 
+    feature_fill = 0x74/0xFF, 0xA5/0xFF, 0x78/0xFF
+
     black = 0x00, 0x00, 0x00
     off_white = 0xFF/0xFF, 0xFC/0xFF, 0xF9/0xFF
-    point_fill = 0x74/0xFF, 0xA5/0xFF, 0x78/0xFF
     water_fill = 0xC7/0xFF, 0xDE/0xFF, 0xF5/0xFF # 0xDD/0xFF, 0xEA/0xFF, 0xF8/0xFF
     road_stroke = 0xC0/0xFF, 0xE0/0xFF, 0xE0/0xFF # 0xE0/0xFF, 0xE3/0xFF, 0xE5/0xFF
     park_fill = 0xDD/0xFF, 0xF6/0xFF, 0xDE/0xFF
@@ -73,35 +67,46 @@ def render(src_filename, png_filename, width, resolution, mapbox_key):
 
     context.set_line_width(.25 * muppx)
 
-    for (x, y) in read_points(points_filename):
-        context.arc(x, y, 15, 0, 2 * pi)
-        context.set_source_rgb(*point_fill)
-        context.fill()
-        context.arc(x, y, 15, 0, 2 * pi)
-        context.set_source_rgb(*black)
-        context.stroke()
+    for geom in iterate_file_geoms(src_filename):
+        if geom.GetGeometryType() == ogr.wkbMultiLineString or geom.GetGeometryType() == ogr.wkbLineString:
+            stroke_geometries(context, [geom])
+        elif geom.GetGeometryType() == ogr.wkbMultiPolygon or geom.GetGeometryType() == ogr.wkbPolygon:
+            fill_geometries(context, [geom], muppx, feature_fill)
+        else:
+            (x, y, e) = geom.PointOnSurface().GetPoint()
+
+            context.arc(x, y, 15, 0, 2 * pi)
+            context.set_source_rgb(*feature_fill)
+            context.fill()
+            context.arc(x, y, 15, 0, 2 * pi)
+            context.set_source_rgb(*black)
+            context.stroke()
+            (x, y, e) = geom.PointOnSurface().GetPoint()
 
-    os.remove(points_filename)
     surface.write_to_png(png_filename)
 
-def iterate_file_lonlats(filename):
-    ''' Stream (lon, lat) coordinates from an input GeoJSON
+def iterate_file_geoms(filename):
+    ''' Stream Geometries from an input GeoJSON+LD File
     '''
 
     with open(filename, 'r') as file:
+        project = get_projection()
+
         for line in file:
             try:
                 line = json.loads(line)
 
-                lon, lat, x = ogr.CreateGeometryFromJson(json.dumps(line['geometry'])).PointOnSurface().GetPoint()
+                geom = ogr.CreateGeometryFromJson(json.dumps(line['geometry']))
+
+                geom.Transform(project)
 
-                if -180 <= lon <= 180 and -90 <= lat <= 90:
-                    yield (lon, lat)
+                yield geom
             except Exception as e:
                 print('ERROR', e)
-            except:
                 continue
 
+    del project, geom
+
 def get_map_features(xmin, ymin, xmax, ymax, resolution, scale, mapbox_key):
     '''
     '''
@@ -188,57 +193,29 @@ def get_projection():
     sref_map = osr.SpatialReference(); sref_map.ImportFromProj4(EPSG900913)
     return osr.CoordinateTransformation(sref_geo, sref_map)
 
-def project_lonlats(lonlats):
-    ''' Stream Mercator (x, y) points from a stream of (lon, lat) coordinates.
-    '''
-    project, geom = get_projection(), ogr.Geometry(ogr.wkbPoint)
-
-    for (lon, lat) in lonlats:
-        geom.SetPoint(0, lon, lat)
-        try:
-            geom.Transform(project)
-        except:
-            pass
-        else:
-            yield (geom.GetX(), geom.GetY())
-
-    del project, geom
-
-def write_points(points, points_filename):
-    ''' Write a stream of (x, y) points into a file of packed values.
+def write_geoms(geoms, geoms_filename):
+    ''' Write a stream of geoms into a file of packed values.
     '''
     count = 0
 
-    with open(points_filename, mode='wb') as file:
-        for (x, y) in points:
-            file.write(struct.pack(FORMAT, x, y))
+    with open(geoms_filename, mode='wb') as file:
+        for geom in geoms:
+            file.write(geom)
             count += 1
 
-    _L.info('Wrote {} points to {}'.format(count, points_filename))
-
-def read_points(points_filename):
-    ''' Read a file of packed values into a stream of (x, y) points.
-    '''
-    _L.debug('Reading from {}'.format(points_filename))
-    chunk_size = struct.calcsize(FORMAT)
-
-    with open(points_filename, mode='rb') as file:
-        while True:
-            chunk = file.read(chunk_size)
-            if chunk:
-                yield struct.unpack(FORMAT, chunk)
-            else:
-                return
+    _L.info('Wrote {} points to {}'.format(count, geoms_filename))
 
-def stats(points_filename):
-    ''' Return means and standard deviations for points in file.
+def stats(geoms_filename):
+    ''' Return means and standard deviations for iterator geoms
 
         Uses Welford's numerically stable algorithm from
         https://en.wikipedia.org/wiki/Algorithms_for_calculating_variance#Online_algorithm
     '''
     n, xmean, xM2, ymean, yM2 = 0, 0, 0, 0, 0
 
-    for (x, y) in read_points(points_filename):
+    for geom in iterate_file_geoms(geoms_filename):
+        (x, y, e) = geom.PointOnSurface().GetPoint()
+
         n += 1
 
         xdelta = x - xmean
@@ -264,10 +241,10 @@ def calculate_zoom(scale, resolution):
 
     return zoom
 
-def calculate_bounds(points_filename):
+def calculate_bounds(geoms_filename):
     '''
     '''
-    xmean, xsdev, ymean, ysdev = stats(points_filename)
+    xmean, xsdev, ymean, ysdev = stats(geoms_filename)
 
     # use standard deviation to avoid far-flung mistakes, and look further
     # horizontally to account for Github comment thread image appearance.
@@ -278,7 +255,9 @@ def calculate_bounds(points_filename):
     left, right = xmax, xmin
     bottom, top = ymax, ymin
 
-    for (x, y) in read_points(points_filename):
+    for geom in iterate_file_geoms(geoms_filename):
+        (x, y, e) = geom.PointOnSurface().GetPoint()
+
         if xmin <= x <= xmax:
             left, right = min(left, x), max(right, x)
         if ymin <= y <= ymax:
@@ -358,6 +337,8 @@ def fill_geometries(ctx, geometries, muppx, rgb):
         else:
             raise NotImplementedError()
 
+        print(parts)
+
         for part in parts:
             for ring in part:
                 points = ring.GetPoints()
@@ -374,7 +355,7 @@ def draw_line(ctx, start, points):
 
 parser = ArgumentParser(description='Draw a map of a single source preview.')
 
-parser.add_argument('src_filename', help='Input GeoJSON')
+parser.add_argument('src_geojson', help='Input GeoJSON')
 parser.add_argument('png_filename', help='Output PNG filename.')
 
 parser.set_defaults(resolution=1, width=668)
@@ -401,7 +382,7 @@ def draw_line(ctx, start, points):
 
 def main():
     args = parser.parse_args()
-    render(args.src_filename, args.png_filename, args.width, args.resolution, args.mapbox_key)
+    render(args.src_geojson, args.png_filename, args.width, args.resolution, args.mapbox_key)
 
 if __name__ == '__main__':
     exit(main())

openaddr/tests/preview.py

@@ -1,6 +1,7 @@
 from __future__ import division
 
 import os
+import json
 import unittest
 import tempfile
 import subprocess
@@ -21,24 +22,43 @@ def tearDown(self):
         rmtree(self.temp_dir)
 
     def test_stats(self):
-        points = [(n, n) for n in range(-1000, 1001)]
-        points_filename = join(self.temp_dir, 'points.bin')
-        preview.write_points(points, points_filename)
+        points = [(-108 + (n * 0.001), -37 + (n * 0.001)) for n in range(0, 1000)]
+        points_filename = join(self.temp_dir, 'points.geojson')
+
+        with open(points_filename, 'w') as file:
+            for point in points:
+                file.write(json.dumps({
+                    "type": "Feature",
+                    "properties": {},
+                    "geometry": {
+                        "type": "Point",
+                        "coordinates": point
+                    }
+                }) + '\n')
 
         xmean, xsdev, ymean, ysdev = preview.stats(points_filename)
-        self.assertAlmostEqual(xmean, 0)
-        self.assertAlmostEqual(xsdev, 577.783263863)
-        self.assertAlmostEqual(ymean, xmean)
-        self.assertAlmostEqual(ysdev, xsdev)
+        self.assertAlmostEqual(xmean, -11966900.920021897)
+        self.assertAlmostEqual(xsdev, 32151.232557143696)
 
     def test_calculate_bounds(self):
-        points = [(-10000, -10000), (10000, 10000)]
-        points += [(-1, -1), (0, 0), (1, 1)] * 100
-        points_filename = join(self.temp_dir, 'points.bin')
-        preview.write_points(points, points_filename)
+        points = [(-108 + (n * 0.001), -37 + (n * 0.001)) for n in range(0, 1000)]
+        points += [(-1, -1), (0, 0), (1, 1)]
+
+        points_filename = join(self.temp_dir, 'points.geojson')
+
+        with open(points_filename, 'w') as file:
+            for point in points:
+                file.write(json.dumps({
+                    "type": "Feature",
+                    "properties": {},
+                    "geometry": {
+                        "type": "Point",
+                        "coordinates": point
+                    }
+                }) + '\n')
 
         bbox = preview.calculate_bounds(points_filename)
-        self.assertEqual(bbox, (-1.04, -1.04, 1.04, 1.04), 'The two outliers are ignored')
+        self.assertEqual(bbox, (-12024729.169099594, -4441873.743568107, -11909072.670945017, -4297992.015057018), 'The two outliers are ignored')
 
     def test_render_geojson(self):
         '''