def test_line_intersection(self):
     line1 = LineString([(0, 0, 0), (1, 1, 1)])
     line2 = LineString([(0, 1, 1), (1, 0, 0)])
     interxn = line1.intersection(line2)
     self.assertTrue(interxn.has_z)
     self.assertEqual(interxn._ndim, 3)
     self.assertTrue(0.0 <= interxn.z <= 1.0)

 def offset(self, distance):
     self.points.append(self.points[0])
     line = LineString(self.getSequence())
     offset = line.parallel_offset(distance, 'right', join_style=1)
     # return list(offset.coords)
     self.setSequence(list(offset.coords))
     return self

def clipLine(lineGeometry, pt1, pt2):
    """Returns a line cliipped to the extent of two given points"""
    # Assumes pt1, pt2 lie on line
    if lineGeometry is None or lineGeometry.isEmpty() or pt1 is None or pt2 is None:
        return QgsGeometry()
    line = LineString(lineGeometry.geometry())
    d1 = line.project(Point(pt1))
    d2 = line.project(Point(pt2))
    if d1 < d2:
        start = pt1
        ds = d1
        end = pt2
        de = d2
    else:
        start = pt2
        ds = d2
        end = pt1
        de = d1
    clip = []
    clip.append(start)
    for coord in line.coords:
        pt = Point(coord)
        dp = line.project(pt)
        if dp > ds and dp < de:
            clip.append(QgsPointV2(pt.x, pt.y))
    clip.append(end)
    return QgsGeometry.fromPolyline(clip)

def acc2latlng(caseno):
    """Calculate the latitude/longitude and update into mongodb

    :param caseno: case number of the accident
    """
    # get the cnty_rte and mile post
    acc = acc_col.find_one({'caseno':caseno})

    # check exists
    if 'lat' in acc.keys() and 'lng' in acc.keys():
        return

    rte_nbr = acc['rte_nbr']
    cnty_rte = acc['cnty_rte']
    milepost = acc['milepost']
    if cnty_rte.index(rte_nbr) != 2:
        raise Exception("Unknown format")

    # get the road information
    rid = rte_nbr + cnty_rte[0:2]
    r_info = cnty_rte2linestring(rid)
    max_mp = r_info['rlist'][-1]['endmp']
    ls = LineString( r_info['plist'] )
    deg_mp = milepost / max_mp * ls.length # mile to degree
    acc_pos = ls.interpolate( deg_mp )
    lng = acc_pos.coords[0][0]
    lat = acc_pos.coords[0][1]
    #print ls
    #print acc_pos
    # TODO, update the lat lng information
    acc_col.update(
            {'caseno':caseno},
            {'$set':{'lat':lat, 'lng':lng}}
            )

    def find_common_path(self, current_adc_trajectory, last_adc_trajectory):
        current_path_points = current_adc_trajectory.trajectory_point
        last_path_points = last_adc_trajectory.trajectory_point
        current_path = []
        for point in current_path_points:
            current_path.append([point.path_point.x, point.path_point.y])
        last_path = []
        for point in last_path_points:
            last_path.append([point.path_point.x, point.path_point.y])
        if len(current_path) == 0 or len(last_path) == 0:
            return [], []

        current_ls = LineString(current_path)
        last_ls = LineString(last_path)
        current_start_point = Point(current_path[0])

        dist = last_ls.project(current_start_point)
        cut_lines = self.cut(last_ls, dist)
        if len(cut_lines) == 1:
            return [], []
        last_ls = cut_lines[1]
        dist = current_ls.project(Point(last_path[-1]))
        if dist <= current_ls.length:
            current_ls = self.cut(current_ls, dist)[0]
        else:
            dist = last_ls.project(Point(current_path[-1]))
            last_ls = self.cut(last_ls, dist)[0]
        return current_ls.coords, last_ls.coords

    def _add_slope(self, bounds, altitude1, altitude2, point1, point2, bottom=False):
        altitude_diff = altitude2-altitude1
        altitude_middle = (altitude1+altitude2)/2
        altitude_halfdiff = altitude_diff/2
        altitude_base = altitude1
        line = LineString([point1, point2])

        minx, miny, maxx, maxy = bounds
        points_2d = [(minx-100, miny-100), (maxx+100, miny-100), (maxx+100, maxy+100), (minx-100, maxy+100)]
        points_3d = []
        for i, (x, y) in enumerate(points_2d):
            point = Point((x, y))
            pos = line.project(point)
            while pos <= 0 or pos >= line.length-1:
                line = scale(line, xfact=2, yfact=2, zfact=2)
                altitude_diff *= 2
                altitude_halfdiff *= 2
                altitude_base = altitude_middle-altitude_halfdiff
                pos = line.project(point)
            z = ((pos/line.length)*altitude_diff)+altitude_base
            points_3d.append((x, y, z/1000))

        extrude = abs(altitude1-altitude2)/1000+100
        if bottom:
            extrude = -extrude
        self._add_python(
            'last_slope = add_polygon_3d(name=%(name)r, coords=%(coords)r, extrude=%(extrude)f)' % {
                'name': 'tmpslope',
                'coords': tuple(points_3d),
                'extrude': extrude,
            }
        )

    def _add_crossroads(new_edges, current_edges):
        new_new_edges = set()
        new_current_edges = set()
        split_current_edge = [False] * len(current_edges)

        for new_edge in new_edges:
            found_crossroad = False
            for current_edge_index, current_edge in enumerate(current_edges):
                if State._different_starts_and_ends(new_edge, current_edge):
                    new_line_segment = LineString([(new_edge[0].x, new_edge[0].y), (new_edge[1].x, new_edge[1].y)])
                    current_line_segment = LineString(([(current_edge[0].x, current_edge[0].y),
                                                        (current_edge[1].x, current_edge[1].y)]))

                    intersection = new_line_segment.intersection(current_line_segment)
                    if type(intersection) is ShapelyPoint:
                        new_crossroad = Point(intersection.x, intersection.y, Point.TYPE_CROSSROAD)
                        new_current_edges.update(State._split_edge(current_edge, new_crossroad))
                        new_new_edges.update(State._split_edge(new_edge, new_crossroad))
                        split_current_edge[current_edge_index] = True
                        found_crossroad = True
                        break
            if not found_crossroad:
                new_new_edges.add(new_edge)

        for current_edge_index, current_edge in enumerate(current_edges):
            if not split_current_edge[current_edge_index]:
                new_current_edges.add(current_edge)

        if set(current_edges) == new_current_edges and set(new_edges) == new_new_edges:
            return list(current_edges) + list(new_edges)
        else:
            return State._add_crossroads(list(new_new_edges), list(new_current_edges))

