def pointinside(lat, lon, shapefile):

    # Verifica todos os pontos que estão dentro do polígono
    # http://streamhacker.com/2010/03/23/python-point-in-polygon-shapely/
    # '/home/rodrigues/AmbientePython27/lib/python2.7/site-packages/PyFuncemeClimateTools/shp/pontos_ce.txt'

    # Ler os pontos do vértice e tranformar em um poligno

    poly = Polygon(shapefile)

    nlons = len(lon)
    nlats = len(lat)

    points_grid = []
    lonlat_grid = []

    array_bool = np.ones((nlats, nlons), dtype="bool")

    for xlon in range(0, nlons):

        for ylat in range(0, nlats):

            point = Point((lon[xlon], lat[ylat]))

            a = poly.contains(point)

            if a == True:

                array_bool[ylat, xlon] = False

                points_grid.append((ylat, xlon))

                lonlat_grid.append((lat[ylat], lon[xlon]))

    return points_grid, lonlat_grid, array_bool

 def nms_discard(self, proposal, accepted_detections, dataframe):
 
   p_idx = proposal[0]
   p_label = proposal[1].index[0]
   p_xmin = dataframe.iloc[p_idx]['xmin']
   p_xmax = dataframe.iloc[p_idx]['xmax']
   p_ymin = dataframe.iloc[p_idx]['ymin']
   p_ymax = dataframe.iloc[p_idx]['ymax']
   p_poly = Polygon([(p_xmin,p_ymin), (p_xmax,p_ymin), (p_xmax,p_ymax), (p_xmin, p_ymax)])
 
   for detection in accepted_detections:
     detection = accepted_detections[detection]
     d_idx = detection[0]
     d_label = detection[1].index[0]
     if d_label != p_label:
       # No point checking if it isn't the same class of object
       continue
     else:
       d_xmin = dataframe.iloc[d_idx]['xmin']
       d_xmax = dataframe.iloc[d_idx]['xmax']
       d_ymin = dataframe.iloc[d_idx]['ymin']
       d_ymax = dataframe.iloc[d_idx]['ymax']
       d_poly = Polygon([(d_xmin,d_ymin), (d_xmax,d_ymin), (d_xmax,d_ymax), (d_xmin, d_ymax)])
  
       intersection = p_poly.intersection(d_poly)
       union = p_poly.union(d_poly)
       if intersection.area / union.area > 0.3:
         return True
         break
 
   return False

def return_sample_list(num, regions, dist, scale):
    """
    Return a list containing the position (x,y) of the samples.
    Output is a list of lists: the ith list contains samples for
    the ith region.
    """
    sample_list = []
    for index, region in enumerate(regions):
        min_x = min([region[i][0] for i in xrange(len(region))])
        max_x = max([region[i][0] for i in xrange(len(region))])
        min_y = min([region[i][1] for i in xrange(len(region))])
        max_y = max([region[i][1] for i in xrange(len(region))])
        poly = Polygon(region)
        region_samples = []
        while len(region_samples) != num:
            if dist == "normal":            
                candidate = (np.random.normal(loc = (min_x + max_x)/2, scale = scale), \
                np.random.normal(loc = (min_y + max_y)/2, scale = scale))            
            if dist == "uniform":           
                candidate = (np.random.uniform(min_x, max_x),np.random.uniform(min_y, max_y))                        
            if poly.contains(Point(candidate)):            
                region_samples.append(candidate)
        sample_list.append(region_samples)
 
    print "-done sampling-"
    return sample_list    

    def pack_shape_scale_linear(self):
        center = self.random_point()
        base = self.base_shapes[0]
        ph = np.random.random() * 2 * np.pi
        R = np.matrix([[np.cos(ph), -np.sin(ph)], [np.sin(ph), np.cos(ph)]])
        rbase = base * R

        # linear search on scale to find best fit
        r = 0
        delta = 2 ** -4
        while True:
            p = Polygon(r * rbase + center)
            intersected = False
            for shape in self.shapes:
                for poly in shape:
                    if p.intersects(poly):
                        intersected = True
                        break
                if intersected:
                    break

            # if any([p.intersects(poly) for poly in polys]):
            if intersected:
                break
            r += delta
        print('  %f' % r)

        self.shapes.append(p)

