Python shapely.geometry.LineString() Examples

def _construct_centerline(self):
        vertices, ridges = self._get_voronoi_vertices_and_ridges()
        linestrings = []
        for ridge in ridges:
            if self._ridge_is_finite(ridge):
                starting_point = self._create_point_with_restored_coordinates(
                    x=vertices[ridge[0]][0], y=vertices[ridge[0]][1]
                )
                ending_point = self._create_point_with_restored_coordinates(
                    x=vertices[ridge[1]][0], y=vertices[ridge[1]][1]
                )
                linestring = LineString((starting_point, ending_point))

                if self._linestring_is_within_input_geometry(linestring):
                    linestrings.append(linestring)

        if len(linestrings) < 2:
            raise exceptions.TooFewRidgesError

        return unary_union(linestrings) 

def wkt_to_shp(wkt_list, shp_file):
    '''Take output of build_graph_wkt() and render the list of linestrings
    into a shapefile
    # https://gis.stackexchange.com/questions/52705/how-to-write-shapely-geometries-to-shapefiles
    '''

    # Define a linestring feature geometry with one attribute
    schema = {
        'geometry': 'LineString',
        'properties': {'id': 'int'},
    }

    # Write a new shapefile
    with fiona.open(shp_file, 'w', 'ESRI Shapefile', schema) as c:
        for i, line in enumerate(wkt_list):
            shape = shapely.wkt.loads(line)
            c.write({
                    'geometry': mapping(shape),
                    'properties': {'id': i},
                    })

    return


############################################################################### 

def WriteHillslopeTracesShp(DataDirectory,FilenamePrefix,ThinningFactor=1, CustomExtent=[-9999]):
    """
    This function writes a shapefile of hillslope traces

    Args:
        DataDirectory: the data directory
        FilenamePrefix: the file name prefix
        ThinningFactor (int): This determines how many of the traces are discarded. Higher numbers mean more traces are discarded
        CustomExtent (list): if this is [-9999] the extent is grabbed from the raster. Otherwise you give it a 4 element list giving the extent of the area of interest.

    Returns: None, but writes a shapefile

    Author: MDH
    """

    #read the raw data to geodataframe
    geo_df = ReadHillslopeTraces(DataDirectory,FilenamePrefix,ThinningFactor, CustomExtent)
    Suffix = '_hillslope_traces'
    WriteFilename = DataDirectory+FilenamePrefix+Suffix+'.shp'

    #aggregate these points with group by
    geo_df = geo_df.groupby(['HilltopID'])['geometry'].apply(lambda x: LineString(x.tolist()))
    geo_df = GeoDataFrame(geo_df, geometry='geometry')
    geo_df.to_file(WriteFilename, driver='ESRI Shapefile') 

def read_terrace_centrelines(DataDirectory, shapefile_name):
    """
    This function reads in a shapefile of terrace centrelines
    using shapely and fiona

    Args:
        DataDirectory (str): the data directory
        shapefile_name (str): the name of the shapefile

    Returns: shapely polygons with terraces

    Author: FJC
    """
    Lines = {}
    with fiona.open(DataDirectory+shapefile_name, 'r') as input:
        for f in input:
            this_line = LineString(shape(f['geometry']))
            this_id = f['properties']['id']
            Lines[this_id] = this_line
    return Lines 

def dummy_constant_bed():
    dx = 1.

    hmax = 3000.
    hmin = 1000.
    nx = 200
    map_dx = 100.
    widths = 3.

    surface_h = np.linspace(hmax, hmin, nx)
    bed_h = surface_h
    widths = surface_h * 0. + widths
    coords = np.arange(0, nx - 0.5, 1)
    line = shpg.LineString(np.vstack([coords, coords * 0.]).T)
    return [flowline.RectangularBedFlowline(line, dx, map_dx, surface_h,
                                            bed_h, widths)] 

