Create WeoGeoTableOfContents.json file for each FIPS code.

mk_county_toc.py

import os,sys
import math
import json
# 3rd party libs
import fiona

def get_s_swaths():
    extent = 20037508.342789244
    s_swaths = []
    for cnt in range(30):
        s_swaths.append(extent*2 / 2**cnt)
    return s_swaths

def get_l_swaths():
    l_swaths = []
    for cnt in range(30):
        l_swaths.append(180./2**cnt)
    return l_swaths
    
def find_optimal_first_scale_geo(swath,smorgeo):
    if smorgeo == "geo":
        l_swaths = get_l_swaths()
    else:
        l_swaths = get_s_swaths()
    
    for cnt in range(len(l_swaths)):
        if (swath) / 2. > l_swaths[cnt]:
            return cnt-2

def get_base(filelist):
    bases = []
    for afile in filelist:
        head, tail = os.path.split(afile)
        base, ext = os.path.splitext(tail)
        bases.append(base)
    return bases

def get_extent(features):
    most_west = most_east = most_south = most_north = None
    for feature in features:
        west = min([x[0] for x in feature["geometry"]["coordinates"][0]])
        east = max([x[0] for x in feature["geometry"]["coordinates"][0]])
        south = min([x[1] for x in feature["geometry"]["coordinates"][0]])
        north = max([x[1] for x in feature["geometry"]["coordinates"][0]])
        if most_west == None or west < most_west: most_west = west
        if most_east == None or east > most_east: most_east = east
        if most_south == None or south < most_south: most_south = south
        if most_north == None or north > most_north: most_north = north
    return most_west, most_south, most_east, most_north

if __name__ == "__main__":

    counties = [x.strip() for x in open("/mnt/tiger/counties_fips_list.txt").readlines()]

    for county in counties:
        print "working on county fips {}".format(county)
        out_path = "/mnt/tiger/TIGER2011/wg_fips/tl_2011_{}/".format(county)
        
        toc = {}
        toc["name"] = "NewFeatureType"
        toc["type"] = "FeatureCollection"
        features = []

        lines = open("/mnt/tiger/TIGER2011/wg_fips/tl_2011_{}/files_list.txt".format(county)).readlines()
        lines = [x.strip() for x in lines]  # remove line breaks from list
        dbfs = [x for x in lines if x[-3:]=="dbf"] # filter list to only .dbf files
        dbfs_base = get_base(dbfs)
        shps = [x for x in lines if x[-3:]=="shp"] # filter list to only .shp files
        shps_base = get_base(shps)

        non_spatial_dbfs = [d for d in dbfs_base if d not in shps_base]
        
        layer_cnt = 0
        for shp in shps:
            base,ext = os.path.splitext(shp)
            j_file_name = base+".json"
            j = json.loads(open(j_file_name).read())

            # add layers to feature
            j["features"][0]["properties"]["LAYERS"] = str(layer_cnt)
            
        # create a GeoJSON feature for the file
            features.append(j["features"][0])
            
        # create a .shp.xml feature
            features.append({
                "type" : "Feature",
                "geometry" : j["features"][0]["geometry"],
                "properties" : {
                    "LAYERS" : str(layer_cnt),
                    "PATH" : j["features"][0]["properties"]["PATH"],
                    "EXTS" : "shp.xml",
                    "WEO_TYPE" : "WEO_FEATURE",
                    "WEO_MISCELLANEOUS_FILE" : "yes"
                }
            })
            layer_cnt += 1
            
        toc["features"] = features
        
        # Add the non-spatial files    
        west, south, east, north = get_extent(features)
        
        for a_non_spatial_dbf in non_spatial_dbfs:
            path =  "./TIGER2011/{}".format(a_non_spatial_dbf)
            features.append(
            {
            "type":"Feature",
            "geometry":{"type": "Polygon",
                "coordinates":[[
                [west,north], #UL  X,Y   c.bounds
                [east,north], #UR  X,Y
                [east,south], #LR  X,Y
                [west,south], #LL  X,Y
                [west,north]  #UL  X,Y
                ]]},
            "properties":{
            "PATH": path,
            "EXTS": 'dbf;dbf.xml',
            "LAYERS" : str(layer_cnt),
            "WEO_MISCELLANEOUS_FILE":"Yes",
            "WEO_TYPE":"WEO_FEATURE"
                }
            })
            layer_cnt += 1

        # LUT Feature
        features.append({
                    "geometry" : None,
                    "type" : "Feature",
                    "properties" : {
                        "WEO_TYPE" : "LOOK_UP_TABLE",}})

        for i in range(int(layer_cnt)):
            features[-1]["properties"][i] = "WEOALL=WEOALL"
            
        toc["features"] = features
        
        outfilename = "WeoGeoTableOfContents.json"
        outname = os.path.join(out_path + outfilename)
        out = open(outname, "w")
        
        # create a JSON object
        e = json.JSONEncoder()
        out.write(e.encode(toc))
        out.close()
        
        # write render.sh
        of = open(os.path.join(out_path + "render.sh"),"w")
        optimal_first_scale = find_optimal_first_scale_geo(max(east-west,north-south),"geo")
        
        of.write("render.py -w {} -s {} -e {} -n {} -c {} -z {} --generate-tilepack=. style.xml".format(
                                west,south,east,north,
                                optimal_first_scale, 
                                13-optimal_first_scale+1))
        of.close()