#!/usr/bin/python2
###############################################################################
# $Id: tigerpoly.py c040b1b6de0393a47c70cd5b8f8b2b692d43a3f4 2018-05-09 02:30:58 +1000 Ben Elliston $
#
# Project:  OGR Python samples
# Purpose:  Assemble TIGER Polygons.
# Author:   Frank Warmerdam, warmerdam@pobox.com
#
###############################################################################
# Copyright (c) 2003, Frank Warmerdam <warmerdam@pobox.com>
# Copyright (c) 2009, Even Rouault <even dot rouault at mines-paris dot org>
#
# Permission is hereby granted, free of charge, to any person obtaining a
# copy of this software and associated documentation files (the "Software"),
# to deal in the Software without restriction, including without limitation
# the rights to use, copy, modify, merge, publish, distribute, sublicense,
# and/or sell copies of the Software, and to permit persons to whom the
# Software is furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included
# in all copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
# FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
# DEALINGS IN THE SOFTWARE.
###############################################################################

import sys

from osgeo import ogr
from osgeo import osr


#############################################################################
class Module(object):

    def __init__(self):
        self.lines = {}
        self.poly_line_links = {}

#############################################################################


def Usage():
    print('Usage: tigerpoly.py infile [outfile].shp')
    print('')
    sys.exit(1)

#############################################################################
# Argument processing.


infile = None
outfile = None

i = 1
while i < len(sys.argv):
    arg = sys.argv[i]

    if infile is None:
        infile = arg

    elif outfile is None:
        outfile = arg

    else:
        Usage()

    i = i + 1

if outfile is None:
    outfile = 'poly.shp'

if infile is None:
    Usage()

#############################################################################
# Open the datasource to operate on.

ds = ogr.Open(infile, update=0)

poly_layer = ds.GetLayerByName('Polygon')

#############################################################################
# Create output file for the composed polygons.

nad83 = osr.SpatialReference()
nad83.SetFromUserInput('NAD83')

shp_driver = ogr.GetDriverByName('ESRI Shapefile')
shp_driver.DeleteDataSource(outfile)

shp_ds = shp_driver.CreateDataSource(outfile)

shp_layer = shp_ds.CreateLayer('out', geom_type=ogr.wkbPolygon,
                               srs=nad83)

src_defn = poly_layer.GetLayerDefn()
poly_field_count = src_defn.GetFieldCount()

for fld_index in range(poly_field_count):
    src_fd = src_defn.GetFieldDefn(fld_index)

    fd = ogr.FieldDefn(src_fd.GetName(), src_fd.GetType())
    fd.SetWidth(src_fd.GetWidth())
    fd.SetPrecision(src_fd.GetPrecision())
    shp_layer.CreateField(fd)

#############################################################################
# Read all features in the line layer, holding just the geometry in a hash
# for fast lookup by TLID.

line_layer = ds.GetLayerByName('CompleteChain')
line_count = 0

modules_hash = {}

feat = line_layer.GetNextFeature()
geom_id_field = feat.GetFieldIndex('TLID')
tile_ref_field = feat.GetFieldIndex('MODULE')
while feat is not None:
    geom_id = feat.GetField(geom_id_field)
    tile_ref = feat.GetField(tile_ref_field)

    try:
        module = modules_hash[tile_ref]
    except KeyError:
        module = Module()
        modules_hash[tile_ref] = module

    module.lines[geom_id] = feat.GetGeometryRef().Clone()
    line_count = line_count + 1

    feat.Destroy()

    feat = line_layer.GetNextFeature()

print('Got %d lines in %d modules.' % (line_count, len(modules_hash)))

#############################################################################
# Read all polygon/chain links and build a hash keyed by POLY_ID listing
# the chains (by TLID) attached to it.

link_layer = ds.GetLayerByName('PolyChainLink')

feat = link_layer.GetNextFeature()
geom_id_field = feat.GetFieldIndex('TLID')
tile_ref_field = feat.GetFieldIndex('MODULE')
lpoly_field = feat.GetFieldIndex('POLYIDL')
rpoly_field = feat.GetFieldIndex('POLYIDR')

link_count = 0

while feat is not None:
    module = modules_hash[feat.GetField(tile_ref_field)]

    tlid = feat.GetField(geom_id_field)

    lpoly_id = feat.GetField(lpoly_field)
    rpoly_id = feat.GetField(rpoly_field)

    if lpoly_id == rpoly_id:
        feat.Destroy()
        feat = link_layer.GetNextFeature()
        continue

    try:
        module.poly_line_links[lpoly_id].append(tlid)
    except KeyError:
        module.poly_line_links[lpoly_id] = [tlid]

    try:
        module.poly_line_links[rpoly_id].append(tlid)
    except KeyError:
        module.poly_line_links[rpoly_id] = [tlid]

    link_count = link_count + 1

    feat.Destroy()

    feat = link_layer.GetNextFeature()

print('Processed %d links.' % link_count)

#############################################################################
# Process all polygon features.

feat = poly_layer.GetNextFeature()
tile_ref_field = feat.GetFieldIndex('MODULE')
polyid_field = feat.GetFieldIndex('POLYID')

poly_count = 0
degenerate_count = 0

while feat is not None:
    module = modules_hash[feat.GetField(tile_ref_field)]
    polyid = feat.GetField(polyid_field)

    tlid_list = module.poly_line_links[polyid]

    link_coll = ogr.Geometry(type=ogr.wkbGeometryCollection)
    for tlid in tlid_list:
        geom = module.lines[tlid]
        link_coll.AddGeometry(geom)

    try:
        poly = ogr.BuildPolygonFromEdges(link_coll)

        if poly.GetGeometryRef(0).GetPointCount() < 4:
            degenerate_count = degenerate_count + 1
            poly.Destroy()
            feat.Destroy()
            feat = poly_layer.GetNextFeature()
            continue

        # print poly.ExportToWkt()
        # feat.SetGeometryDirectly( poly )

        feat2 = ogr.Feature(feature_def=shp_layer.GetLayerDefn())

        for fld_index in range(poly_field_count):
            feat2.SetField(fld_index, feat.GetField(fld_index))

        feat2.SetGeometryDirectly(poly)

        shp_layer.CreateFeature(feat2)
        feat2.Destroy()

        poly_count = poly_count + 1
    except:
        print('BuildPolygonFromEdges failed.')

    feat.Destroy()

    feat = poly_layer.GetNextFeature()

if degenerate_count:
    print('Discarded %d degenerate polygons.' % degenerate_count)

print('Built %d polygons.' % poly_count)

#############################################################################
# Cleanup

shp_ds.Destroy()
ds.Destroy()