def intersectNodes(path, srs, projName, projSRS, isGridProject, gridResolution):
    j = []
    isGridLine = False
    sf = shapefile.Reader(path)
    nodes = list(getShapelyNodes(projName))
    for shape in sf.shapes():
        shType = shape.shapeType
        # http://en.wikipedia.org/wiki/Shapefile#Shapefile_shape_format_.28.shp.29
        if shType == 5: # Polygon
            sh = Polygon(shape.points)
        elif shType == 3: # Line
            if isGridProject:
                sh = LineString(shape.points)
                isGridLine = True
            else:
                pass
        else:
            consoleAppend('Unknown shape type %s. Continue without access' %shType)

        if srs != projSRS:
            sh = shapelyReproject(sh, srs, projSRS)

        if isGridLine:
            sh = sh.buffer(gridResolution)

        for node in nodes:
            if sh.contains(node[1]): # node.geom
                j.append(node[0]) # node.node_id

    return j if j else None

def compare_location_results(expected, found):
    true_positive = 0
    false_positive = 0
    false_negative = 0

    ambiguous = 0

    paired = [False]*len(found)

    for e in expected:
        p1=Polygon(reshape_list(e))
        total_matched = 0
        for idx,f in enumerate(found):
            p2=Polygon(reshape_list(f))
            try:
                x = p1.intersection(p2)
                if x.area/p1.area > 0.1:
                    paired[idx] = True
                    total_matched += 1
                    true_positive += 1
            except:
                pass # not sure what to do here
        if total_matched == 0:
            false_negative += 1
        elif total_matched > 1:
            ambiguous += 1

    for idx in range(len(found)):
        if not paired[idx]:
            false_positive += 1

    return {"tp":true_positive,"fp":false_positive,"fn":false_negative,"ambiguous":ambiguous}

    def get_reachable(self, coord, extra_coords):
        res = []
        for c in [(node.x, node.y) for node in self.nodes] + extra_coords:
            if c == coord:
                continue

            dirvec = (c[0]-coord[0], c[1]-coord[1])
            norm = (dirvec[0]**2 + dirvec[1]**2)**.5
            scl = self.uav_radius / norm
            norvecs = [(v[0]*scl, v[1]*scl) for v in [(-dirvec[1], dirvec[0]), (dirvec[1], -dirvec[0])]]
            corners = [
                (c[0]+norvecs[0][0], c[1]+norvecs[0][1]),
                (c[0]+norvecs[1][0], c[1]+norvecs[1][1]),
                (coord[0]+norvecs[0][0], coord[1]+norvecs[0][1]),
                (coord[0]+norvecs[1][0], coord[1]+norvecs[1][1])
            ]
            path_poly = Polygon(corners)

            # ls = LineString([coord, c])
            canReach = True
            for o in self.obstacles:
                # if ls.intersects(o):
                    # canReach = False
                if path_poly.intersects(o):
                    canReach = False
            #if not self.fly_zone.contains(ls):
                #canReach = False
            if canReach:
                res.append(c)
        return res

class Polygon(Geometry):
    """ class for convex polygon intersection test """
    def __init__(self, vertices):
        """ constructor for polygon, vertices must be specified
        in counter-clockwise order """

        # shapely should not be used if the students are implementing this
        self.poly = ShapelyPolygon(vertices)

    def intersects(self, geometry):
        if isinstance(geometry, Collection):
            return geometry.intersects(self)
        elif isinstance(geometry, Point):
            return self.point_poly_test(geometry)
        else:
            return self.poly_poly_test(geometry)

    def point_poly_test(self, p):
        """ This method should be implemented by the students but for
        demo purposes, shapely is used """
        return self.poly.intersects(p.point)

    def poly_poly_test(self, p):
        """ This method should be implemented by the students but for
        demo purposes, shapely is used """
        return self.poly.intersects(p.poly)

    @property
    def vertices(self):
        return list(self.poly.exterior.coords)

