15

A common requirement in GIS is to apply a processing tool to a number of files or apply a process for a number of features in one file to another file.

Much of these operations are embarrassingly parallel in that the results of the calculations in no way influence any other operation in the loop. Not only that but often the input files are each distinct.

A classic case in point is the tiling out of shape files against files containing polygons to clip them against.

Here is a (tested) classical procedural method to achieve this in a python script for QGIS. (fyi the output of temporary memory files to real files more than halved the time to process my test files)

import processing
import os
input_file="/path/to/input_file.shp"
clip_polygons_file="/path/to/polygon_file.shp"
output_folder="/tmp/test/"
input_layer = QgsVectorLayer(input_file, "input file", "ogr")
QgsMapLayerRegistry.instance().addMapLayer(input_layer)
tile_layer  = QgsVectorLayer(clip_polygons_file, "clip_polys", "ogr")
QgsMapLayerRegistry.instance().addMapLayer(tile_layer)
tile_layer_dp=input_layer.dataProvider()
EPSG_code=int(tile_layer_dp.crs().authid().split(":")[1])
tile_no=0
clipping_polygons = tile_layer.getFeatures()
for clipping_polygon in clipping_polygons:
    print "Tile no: "+str(tile_no)
    tile_no+=1
    geom = clipping_polygon.geometry()
    clip_layer=QgsVectorLayer("Polygon?crs=epsg:"+str(EPSG_code)+\
    "&field=id:integer&index=yes","clip_polygon", "memory")
    clip_layer_dp = clip_layer.dataProvider()
    clip_layer.startEditing()
    clip_layer_feature = QgsFeature()
    clip_layer_feature.setGeometry(geom)
    (res, outFeats) = clip_layer_dp.addFeatures([clip_layer_feature])
    clip_layer.commitChanges()
    clip_file = os.path.join(output_folder,"tile_"+str(tile_no)+".shp")
    write_error = QgsVectorFileWriter.writeAsVectorFormat(clip_layer, \
    clip_file, "system", \
    QgsCoordinateReferenceSystem(EPSG_code), "ESRI Shapefile")
    QgsMapLayerRegistry.instance().addMapLayer(clip_layer)
    output_file = os.path.join(output_folder,str(tile_no)+".shp")
    processing.runalg("qgis:clip", input_file, clip_file, output_file)
    QgsMapLayerRegistry.instance().removeMapLayer(clip_layer.id())

This would be fine except that my input file is 2GB and the polygon clipping file contains 400+ polygons. The resulting process takes over a week on my quad core machine. All the while three cores are just idling.

The solution I have in my head is to export the process out to script files and run them asynchronously using gnu parallel for example. However it seems a shame to have to drop out of QGIS into an OS specific solution rather than use something native to QGIS python. So my question is:

Can I parallelise embarrassingly parallel geographic operations natively inside python QGIS?

If not, then perhaps someone already has the code to send this sort of work off to asynchronous shell scripts?

  • Not familiar with multiprocessing in QGIS, but this ArcGIS-specific example may be of some use: gis.stackexchange.com/a/20352/753 – blah238 Oct 28 '14 at 3:09
  • Looks interesting. I'll see what I can do with it. – Mr Purple Oct 28 '14 at 6:13
11

If you change your program to read the file name from the command line and split up your input file in smaller chunks, you can do something like this using GNU Parallel:

parallel my_processing.py {} /path/to/polygon_file.shp ::: input_files*.shp

This will run 1 job per core.

All new computers have multiple cores, but most programs are serial in nature and will therefore not use the multiple cores. However, many tasks are extremely parallelizeable:

  • Run the same program on many files
  • Run the same program for every line in a file
  • Run the same program for every block in a file

GNU Parallel is a general parallelizer and makes is easy to run jobs in parallel on the same machine or on multiple machines you have ssh access to.

If you have 32 different jobs you want to run on 4 CPUs, a straight forward way to parallelize is to run 8 jobs on each CPU:

Simple scheduling

GNU Parallel instead spawns a new process when one finishes - keeping the CPUs active and thus saving time:

GNU Parallel scheduling

Installation

If GNU Parallel is not packaged for your distribution, you can do a personal installation, which does not require root access. It can be done in 10 seconds by doing this:

