-1

This question already has an answer here:

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.

marked as duplicate by PolyGeo Nov 5 '15 at 6:45

This question has been asked before and already has an answer. If those answers do not fully address your question, please ask a new question.

2

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

  • Thank you for reply, as mentioned above it helped me a great deal. However, as you will see although I am able to implement this - the solution is not feasible: in arcGIS the intersect takes seconds whereas here it takes hours. Do you know if I could install a qGIS library and use that within my script to perform the intersection? – Ilia Oct 21 '15 at 9:48
  • You cannot compare a program written in C/C++ (ArcGIS) and a Python script – gene Oct 21 '15 at 18:58
0

@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?

Not the answer you're looking for? Browse other questions tagged or ask your own question.