def pg_pix2latlon_strdf(im_name):
    # load the polygon data for each region
    poly_path = './shpres/%s.json'%(im_name)
    if not os.path.isfile(poly_path):
        #print 'no shape files', im_name
        exit()

    #print 'shift', shift_dict[im_name]
    rs, cs = shift_dict[im_name]

    pg = simplejson.load(open(poly_path))
    #print pg

    exter = np.array(pg['ext'])
    exter[:,0] += cs
    exter[:,1] += rs
    ex_lonlat = pix2ll(exter,t)

    inters = pg['intlist']
    
    inter_list = []
    for inter in inters:
        np_inter = np.array(inter)
        np_inter[:,0] += cs
        np_inter[:,1] += rs
        np_inter = pix2ll(np_inter,t)
        inter_list.append(np_inter)
    
    pg_obj = Polygon(ex_lonlat, inter_list)
    print pg_obj.contains( Point(69.178746, 35.813774) )

def box3d_intersection(box_a, box_b):
    """
    A simplified calculation of 3d bounding box intersection.
    It is assumed that the bounding box is only rotated
    around Z axis (yaw) from an axis-aligned box.
    :param box_a, box_b: obstacle bounding boxes for comparison
    :return: intersection volume (float)
    """
    # height (Z) overlap
    min_h_a = np.min(box_a[2])
    max_h_a = np.max(box_a[2])
    min_h_b = np.min(box_b[2])
    max_h_b = np.max(box_b[2])
    max_of_min = np.max([min_h_a, min_h_b])
    min_of_max = np.min([max_h_a, max_h_b])
    z_intersection = np.max([0, min_of_max - max_of_min])
    if z_intersection == 0:
        return 0.

    # oriented XY overlap
    xy_poly_a = Polygon(zip(*box_a[0:2, 0:4]))
    xy_poly_b = Polygon(zip(*box_b[0:2, 0:4]))
    xy_intersection = xy_poly_a.intersection(xy_poly_b).area
    if xy_intersection == 0:
        return 0.

    return z_intersection * xy_intersection

    def _center_pts(pts):
        '''Fancy label position generator, using erosion to get label coordinate'''
        min = pts.min(0)
        pts -= min
        max = pts.max(0)
        pts /= max

        #probably don't need more than 20 points, reduce detail of the polys
        if len(pts) > 20:
            pts = pts[::len(pts)//20]

        try:
            poly = Polygon([tuple(p) for p in pts])
            for i in np.linspace(0,1,100):
                if poly.buffer(-i).is_empty:
                    return list(poly.buffer(-last_i).centroid.coords)[0] * max + min
                last_i = i

            print("unable to find zero centroid...")
            return list(poly.buffer(-100).centroid.coords)[0] * max + min
        except:
            # This may not be worth being so verbose about... I think this is only for label positions.
            import warnings
            warnings.warn("Shapely error - computing mean of points instead of geometric center")
            return np.nanmean(pts, 0)

def fix_geometry(geometry):
    """Attempts to fix an invalid geometry (from https://goo.gl/nfivMh)"""
    try:
        return geometry.buffer(0)
    except ValueError:
        pass

    polygons = geom_as_list(geometry)

    fixed_polygons = list()
    for i, polygon in enumerate(polygons):
        if not linear_ring_is_valid(polygon.exterior):
            continue

        interiors = []
        for ring in polygon.interiors:
            if linear_ring_is_valid(ring):
                interiors.append(ring)

        fixed_polygon = Polygon(polygon.exterior, interiors)

        try:
            fixed_polygon = fixed_polygon.buffer(0)
        except ValueError:
            continue

        fixed_polygons.extend(geom_as_list(fixed_polygon))

    if len(fixed_polygons) > 0:
        return MultiPolygon(fixed_polygons)
    else:
        return None

def replace_lines_to_point(line_lyrs, old_point, new_point):
    x = old_point[0].coords[0][0]
    y = old_point[0].coords[0][1]
    buff = Polygon([(x-1, y-1), (x-1, y+1), (x+1, y+1), (x+1, y-1), (x-1, y-1)])
    buff.srid = old_point.srid

    for line_lyr in line_lyrs:
        # get intersection lines
        query = line_lyr.feature_query()
        #query.intersects(buff)
        query.geom()

        # replace point in lines
        for f in query():
            line = f.geom[0]
            new_line_points = []
            need_reconstruct = False
            for vertex in line.coords:
                if vertex == old_point[0].coords[0]:
                    new_line_points.append(new_point[0].coords[0])
                    need_reconstruct = True
                else:
                    new_line_points.append(vertex)

            if need_reconstruct:
                new_geom = MultiLineString([new_line_points,])
                f.geom = new_geom
                line_lyr.feature_put(f)

    def __init__(self, robot_name, data_polygons):
        self.robot_name = robot_name
        self.open_area = Polygon()
        self.full_area = Polygon()


        for id_robot, pol_data in data_polygons.items():
            # pol_data = [polygon, closed, time]

            try:
                available = not pol_data[1]
                pol = Polygon(pol_data[0])

                # Invalid polygons cannot be join.
                if not pol.is_valid:
                    pol = MultiPoint(pol_data[0]).convex_hull

                if available or robot_name == id_robot:

                    self.open_area = self.open_area.union(pol)
                else:
                    self.full_area = self.full_area.union(pol)

            except Exception:
                print "Error Joining polygons", SystemError.message

def parse_modis_coordinates(url_xml, coordinates, verbose):
    upperleft = (float(coordinates.split(',')[0]), float(coordinates.split(',')[1]))  # lat, lon
    downright = (float(coordinates.split(',')[2]), float(coordinates.split(',')[3]))  # lat, lon
    upperright = (upperleft[0], downright[1])
    downleft = (downright[0], upperleft[1])
    requested_bbox = Polygon((upperleft, upperright, downright, downleft))
    if verbose:
        LOG.info("UL: LAT -> %s, LON -> %s" % upperleft)
        LOG.info("DR: LAT -> %s, LON -> %s" % downright)
    req = urllib2.Request("%s" % url_xml, None, HEADERS)
    print(url_xml)
    root = etree.parse(urllib2.urlopen(req))
    bbox = []
    for point in root.xpath('/GranuleMetaDataFile/GranuleURMetaData/SpatialDomainContainer/'
                            'HorizontalSpatialDomainContainer/GPolygon/Boundary/Point'):
        lon = point.xpath('./PointLongitude')
        lat = point.xpath('./PointLatitude')
        bbox.append((float(lat[0].text), float(lon[0].text)))
    product_bbox = MultiPoint(bbox).convex_hull
    if verbose:
        for point in bbox:
            (lat, lon) = point
            LOG.info("Point: LAT -> %s LON -> %s" % (lat, lon))

    if requested_bbox.intersects(product_bbox):
        LOG.info("Compatible")
        return True
    else:
        LOG.info("Not Compatible")
        return False

    def test_attribute_chains(self):

        # Attribute Chaining
        # See also ticket #151.
        p = Polygon(((0.0, 0.0), (0.0, 1.0), (-1.0, 1.0), (-1.0, 0.0)))
        self.assertEqual(
            list(p.boundary.coords),
            [(0.0, 0.0), (0.0, 1.0), (-1.0, 1.0), (-1.0, 0.0), (0.0, 0.0)])

        ec = list(Point(0.0, 0.0).buffer(1.0, 1).exterior.coords)
        self.assertIsInstance(ec, list)  # TODO: this is a poor test

        # Test chained access to interiors
        p = Polygon(
            ((0.0, 0.0), (0.0, 1.0), (-1.0, 1.0), (-1.0, 0.0)),
            [((-0.25, 0.25), (-0.25, 0.75), (-0.75, 0.75), (-0.75, 0.25))]
        )
        self.assertEqual(p.area, 0.75)

        """Not so much testing the exact values here, which are the
        responsibility of the geometry engine (GEOS), but that we can get
        chain functions and properties using anonymous references.
        """
        self.assertEqual(
            list(p.interiors[0].coords),
            [(-0.25, 0.25), (-0.25, 0.75), (-0.75, 0.75), (-0.75, 0.25),
             (-0.25, 0.25)])
        xy = list(p.interiors[0].buffer(1).exterior.coords)[0]
        self.assertEqual(len(xy), 2)

        # Test multiple operators, boundary of a buffer
        ec = list(p.buffer(1).boundary.coords)
        self.assertIsInstance(ec, list)  # TODO: this is a poor test

def doPolygonize():
  blocks = polygonize(lines)
  writeBlocks(blocks, args[0] + '-blocks.geojson')

  blocks = polygonize(lines)
  bounds = Polygon([
    [minlng, minlat],
    [minlng, maxlat],
    [maxlng, maxlat],
    [maxlng, minlat],
    [minlng, minlat]
  ])
  # Geometry transform function based on pyproj.transform
  project = partial(
    pyproj.transform,
    pyproj.Proj(init='EPSG:3785'),
    pyproj.Proj(init='EPSG:4326'))
  print bounds
  print transform(project, bounds)

  print 'finding holes'
  for index, block in enumerate(blocks):
    if index % 1000 == 0:
      print "diff'd  %s" % (index)
    if not block.is_valid:
      print explain_validity(block)
      print transform(project, block)
    else:
      bounds = bounds.difference(block)
  print bounds

	def setStInShape(self,shpfile):
		"""
		Функция возвращает список станций попадающий в полигон(ы) из шэйпфайла файла
		"""
		import shapefile as shp
		import geocalc
		from shapely.geometry import Polygon,Point
		res=[]
		sf = shp.Reader(shpfile)
		for sp in sf.shapes():
			res_tmp=[]
			lonmin,latmin,lonmax,latmax=sp.bbox
			lonmin,lonmax=geocalc.cLon(lonmin),geocalc.cLon(lonmax)
			if lonmin<0 or lonmax<0:
				polygonPoints=[[geocalc.cLon(cors[0]),cors[1]] for cors in sp.points]
			else:
				polygonPoints=sp.points
			poly=Polygon(polygonPoints)
			indsInBox=[ind for ind in self.stInds if lonmin<=geocalc.cLon(self.stMeta[ind]['lon'])<=lonmax and latmin<=self.stMeta[ind]['lat']<=latmax]
			for ind in indsInBox:
				lat,lon=self.stMeta[ind]['lat'], geocalc.cLon(self.stMeta[ind]['lon'])
				pnt=Point(lon,lat)
				if poly.contains(pnt): res_tmp.append(ind)
			res=res+res_tmp
		return list(set(res))

    def partition(self):
        partition_by_station_dict = dict()
        population_by_station_dict = dict()
        station_coordinates = []
        for station in self.stations.values():
            station_coordinates.append([station.lon, station.lat])
        points = np.array(station_coordinates)
        partitions = Voronoi(points)
        regions, vertices = LoadEstimator.voronoi_finite_polygons_2d(partitions)

        polygons = []
        for region in regions:
            vertices_coordinates = []
            for vertice_index in region:
                vertices_coordinates.append((vertices[vertice_index][0], vertices[vertice_index][1]))
            vertices_coordinates.append(vertices_coordinates[0])
            partition_polygon = Polygon(vertices_coordinates)
            polygons.append(partition_polygon)
            for station in self.stations.values():
                if Point(station.lon, station.lat).within(partition_polygon):
                    partition_by_station_dict[str(station.id)] = partition_polygon.intersection(self.boundary)
                    population_of_region = self.population_of_region(partition_polygon)
                    self.root.info('Region of station %s has %s people', str(station.id), str(population_of_region))
                    population_by_station_dict[str(station.id)] = population_of_region
                    break
        return partition_by_station_dict, population_by_station_dict

def normalize_footprint(footprint: List[Tuple[float, float]]) -> MultiPolygon:
    """Split footprints which cross the anti-meridian

    Most applications, including RasterFoundry, cannot handle a single polygon representation
    of anti-meridian crossing footprint.
    Normalizing footprints by splitting them over the anti-meridian fixes cases where
    scenes appear to span all longitudes outside their actual footprint.
    If a footprint covers the anti-meridian, the shape is shifted 360 degrees, split,
    then the split section is moved back.
    """
    intersects = intersects_antimeridian(footprint)
    if not intersects:
        return MultiPolygon([Polygon(footprint)])
    else:
        shifted_footprint = Polygon([shift_point(p, 0, Side.RIGHT, False, 360) for p in footprint])
        left_hemisphere_mask = Polygon([(0, -90), (180, -90), (180, 90), (0, 90), (0, -90)])
        left_split = shifted_footprint.intersection(left_hemisphere_mask)
        right_split = Polygon([
            shift_point((x, y), 180, Side.LEFT, True, -360)
            for x, y
            in shifted_footprint.difference(left_hemisphere_mask).exterior.coords
        ])
        if left_split.area > 0 and right_split.area > 0:
            return MultiPolygon([left_split, right_split])
        elif left_split.area > 0:
            return MultiPolygon([left_split])
        elif right_split.area > 0:
            return MultiPolygon([right_split])
        else:
            return MultiPolygon([])