def clip_lines_by_polygon(lines, polygon,
                          closed=True,
                          check_input=True):
    """Clip multiple lines by polygon

    Input:
       lines: Sequence of polylines: [[p0, p1, ...], [q0, q1, ...], ...]
              where pi and qi are point coordinates (x, y).

       polygon: list or Nx2 array of polygon vertices

       closed: (optional) determine whether points on boundary should be
       regarded as belonging to the polygon (closed = True)
       or not (closed = False). False is not recommended here.
       This parameter can also be None in which case it is undefined
       whether a line on the boundary will fall inside or outside.
       This will make the algorithm about 20% faster.

       check_input: Allows faster execution if set to False

    Output:
       inside_lines: Dictionary of lines that are inside polygon
       outside_lines: Dictionary of lines that are outside polygon

       Elements in output dictionaries can be a list of multiple lines.
       One line is a numpy array of vertices.

       Both output dictionaries use the indices of the original line as keys.
       This makes it possible to track which line the new clipped lines
       come from, if one e.g. wants to assign the original attribute values
       to clipped lines.

    This is a wrapper around clip_line_by_polygon
    """

    if check_input:
        # Input checks
        msg = 'Keyword argument "closed" must be boolean'
        if not isinstance(closed, bool):
            raise RuntimeError(msg)

        for i in range(len(lines)):
            try:
                lines[i] = ensure_numeric(lines[i], numpy.float)
            except Exception, e:
                msg = ('Line could not be converted to numeric array: %s'
                       % str(e))
                raise Exception(msg)

            msg = ('Lines must be 2d array of vertices: '
                   'I got %d dimensions' % len(lines[i].shape))
            if not len(lines[i].shape) == 2:
                raise RuntimeError(msg)

        try:
            polygon = ensure_numeric(polygon, numpy.float)
        except Exception, e:
            msg = ('Polygon could not be converted to numeric array: %s'
                   % str(e))
            raise Exception(msg)

def rings_equal(x, y, rtol=1.0e-6, atol=1.0e-8):
    """Compares to linear rings as numpy arrays

    Args
        * x, y: Nx2 numpy arrays

    Returns:
        * True if x == y or x' == y (up to the specified tolerance)

        where x' is x reversed in the first dimension. This corresponds to
        linear rings being seen as equal irrespective of whether they are
        organised in clock wise or counter clock wise order
    """

    x = ensure_numeric(x, numpy.float)
    y = ensure_numeric(y, numpy.float)

    msg = 'Arrays must a 2d arrays of vertices. I got %s and %s' % (x, y)
    verify(len(x.shape) == 2 and len(y.shape) == 2, msg)

    msg = 'Arrays must have two columns. I got %s and %s' % (x, y)
    verify(x.shape[1] == 2 and y.shape[1] == 2, msg)

    if (numpy.allclose(x, y, rtol=rtol, atol=atol) or
        numpy.allclose(x, y[::-1], rtol=rtol, atol=atol)):
        return True
    else:
        return False

def array_to_line(A, geometry_type=ogr.wkbLinearRing):
    """Convert coordinates to linear_ring

    :param A: Nx2 Array of coordinates representing either a polygon or a line.
        A can be either a numpy array or a list of coordinates.
    :type A: numpy.ndarray, list

    :param geometry_type: A valid OGR geometry type.
        Default type ogr.wkbLinearRing
    :type geometry_type: ogr.wkbLinearRing, include ogr.wkbLineString

    Returns:
        * ring: OGR line geometry

    Note:
    Based on http://www.packtpub.com/article/working-geospatial-data-python
    """

    try:
        A = ensure_numeric(A, numpy.float)
    except Exception, e:
        msg = ('Array (%s) could not be converted to numeric array. '
               'I got type %s. Error message: %s'
               % (A, str(type(A)), e))
        raise Exception(msg)

