Merge pull request #80 from openaddresses/add-us-map-render-#77

Add USA map render

openaddr/geodata/cb_2013_us_county_20m-2163.prj

@@ -0,0 +1 @@
+PROJCS["US_National_Atlas_Equal_Area",GEOGCS["GCS_Unspecified datum based upon the Clarke 1866 Authalic Sphere",DATUM["D_Sphere_Clarke_1866_Authalic",SPHEROID["Clarke_1866_Authalic_Sphere",6370997,0]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],PROJECTION["Lambert_Azimuthal_Equal_Area"],PARAMETER["latitude_of_origin",45],PARAMETER["central_meridian",-100],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["Meter",1]]

openaddr/geodata/cb_2013_us_nation_20m-2163.prj

@@ -0,0 +1 @@
+PROJCS["US_National_Atlas_Equal_Area",GEOGCS["GCS_Unspecified datum based upon the Clarke 1866 Authalic Sphere",DATUM["D_Sphere_Clarke_1866_Authalic",SPHEROID["Clarke_1866_Authalic_Sphere",6370997,0]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],PROJECTION["Lambert_Azimuthal_Equal_Area"],PARAMETER["latitude_of_origin",45],PARAMETER["central_meridian",-100],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["Meter",1]]

openaddr/geodata/cb_2013_us_state_20m-2163.prj

@@ -0,0 +1 @@
+PROJCS["US_National_Atlas_Equal_Area",GEOGCS["GCS_Unspecified datum based upon the Clarke 1866 Authalic Sphere",DATUM["D_Sphere_Clarke_1866_Authalic",SPHEROID["Clarke_1866_Authalic_Sphere",6370997,0]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],PROJECTION["Lambert_Azimuthal_Equal_Area"],PARAMETER["latitude_of_origin",45],PARAMETER["central_meridian",-100],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["Meter",1]]

openaddr/process_all.py

@@ -65,10 +65,19 @@ def main():
     # Talk about the work
     #
     _L.info("Talking about the work")
-    png_filename = 'render-{}.png'.format(run_name)
-    render.render(paths.sources, good_sources, 960, 2, png_filename)
-    render_data = open(png_filename).read()
-    render_key = s3.new_key(join('runs', run_name, 'render.png'))
+    for (area, area_name) in [(render.WORLD, 'world'), (render.USA, 'usa')]:
+        png_filename = 'render-{}-{}.png'.format(run_name, area_name)
+        render.render(paths.sources, good_sources, 960, 2, png_filename, area)
+        render_data = open(png_filename).read()
+        render_path = 'render-{}.png'.format(area_name)
+        render_key = s3.new_key(join('runs', run_name, render_path))
+        render_key.set_contents_from_string(render_data, policy='public-read',
+                                            headers={'Content-Type': 'image/png'})
+
+    png_usa_filename = 'render-{}-usa.png'.format(run_name)
+    render.render(paths.sources, good_sources, 960, 2, png_usa_filename, render.USA)
+    render_data = open(png_usa_filename).read()
+    render_key = s3.new_key(join('runs', run_name, 'render-usa.png'))
     render_key.set_contents_from_string(render_data, policy='public-read',
                                         headers={'Content-Type': 'image/png'})
 

openaddr/render.py

@@ -16,13 +16,29 @@
 
 from . import paths
 
-def make_context(width=960, resolution=1):
+# Areas
+WORLD, USA = 54029, 2163
+
+def make_context(width=960, resolution=1, area=WORLD):
     ''' Get Cairo surface, context, and drawing scale.
     
-        World extent: (-19918964.35, -8269767.91) - (19918964.18, 14041770.96)
+        Global extent, World Van der Grinten I:
+        (-19918964.35, -8269767.91) - (19918964.18, 14041770.96)
+        http://spatialreference.org/ref/esri/54029/
+
+        U.S. extent, National Atlas Equal Area:
+        (-2031905.05, -2114924.96) - (2516373.83, 732103.34)
+        http://spatialreference.org/ref/epsg/2163/
     '''
-    left, top = -18000000, 14050000
-    right, bottom = 19500000, -7500000
+    if area == WORLD:
+        left, top = -18000000, 14050000
+        right, bottom = 19500000, -7500000
+    elif area == USA:
+        left, top = -2040000, 740000
+        right, bottom = 2525000, -2130000
+    else:
+        raise RuntimeError('Unknown area "{}"'.format(area))
+
     aspect = (right - left) / (top - bottom)
 
     hsize = int(resolution * width)
@@ -196,7 +212,7 @@ def draw_line(ctx, start, points):
 
 parser = ArgumentParser(description='Draw a map of worldwide address coverage.')
 
-parser.set_defaults(resolution=1, width=960)
+parser.set_defaults(resolution=1, width=960, area=WORLD)
 
 parser.add_argument('--2x', dest='resolution', action='store_const', const=2,
                     help='Draw at double resolution.')
@@ -212,12 +228,19 @@ def draw_line(ctx, start, points):
 
 parser.add_argument('filename', help='Output PNG filename.')
 
+parser.add_argument('--world', dest='area', action='store_const', const=WORLD,
+                    help='Render the whole world.')
+
+parser.add_argument('--usa', dest='area', action='store_const', const=USA,
+                    help='Render the United States.')
+
 def main():
     args = parser.parse_args()
     good_sources = load_live_state() if args.use_state else load_fake_state(paths.sources)
-    return render(paths.sources, good_sources, args.width, args.resolution, args.filename)
+    return render(paths.sources, good_sources, args.width, args.resolution,
+                  args.filename, args.area)
 
-def render(sources_dir, good_sources, width, resolution, filename=None):
+def render(sources_dir, good_sources, width, resolution, filename=None, area=WORLD):
     ''' Resolution: 1 for 100%, 2 for 200%, etc.
     '''
     if filename is None:
@@ -229,36 +252,56 @@ def render(sources_dir, good_sources, width, resolution, filename=None):
         # Use fake sources
         good_sources = load_fake_state(sources_dir)
 
-    return _render_state(sources_dir, good_sources, width, resolution, filename)
+    return _render_state(sources_dir, good_sources, width, resolution, filename, area)
+
+def first_layer_list(datasource):
+    ''' Return features as a list, or an empty list.
+    '''
+    return list(datasource.GetLayer(0) if hasattr(datasource, 'GetLayer') else [])
 
-def _render_state(sources_dir, good_sources, width, resolution, filename):
+def _render_state(sources_dir, good_sources, width, resolution, filename, area):
     ''' Resolution: 1 for 100%, 2 for 200%, etc.
     '''
     # Prepare output surface
-    surface, context, scale = make_context(width, resolution)
+    surface, context, scale = make_context(width, resolution, area)
 
     # Load data
     good_geoids, bad_geoids = load_geoids(sources_dir, good_sources)
     good_iso3166s, bad_iso3166s = load_iso3166s(sources_dir, good_sources)
     good_geometries, bad_geometries = load_geometries(sources_dir, good_sources)
 
     geodata = join(dirname(__file__), 'geodata')
-    coastline_ds = ogr.Open(join(geodata, 'ne_50m_coastline-54029.shp'))
-    lakes_ds = ogr.Open(join(geodata, 'ne_50m_lakes-54029.shp'))
-    countries_ds = ogr.Open(join(geodata, 'ne_50m_admin_0_countries-54029.shp'))
-    countries_borders_ds = ogr.Open(join(geodata, 'ne_50m_admin_0_boundary_lines_land-54029.shp'))
-    admin1s_ds = ogr.Open(join(geodata, 'ne_10m_admin_1_states_provinces-54029.shp'))
-    us_state_ds = ogr.Open(join(geodata, 'cb_2013_us_state_20m-54029.shp'))
-    us_county_ds = ogr.Open(join(geodata, 'cb_2013_us_county_20m-54029.shp'))
+
+    if area == WORLD:
+        landarea_ds = ogr.Open(join(geodata, 'ne_50m_admin_0_countries-54029.shp'))
+        coastline_ds = ogr.Open(join(geodata, 'ne_50m_coastline-54029.shp'))
+        lakes_ds = ogr.Open(join(geodata, 'ne_50m_lakes-54029.shp'))
+        countries_ds = ogr.Open(join(geodata, 'ne_50m_admin_0_countries-54029.shp'))
+        countries_borders_ds = ogr.Open(join(geodata, 'ne_50m_admin_0_boundary_lines_land-54029.shp'))
+        admin1s_ds = ogr.Open(join(geodata, 'ne_10m_admin_1_states_provinces-54029.shp'))
+        us_state_ds = ogr.Open(join(geodata, 'cb_2013_us_state_20m-54029.shp'))
+        us_county_ds = ogr.Open(join(geodata, 'cb_2013_us_county_20m-54029.shp'))
+    elif area == USA:
+        landarea_ds = ogr.Open(join(geodata, 'cb_2013_us_nation_20m-2163.shp'))
+        coastline_ds = ogr.Open(join(geodata, 'cb_2013_us_nation_20m-2163.shp'))
+        lakes_ds = None
+        countries_ds = None
+        countries_borders_ds = None
+        admin1s_ds = None
+        us_state_ds = ogr.Open(join(geodata, 'cb_2013_us_state_20m-2163.shp'))
+        us_county_ds = ogr.Open(join(geodata, 'cb_2013_us_county_20m-2163.shp'))
+    else:
+        raise RuntimeError('Unknown area "{}"'.format(area))
 
     # Pick out features
-    coastline_features = list(coastline_ds.GetLayer(0))
-    lakes_features = [f for f in lakes_ds.GetLayer(0) if f.GetField('scalerank') == 0]
-    countries_features = list(countries_ds.GetLayer(0))
-    countries_borders_features = list(countries_borders_ds.GetLayer(0))
-    admin1s_features = list(admin1s_ds.GetLayer(0))
-    us_state_features = list(us_state_ds.GetLayer(0))
-    us_county_features = list(us_county_ds.GetLayer(0))
+    landarea_features = first_layer_list(landarea_ds)
+    coastline_features = first_layer_list(coastline_ds)
+    lakes_features = [f for f in first_layer_list(lakes_ds) if f.GetField('scalerank') == 0]
+    countries_features = first_layer_list(countries_ds)
+    countries_borders_features = first_layer_list(countries_borders_ds)
+    admin1s_features = first_layer_list(admin1s_ds)
+    us_state_features = first_layer_list(us_state_ds)
+    us_county_features = first_layer_list(us_county_ds)
 
     # Assign features to good or bad lists
     good_data_states = [f for f in us_state_features if f.GetFieldAsString('GEOID') in good_geoids]
@@ -284,8 +327,8 @@ def _render_state(sources_dir, good_sources, width, resolution, filename):
     light_green = 0x74/0xff, 0xA5/0xff, 0x78/0xff
     dark_green = 0x1C/0xff, 0x89/0xff, 0x3F/0xff
 
-    # Fill countries background
-    fill_features(context, countries_features, silver)
+    # Fill land area background
+    fill_features(context, landarea_features, silver)
 
     # Fill populated countries
     fill_features(context, bad_data_countries, light_red)

openaddr/templates/state.html

@@ -161,6 +161,18 @@
 
         return false;
     }
+
+    function choosemap(src)
+    {
+        var img = document.getElementById('render-img');
+
+        if(img)
+        {
+            img.src = src;
+        }
+
+        return false;
+    }
 
 </script>
 </head>
@@ -184,7 +196,12 @@
     A plain text version of data on this page can be found at
     <a href="state.txt">state.txt</a>.
 </p>
-<p><img src="render.png" width="100%"></p>
+<p>
+    Show map:
+    <button onclick="choosemap('render-world.png')">World</button>
+    <button onclick="choosemap('render-usa.png')">United States</button>
+</p>
+<p><img src="render-world.png" width="100%" id="render-img"></p>
 <table>
     <tr>
         <th class="name">Name</th>