python - Clipping vs Intersecting

I was curious why clipping shapefiles was just so difficult and whether it was different to intersecting shapefiles?

You will see from one my questions that I gave up trying to clip a polyline with a polygon (here)

I tried various different methods:

1. 
   

   Conversion ogr2ogr (too slow)

   
   

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



2. 
   

   SQL ogr (too slow)

   
   

   from osgeo import ogr
   ogr.UseExceptions()
   ogr_ds = ogr.Open('K:/compdata/OS_Shape_Files/os_open_roads/trim', True)  # Windows: r'C:\path\to\data'
   start = time.clock()
   print('Starting:')
   print(ogr_ds)
   SQL = """\
   
   SELECT ST_Intersection(A.geometry, B.geometry) AS geometry, A.*, B.*
   FROM ROADLINK_Project A, cut_se59ef_poly_60 B
   WHERE ST_Intersects(A.geometry, B.geometry);"""
    layer = ogr_ds.ExecuteSQL(SQL, dialect='SQLITE')
   # copy result back to datasource as a new shapefile
   layer2 = ogr_ds.CopyLayer(layer, 'h1_buf_int_ct3')
   # save, close
   layer = layer2 = ogr_ds = None
   print("Finished in %.0f seconds" % time.clock() - start)
   




3. 
   

   pyclipper (couldn't get this to work)

   
   

   solution = pc.Execute2(pyclipper.CT_INTERSECTION,     pyclipper.PFT_NONZERO, pyclipper.PFT_NONZERO)
   



4. 
   

   rTree (couldn't get this work as I have python 3.4)

   
   

   import fiona
   from shapely.geometry import shape, mapping
   import rtree
   
   
   bufSHP = produce_poly_shape
   intSHP = 'K:/compdata/OS_Shape_Files/os_open_roads    /trim/ROADLINK_Project.shp'
   ctSHP  = 'output.shp'
   
   start = time.clock()
   print('Starting')
   
   with fiona.open(bufSHP, 'r') as layer1:
   with fiona.open(ctSHP, 'r') as layer2:
       # We copy schema and add the  new property for the new resulting shp
   
       schema = layer2.schema.copy()
       schema['properties']['uid'] = 'int:10'
       # We open a first empty shp to write new content from both others shp
       with fiona.open(intSHP, 'w', 'ESRI Shapefile', schema) as layer3:
           index = rtree.index.Index()
           for feat1 in layer1:
               fid = int(feat1['id'])
               geom1 = shape(feat1['geometry'])
               index.insert(fid, geom1.bounds)
   
   
           for feat2 in layer2:
               geom2 = shape(feat2['geometry'])
               for fid in list(index.intersection(geom2.bounds)):
                   if fid != int(feat2['id']):
                       feat1 = layer1[fid]
                       geom1 = shape(feat1['geometry'])
                       if geom1.intersects(geom2):
                           # We take attributes from ctSHP
                           props = feat2['properties']
                           # Then append the uid attribute we want from the other shp
   
                           props['uid'] = feat1['properties']['uid']
                           # Add the content to the right schema in the new shp
                           layer3.write({
                               'properties': props,
                               'geometry': mapping(geom1.intersection(geom2))
                           })
   
   print("Finished in %.0f seconds" % time.clock() - start)
   

Finally, I looked at this but couldn't get it to work with python 3.4.

Then in R I tried:

1. 
   

   gIntersection(x,y, byid=TRUE)

   
   


2. 
   

   gDifference

   
   


3. 
   
   

   custom gClip command I found online

   
   

However with a large polylines shapefile (around 500mb) and a very small polygon shapefile (1mb) the intersection would take a day. This led me to think perhaps I am doing something incredibly stupid and there is a quick clip command?

For example, in arcGIS the intersection command takes a few seconds so surely clipping is even easier? I know very little about GIS but to me it seems as simple as aligning the two layers by one coordinate (assuming same projection) and then selecting the outside of the clipper shape (e.g. in paint) and just deleting it from the other layer. I guess obviously this isn't the case ..

Thus my goal: I create a polygon in my code and I want to clip the UK roads network I load in to the extent of that polygon, buffer the points, dissolve and then output as an array containing the polygon coordinates.

I don't need to retain any features that an intersection would provide just a straight clip. I guess I also don't strictly need shape files and could convert my 500mb UK road network into a geodatabase?