def array2wkt(A, geom_type="POLYGON"):
    """Convert coordinates to wkt format

    Args:
        * A: Nx2 Array of coordinates representing either a polygon or a line.
             A can be either a numpy array or a list of coordinates.
        * geom_type: Determines output keyword 'POLYGON' or 'LINESTRING'

    Returns:
        * wkt: geometry in the format known to ogr: Examples

    Note:
        POLYGON((1020 1030,1020 1045,1050 1045,1050 1030,1020 1030))
        LINESTRING(1000 1000, 1100 1050)

    """

    try:
        A = ensure_numeric(A, numpy.float)
    except Exception, e:
        msg = "Array (%s) could not be converted to numeric array. " "I got type %s. Error message: %s" % (
            geom_type,
            str(type(A)),
            e,
        )
        raise Exception(msg)

def clip_lines_by_polygons(lines, polygons, check_input=True, closed=True):
    """Clip multiple lines by multiple polygons

    Args:
        * lines: Sequence of polylines: [[p0, p1, ...], [q0, q1, ...], ...]
            where pi and qi are point coordinates (x, y).
        * polygons: list of polygons, each an array of vertices
        * closed: optional parameter to determine whether lines that fall on
            an polygon boundary should be considered to be inside
            (closed=True), outside (closed=False) or
            undetermined (closed=None). The latter case will speed the
            algorithm up but lines on boundaries may or may not be
            deemed to fall inside the polygon and so will be
            indeterministic.

    Returns:
        lines_covered: List of polylines inside a polygon -o ne per input
        polygon.


    .. note:: If multiple polygons overlap, the one first encountered will be
        used.
    """

    if check_input:
        for i in range(len(lines)):
            try:
                lines[i] = ensure_numeric(lines[i], numpy.float)
            except Exception, e:
                msg = ('Line could not be converted to numeric array: %s'
                       % str(e))
                raise Exception(msg)

            msg = 'Lines must be 2d array of vertices'
            if not len(lines[i].shape) == 2:
                raise RuntimeError(msg)

        for i in range(len(polygons)):
            try:
                polygons[i] = ensure_numeric(polygons[i], numpy.float)
            except Exception, e:
                msg = ('Polygon could not be converted to numeric array: %s'
                       % str(e))
                raise Exception(msg)

def populate_polygon(polygon, number_of_points, seed=None, exclude=None):
    """Populate given polygon with uniformly distributed points.

    Input:
       polygon - list of vertices of polygon

       number_of_points - (optional) number of points

       seed - seed for random number generator (default=None)

       exclude - list of polygons (inside main polygon) from where points
       should be excluded

    Output:
       points - list of points inside polygon

    Examples:
       populate_polygon( [[0,0], [1,0], [1,1], [0,1]], 5 )
       will return five randomly selected points inside the unit square
    """

    polygon = ensure_numeric(polygon)

    # Find outer extent of polygon
    minpx = min(polygon[:, 0])
    maxpx = max(polygon[:, 0])
    minpy = min(polygon[:, 1])
    maxpy = max(polygon[:, 1])

    # Generate random points until enough are in polygon
    seed_function(seed)
    points = []
    while len(points) < number_of_points:
        x = uniform(minpx, maxpx)
        y = uniform(minpy, maxpy)

        append = False
        if is_inside_polygon([x, y], polygon):
            append = True

            #Check exclusions
            if exclude is not None:
                for ex_poly in exclude:
                    if is_inside_polygon([x, y], ex_poly):
                        append = False

        if append is True:
            points.append([x, y])

    return points

def polygon2segments(polygon):
    """Convert polygon to segments structure suitable for use in intersection

    Args:
        polygon: Nx2 array of polygon vertices

    Returns:
        A collection of line segments (x0, y0) -> (x1, y1) vectorised
        following the format
               line[0, 0, :] = x0
               line[0, 1, :] = y0
               line[1, 0, :] = x1
               line[1, 1, :] = y1

    """

    try:
        polygon = ensure_numeric(polygon, numpy.float)
    except Exception, e:
        msg = "Polygon could not be converted to numeric array: %s" % str(e)
        raise Exception(msg)

