Merge pull request #14 from maxhully/master

Version used to create vtd adjacency graphs

Author: maxhully

Dockerfile

@@ -0,0 +1,13 @@
+FROM python:3.6.5
+
+RUN pip install pipenv
+
+COPY ./Pipfile ./Pipfile.lock /code/
+
+WORKDIR /code
+
+RUN pipenv install
+
+COPY ./graphmaker .
+
+ENTRYPOINT [ "bash" ]

docker-compose.yml

@@ -0,0 +1,18 @@
+version: '3.3'
+
+services: 
+  graphmaker:
+    build: .
+    volumes:
+    - ./graphmaker:/code
+    links:
+    - mongo
+    stdin_open: true
+    tty: true
+  mongo:
+    image: mongo
+    environment:
+      MONGO_INITDB_ROOT_USERNAME: root
+      MONGO_INITDB_ROOT_PASSWORD: example
+    volumes:
+    - ./mongodata:/data/db

graphmaker/Pipfile

@@ -0,0 +1,12 @@
+[[source]]
+verify_ssl = true
+name = "pypi"
+url = "https://pypi.org/simple"
+
+[packages]
+utm = "*"
+
+[requires]
+python_version = "3.5"
+
+[dev-packages]

graphmaker/Pipfile.lock

@@ -0,0 +1,21 @@
+from collections import Counter
+
+import utm
+
+
+def utm_of_point(point):
+    print(point.y, point.x)
+    return utm.from_latlon(point.y, point.x)[2]
+
+
+def identify_utm_zone(df):
+    utms = map(utm_of_point, df['geometry'].centroid)
+    utm_counts = Counter(utms)
+    # most_common returns a list of tuples, and we want the 0,0th entry
+    most_common = utm_counts.most_common(1)[0][0]
+    return most_common
+
+
+def reprojected(df):
+    utm = identify_utm_zone(df)
+    return df.to_crs(f"+proj=utm +zone={utm} +ellps=WGS84 +datum=WGS84 +units=m +no_defs")

graphmaker/geospatial.py

@@ -1,19 +1,15 @@
 import io
-import json
 import os
 import zipfile
 
-import networkx
 import requests
 
-from make_graph import construct_graph_from_file
-
 
 def zfill2(x):
     return str(x).zfill(2)
 
 
-fips_codes = ['01', '02'] + list(map(zfill2, range(4, 57)))
+fips_codes = ['01', '02', '04', '05', '06'] + list(map(zfill2, range(8, 57)))
 
 
 def state_zip_url(fips, year='2012'):
@@ -45,10 +41,6 @@ def main():
             continue
         try:
             download_state_vtds(fips)
-            graph = construct_graph_from_file(shp_location(
-                fips), 'GEOID10', ['ALAND10', 'COUNTYFP10'])
-            with open(output_location(fips), "w+") as f:
-                json.dump(networkx.json_graph.adjacency_data(graph), f)
         except Exception:
             print("An error occurred for FIPS code " + fips)
 

graphmaker/get_vtds.py

@@ -0,0 +1,13 @@
+from main import main
+import os
+
+def path_to_shp(fips):
+    return os.path.join('/code/tiger_data', fips, 'tl_2012_' + fips + '_vtd10.shp')
+
+states = os.listdir('/code/tiger_data/')
+states = [state for state in states if len(state) == 2]
+
+print(list(path_to_shp(fips) for fips in states))
+
+for fips in states:
+    main([path_to_shp(fips)])

graphmaker/make_all.py

@@ -7,6 +7,8 @@
 from networkx.readwrite import json_graph
 from shapely.ops import cascaded_union
 
+from geospatial import reprojected
+
 
 def get_list_of_data(filepath, col_name, geoid=None):
     """Pull a column data from a shape or CSV file.
@@ -66,6 +68,9 @@ def construct_graph_from_df(df, geoid_col=None, cols_to_add=None, queen=False):
     :returns: NetworkX Graph.
 
     """
+    # reproject to a UTM projection for accurate areas and perimeters in meters
+    df = reprojected(df)
+
     if geoid_col is not None:
         df = df.set_index(geoid_col)
 

graphmaker/make_graph.py

@@ -5,7 +5,6 @@
 import networkx
 import numpy
 import scipy
-
 from constants import round_to
 
 # import planarity

graphmaker/reports.py

@@ -0,0 +1,31 @@
+import geopandas
+from graphmaker.geospatial import identify_utm_zone, reprojected, utm_of_point
+from shapely.geometry import Point
+
+
+def example_point():
+    return Point(42.3, -83.7)
+
+
+def test_identify_utm_zone():
+    point = example_point()
+    data = [{'name': 'example point', 'geometry': point}]
+    df = geopandas.GeoDataFrame(data, crs="+init=epsg:4326")
+
+    assert identify_utm_zone(df) == 17
+
+
+def test_reprojected():
+    point = example_point()
+    data = [{'name': 'example point', 'geometry': point}]
+    df = geopandas.GeoDataFrame(data, crs="+init=epsg:4326")
+
+    reprojected_df = reprojected(df)
+    reprojected_point = reprojected_df['geometry'][0]
+
+    assert reprojected_point.x != point.x and reprojected_point.y != point.y
+
+
+def test_utm_of_point():
+    point = example_point()
+    assert utm_of_point(point) == 17