Python shapely.ops.split() Examples

def check_track(self, track):
        """
        Sometime track GPX files are effectively doubled -- tracks are there and back,
        rather than just 1-way segments.  This method returns only the 1-way segment
        for the track
        """
        mp        = track.length/2
        mp        = Point(track.interpolate(mp))
        new_track = ops.snap(track, mp, snap_tolerance)
        if not new_track.contains(mp):
            print("Unable to snap track")
            return track
        result    = ops.split(new_track, mp)  
              
        result[0].intersects(result[1])
        isect = result[0].intersection(result[1])
        
        if isect.type == "Point":
            return track
            
        if len(isect)/(len(result[0].coords)-1) > .5:
            return result[0]
        else:
            return track 

def split_track(self, track, point):
        # Snap Point
        new_track  = ops.snap(track, point, snap_tolerance)
        
        while not new_track.contains(point):
            if new_track.project(point) == 0:
                # The split has a tolerance issue with the prior track.
                # We therefore create a new mini-trail, and let the old
                # trail be carried forward.
                print("Making mini segment at %s" % str(point))
                point = new_track.interpolate(snap_tolerance*1.1)
            # Use the nearest point functionality to try and estimate a suitable node point on the line
            point, __ = ops.nearest_points(track, point)
            new_track = ops.snap(track, point, snap_tolerance)
            if not new_track.contains(point):
                raise Exception("Unable to snap %s to track at %f" % (str(point),snap_tolerance))
            
        result = ops.split(new_track, point)
        
        if len(result) == 1:
            return(result[0], None)
            
        return (result[0], result[1]) 

def split_line_by_nearest_points(gdf_line, gdf_points, tolerance_grid_snap, crs):
    """
    Split the union of lines with the union of points resulting

    :param GeoDataFrame gdf_line:  GeoDataFrame with multiple rows of connecting line segments
    :param GeoDataFrame gdf_points: geodataframe with multiple rows of single points

    :returns: ``gdf_segments`` (GeoDataFrame of segments)
    :rtype: GeoDataFrame

    https://github.com/ojdo/python-tools/blob/master/shapelytools.py#L144
    """

    # union all geometries
    line = gdf_line.geometry.unary_union
    line._crs = crs
    snap_points = gdf_points.geometry.unary_union
    snap_points._crs = crs

    # snap and split coords on line
    # returns GeometryCollection
    # snap_points = snap(coords, line, tolerance)
    # snap_points._crs = crs
    split_line = split(line, snap(snap_points, line, tolerance_grid_snap))
    split_line._crs = crs
    segments = [feature for feature in split_line if feature.length > 0.01]

    gdf_segments = gdf(geometry=segments, crs=crs)
    # gdf_segments.columns = ['index', 'geometry']

    return gdf_segments 

def snap_points(points, lines, tolerance, crs):
    length = lines.shape[0]
    for i in range(length):
        for point in points.geometry:
            line = lines.loc[i, "geometry"]
            line._crs = crs
            point._crs = crs
            point_inline_projection = line.interpolate(line.project(point))
            point_inline_projection._crs = crs
            distance_to_line = point.distance(point_inline_projection)
            if (point.x, point.y) in line.coords:
                x = "nothing"
            else:
                if distance_to_line < tolerance:
                    buff = point.buffer(0.1)
                    ### Split the line on the buffer
                    geometry = split(line, buff)
                    geometry._crs = crs
                    line_1_points = [tuple(xy) for xy in geometry[0].coords[:-1]]
                    line_1_points.append((point.x, point.y))
                    line_2_points = []
                    line_2_points.append((point.x, point.y))
                    line_2_points.extend([x for x in geometry[-1].coords[1:]])
                    ### Stitch together the first segment, the interpolated point, and the last segment
                    new_line = linemerge((LineString(line_1_points), LineString(line_2_points)))
                    lines.loc[i, "geometry"] = new_line

    G = points["geometry"].apply(lambda geom: geom.wkb)
    points = points.loc[G.drop_duplicates().index]

    G = lines["geometry"].apply(lambda geom: geom.wkb)
    lines = lines.loc[G.drop_duplicates().index]

    return points, lines 