def array_to_wkt(A, geom_type='POLYGON'):
    """Convert coordinates to wkt format

    :param A: Nx2 Array of coordinates representing either a polygon or a line.
        A can be either a numpy array or a list of coordinates.
    :type A: numpy.array

    :param geom_type: Determines output keyword 'POLYGON' or 'LINESTRING'
    :type geom_type: str

    :returns: wkt: geometry in the format known to ogr: Examples

    Note:
        POLYGON((1020 1030,1020 1045,1050 1045,1050 1030,1020 1030))
        LINESTRING(1000 1000, 1100 1050)
    """

    try:
        A = ensure_numeric(A, numpy.float)
    except Exception, e:
        msg = ('Array (%s) could not be converted to numeric array. '
               'I got type %s. Error message: %s'
               % (geom_type, str(type(A)), e))
        raise Exception(msg)

def clip_line_by_polygon(line, polygon,
                         closed=True,
                         polygon_bbox=None,
                         check_input=True):
    """Clip line segments by polygon

    Input
       line: Sequence of line nodes: [[x0, y0], [x1, y1], ...] or
             the equivalent Nx2 numpy array
       polygon: list or Nx2 array of polygon vertices
       closed: (optional) determine whether points on boundary should be
       regarded as belonging to the polygon (closed = True)
       or not (closed = False) - False is not recommended here
       polygon_bbox: Provide bounding box around polygon if known.
           This is a small optimisation
       check_input: Allows faster execution if set to False

    Outputs
       inside_lines: Clipped lines that are inside polygon
       outside_lines: Clipped lines that are outside polygon

       Both outputs lines take the form of lists of Nx2 numpy arrays,
       i.e. each line is an array of multiple segments

    Example:

        U = [[0,0], [1,0], [1,1], [0,1]]  # Unit square

        # Simple horizontal fully intersecting line
        line = [[-1, 0.5], [2, 0.5]]

        inside_line_segments, outside_line_segments = \
            clip_line_by_polygon(line, polygon)

        print numpy.allclose(inside_line_segments,
                              [[[0, 0.5], [1, 0.5]]])

        print numpy.allclose(outside_line_segments,
                              [[[-1, 0.5], [0, 0.5]],
                               [[1, 0.5], [2, 0.5]]])

    Remarks:
       The assumptions listed in separate_points_by_polygon apply

       Output line segments are listed as separate lines i.e. not joined
    """

    if check_input:
        # Input checks
        msg = 'Keyword argument "closed" must be boolean'
        if not isinstance(closed, bool):
            raise RuntimeError(msg)

        try:
            line = ensure_numeric(line, numpy.float)
        except Exception, e:
            msg = ('Line could not be converted to numeric array: %s'
                   % str(e))
            raise Exception(msg)

        msg = 'Line segment array must be a 2d array of vertices'
        if not len(line.shape) == 2:
            raise RuntimeError(msg)

        msg = 'Line array must have two columns'
        if not line.shape[1] == 2:
            raise RuntimeError(msg)

        try:
            polygon = ensure_numeric(polygon, numpy.float)
        except Exception, e:
            msg = ('Polygon could not be converted to numeric array: %s'
                   % str(e))
            raise Exception(msg)