def dummy_parabolic_bed(hmax=3000., hmin=1000., nx=200, map_dx=100.,
                        default_shape=5.e-3,
                        from_other_shape=None, from_other_bed=None):
    dx = 1.

    surface_h = np.linspace(hmax, hmin, nx)
    bed_h = surface_h * 1
    shape = surface_h * 0. + default_shape
    if from_other_shape is not None:
        shape[0:len(from_other_shape)] = from_other_shape

    if from_other_bed is not None:
        bed_h[0:len(from_other_bed)] = from_other_bed

    coords = np.arange(0, nx - 0.5, 1)
    line = shpg.LineString(np.vstack([coords, coords * 0.]).T)
    return [flowline.ParabolicBedFlowline(line, dx, map_dx, surface_h,
                                          bed_h, shape)] 

def gdf_to_nx(gdf_network, approach="primal", length="mm_len"):
    """
    Convert LineString GeoDataFrame to networkx.MultiGraph

    Parameters
    ----------
    gdf_network : GeoDataFrame
        GeoDataFrame containing objects to convert
    approach : str, default 'primal'
        Decide wheter genereate ``'primal'`` or ``'dual'`` graph.
    length : str, default mm_len
        name of attribute of segment length (geographical) which will be saved to graph

    Returns
    -------
    networkx.Graph
        Graph

    """
    gdf_network = gdf_network.copy()
    if "key" in gdf_network.columns:
        gdf_network.rename(columns={"key": "__key"}, inplace=True)
    # generate graph from GeoDataFrame of LineStrings
    net = nx.MultiGraph()
    net.graph["crs"] = gdf_network.crs
    gdf_network[length] = gdf_network.geometry.length
    fields = list(gdf_network.columns)

    if approach == "primal":
        _generate_primal(net, gdf_network, fields)

    elif approach == "dual":
        _generate_dual(net, gdf_network, fields)

    else:
        raise ValueError(
            "Approach {} is not supported. Use 'primal' or 'dual'.".format(approach)
        )

    return net 

def convert_pix_lstring_to_geo(wkt_lstring, im_file):
    '''Convert linestring in pixel coords to geo coords'''
    shape = wkt_lstring  # shapely.wkt.loads(lstring)
    x_pixs, y_pixs = shape.coords.xy
    coords_latlon = []
    coords_utm = []
    for (x, y) in zip(x_pixs, y_pixs):
        lon, lat = apls_utils.pixelToGeoCoord(x, y, im_file)
        [utm_east, utm_north, utm_zone, utm_letter] = utm.from_latlon(lat, lon)
        coords_utm.append([utm_east, utm_north])
        coords_latlon.append([lon, lat])

    lstring_latlon = LineString([Point(z) for z in coords_latlon])
    lstring_utm = LineString([Point(z) for z in coords_utm])

    return lstring_latlon, lstring_utm, utm_zone, utm_letter


############################################################################### 

def create_frames(mapfile):

    # increment in steps of 3 units from 1 to 35
    min_buffer = 1
    max_buffer = 35
    step = 3 

    # create a line
    line = LineString([(0, 0), (100, 100), (0, 200), (200, 200), (300, 100), (100, 0)])
    # set the map extent to this line
    mapfile["extent"] = " ".join(map(str, line.buffer(max_buffer*2).bounds))

    all_images = []

    # create expanding buffers
    for dist in range(min_buffer, max_buffer, step):
        all_images.append(create_frame(mapfile, line, dist))

    # create shrinking buffers
    for dist in range(max_buffer, min_buffer, -step):
        all_images.append(create_frame(mapfile, line, dist))

    return all_images 

def dummy_constant_bed_cliff(hmax=3000., hmin=1000., nx=200, map_dx=100.,
                             cliff_height=250.):
    """
    I introduce a cliff in the bed to test the mass conservation of the models
    Such a cliff could be real or a DEM error/artifact
    """
    dx = 1.

    surface_h = np.linspace(hmax, hmin, nx)

    surface_h[50:] = surface_h[50:] - cliff_height

    bed_h = surface_h
    widths = surface_h * 0. + 1.

    coords = np.arange(0, nx - 0.5, 1)
    line = shpg.LineString(np.vstack([coords, coords * 0.]).T)
    return [flowline.RectangularBedFlowline(line, dx, map_dx, surface_h,
                                            bed_h, widths)] 