def _quadrat_cut_geometry(geometry, quadrat_width, min_num=3):
    """
    Split a Polygon or MultiPolygon up into sub-polygons of a specified size.

    Parameters
    ----------
    geometry : shapely.geometry.Polygon or shapely.geometry.MultiPolygon
        the geometry to split up into smaller sub-polygons
    quadrat_width : numeric
        the linear width of the quadrats with which to cut up the geometry (in
        the units the geometry is in)
    min_num : int
        the minimum number of linear quadrat lines (e.g., min_num=3 would
        produce a quadrat grid of 4 squares)

    Returns
    -------
    geometry : shapely.geometry.MultiPolygon
    """
    # create n evenly spaced points between the min and max x and y bounds
    west, south, east, north = geometry.bounds
    x_num = math.ceil((east - west) / quadrat_width) + 1
    y_num = math.ceil((north - south) / quadrat_width) + 1
    x_points = np.linspace(west, east, num=max(x_num, min_num))
    y_points = np.linspace(south, north, num=max(y_num, min_num))

    # create a quadrat grid of lines at each of the evenly spaced points
    vertical_lines = [LineString([(x, y_points[0]), (x, y_points[-1])]) for x in x_points]
    horizont_lines = [LineString([(x_points[0], y), (x_points[-1], y)]) for y in y_points]
    lines = vertical_lines + horizont_lines

    # recursively split the geometry by each quadrat line
    for line in lines:
        geometry = MultiPolygon(split(geometry, line))

    return geometry 

def setup_paths(self):
        """
        Splits the track at each node to generate
        a set of line segments that make up the path
        """
        working_path = self.track
        nodes        = self.generate_nodes()
        node_place   = [x for x in nodes]
        node_place.sort()
        
        for i, dist in enumerate(node_place):
            node = nodes[dist]
            if i == 0:
                continue
            origin      = nodes[node_place[i-1]]
            destination = nodes[node_place[i]]    
            path_name = "%i_%i_%s" % (i-1,i, self.name)
            if i < len(node_place)-1:
                # The last point in the track, therefore it can't be split
                try:
                    path_pts, working_path = self.split_track(working_path, node)
                except:
                    # This was previously an exception, but sometimes paths are dumb.
                    # I don't know the best way to kick the path forward, or to remove it.
                    print("Unable to split track for %s" % path_name)
                    path_pts = working_path
            else:
                path_pts = working_path

            try:
                path = Path(path_name, path_pts, origin.coords[0], destination.coords[0])
                self.paths[path_name] = path
            except:
                raise Exception("Unable to add path to class on:%s" % self.name, i,"of ", len(node_place)-1) 

def calc_connectivity_network(path_streets_shp, building_centroids_df, temp_path_potential_network_shp):
    """
    This script outputs a potential network connecting a series of building points to the closest street network
    the street network is assumed to be a good path to the district heating or cooling network

    :param path_streets_shp: path to street shapefile
    :param building_centroids_df: path to substations in buildings (or close by)
    :param path_potential_network: output path shapefile
    :return:
    """
    # first get the street network
    street_network = gdf.from_file(path_streets_shp)

    # check coordinate system
    street_network = street_network.to_crs(get_geographic_coordinate_system())
    lon = street_network.geometry[0].centroid.coords.xy[0][0]
    lat = street_network.geometry[0].centroid.coords.xy[1][0]
    street_network = street_network.to_crs(get_projected_coordinate_system(lat, lon))
    crs = street_network.crs

    street_network = simplify_liness_accurracy(street_network.geometry.values, SHAPEFILE_TOLERANCE, crs)

    # create terminals/branches form street to buildings
    prototype_network = create_terminals(building_centroids_df, crs, street_network)
    config = cea.config.Configuration()
    locator = cea.inputlocator.InputLocator(scenario=config.scenario)

    # first split in intersections
    prototype_network = one_linestring_per_intersection(prototype_network.geometry.values,
                                                        crs)
    # snap endings of all vectors to ending of all other vectors
    prototype_network = snappy_endings(prototype_network.geometry.values, SNAP_TOLERANCE, crs)

    # calculate intersections
    gdf_points_snapped = calculate_end_points_intersections(prototype_network, crs)

    # snap these points to the lines and transform lines
    gdf_points_snapped, prototype_network = snap_points(gdf_points_snapped, prototype_network, SNAP_TOLERANCE, crs)

    # save for verification purposes
    prototype_network.to_file(locator.get_temporary_file("prototype_network.shp"), driver='ESRI Shapefile')
    # get segments
    potential_network_df = split_line_by_nearest_points(prototype_network, gdf_points_snapped, 1.0, crs)

    # calculate Shape_len field
    potential_network_df["Shape_Leng"] = potential_network_df["geometry"].apply(lambda x: x.length)

    potential_network_df.to_file(temp_path_potential_network_shp, driver='ESRI Shapefile')

    return crs