def point_on_line(points, line, rtol=1.0e-5, atol=1.0e-8,
                  check_input=True):
    """Determine if a point is on a line segment

    Input
        points: Coordinates of either
                * one point given by sequence [x, y]
                * multiple points given by list of points or Nx2 array
        line: Endpoint coordinates [[x0, y0], [x1, y1]] or
              the equivalent 2x2 numeric array with each row corresponding
              to a point.
        rtol: Relative error for how close a point must be to be accepted
        atol: Absolute error for how close a point must be to be accepted

    Output
        True or False

    Notes

    Line can be degenerate and function still works to discern coinciding
    points from non-coinciding.

    Tolerances rtol and atol are used with numpy.allclose()
    """

    one_point = False
    if check_input:
        # Prepare input data
        points = ensure_numeric(points)
        line = ensure_numeric(line)

        if len(points.shape) == 1:
            # One point only - make into 1 x 2 array
            points = points[numpy.newaxis, :]
            one_point = True
        else:
            one_point = False

        msg = 'Argument points must be either [x, y] or an Nx2 array of points'
        if len(points.shape) != 2:
            raise Exception(msg)
        if not points.shape[0] > 0:
            raise Exception(msg)
        if points.shape[1] != 2:
            raise Exception(msg)

    N = points.shape[0]  # Number of points

    x = points[:, 0]
    y = points[:, 1]
    x0, y0 = line[0]
    x1, y1 = line[1]

    # Vector from beginning of line to point
    a0 = x - x0
    a1 = y - y0

    # It's normal vector
    a_normal0 = a1
    a_normal1 = -a0

    # Vector parallel to line
    b0 = x1 - x0
    b1 = y1 - y0

    # Dot product
    nominator = abs(a_normal0 * b0 + a_normal1 * b1)
    denominator = b0 * b0 + b1 * b1

    # Determine if point vector is parallel to line up to a tolerance
    is_parallel = numpy.zeros(N, dtype=numpy.bool)  # All False
    is_parallel[nominator <= atol + rtol * denominator] = True

    # Determine for points parallel to line if they are within end points
    a0p = a0[is_parallel]
    a1p = a1[is_parallel]

    len_a = numpy.sqrt(a0p * a0p + a1p * a1p)
    len_b = numpy.sqrt(b0 * b0 + b1 * b1)
    cross = a0p * b0 + a1p * b1

    # Initialise result to all False
    result = numpy.zeros(N, dtype=numpy.bool)

    # Result is True only if a0 * b0 + a1 * b1 >= 0 and len_a <= len_b
    result[is_parallel] = (cross >= 0) * (len_a <= len_b)

    # Return either boolean scalar or boolean vector
    if one_point:
        return result[0]
    else:
        return result

def separate_points_by_polygon(points, polygon,
                               polygon_bbox=None,
                               closed=True,
                               check_input=True,
                               use_numpy=True):
    """Determine whether points are inside or outside a polygon.

    Args:
        * points: Tuple of (x, y) coordinates, or list of tuples
        * polygon: list or Nx2 array of polygon vertices
        * polygon_bbox: (optional) bounding box for polygon
        * closed: (optional) determine whether points on boundary should be
              regarded as belonging to the polygon (closed = True)
              or not (closed = False). If None, boundary is left undefined,
              i.e. some points on boundary may be deemed to be inside while
              others may be deemed to be outside. This options makes
              the code faster.
        * check_input: Allows faster execution if set to False
        * use_numpy: Use the fast numpy implementation

    Returns:
        * indices_inside_polygon: array of indices of points
              falling inside the polygon
        * indices_outside_polygon: array of indices of points
              falling outside the polygon

    Raises: A generic Exception is raised for unexpected input.

    Example:

        U = [[0,0], [1,0], [1,1], [0,1]]  # Unit square
        separate_points_by_polygon( [[0.5, 0.5], [1, -0.5], [0.3, 0.2]], U)

        will return the indices [0, 2, 1] and count == 2 as only the first
        and the last point are inside the unit square

    Remarks:
        The vertices may be listed clockwise or counterclockwise and
        the first point may optionally be repeated.
        Polygons do not need to be convex.
        Polygons can have holes in them and points inside a hole is
        regarded as being outside the polygon.

    Algorithm is based on work by Darel Finley,
    http://www.alienryderflex.com/polygon/
    """

    # FIXME (Ole): Make sure bounding box here follows same format as
    #              those returned by layers. Methinks they don't at the moment
    if check_input:
        # Input checks
        msg = 'Keyword argument "closed" must be boolean or None'
        if not (isinstance(closed, bool) or closed is None):
            raise PolygonInputError(msg)

        try:
            points = ensure_numeric(points, numpy.float)
        except Exception, e:
            msg = ('Points could not be converted to numeric array: %s'
                   % str(e))
            raise PolygonInputError(msg)

        try:
            polygon = ensure_numeric(polygon, numpy.float)
        except Exception, e:
            msg = ('Polygon could not be converted to numeric array: %s'
                   % str(e))
            raise PolygonInputError(msg)