def WriteHillslopeTracesShp(DataDirectory,FilenamePrefix):
    """
    This function writes a shapefile of hillslope traces

    Args:
        DataDirectory: the data directory
        FilenamePrefix: the file name prefix

    Author: MDH
    """
    
    #read the raw data to geodataframe
    geo_df = ReadHillslopeTraces(DataDirectory,FilenamePrefix)
    Suffix = '_hillslope_traces'
    WriteFilename = DataDirectory+FilenamePrefix+Suffix+'.shp'
    
    #aggregate these points with group by
    geo_df = geo_df.groupby(['HilltopID'])['geometry'].apply(lambda x: LineString(x.tolist()))
    geo_df = GeoDataFrame(geo_df, geometry='geometry')
    geo_df.to_file(WriteFilename, driver='ESRI Shapefile') 

def _line_iterator(obj):
    if isinstance(obj, (sg.LineString)):
        yield obj
    elif isinstance(obj, (sg.MultiLineString)):
        for obj2 in obj.geoms:
            yield obj2
    elif isinstance(obj, (sg.Polygon)):
        yield sg.LineString(obj.exterior)
        for obj2 in obj.interiors:
            yield sg.LineString(obj2)
    elif isinstance(obj, (sg.MultiPolygon)):
        for obj2 in obj.geoms:
            yield sg.LineString(obj2.exterior)
            for obj3 in obj2.interiors:
                yield sg.LineString(obj3)
    else:
        try:
            tup = tuple(obj)
        except TypeError:
            raise TypeError('Could not use type %s' % type(obj))
        else:
            for obj2 in tup:
                for line in _line_iterator(obj2):
                    yield line 

def __init__(self, line=None, dx=None, map_dx=None,
                 surface_h=None, bed_h=None, widths=None, rgi_id=None,
                 water_level=None):
        """ Instanciate.

        Parameters
        ----------
        line : :py:class:`shapely.geometry.LineString`
            the geometrical line of a :py:class:`oggm.Centerline`

        Properties
        ----------
        #TODO: document properties
        """
        super(RectangularBedFlowline, self).__init__(line, dx, map_dx,
                                                     surface_h, bed_h,
                                                     rgi_id=rgi_id,
                                                     water_level=water_level)

        self._widths = widths 

def __init__(self, line=None, dx=None, map_dx=None,
                 surface_h=None, bed_h=None, bed_shape=None, rgi_id=None,
                 water_level=None):
        """ Instanciate.

        Parameters
        ----------
        line : :py:class:`shapely.geometry.LineString`
            the geometrical line of a :py:class:`oggm.Centerline`

        Properties
        ----------
        #TODO: document properties
        """
        super(ParabolicBedFlowline, self).__init__(line, dx, map_dx,
                                                   surface_h, bed_h,
                                                   rgi_id=rgi_id,
                                                   water_level=water_level)

        assert np.all(np.isfinite(bed_shape))
        self.bed_shape = bed_shape 