def append_to_included_groups(locs, d):
    for group_key in d.keys():
        ref_gene_poly = LineString([(0.0, float(group_key[0])),(0.0, float(group_key[1]))])
        locs_poly = LineString([(0, locs[2]), (0, locs[3])])
        if ref_gene_poly.intersects(locs_poly):
            d[group_key].append(tuple(locs))
    return d

def cut_line(cut_point, line, eps_mult=1e2):
    dist = line.project(cut_point)
    point = line.interpolate(dist)
    eps = line.distance(point) * eps_mult

    coords = list(line.coords)

    if point.coords[0] in coords:
        i = coords.index(point.coords[0])

        if i == 0:
            return LineString(), line
        if i == len(coords) - 1:
            return line, LineString()

        start_segment = LineString(coords[:i + 1])
        end_segment = LineString(coords[i:])

        return start_segment, end_segment


    for i, p in enumerate(coords[:-1]):
        line_segment = LineString([coords[i], coords[i + 1]])
        line_segment_buffer = line_segment.buffer(eps, resolution=1)

        if line_segment_buffer.contains(point):
            start_segment = LineString(coords[:i + 1] + [point])
            end_segment = LineString([point] + coords[i + 1:])

            return start_segment, end_segment

    raise Exception('point not found in line, consider raising eps_mult')

def getShortestDisLinePoint(pointLon, pointLat, lineCoords, earthR):
	#get shortest distance between point and line
	line = LineString(lineCoords) 
	point = Point(pointLon, pointLat)
	closestCoord = line.interpolate(line.project(point))
	distance = haversine(pointLon, pointLat, closestCoord.x, closestCoord.y, earthR)
	return distance