def intersection(line0, line1):
    """Returns intersecting point between two line segments.

    If the lines are parallel or coincide partly (i.e. share a common segment),
    they are considered to not intersect.

    Inputs:
        line0: A simple line segment defined by two end points:
              [[x0, y0], [x1, y1]]

        line1: A collection of line segments vectorised following the format
               line[0, 0, :] = x2
               line[0, 1, :] = y2
               line[1, 0, :] = x3
               line[1, 1, :] = y3

    Output:
        intersections: Nx2 array with intersection points or nan
                       (in case of no intersection)
                       If line1 consisted of just one line,
                       None is returned for backwards compatibility


    Notes

    A vectorised input line can be constructed either as list:
    line1 = [[[0, 24, 0, 15],    # x2
              [12, 0, 24, 0]],   # y2
             [[24, 0, 0, 5],     # x3
              [0, 12, 12, 15]]]  # y3

    or as an array

    line1 = numpy.zeros(16).reshape(2, 2, 4)  # Four segments
    line1[0, 0, :] = [0, 24, 0, 15]   # x2
    line1[0, 1, :] = [12, 0, 24, 0]   # y2
    line1[1, 0, :] = [24, 0, 0, 5]    # x3
    line1[1, 1, :] = [0, 12, 12, 15]  # y3


    To select array of intersections from result, use the following idiom:

    value = intersection(line0, line1)
    mask = -numpy.isnan(value[:, 0])
    v = value[mask]
    """

    line0 = ensure_numeric(line0, numpy.float)
    line1 = ensure_numeric(line1, numpy.float)

    x0, y0 = line0[0, :]
    x1, y1 = line0[1, :]

    # Special treatment of return value if line1 was non vectorised
    if len(line1.shape) == 2:
        one_point = True
        # Broadcast to vectorised version
        line1 = line1.reshape(2, 2, 1)
    else:
        one_point = False

    # Extract vectorised coordinates
    x2 = line1[0, 0, :]
    y2 = line1[0, 1, :]
    x3 = line1[1, 0, :]
    y3 = line1[1, 1, :]

    # Calculate denominator (lines are parallel if it is 0)
    y3y2 = y3 - y2
    x3x2 = x3 - x2
    x1x0 = x1 - x0
    y1y0 = y1 - y0
    x2x0 = x2 - x0
    y2y0 = y2 - y0
    denominator = y3y2 * x1x0 - x3x2 * y1y0

    # Suppress numpy warnings (as we'll be dividing by zero)
    original_numpy_settings = numpy.seterr(invalid='ignore', divide='ignore')

    u0 = (y3y2 * x2x0 - x3x2 * y2y0) / denominator
    u1 = (x2x0 * y1y0 - y2y0 * x1x0) / denominator

    # Restore numpy warnings
    numpy.seterr(**original_numpy_settings)

     # Only points that lie within given line segments are true intersections
    mask = (0.0 <= u0) * (u0 <= 1.0) * (0.0 <= u1) * (u1 <= 1.0)

    # Calculate intersection points
    x = x0 + u0[mask] * x1x0
    y = y0 + u0[mask] * y1y0

    # Return intersection points as N x 2 array
    N = line1.shape[2]
    result = numpy.zeros((N, 2)) * numpy.nan
    result[mask, 0] = x
    result[mask, 1] = y

    # Special treatment of return value if line1 was non vectorised
    if one_point:
#.........这里部分代码省略.........

