Python buffer points and clip with line-segment (OS Open Roads in BNG) [duplicate]

Question

I have an array of points (in latitude, longitude coordinates) which I would like to buffer and dissolve (into a polygon) and then intersect the Ordnance Survey Open Roads (projected in BNG) with that polygon and export the results.

Note: I am not interested in preserving any of the shapefile information -> I just want to effectively clip the road-network to show a line only within my polygon as efficiently as possible.

Answer 1

What are you trying to do exactly ?

The first part of your script gives you a MultiPoint geometry with EPSG:3857 (Google Mercator and not WGS84) as projection

1) read the original csv file and write the resulting shapefile

from csv import DictReader
from shapely.geometry import Point, mapping
# unary_union is preferable now
from shapely.ops import unary_union
import fiona 
from fiona.crs import from_epsg
# open the csv file as a simple list of dictionaries
pts = [Point(float(pt['X']),float(pt['Y'])) for pt in DictReader(open("test.csv"), delimiter=',')]
buffer = [point.buffer(0.001) for point in pts]
merged = unary_union(buffer)
# save the shapefile with Fiona
with fiona.open("some_union.shp", "w", driver='ESRI Shapefile', crs=from_epsg(3857), schema=schema) as output:
    output.write({'geometry': mapping(merged), 'properties': {'id': 123}})

2) Change projection

The projection of UK roads network is OSGB36:EPSG7405.

You need to change the projection of your polygons to the the projection of the UK roads network or easier, the projection of the points before buffering. For that you can use the pyproj module (same as the R proj4 package) or the osgeo.GDAL module (same as rgdal) -> look at Convert shapely polygon coordinates with pyproj, shapely and Fiona, for example).

3) intersect the lines and the polygons

Look at Line vs. Polygon Intersection Coordinates with Fiona and Shapely

Answer 2

@Gene - thanks so much for the help! It is incredibly useful. I feel close but still struggling a bit.

To answer your question -> I extract from Google Distance Matrix API a list of points which are within 10 minutes drive of a location (by querying all points in a structured grid in a bounding circle based on 60mph drive). I then buffer those points and join into to a polygon.

What I want to do is load up the UK road network and then clip the roads that fall within the extent of the polygon (I think I worded this poorly since intersection of polygon and line will be a point but I want to clip the line).

The following code works however is HORRIBLY slow. I can clip the full UK road network in seconds in arcGIS. However, here a tiny sub-sample took around 20 minutes? My backup option is installing an R-wrapper and calling gIntersection, however that is also much slower:

# Create shape-file of Intersected
pts = []

# Note try using dictreader
with open(path_to_data) as f:
    for x in csv.reader(f):
        # Lat, Lng, Minutes and keep those points within 10 minutes
        if float(x[2]) <= 10:
                # Transform points at this stage
                pts.append(Point(
                        float(x[1]), float(x[0])
                    ))

# Buffer
buffer = [point.buffer(0.001) for point in pts]
merged = unary_union(buffer)

# Save the Shapefile with Fiona
schema = {
    'geometry': 'Polygon',
    'properties': {'id': 'int'},
}

with fiona.open("H:/polygon_out.shp", "w", driver='ESRI Shapefile', crs=from_epsg(4326), schema=schema) as output:
    output.write({'geometry': mapping(merged), 'properties': {'id': 123}})

input_shp = my_roads_shape
clipping_shp = produce_poly_shape
output_shp = "H:/my_cropped_output.shp"

import subprocess
subprocess.call(["C:\Program Files\GDAL\ogr2ogr", "-f", "ESRI Shapefile", "-clipsrc", clipping_shp, output_shp, input_shp])
print('Done!')

Edit:

I've found a wrapper for the C clipper library (pyclipper):

import pyclipper

pc = pyclipper.Pyclipper()

# Add a single line as the subject.
pc.AddPath([(-1, -1), (2, 1)], pyclipper.PT_SUBJECT, False)

# Add a square as the clipping region.
pc.AddPath([(0, 0), (1, 0), (1, 1), (0, 1)], pyclipper.PT_CLIP, True)

# Clip the line using the rectangle.
solution = pc.Execute2(pyclipper.CT_INTERSECTION, pyclipper.PFT_NONZERO, pyclipper.PFT_NONZERO)

print(pyclipper.PolyTreeToPaths(solution))

I was thinking perhaps I can get this to work?