def test_inversion_tributary_sf_adhikari(self, inversion_gdir):

        old_model_sf = cfg.PARAMS['use_shape_factor_for_fluxbasedmodel']
        old_inversion_sf = cfg.PARAMS['use_shape_factor_for_inversion']
        cfg.PARAMS['use_shape_factor_for_fluxbasedmodel'] = 'Adhikari'
        cfg.PARAMS['use_shape_factor_for_inversion'] = 'Adhikari'

        fls = dummy_width_bed_tributary(map_dx=inversion_gdir.grid.dx)
        for fl in fls:
            fl.is_rectangular = np.ones(fl.nx).astype(np.bool)
        mb = LinearMassBalance(2600.)

        model = FluxBasedModel(fls, mb_model=mb, y0=0.)
        model.run_until_equilibrium()

        fls = []
        for fl in model.fls:
            pg = np.where(fl.thick > 0)
            line = shpg.LineString([fl.line.coords[int(p)] for p in pg[0]])
            sh = fl.surface_h[pg]
            flo = centerlines.Centerline(line, dx=fl.dx,
                                         surface_h=sh)
            flo.widths = fl.widths[pg]
            flo.is_rectangular = np.ones(flo.nx).astype(np.bool)
            fls.append(flo)

        fls[0].set_flows_to(fls[1])

        inversion_gdir.write_pickle(copy.deepcopy(fls), 'inversion_flowlines')

        climate.apparent_mb_from_linear_mb(inversion_gdir)
        inversion.prepare_for_inversion(inversion_gdir)
        v, _ = inversion.mass_conservation_inversion(inversion_gdir)

        assert_allclose(v, model.volume_m3, rtol=0.02)
        if do_plot:  # pragma: no cover
            self.double_plot(model, inversion_gdir)

        cfg.PARAMS['use_shape_factor_for_fluxbasedmodel'] = old_model_sf
        cfg.PARAMS['use_shape_factor_for_inversion'] = old_inversion_sf 

def test_inversion_cliff_sf_huss(self, inversion_gdir):
        old_model_sf = cfg.PARAMS['use_shape_factor_for_fluxbasedmodel']
        old_inversion_sf = cfg.PARAMS['use_shape_factor_for_inversion']
        cfg.PARAMS['use_shape_factor_for_fluxbasedmodel'] = 'Huss'
        cfg.PARAMS['use_shape_factor_for_inversion'] = 'Huss'

        fls = dummy_constant_bed_cliff(map_dx=inversion_gdir.grid.dx,
                                       cliff_height=100)
        for fl in fls:
            fl.is_rectangular = np.ones(fl.nx).astype(np.bool)
        mb = LinearMassBalance(2600.)

        model = FluxBasedModel(fls, mb_model=mb, y0=0.)
        model.run_until_equilibrium()
        fls = []
        for fl in model.fls:
            pg = np.where(fl.thick > 0)
            line = shpg.LineString([fl.line.coords[int(p)] for p in pg[0]])
            sh = fl.surface_h[pg]
            flo = centerlines.Centerline(line, dx=fl.dx,
                                         surface_h=sh)
            flo.widths = fl.widths[pg]
            flo.is_rectangular = np.ones(flo.nx).astype(np.bool)
            fls.append(flo)
        inversion_gdir.write_pickle(copy.deepcopy(fls), 'inversion_flowlines')

        climate.apparent_mb_from_linear_mb(inversion_gdir)
        inversion.prepare_for_inversion(inversion_gdir)
        v, _ = inversion.mass_conservation_inversion(inversion_gdir)

        assert_allclose(v, model.volume_m3, rtol=0.05)
        if do_plot:  # pragma: no cover
            self.simple_plot(model, inversion_gdir)

        cfg.PARAMS['use_shape_factor_for_fluxbasedmodel'] = old_model_sf
        cfg.PARAMS['use_shape_factor_for_inversion'] = old_inversion_sf 