def interpolate_polygon_points(source, target,
                               layer_name=None):
    """Interpolate from polygon vector layer to point vector data

    Args:
        * source: Vector data set (polygon)
        * target: Vector data set (points)
        * layer_name: Optional name of returned interpolated layer.
              If None the name of target is used for the returned layer.

    Output
        I: Vector data set; points located as target with values interpolated
        from source

    Note
        All attribute names from polygons are transferred to the points
        that are inside them.
    """

    msg = ('Vector layer to interpolate to must be point geometry. '
           'I got OGR geometry type %s'
           % geometry_type_to_string(target.geometry_type))
    verify(target.is_point_data, msg)

    msg = ('Name must be either a string or None. I got %s'
           % (str(type(target)))[1:-1])
    verify(layer_name is None or
           isinstance(layer_name, basestring), msg)

    attribute_names = source.get_attribute_names()

    #----------------
    # Start algorithm
    #----------------

    # Extract point features
    points = ensure_numeric(target.get_geometry())
    attributes = target.get_data()
    original_geometry = target.get_geometry()  # Geometry for returned data

    # Extract polygon features
    geom = source.get_geometry(as_geometry_objects=True)
    data = source.get_data()
    verify(len(geom) == len(data))

    # Include polygon_id as attribute
    attribute_names.append('polygon_id')
    attribute_names.append(DEFAULT_ATTRIBUTE)

    # Augment point features with empty attributes from polygon
    for a in attributes:
        # Create all attributes that exist in source
        for key in attribute_names:
            a[key] = None

    # Traverse polygons and assign attributes to points that fall inside
    for i, polygon in enumerate(geom):
        # Carry all attributes across from source
        poly_attr = data[i]

        # Assign default attribute to indicate points inside
        poly_attr[DEFAULT_ATTRIBUTE] = True

        # Clip data points by polygons and add polygon attributes
        indices = inside_polygon(points, polygon.outer_ring,
                                 holes=polygon.inner_rings)

        for k in indices:
            for key in poly_attr:
                # Assign attributes from polygon to points
                attributes[k][key] = poly_attr[key]
            attributes[k]['polygon_id'] = i  # Store id for associated polygon

    # Create new Vector instance and return
    V = Vector(data=attributes,
               projection=target.get_projection(),
               geometry=original_geometry,
               name=layer_name)
    return V

def interpolate_polygon_points(source, target,
                               layer_name=None,
                               attribute_name=None):
    """Interpolate from polygon vector layer to point vector data

    Args:
        * source: Vector data set (polygon)
        * target: Vector data set (points)
        * layer_name: Optional name of returned interpolated layer.
              If None the name of target is used for the returned layer.
        * attribute_name: Name for new attribute.
              If None (default) the name of source is used

    Output
        I: Vector data set; points located as target with values interpolated
        from source
    """

    msg = ('Vector layer to interpolate to must be point geometry. '
           'I got OGR geometry type %s'
           % geometrytype2string(target.geometry_type))
    verify(target.is_point_data, msg)

    msg = ('Name must be either a string or None. I got %s'
           % (str(type(target)))[1:-1])
    verify(layer_name is None or
           isinstance(layer_name, basestring), msg)

    msg = ('Attribute must be either a string or None. I got %s'
           % (str(type(target)))[1:-1])
    verify(attribute_name is None or
           isinstance(attribute_name, basestring), msg)

    attribute_names = source.get_attribute_names()
    if attribute_name is not None:
        msg = ('Requested attribute "%s" did not exist in %s'
               % (attribute_name, attribute_names))
        verify(attribute_name in attribute_names, msg)

    #----------------
    # Start algorithm
    #----------------

    # Extract point features
    points = ensure_numeric(target.get_geometry())
    attributes = target.get_data()
    original_geometry = target.get_geometry()  # Geometry for returned data

    # Extract polygon features
    geom = source.get_geometry()
    data = source.get_data()
    verify(len(geom) == len(data))

    # Augment point features with empty attributes from polygon
    for a in attributes:
        if attribute_name is None:
            # Use all attributes
            for key in attribute_names:
                a[key] = None
        else:
            # Use only requested attribute
            # FIXME (Ole): Test for this is not finished
            a[attribute_name] = None

        # Always create default attribute flagging if point was
        # inside any of the polygons
        a[DEFAULT_ATTRIBUTE] = None

    # Traverse polygons and assign attributes to points that fall inside
    for i, polygon in enumerate(geom):
        if attribute_name is None:
            # Use all attributes
            poly_attr = data[i]
        else:
            # Use only requested attribute
            poly_attr = {attribute_name: data[i][attribute_name]}

        # Assign default attribute to indicate points inside
        poly_attr[DEFAULT_ATTRIBUTE] = True

        # Clip data points by polygons and add polygon attributes
        indices = inside_polygon(points, polygon)
        for k in indices:
            for key in poly_attr:
                # Assign attributes from polygon to points
                attributes[k][key] = poly_attr[key]

    # Create new Vector instance and return
    V = Vector(data=attributes,
               projection=target.get_projection(),
               geometry=original_geometry,
               name=layer_name)
    return V

