Update#1 I used the code provided by @GeospatialPython.com in the answer:
import os
import shapefile
from shapely.geometry import *
def cut(line_file, point_file, cut_file):
"""Cuts a line shapefile using a point shapefile."""
# Line shapefile we are going to cut
line_sf = shapefile.Reader(line_file)
# Point shapefile containing the cut points
point_sf = shapefile.Reader(point_file)
# Prefix for output shapefile file name
cut_fn = cut_file
# Create the shapefile writer for the output line shapefile
cut_sf = shapefile.Writer(shapefile.POLYLINE)
cut_sf.fields = list(line_sf.fields)
for sr in line_sf.shapeRecords():
line = LineString(sr.shape.points)
envelope = box(*line.bounds)
pt = None
for shape in point_sf.shapes():
pt2 = Point(*shape.points[0])
if envelope.contains(pt2):
if not pt:
pt = pt2
continue
if line.distance(pt2) < line.distance(pt):
pt = pt2
cut_line = None
if not pt:
cut_line = [line]
else:
coords = list(line.coords)
distance = line.project(pt)
if distance <= 0.0 or distance >= line.length:
return [LineString(line)]
for i, p in enumerate(coords):
pd = line.project(Point(p))
if pd == distance:
return [LineString(coords[:i+1]), LineString(coords[i:])]
if pd > distance:
cp = line.interpolate(distance)
return [
LineString(coords[:i] + [(cp.x, cp.y)]),
LineString([(cp.x, cp.y)] + coords[i:])]
for part in cut_line:
coords = list(part.coords)
cut_sf.line([coords])
cut_sf.record(*sr.record)
cut_sf.save("{}.shp".format(cut_fn))
# If the line shapefile has a cry file, copy it for the new file.
shp_name = os.path.splitext(line_sf.shp.name)[0]
line_prj = "{}.prj".format(shp_name)
if os.path.exists(line_prj):
with open(line_prj, "r") as line_crs:
crs = line_crs.read()
with open("{}.prj".format(cut_fn), "w") as cut_crs:
cut_crs.write(crs)
# Call the function with some shapefile names
# cut("river", "points", "cut")
# cut('tc_500.shp', 'tc_stations.shp', 'cut.shp')
if __name__ == "__main__":
cut('tc_line.shp', 'tc_stations.shp', 'tc_cut')
But returned nothing, is there some misses in this code?
I want to write a python function to split line shapefile with point shapefile (point may be near the line or falls exactly on the line), and the following is the schematic function I expect:
cut line = cut(line shapefile, point shapefile)
the function cut input line based on input point, and if the point is not one the line, it cuts the line at the nearest point (projection point) to the input point.
I now know that pyshp, shapely, etc., could handle shapefiles, and I used and modified code below, which was originally posted in the shapely manual to test the ability of pyshp and shapely.
import shapefile
from shapely.geometry import *
from pprint import pprint
def cut(line, distance):
# Cuts a line in two at a distance from its starting point
if distance <= 0.0 or distance >= line.length:
return [LineString(line)]
coords = list(line.coords)
for i, p in enumerate(coords):
pd = line.project(Point(p))
if pd == distance:
return [
LineString(coords[:i+1]),
LineString(coords[i:])]
if pd > distance:
cp = line.interpolate(distance)
return [
LineString(coords[:i] + [(cp.x, cp.y)]),
LineString([(cp.x, cp.y)] + coords[i:])]
if __name__ == "__main__":
sf = shapefile.Reader('tc_2000.shp')
pprint([list(x.coords) for x in cut(sf, 0.5)])
and an error occurred:
AttributeError: Reader instance has no attribute 'length'
I have a polyline shapefile, which is the input shapefile 'tc_2000.shp' in the test code above, and another point shapefile loaded in QGIS, as the figure below, and obviously the points are not on the line. (The second figure is part of figure 1)
I want to write a python function which inputs are point and polyline shapefile, and the output is a new polyline shapefile which is the original input polyline shapefile being split at position nearest to the point.
The schematic function may look like:
splited polyline = cut(line, point)
and here's the code I have written:
import shapefile
from shapely.geometry import *
from pprint import pprint
def cut(line, point):
coords = list(line.coords)
distance = line.project(point)
if distance <= 0.0 or distance >= line.length:
return [LineString(line)]
for i, p in enumerate(coords):
pd = line.project(Point(p))
print pd
if pd == distance:
print 'pd=distance'
return [LineString(coords[:i+1]), LineString(coords[i:])]
if pd > distance:
print 'pd>distance'
cp = line.interpolate(distance)
return [
LineString(coords[:i] + [(cp.x, cp.y)]),
LineString([(cp.x, cp.y)] + coords[i:])]
if __name__ == "__main__":
line = LineString([(0, 0), (0, 1), (1, 1)])
p = Point([(0.5, 0.5)])
pprint([list(x.coords) for x in cut(line, p)])
and its result is correct:
0.0
1.0
pd>distance
[[(0.0, 0.0), (0.0, 0.5)], [(0.0, 0.5), (0.0, 1.0), (1.0, 1.0)]]
So now I have 2 major questions:
- How to make the input shapefile directly readable by python function
cut(line, point)? - Could this kind of function
cut(line, point)I wrote handle point shapefile with many points and complicated line shapefile with tributaries (as the figure above)?
How to achieve the ultimate goal of getting a splited line?
I am a newbie trying to use python to handle shapefiles.
