import sys
import matplotlib.pyplot as plt
import matplotlib.tri as tri
from matplotlib import path
from shapely.geometry import geo, Polygon, Point
import numpy

# create contours for test case
def create_test_contours():

  x = numpy.array( [-76.34244,-76.36786,-76.35223,-76.40118,-76.3905,
       -76.41039,-76.42416,-76.43492,-76.4437,-76.44031,
       -76.47977,-76.47267,-76.45914,-76.50938,-76.50371,
       -76.50534,-76.53763,-76.53808,-76.53674,-76.57217,
       -76.57292,-76.56992,-76.61163,-76.60898,-76.59344,
       -76.65379,-76.65218,-76.62014] )
  y = numpy.array( [37.00411,36.94697,36.91855,36.97744,36.89212,36.94417,
       36.9803,36.90823,36.95995,36.99866,36.9427,36.98559,
       37.02776,36.97269,37.01178,37.05462,36.99284,37.03721,
       37.08237,37.02015,37.06149,37.09983,37.04387,37.08551,
       37.1161,37.08103,37.11076,37.15342] )
  z = numpy.ma.masked_equal( [0.57833663768,0.59585130282,0.57725037291,0.61816105718,
       -99999,0.62808432966,0.62889707758,0.63214800973,0.62361256066,
       0.63571681755,0.63464181838,0.63445566205,-99999,0.63305024240,
       0.63890298185,-99999,-99999,0.0014958446745,-99999,
       0.0035448925787,0.00056665199692,-99999,0.0022350042433,0.0013006928908,
       -99999,0.0035717841844,0.0036461797716,-99999], -99999 )

  ele = [(1,0,3),(2,1,4),(4,1,5),(5,1,3),(5,3,6),(7,4,5),(7,5,8),(8,5,6),
         (6,9,8),(7,8,10),(10,8,11),(8,9,11),(9,12,11),(10,11,13),(13,11,14),(11,12,14),
         (14,12,15),(13,14,16),(16,14,17),(14,15,17),(17,15,18),(16,17,19),(19,17,20),(17,18,20),
         (20,18,21),(19,20,22),(22,20,23),(20,21,23),(23,21,24),(22,23,25),(25,23,26),(23,24,26),
         (24,27,26)]

  triang = tri.Triangulation(x, y, triangles=ele)
  levels = numpy.linspace(0, numpy.max(z), num=31)

  contours = plt.tricontourf(triang, z, levels=levels)
  plt.tricontour(triang, z, levels=levels)

  return contours

# classify polygons based on their area
def classify_polygons(polys):

  def signed_area(ring):
    v2 = numpy.roll(ring, -1, axis=0)
    return numpy.cross(ring, v2).sum() / 2.0

  exteriors = []
  interiors = []
  for p in polys:
    if signed_area(p) >= 0:
      exteriors.append(p)
    else:
      interiors.append(p)

  return exteriors, interiors


# returns mask for points which are inside given polygon
def points_inside_poly(points, polygon):
    p = path.Path(polygon)
    return p.contains_points(points)


# extract shapely Polygon objects from the contours
def extract_geometries(contours):

  geoms = []
  for coll in contours.collections:  

    for p in coll.get_paths():

      polys = [p for p in p.to_polygons() if p.shape[0] >= 3]
      exteriors, interiors = classify_polygons(polys)
      if len(interiors) > 0:
        interior_points = [pts[0] for pts, a in interiors]

      for exterior in exteriors:
        if len(interiors) > 0:
          inout = points_inside_poly(interior_points, exterior)
          my_interiors = [g for i, g in enumerate(interiors) if inout[i]]
          poly = Polygon(exterior, my_interiors)

        else:
          poly = Polygon(exterior)

        geoms.append(poly)

  return geoms


def main(argv = None):
    if argv is None:
        argv = sys.argv
    contours = create_test_contours()
    geoms = extract_geometries(contours)

    plt.show()

if __name__ == "__main__":
    sys.exit(main())