def test_inversion_noisy_sf_adhikari(self, inversion_gdir):
        old_model_sf = cfg.PARAMS['use_shape_factor_for_fluxbasedmodel']
        old_inversion_sf = cfg.PARAMS['use_shape_factor_for_inversion']
        cfg.PARAMS['use_shape_factor_for_fluxbasedmodel'] = 'Adhikari'
        cfg.PARAMS['use_shape_factor_for_inversion'] = 'Adhikari'

        fls = dummy_noisy_bed(map_dx=inversion_gdir.grid.dx)
        for fl in fls:
            fl.is_rectangular = np.ones(fl.nx).astype(np.bool)
        mb = LinearMassBalance(2600.)

        model = FluxBasedModel(fls, mb_model=mb, y0=0.)
        model.run_until_equilibrium()
        fls = []
        for fl in model.fls:
            pg = np.where(fl.thick > 0)
            line = shpg.LineString([fl.line.coords[int(p)] for p in pg[0]])
            sh = fl.surface_h[pg]
            flo = centerlines.Centerline(line, dx=fl.dx,
                                         surface_h=sh)
            flo.widths = fl.widths[pg]
            flo.is_rectangular = np.ones(flo.nx).astype(np.bool)
            fls.append(flo)
        inversion_gdir.write_pickle(copy.deepcopy(fls), 'inversion_flowlines')

        climate.apparent_mb_from_linear_mb(inversion_gdir)
        inversion.prepare_for_inversion(inversion_gdir)
        v, _ = inversion.mass_conservation_inversion(inversion_gdir)

        assert_allclose(v, model.volume_m3, rtol=0.05)
        if do_plot:  # pragma: no cover
            self.simple_plot(model, inversion_gdir)

        cfg.PARAMS['use_shape_factor_for_fluxbasedmodel'] = old_model_sf
        cfg.PARAMS['use_shape_factor_for_inversion'] = old_inversion_sf 

def test_inversion_tributary(self, inversion_gdir):

        fls = dummy_width_bed_tributary(map_dx=inversion_gdir.grid.dx)
        mb = LinearMassBalance(2600.)

        model = FluxBasedModel(fls, mb_model=mb, y0=0.)
        model.run_until_equilibrium()

        fls = []
        for fl in model.fls:
            pg = np.where(fl.thick > 0)
            line = shpg.LineString([fl.line.coords[int(p)] for p in pg[0]])
            sh = fl.surface_h[pg]
            flo = centerlines.Centerline(line, dx=fl.dx,
                                         surface_h=sh)
            flo.widths = fl.widths[pg]
            flo.is_rectangular = np.ones(flo.nx).astype(np.bool)
            fls.append(flo)

        fls[0].set_flows_to(fls[1])

        inversion_gdir.write_pickle(copy.deepcopy(fls), 'inversion_flowlines')

        climate.apparent_mb_from_linear_mb(inversion_gdir)
        inversion.prepare_for_inversion(inversion_gdir)
        v, _ = inversion.mass_conservation_inversion(inversion_gdir)

        assert_allclose(v, model.volume_m3, rtol=0.02)
        if do_plot:  # pragma: no cover
            self.double_plot(model, inversion_gdir) 

def test_inversion_noisy_sf_huss(self, inversion_gdir):
        old_model_sf = cfg.PARAMS['use_shape_factor_for_fluxbasedmodel']
        old_inversion_sf = cfg.PARAMS['use_shape_factor_for_inversion']
        cfg.PARAMS['use_shape_factor_for_fluxbasedmodel'] = 'Huss'
        cfg.PARAMS['use_shape_factor_for_inversion'] = 'Huss'

        fls = dummy_noisy_bed(map_dx=inversion_gdir.grid.dx)
        for fl in fls:
            fl.is_rectangular = np.ones(fl.nx).astype(np.bool)
        mb = LinearMassBalance(2600.)

        model = FluxBasedModel(fls, mb_model=mb, y0=0.)
        model.run_until_equilibrium()
        fls = []
        for fl in model.fls:
            pg = np.where(fl.thick > 0)
            line = shpg.LineString([fl.line.coords[int(p)] for p in pg[0]])
            sh = fl.surface_h[pg]
            flo = centerlines.Centerline(line, dx=fl.dx,
                                         surface_h=sh)
            flo.widths = fl.widths[pg]
            flo.is_rectangular = np.ones(flo.nx).astype(np.bool)
            fls.append(flo)
        inversion_gdir.write_pickle(copy.deepcopy(fls), 'inversion_flowlines')

        climate.apparent_mb_from_linear_mb(inversion_gdir)
        inversion.prepare_for_inversion(inversion_gdir)
        v, _ = inversion.mass_conservation_inversion(inversion_gdir)

        assert_allclose(v, model.volume_m3, rtol=0.05)
        if do_plot:  # pragma: no cover
            self.simple_plot(model, inversion_gdir)

        cfg.PARAMS['use_shape_factor_for_fluxbasedmodel'] = old_model_sf
        cfg.PARAMS['use_shape_factor_for_inversion'] = old_inversion_sf 