def identify_first_point_in_polygon(points, ddd=0.5):
    """
    Identify the cycle beginning for the last point
    :param points: array of objects of PointX
    :param ddd: distance to the fist point ddd>0.
    :return -1 if do not close, or i if close in point i.
    """
    if len(points) < 3:
        return -1

    # last point
    lp = Point(points[-1])

    first_point = len(points) - 2

    # Distance between last point and a line segment
    distance_point_line = 2 * ddd

    # Distance is not less than ddd or less than 3 points.
    while not (distance_point_line < ddd and polyline_length(points[first_point:]) > 2):
        # next segment
        first_point -= 1

        if first_point == -1:
            break

        pt1 = points[first_point + 1]
        pt2 = points[first_point]

        # Distance to the first point to the line segment
        line_segment = LineString([pt1, pt2])
        distance_point_line = line_segment.distance(lp)

    return first_point

    def divide_polygon_for_intersection(self, segments):
        """ Generates multiple polygons based on cutting the 
        fracture faces by line segments. 
        """
        R = self.build_rotation_matrix()
        fracture_poly_list = []
        # face_poly_list = []

        for point in self.points:
            rot_point = np.linalg.solve(R, point - self.center)
            fracture_poly_list.append(rot_point[:2])

        seg_rot = []
        for seg in segments:
            p1 = seg[0]
            p2 = seg[1]

            vec = p2 - p1
            p1 = p1 - 100.0 * vec
            p2 = p2 + 100.0 * vec

            p1_rot = np.linalg.solve(R, p1 - self.center)
            p2_rot = np.linalg.solve(R, p2 - self.center)

            line = LineString((p1_rot[:2], p2_rot[:2]))
            dilated = line.buffer(1.0e-10)
            seg_rot.append(dilated)

        fracture_poly = Polygon(fracture_poly_list)

        divided_polygons = fracture_poly.difference(seg_rot[0])

        return (divided_polygons, fracture_poly)

 def test_line_intersection(self):
     line1 = LineString([(0, 0, 0), (1, 1, 1)])
     line2 = LineString([(0, 1, 1), (1, 0, 0)])
     interxn = line1.intersection(line2)
     self.failUnless(interxn.has_z)
     self.failUnless(interxn._ndim == 3)
     self.failUnless(0.0 <= interxn.z <= 1.0)

 def setUp(self):
     self.point = Point(1, 1)
     self.line1 = LineString(([0, 0], [2, 0]))
     self.line2 = LineString(([3, 0], [3, 6]))
     self.multiline = MultiLineString([
         list(self.line1.coords), list(self.line2.coords)
     ])

def lines_fusion(line1, line2):
    """
    Validate each line segment for intersection. Intersection is checked as
    http://en.wikipedia.org/wiki/Line_segment_intersection

    :param line1: main line
    :param line2: line to be fused
    :return same line1 if there is no intersections with line2. else fused lines.
    """
    # TODO use bounding boxes to speed up
    # Cross validation for each line segment in polylines.
    for i in range(len(line1) - 1, 0, -1):
        p1 = line1[i]
        p2 = line1[i - 1]
        # line segment
        l1 = LineString([p1, p2])

        for j in range(len(line2) - 1, 0, -1):
            p3 = line2[j]
            p4 = line2[j - 1]

            l2 = LineString([p3, p4])
            # Check Line intersection
            intersection = l1.intersection(l2)

            # if intersected
            if not intersection.is_empty:
                # take the line 1 from first point until the intersection and line2
                inter_p = (intersection.coords.xy[0][0], intersection.coords.xy[1][0])
                new_line = line2[:j] + [inter_p] + line1[i:]
                return new_line

    # no fusion
    return line1

    def create_centerline(self):
        """
        Calculates the centerline of a polygon.
        Densifies the border of a polygon which is then represented
        by a Numpy array of points necessary for creating the
        Voronoi diagram. Once the diagram is created, the ridges
        located within the polygon are joined and returned.
        Returns:
            a MultiLinestring located within the polygon.
        """
        minx = int(min(self.inputGEOM.envelope.exterior.xy[0]))
        miny = int(min(self.inputGEOM.envelope.exterior.xy[1]))
        border = np.array(self.densify_border(self.inputGEOM, minx, miny))
        #print(border, minx, miny, self.dist)
        vor = Voronoi(border)
        vertex = vor.vertices

        lst_lines = []
        for j, ridge in enumerate(vor.ridge_vertices):
            if -1 not in ridge:
                line = LineString([
                    (vertex[ridge[0]][0] + minx, vertex[ridge[0]][1] + miny),
                    (vertex[ridge[1]][0] + minx, vertex[ridge[1]][1] + miny)])
                if line.within(self.inputGEOM) and len(line.coords[0]) > 1:
                    lst_lines.append(line)
        return MultiLineString(lst_lines)

    def test_linestring_geojson(self):
        '''Create a line that goes from west to east (clip on)'''
        
        self.defineGeometry('LINESTRING')

        geom = LineString( [(-180, 32), (180, 32)] )

        self.insertTestRow(geom.wkt)

        # we should have a line that clips at 0...

        # for western hemisphere....
        tile_mimetype, tile_content = utils.request(self.config_file_content, "vector_test", "geojson", 0, 0, 0)
        self.assertEqual(tile_mimetype, "text/json")
        geojson_result = json.loads(tile_content)
        west_hemisphere_geometry = asShape(geojson_result['features'][0]['geometry'])
        expected_geometry = LineString([(-180, 32), (0, 32)])
        self.assertTrue(expected_geometry.almost_equals(west_hemisphere_geometry))

        # for eastern hemisphere....
        tile_mimetype, tile_content = utils.request(self.config_file_content, "vector_test", "geojson", 0, 1, 0)
        self.assertEqual(tile_mimetype, "text/json")
        geojson_result = json.loads(tile_content)
        east_hemisphere_geometry = asShape(geojson_result['features'][0]['geometry'])
        expected_geometry = LineString([(0, 32), (180, 32)])
        self.assertTrue(expected_geometry.almost_equals(east_hemisphere_geometry))

