Clipping vs Intersecting

Question

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?

Answer 1

I wanted to post the below in-case it helps anyone. I managed to improve my clipping time considerably to 35 seconds by removing any holes in my polygon (fortunately I didn't need these holes anyway so a double-win).

The code below creates a rough polygon from points which are reachable within X minutes (by buffering); intersects them with a road network; then buffers the road-network slightly to create a nice isochrone. (Only thing I would like to add is a way to join the polygons along the road)

These match arcGIS generated isochrones quite well:

import subprocess
import fiona
from shapely.geometry import Point, mapping, shape, MultiLineString, MultiPolygon, Polygon, LineString
from shapely.ops import unary_union
from fiona.crs import from_epsg

# Projections:
# from pyproj import Proj, transform
# wgs84=Proj("+init=EPSG:4326")  # LatLon with WGS84 datum used by GPS units and Google Earth
# osgb36=Proj("+init=EPSG:27700")  # UK Ordnance Survey, 1936 datum
# original = Proj(init='EPSG:4326')
# destination = Proj(init='EPSG:4326')

# Create shape-file of Intersected
# Note try using dictreader
pts = []
with open(produce_csv_file) as f:
    for x in csv.reader(f):
        # Lat, Lng, Minutes and keep those points within
        if float(x[2]) <= float(duration):
            # Shapely:
            pts.append(Point(
                # transform(original,destination,float(x[1]), float(x[0]))
                float(x[1]), float(x[0])
            ))

# Buffer (divide spacing by 100)
# Buffer + 10%
buffer = [point.buffer(1.1*float(RADIUS_KM)/100) for point in pts]

# Union
merged = unary_union(buffer)

# Remove holes
exterior_polys = []
for poly_l in merged:
    exterior_polys.append(Polygon(poly_l.exterior))
merged = MultiPolygon(exterior_polys)

# Save the interim shape-file
schema = {
    'geometry': 'Polygon',
    'properties': {'id': 'int'},
}

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

# Bounds
bounds = merged.bounds
x_min = '%.5f' % bounds[0]
x_max = '%.5f' % bounds[2]
y_min = '%.5f' % bounds[1]
y_max = '%.5f' % bounds[3]
print(x_min, x_max, y_min, y_max)

# Clip using GDAL
# http://www.gisinternals.com/release.php
print("Generation Time to Rough Shapefile: %.0f seconds" % (time.clock() - start))
print("Beginning Clip")
start = time.clock()

#Make sure shape-file has an index
#ogrinfo "../ROADLINK_Project.shp" -sql "CREATE SPATIAL INDEX ON ROADLINK_Project"
subprocess.call(["C:\Program Files\GDAL\ogr2ogr", "-f", "ESRI Shapefile", "-spat", x_min, y_min, x_max, y_max, "-clipsrc", produce_poly_shape, interim_poly_shape, my_roads_shape])

print('Clipped in %.0f seconds' % (time.clock()-start))

road_lines = ([shape(pol['geometry']) for pol in fiona.open(interim_poly_shape)])
buffer = []
# Buffer + 10%
for road in road_lines:
    buffer.append(road.buffer(1.1*float(RADIUS_KM)/100))

# Union
merged = unary_union(buffer)

# Remove holes
exterior_polys = []
for poly_l in merged:
    exterior_polys.append(Polygon(poly_l.exterior))
merged = MultiPolygon(exterior_polys)

# Save the final shape-file
schema = {
    'geometry': 'Polygon',
    'properties': {'id': 'int'},
}

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

print('Complete! Total time-taken: %.0f' % (time.clock()- overall_start))

Answer 2

@user30184 helped me to make pretty good progress and I want to share my results.

I save the bounds of my small clipper polygon:

# Bounds
bounds = merged.bounds
x_min = '%.5f' % bounds[0]
x_max = '%.5f' % bounds[2]
y_min = '%.5f' % bounds[1]
y_max = '%.5f' % bounds[3]

And then run the following command:

subprocess.call(["C:\Program Files\GDAL\ogr2ogr", "-f", "ESRI Shapefile", "-spat", x_min, y_min, x_max, y_max, "-clipsrc", clipping_shp, output_shp, input_shp])

Specifying --spat -> this makes my command much faster (more details below).

I tried to also create an index:

ogrinfo "../ROADLINK_Project.shp" -sql "CREATE SPATIAL INDEX ON ROADLINK_Project"

However, didn't notice any speed improvement. I think this is because I create the shape using arcGIS and thus ogr2ogr uses the arcGIS index. So I ran a few experiments (my .shp is 504MB, my arcGIS index .shx is 24MB and my crated quadtree .qix index is also around 25MB ... I wasn't sure what to set the recursion level to and so used the default)

1) Removing all index files and adding the default index using ogrinfo it takes 592 seconds to intersect a 500mb UK Roads shapefile with a 5mb polygon using the ogr2ogr command

2) Using arcGIS it takes 56 seconds to perform an intersection

3) Using arcGIS it takes 105 seconds to perform a clip

4) If I remove the index file and re-run the ogr2ogr command it takes 621 seconds (so it appears that the index file saves me 30 seconds).