(wget -O - pi.dk/3 || curl pi.dk/3/ || fetch -o - http://pi.dk/3) | bash

For other installation options see http://git.savannah.gnu.org/cgit/parallel.git/tree/README

Learn more

See more examples: http://www.gnu.org/software/parallel/man.html

Watch the intro videos: https://www.youtube.com/playlist?list=PL284C9FF2488BC6D1

Walk through the tutorial: http://www.gnu.org/software/parallel/parallel_tutorial.html

Sign up for the email list to get support: https://lists.gnu.org/mailman/listinfo/parallel

  • This is something like I was going to try and attempt, but I need it to all remain inside python. The line needs rewriting to use say Popen for example... Something like: from subprocess import Popen, PIPE p = Popen(["parallel", "ogr2ogr","-clipsrc","clip_file*.shp","output*.shp" input.shp"], stdin=PIPE, stdout=PIPE, stderr=PIPE) Trouble is I dont yet know how to prepare the syntax properly – Mr Purple Oct 28 '14 at 7:58
  • Awesome answer. I hadn't come across triple (or quadruple) colon operators before (although I'm currently doing a Haskell mooc on edX, so no doubt something similar will come). I agree with you about santa, ghosts, fairies and gods, but definitely not goblins :D – John Powell Oct 28 '14 at 8:30
  • @MrPurple I think that comment warrants a question on its own. The answer is definitely too long to put in a comment. – Ole Tange Oct 28 '14 at 8:55
  • OK, thanks for the links. If I formulate an answer using gnu parallel I'll post it here. – Mr Purple Oct 29 '14 at 19:20
  • A good formulation for your my_processing.py can be found at gis.stackexchange.com/a/130337/26897 – Mr Purple Apr 8 '16 at 22:35
4

Rather than using the GNU Parallel method you could use the python mutliprocess module to create a pool of tasks and execute them. I don't have access to a QGIS setup to test it on but multiprocess was added in Python 2.6 so provided that you are using 2.6 or later it should be available. There are a lot of examples online on using this module.

  • 2
    I gave multiprocess a go but I have yet to see it successfully implimented in the embedded python of QGIS. I hit a number of issues as I tried it out. I may post them as seperate questions. As far as I can tell there are no public examples accessible to someone starting out with this. – Mr Purple Nov 4 '14 at 22:02
  • It's a real shame. If someone could write an example of the multiprocess module wrapping a single pyQGIS function as I did with the gnu parallel then we could all go off and parallelise whatever we chose. – Mr Purple Nov 4 '14 at 22:06
  • I agree but as I said I don't have access to a QGIS at the moment. – Steve Barnes Nov 4 '14 at 23:02
  • This question & answer may be of some help if you are running under windows, gis.stackexchange.com/questions/35279/… – Steve Barnes Nov 5 '14 at 8:39
  • @MrPurple and this one gis.stackexchange.com/questions/114260/… gives an example – Steve Barnes Nov 5 '14 at 8:43
3

Here is the gnu parallel solution. With some care most emabrrassingly parallel linux based ogr or saga algorithms could be made to run with it inside your QGIS installation.

Obviously this solution requires the installation of gnu parallel. To install gnu parallel in Ubuntu, for example, go to your terminal and type

sudo apt-get -y install parallel

NB: I couldn't get the parallel shell command to work in Popen or subprocess, which I would have preferred, so I hacked together an export to a bash script and ran that with Popen instead.

Here is the specific shell command using parallel that I wrapped in python

parallel ogr2ogr -skipfailures -clipsrc tile_{1}.shp output_{1}.shp input.shp ::: {1..400}

Each {1} gets swapped out for a number from the {1..400} range and then the four hundred shell commands get managed by gnu parallel to concurrently use all the cores of my i7 :).

Here is the actual python code I wrote to solve the example problem I posted. One could paste it in directly after the end of the code in the question.

import stat
from subprocess import Popen
from subprocess import PIPE
feature_count=tile_layer.dataProvider().featureCount()
subprocess_args=["parallel", \
"ogr2ogr","-skipfailures","-clipsrc",\
os.path.join(output_folder,"tile_"+"{1}"+".shp"),\
os.path.join(output_folder,"output_"+"{1}"+".shp"),\
input_file,\
" ::: ","{1.."+str(feature_count)+"}"]
#Hacky part where I write the shell command to a script file
temp_script=os.path.join(output_folder,"parallelclip.sh")
f = open(temp_script,'w')
f.write("#!/bin/bash\n")
f.write(" ".join(subprocess_args)+'\n')
f.close()
st = os.stat(temp_script)
os.chmod(temp_script, st.st_mode | stat.S_IEXEC)
#End of hacky bash script export
p = Popen([os.path.join(output_folder,"parallelclip.sh")],\
stdin=PIPE, stdout=PIPE, stderr=PIPE)
#Below is the commented out Popen line I couldn't get to work
#p = Popen(subprocess_args, stdin=PIPE, stdout=PIPE, stderr=PIPE)
output, err = p.communicate(b"input data that is passed to subprocess' stdin")
rc = p.returncode
print output
print err

#Delete script and old clip files
os.remove(os.path.join(output_folder,"parallelclip.sh"))
for i in range(feature_count):
    delete_file = os.path.join(output_folder,"tile_"+str(i+1)+".shp")
    nosuff=os.path.splitext(delete_file)[0]
    suffix_list=[]
    suffix_list.append('.shx')
    suffix_list.append('.dbf')
    suffix_list.append('.qpj')
    suffix_list.append('.prj')
    suffix_list.append('.shp')
    suffix_list.append('.cpg')
    for suffix in suffix_list:
        try:
            os.remove(nosuff+suffix)
        except:
            pass

Let me tell you it's really something when you see all the cores fire up to full noise :). Special thanks to Ole and the team that built Gnu Parallel.

It would be nice to have a cross platform solution and it would be nice if I could have figured out the multiprocessing python module for the qgis embedded python but alas it was not to be.

Regardless this solution will serve me and maybe you nicely.

  • Obviously one should comment out the "processing.runalg" line in the first piece of code so the clip isn't run sequentially first before it's run in parallel. Other than that it's simply a matter of copying and pasting the code from the answer underneath the code in the question. – Mr Purple Nov 11 '14 at 23:06
  • If you just want to run lots of processing commands like a set of "qgis:dissolve" applied to different files in parallel then you can see my process for this at purplelinux.co.nz/?p=190 – Mr Purple Apr 8 '16 at 8:14

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.