def test_inversion_cliff_sf_adhikari(self, inversion_gdir):
        old_model_sf = cfg.PARAMS['use_shape_factor_for_fluxbasedmodel']
        old_inversion_sf = cfg.PARAMS['use_shape_factor_for_inversion']
        cfg.PARAMS['use_shape_factor_for_fluxbasedmodel'] = 'Adhikari'
        cfg.PARAMS['use_shape_factor_for_inversion'] = 'Adhikari'

        fls = dummy_constant_bed_cliff(map_dx=inversion_gdir.grid.dx,
                                       cliff_height=100)
        for fl in fls:
            fl.is_rectangular = np.ones(fl.nx).astype(np.bool)
        mb = LinearMassBalance(2600.)

        model = FluxBasedModel(fls, mb_model=mb, y0=0.)
        model.run_until_equilibrium()
        fls = []
        for fl in model.fls:
            pg = np.where(fl.thick > 0)
            line = shpg.LineString([fl.line.coords[int(p)] for p in pg[0]])
            sh = fl.surface_h[pg]
            flo = centerlines.Centerline(line, dx=fl.dx,
                                         surface_h=sh)
            flo.widths = fl.widths[pg]
            flo.is_rectangular = np.ones(flo.nx).astype(np.bool)
            fls.append(flo)
        inversion_gdir.write_pickle(copy.deepcopy(fls), 'inversion_flowlines')

        climate.apparent_mb_from_linear_mb(inversion_gdir)
        inversion.prepare_for_inversion(inversion_gdir)
        v, _ = inversion.mass_conservation_inversion(inversion_gdir)

        assert_allclose(v, model.volume_m3, rtol=0.05)
        if do_plot:  # pragma: no cover
            self.simple_plot(model, inversion_gdir)

        cfg.PARAMS['use_shape_factor_for_fluxbasedmodel'] = old_model_sf
        cfg.PARAMS['use_shape_factor_for_inversion'] = old_inversion_sf 

def linestring():
    return geometry.LineString([(0, 0), (0.8, 0.8), (1.8, 0.95), (2.6, 0.5)]) 

def dummy_trapezoidal_bed(hmax=3000., hmin=1000., nx=200):
    map_dx = 100.
    dx = 1.

    surface_h = np.linspace(hmax, hmin, nx)
    bed_h = surface_h
    widths = surface_h * 0. + 1.6

    lambdas = surface_h * 0. + 2

    coords = np.arange(0, nx - 0.5, 1)
    line = shpg.LineString(np.vstack([coords, coords * 0.]).T)

    return [flowline.TrapezoidalBedFlowline(line, dx, map_dx, surface_h,
                                            bed_h, widths, lambdas)] 

def dummy_mixed_bed(deflambdas=3.5, map_dx=100., mixslice=None):
    dx = 1.
    nx = 200

    surface_h = np.linspace(3000, 1000, nx)
    bed_h = surface_h
    shape = surface_h * 0. + 3.e-03
    if mixslice:
        shape[mixslice] = np.NaN
    else:
        shape[10:20] = np.NaN
    is_trapezoid = ~np.isfinite(shape)
    lambdas = shape * 0.
    lambdas[is_trapezoid] = deflambdas

    widths_m = bed_h * 0. + 10
    section = bed_h * 0.

    coords = np.arange(0, nx - 0.5, 1)
    line = shpg.LineString(np.vstack([coords, coords * 0.]).T)

    fls = flowline.MixedBedFlowline(line=line, dx=dx, map_dx=map_dx,
                                    surface_h=surface_h, bed_h=bed_h,
                                    section=section, bed_shape=shape,
                                    is_trapezoid=is_trapezoid,
                                    lambdas=lambdas, widths_m=widths_m)

    return [fls] 

