I have a shapefile made up of 127 separate polygons that represent different districts. Overlaid on this layer is a polyline layer of 2800 streets, some which are completely bound by the polygons and some which crisscross two or more polygons.

My task is two find the total length of streets within each polygon (district). I've already done this with one district using Vector> Geoprocessing Tools > Intersect and then performing summary statsitics. However, in order to avoid doing this manually 127 times, I would like to write a script that will iterate through each polygon - performing the intersect with the polyline layer and leaving me with 127 separate shapefiles.

I not very familiar with Python - I only know some very basic commands. I've attached images to help you visualize and conceptualize the problem.

enter image description here


If you want to create separate shapefiles for each district, first make sure your polygon layer contains unique attributes (if not, you can create one using the Field Calculator with the expression $rownum). You can then run the Split Vector Layer from Vector > Data Management Tools > Split Vector Layer. Choose a directory (for example: C:/Users/you/Desktop/Test/Shapefiles//) and each polygon district should have its own shapefile.

The following script should intersect each polygon district with the line layer:

import glob, os, processing

line_layer = "path/to/line_layer//"    # This is the path to your line layer
poly_layers = "C:/Users/you/Desktop/Test/Shapefiles//"  # This is the path to your district shapefiles
output_dir = "path/to/output_dir//"    # This is the path for the shapefiles to be saved

for fname in glob.glob("*.shp"):
    processing.runalg("qgis:intersection", line_layer, fname, output_dir + fname)
  • I tried this -saved my shapefiles in a folder called "FINALDELIVERY" but I got the following error: Traceback (most recent call last): File "<input>", line 1, in <module> OSError: [Errno 2] No such file or directory: '/Users/alexander/Desktop/Strelka/FINALDELIVERY' Any ideas? – iskandarblue Oct 27 '15 at 19:23
  • Strange, try adding double slashes at the end of the folder path. – Joseph Oct 28 '15 at 10:30
  • That also didn't work. I set the cwd to my Desktop, but for some reason when I ask the list the files in cwd, none of these files are actually on my Desktop. os.listdir(cwd) ['.DocumentRevisions-V100', '.file', '.fseventsd', '.PKInstallSandboxManager', '.Spotlight-V100', '.Trashes', '.vol', 'Applications', 'bin', 'cores', 'dev', 'etc', 'home', 'installer.failurerequests', 'Library', 'net', 'Network', 'private', 'sbin', 'System', 'tmp', 'Users', 'usr', 'var', 'Volumes'] – iskandarblue Oct 28 '15 at 12:19
  • As soon as I try to change the working directory to any folder that is on my Desktop, I get the following error (however, I can change to cwd to Desktop, just none of the folders that are in Desktop) ------os.chdir("/Users/alexander/Documents/Strelka") Traceback (most recent call last): File "<input>", line 1, in <module> OSError: [Errno 2] No such file or directory: '/Users/alexander/Documents/Strelka' – iskandarblue Oct 28 '15 at 12:21
  • I'm assuming you're using mac? Not quite sure what the problem is but the answer by @snaileater should also work. – Joseph Oct 28 '15 at 12:37

I think there's an integrated tool in QGIS doing what you are expecting :
Go to Vector/Analysis Tools/Sum Line Lengths ...

Unless you want to do the task programmatically, of course ... But in that case you can call this tool from inside Python.

  • Thank you for the tip but that did not work - the output shapefile is the same as the input polygon vector layer however the added column does not correspond to the sum of the lengths of the streets within that polygon. – iskandarblue Oct 27 '15 at 15:44
  • Huh ? working nicely on my side ... did u check the CRS of the layers ? – snaileater Oct 27 '15 at 15:47
  • Yes, the CRS is the same for both layers - EPSG: 4326, WGS84. The length column that is added to the output shapefile has all values that are between 0 and 1, which is impossible because the lengths should be in thousands of meters. – iskandarblue Oct 27 '15 at 15:57
  • I think u won't get meters if your layers are in WGS84 ... – snaileater Oct 27 '15 at 16:02
  • What should I do to get meters? – iskandarblue Oct 27 '15 at 16:05

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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