def intersection(line0, line1, rtol=1.0e-12, atol=1.0e-12):
    """Returns intersecting point between two line segments.

    However, if parallel lines coincide partly (i.e. share a common segment),
    the line segment where lines coincide is returned

    Inputs:
        line0, line1: Each defined by two end points as in:
                      [[x0, y0], [x1, y1]]
                      A line can also be a 2x2 numpy array with each row
                      corresponding to a point.
        rtol, atol: Tolerances passed onto numpy.allclose

    Output:
        status, value - where status and value is interpreted as follows:
        status == 0: no intersection, value set to None.
        status == 1: intersection point found and returned in value as [x,y].
        status == 2: Collinear overlapping lines found.
                     Value takes the form [[x0,y0], [x1,y1]] which is the
                     segment common to both lines.
        status == 3: Collinear non-overlapping lines. Value set to None.
        status == 4: Lines are parallel. Value set to None.
    """

    line0 = ensure_numeric(line0, numpy.float)
    line1 = ensure_numeric(line1, numpy.float)

    x0, y0 = line0[0, :]
    x1, y1 = line0[1, :]

    x2, y2 = line1[0, :]
    x3, y3 = line1[1, :]

    denom = (y3 - y2) * (x1 - x0) - (x3 - x2) * (y1 - y0)
    u0 = (x3 - x2) * (y0 - y2) - (y3 - y2) * (x0 - x2)
    u1 = (x2 - x0) * (y1 - y0) - (y2 - y0) * (x1 - x0)

    if numpy.allclose(denom, 0.0, rtol=rtol, atol=atol):
        # Lines are parallel - check if they are collinear
        if numpy.allclose([u0, u1], 0.0, rtol=rtol, atol=atol):
            # We now know that the lines are collinear
            state = (point_on_line([x0, y0], line1, rtol=rtol, atol=atol),
                     point_on_line([x1, y1], line1, rtol=rtol, atol=atol),
                     point_on_line([x2, y2], line0, rtol=rtol, atol=atol),
                     point_on_line([x3, y3], line0, rtol=rtol, atol=atol))

            return collinearmap[state]([x0, y0], [x1, y1],
                                       [x2, y2], [x3, y3])
        else:
            # Lines are parallel but aren't collinear
            return 4, None
    else:
        # Lines are not parallel, check if they intersect
        u0 = u0 / denom
        u1 = u1 / denom

        x = x0 + u0 * (x1 - x0)
        y = y0 + u0 * (y1 - y0)

        # Sanity check - can be removed to speed up if needed
        if not numpy.allclose(x, x2 + u1 * (x3 - x2), rtol=rtol, atol=atol):
            raise Exception
        if not numpy.allclose(y, y2 + u1 * (y3 - y2), rtol=rtol, atol=atol):
            raise Exception

        # Check if point found lies within given line segments
        if 0.0 <= u0 <= 1.0 and 0.0 <= u1 <= 1.0:
            # We have intersection
            return 1, numpy.array([x, y])
        else:
            # No intersection
            return 0, None