def dummy_noisy_bed(map_dx=100.):
    dx = 1.
    nx = 200
    np.random.seed(42)
    coords = np.arange(0, nx - 0.5, 1)
    surface_h = np.linspace(3000, 1000, nx)
    surface_h += 100 * np.random.rand(nx) - 50.

    bed_h = surface_h
    widths = surface_h * 0. + 3.
    line = shpg.LineString(np.vstack([coords, coords * 0.]).T)
    return [flowline.RectangularBedFlowline(line, dx, map_dx, surface_h,
                                            bed_h, widths)] 

def dummy_bumpy_bed():
    map_dx = 100.
    dx = 1.
    nx = 200

    coords = np.arange(0, nx - 0.5, 1)
    surface_h = np.linspace(3000, 1000, nx)
    surface_h += 170. * np.exp(-((coords - 30) / 5) ** 2)

    bed_h = surface_h
    widths = surface_h * 0. + 3.
    line = shpg.LineString(np.vstack([coords, coords * 0.]).T)
    return [flowline.RectangularBedFlowline(line, dx, map_dx, surface_h,
                                            bed_h, widths)] 

def dummy_constant_bed(hmax=3000., hmin=1000., nx=200, map_dx=100.,
                       widths=3.):

    dx = 1.

    surface_h = np.linspace(hmax, hmin, nx)
    bed_h = surface_h
    widths = surface_h * 0. + widths
    coords = np.arange(0, nx - 0.5, 1)
    line = shpg.LineString(np.vstack([coords, coords * 0.]).T)
    return [flowline.RectangularBedFlowline(line, dx, map_dx, surface_h,
                                            bed_h, widths)] 

def set_line(self, line):
        """Update the Shapely LineString coordinate.

        Parameters
        ----------
        line : :py:class`shapely.geometry.LineString`
        """

        self.nx = len(line.coords)
        self.line = line
        dis = [line.project(shpg.Point(co)) for co in line.coords]
        self.dis_on_line = np.array(dis)
        xx, yy = line.xy
        self.head = shpg.Point(xx[0], yy[0])
        self.tail = shpg.Point(xx[-1], yy[-1]) 

def _line_extend(uline, dline, dx):
    """Adds a downstream line to a flowline

    Parameters
    ----------
    uline: a shapely.geometry.LineString instance
    dline: a shapely.geometry.LineString instance
    dx: the spacing

    Returns
    -------
    (line, line) : two shapely.geometry.LineString instances. The first
    contains the newly created (long) line, the second only the interpolated
    downstream part (useful for other analyses)
    """

    # First points is easy
    points = [shpg.Point(c) for c in uline.coords]
    dpoints = []

    # Continue as long as line is not finished
    while True:
        pref = points[-1]
        pbs = pref.buffer(dx).boundary.intersection(dline)
        if pbs.type == 'Point':
            pbs = [pbs]
        # Out of the point(s) that we get, take the one farthest from the top
        refdis = dline.project(pref)
        tdis = np.array([dline.project(pb) for pb in pbs])
        p = np.where(tdis > refdis)[0]
        if len(p) == 0:
            break
        points.append(pbs[int(p[0])])
        dpoints.append(pbs[int(p[0])])

    return shpg.LineString(points), shpg.LineString(dpoints) 

def geometry_collection():
    return geometry.GeometryCollection(
        (
            geometry.Point(0, 0),
            geometry.LineString([(0, 0), (0.8, 0.8), (1.8, 0.95), (2.6, 0.5)]),
            geometry.Polygon([[0, 0], [0, 4], [4, 4], [4, 0]]),
        )
    ) 

def two_line_gdf():
    """ Create Line Objects For Testing """
    linea = LineString([(1, 1), (2, 2), (3, 2), (5, 3)])
    lineb = LineString([(3, 4), (5, 7), (12, 2), (10, 5), (9, 7.5)])
    gdf = gpd.GeoDataFrame([1, 2], geometry=[linea, lineb], crs="epsg:4326")
    return gdf