def fitline_plot(df):
	## best fit straight line of points
	slope, intercept = hcf.fitline(df)

	## define Shapely linestring using westmost and eastmost endpoints
	point1 = Point(df['x_working'].min(), slope * df['x_working'].min() + intercept)
	point2 = Point(df['x_working'].max(), slope * df['x_working'].max() + intercept)
	ls = LineString([point1, point2])

	## iterate through dfframe, project each TS point on Shapely line.  Determine position on line and absolute x,y value
	for index, row in df.iterrows():
		pointN = Point(row['x_working'], row['y_working'])
		projectN = ls.project(pointN)
		df.set_value(index, 'position_on_line', projectN)
		globalN = ls.interpolate(projectN)
		df.set_value(index, 'x_project', globalN.x)
		df.set_value(index, 'y_project', globalN.y)
		offsetN = globalN.distance(pointN)
		df.set_value(index, 'offset_from_line', offsetN)

	df = df.sort_values(by='position_on_line', ascending='true')
	print df

	## plot points projected on best fit line
	plt.scatter(df['x_project'], df['y_project'], color='pink')

	## init Figure 2
	plt.figure(num=2, figsize=(13, 3), dpi=80)
	plt.grid(True)
	# plt.axes().set_aspect('equal', 'dflim')
	plt.axis([df['position_on_line'].min() - 20, df['position_on_line'].max() + 20, df['z_working'].min() - 4,
	          df['z_working'].max() + 4])
	plt.scatter(df['position_on_line'], df['z_working'], color='black')
	plt.plot(df['position_on_line'], df['z_working'], color='black')

def pymol_select_memb_old(pdb: MyPDB) -> set():
    """
    print a pymol selection line for all residues that are in the membrane
    !!! this assumes that cntr is at 0 0 0 and norm at 15 0 0 !!!
    """
    from shapely.geometry import LineString, Point
    # create Points from center & thickness
    cntr_pnt = Point(pdb.memb_res.cntr.x, pdb.memb_res.cntr.y, pdb.memb_res.cntr.z)
    thkn_m_pnt = Point(-pdb.memb_res.thkn.x, pdb.memb_res.thkn.y, pdb.memb_res.thkn.z)
    thkn_pnt = Point(pdb.memb_res.thkn.x, pdb.memb_res.thkn.y, pdb.memb_res.thkn.z)

    # define the line between center and thickness
    line = LineString([thkn_m_pnt, thkn_pnt])
    thickness = cntr_pnt.distance(thkn_pnt)

    result = set()
    # iterate over all CAs in the pdb
    for cid in sorted(pdb.chains.keys()):
        for rid in sorted(pdb[cid].residues.keys()):
            atom = pdb[cid][rid]['CA']

            # the atom as a Point
            p = Point(atom.xyz.x, atom.xyz.y, atom.xyz.z)

            # projection of the CA atom on the center-thickness line
            np = line.interpolate(line.project(p))

            # if the distance on the center-thickness line is smaller than 15, than this is in the membrane
            if cntr_pnt.distance(np) < thickness-0.1:
                result.add(atom.res_seq_num)

    return result

 def appendNewWay(coords, intersects, osmXml):
     way = etree.Element('way', visible='true', id=str(newOsmId('way')))
     firstNid = 0
     for i, coord in enumerate(coords):
         if i == 0: continue # the first and last coordinate are the same
         node = appendNewNode(coord, osmXml)
         if i == 1: firstNid = node.get('id')
         way.append(etree.Element('nd', ref=node.get('id')))
         
         # Check each way segment for intersecting nodes
         int_nodes = {}
         try:
             line = LineString([coord, coords[i+1]])
         except IndexError:
             line = LineString([coord, coords[1]])
         for idx, c in enumerate(intersects):
             if line.buffer(0.0000015).contains(Point(c[0], c[1])) and c not in coords:
                 t_node = appendNewNode(c, osmXml)
                 for n in way.iter('nd'):
                     if n.get('ref') == t_node.get('id'):
                         break
                 else:
                     int_nodes[t_node.get('id')] = Point(c).distance(Point(coord))
         for n in sorted(int_nodes, key=lambda key: int_nodes[key]): # add intersecting nodes in order
             way.append(etree.Element('nd', ref=n))
         
     way.append(etree.Element('nd', ref=firstNid)) # close way
     osmXml.append(way)
     return way