%matplotlib inline
import shapefile
from functools import reduce

from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
import numpy as np

sf = shapefile.Reader("Data/cb_2015_us_zcta510_500k/cb_2015_us_zcta510_500k")

def find_max_and_min(shape_records):
    points = map(lambda x: x.shape.points, shape_records)
    points = reduce(lambda x, y: x + y, points)
    points = np.array(points)
    lon_min, lon_max = points[:, 0].min() - .05, points[:, 0].max() + .05
    lat_min, lat_max = points[:, 1].min() - .05, points[:, 1].max() + .05
    return lon_min, lon_max, lat_min, lat_max

def sub_polys_in_poly(poly_points):
    polys = []
    poly = []
    for index, point in enumerate(poly_points):
        poly.append(point)
        if len(poly) >= 2 and poly[0] == poly[len(poly)-1]:
            polys.append(poly)
            poly = []
    return polys

def plotShapes(plt, shape_records, colors):
    lon_min, lon_max, lat_min, lat_max = find_max_and_min(shape_records)
    
    map = Basemap(projection='merc', 
    lat_0 = np.average([lat_min, lat_max]), lon_0 = np.average([lon_min, lon_max]),
    resolution = 'c', area_thresh = 0.1,
    llcrnrlon=lon_min, llcrnrlat=lat_min,
    urcrnrlon=lon_max, urcrnrlat=lat_max)
    
    map.drawcoastlines()
    map.drawstates()
    map.fillcontinents(color = 'coral')
    map.drawmapboundary()
    
    for index, shape_record in enumerate(shape_records):
        polys = sub_polys_in_poly(shape_record.shape.points)
        for poly in polys:
            points = np.array(poly)
            x, y = map(points[:, 0], points[:,1])
            map.plot(x, y, 'o-', color=colors[index])

"""
Explanation: GIS Shapes -- An Experiment
Data source: Cartographic Boundary Shapefiles - ZIP Code Tabulation Areas (ZCTAs)
Note: Zip Codes are Not Areas
End of explanation
"""


shape_records = {}
for shape_record in sf.shapeRecords():
    zip_code = shape_record.record[0]
    shape_records[zip_code] = shape_record

"""
Explanation: Reading shape & records
End of explanation
"""


# def of record fields
sf.fields

# number of shape records
len(shape_records)

"""
Explanation: Metadata
End of explanation
"""


# object toomsuba & cuba are shape records, each of them has a shape and a record property
toomsuba = shape_records['39364']
cuba = shape_records['36907']
ithaca = shape_records['14850']

toomsuba.record

plt.figure(figsize=(12, 12))
plotShapes(plt, [toomsuba, cuba], ['b', 'y'])
plt.show()

plt.figure(figsize=(12, 12))
plotShapes(plt, [ithaca], ['b'])
plt.show()

"""
Explanation: Read zip codes
Google Maps: 
* Toomsuba, MS
* Cuba, AL
* Ithaca, NY (have excludes inside region)
End of explanation
"""


from shapely.geometry import Polygon
from shapely.ops import cascaded_union

# A Simple union example
triangle0 = Polygon([(0,0), (2,0), (0,2)])
triangle1 = Polygon([(1,0), (1,2), (3,0)])
triangle0.union(triangle1)

polygon0 = Polygon(toomsuba.shape.points)
polygon1 = Polygon(cuba.shape.points)
union_poly = polygon0.union(polygon1)

polygon0

polygon1

union_poly

# get coordinates of unoined shape
union_coords = union_poly.exterior.coords.xy
union_coords = list(zip(union_coords[0], union_coords[1]))
union_coords # [(lon, lat)]

"""
Explanation: Shapely Union - Works
End of explanation
"""


# from osgeo import ogr

# dataSource = ogr.Open("cb_2015_us_zcta510_500k/cb_2015_us_zcta510_500k.shp")
# layer = dataSource.GetLayer()
# print(layer.GetGeomType()) # 3 means polygons

# print(layer.GetFeatureCount()) 

# geometries = []
# for feature in layer:
#     geom = feature.GetGeometryRef()
#     geometries.append(geom)

# geo0 = geometries[0]
# geo1 = geometries[1]
# union_poly = geo0.Union(geo1)

"""
Explanation: OSGEO Unions - Doesn't work
Python GDAL/OGR Cookbook!
From lib osgeo doc: Calculate union between two Geometries
End of explanation
"""


