#!/bin/bash

# EPSG:3310 California Albers
PROJECTION='d3.geoAlbers().parallels([34, 40.5]).rotate([120, 0])'

# The state FIPS code.
STATE=06

# The ACS 5-Year Estimate vintage.
YEAR=2014

# The display size.
WIDTH=960
HEIGHT=1100

# Download the census tract boundaries.
# Extract the shapefile (.shp) and dBASE (.dbf).
if [ ! -f cb_${YEAR}_${STATE}_tract_500k.shp ]; then
  curl -o cb_${YEAR}_${STATE}_tract_500k.zip \
    "http://www2.census.gov/geo/tiger/GENZ${YEAR}/shp/cb_${YEAR}_${STATE}_tract_500k.zip"
  unzip -o \
    cb_${YEAR}_${STATE}_tract_500k.zip \
    cb_${YEAR}_${STATE}_tract_500k.shp \
    cb_${YEAR}_${STATE}_tract_500k.dbf
fi

# Download the census tract population estimates.
if [ ! -f cb_${YEAR}_${STATE}_tract_B01003.json ]; then
  curl -o cb_${YEAR}_${STATE}_tract_B01003.json \
    "http://api.census.gov/data/${YEAR}/acs5?get=B01003_001E&for=tract:*&in=state:${STATE}"
fi

# 1. Convert to GeoJSON.
# 2. Project.
# 3. Join with the census data.
# 4. Compute the population density.
# 5. Simplify.
# 6. Compute the county borders.
geo2topo -n \
  tracts=<(ndjson-join 'd.id' \
    <(shp2json cb_${YEAR}_${STATE}_tract_500k.shp \
      | geoproject "${PROJECTION}.fitExtent([[10, 10], [${WIDTH} - 10, ${HEIGHT} - 10]], d)" \
      | ndjson-split 'd.features' \
      | ndjson-map 'd.id = d.properties.GEOID.slice(2), d') \
    <(ndjson-cat cb_${YEAR}_${STATE}_tract_B01003.json \
      | ndjson-split 'd.slice(1)' \
      | ndjson-map '{id: d[2] + d[3], B01003: +d[0]}') \
    | ndjson-map -r d3=d3-array 'd[0].properties = {density: d3.bisect([1, 10, 50, 200, 500, 1000, 2000, 4000], (d[1].B01003 / d[0].properties.ALAND || 0) * 2589975.2356)}, d[0]') \
  | topomerge -k 'd.id.slice(0, 3)' counties=tracts \
  | topomerge --mesh -f 'a !== b' counties=counties \
  | topomerge -k 'd.properties.density' tracts=tracts \
  | toposimplify -p 1 -f \
  > topo.json

# Re-compute the topology as a further optimization.
# This consolidates unique sequences of arcs.
# https://github.com/topojson/topojson-simplify/issues/4
topo2geo \
  < topo.json \
  tracts=tracts.json \
  counties=counties.json

geo2topo \
  tracts=tracts.json \
  counties=counties.json \
  | topoquantize 1e5 \
  > topo.json

rm tracts.json counties.json