6

I have a script which performs various functions on shapefiles, most of which are polygon-based. However, some shapefiles are point-based and I have to specify where these files are by editing the script. Is there a technique in python which allows the script to identify which shapefiles are polygons, which are points, and run the necessary functions. I have provided a snippet (tested) of my script below which focuses primarily on polygon shapefiles and I have commented out the last line of code which includes the function I use for point shapefiles.

I am using QGIS 2.6.0.

##Test=name
import os
import glob
from os.path import expanduser
home = expanduser("~")

#   Folder path of the directory and the "Results" folder
path_dir = home + "\Desktop\Test\\"
path_res = path_dir + "Results\\"

def run():
#   Set directory and search for all .shp files
    os.chdir(path_dir + "Shapefiles\\")
    for fname in glob.glob("*.shp"): 
#   Clip .shp files with the polygon shapefile "Grid.shp" and save files to Result folder (How to specify this for polygons only?)
        polygon_output=processing.runalg("qgis:clip", path_dir + "Grid.shp", fname, path_res  + "/"+ fname)
#   (How to identify point shapefiles and run the following command?)
#        point_output=processing.runalg("qgis:distancetonearesthub", path_dir + "Grid.shp", fname, 'HubName', 0, 0, path_res  + "/"+ fname)
run()


I managed to get a working script (even though the entire code could be vastly improved!).

##Test=name
import os
import glob
from qgis.core import *
from os.path import expanduser
home = expanduser("~")

#   Folder path of the Results for shapefiles
path_dir = home + "\Desktop\Test\\"
path_res = path_dir + "Results\\"

def run():
#   Set directory and search for all polygon .shp files
    os.chdir(path_dir + "Shapefiles\\")
    for fname in glob.glob("*.shp"): 
        shapefile = QgsVectorLayer( os.chdir(path_dir + "Shapefiles\\"), fname, "ogr" )
        print shapefile.geometryType()

#   Clip .shp files with the polygon shapefile "Grid.shp" and save files to Result folder
        if shapefile.geometryType() == 3:
            polygon_output=processing.runalg("qgis:clip", path_dir + "Grid.shp", fname, path_res  + "/"+ fname)

#   Use Distance to nearest hub function on point.shp files with the polygon shapefile "Grid.shp" as Hub layer and save files to Result folder  
        if shapefile.geometryType() == 3:
            point_output=processing.runalg("qgis:distancetonearesthub", path_dir + "Grid.shp", fname, 'HubName', 0, 0, path_res  + "/"+ fname)
run()
4

You need to open it as a QgsVectorLayer

 layer = QgsVectorLayer('/path/to/shapefile_folder/test.shp', 'test', 'ogr')

QgsVectorLayer objects have a method called geometryType. If you call that method for your newly created layer :

>>> print layer.geometryType()
2

Where 0 is points, 1 is lines and 2 is polygons

  • Many thanks mate, I will attempt to incorporate this into the script and see how far I can get. Question: would the script open the entire file to see the geometryType or just the header as @GIS-Jonathan described? I do have a lot of shapefiles so I wonder if speed would be an issue here. – Joseph Dec 8 '14 at 14:14
3

I guess the "proper" way is to use QGIS' API which I'd hope would expose that. But I don't know how to do that.


An alternative is to look without the shapefile itself. The basics of the specification are on Wikipedia: https://en.wikipedia.org/wiki/Shapefile#Shapefile_shape_format_.28.shp.29

What you want is: bytes 32–35 of the main header which dictates what the geometry type is. The values of this are integers and should be:

Value Shape Type
0 Null Shape
1 Point
3 PolyLine
5 Polygon
8 MultiPoint
11 PointZ
13 PolyLineZ
15 PolygonZ
18 MultiPointZ
21 PointM
23 PolyLineM
25 PolygonM
28 MultiPointM
31 MultiPatch

So all you need to do is parse the header of the file in Python to determine those.

For reference, the official Shapefile "Specification" (actually a white paper) is: http://www.esri.com/library/whitepapers/pdfs/shapefile.pdf

  • Many thanks buddy, I will look into this and see how 'easy' it is for me to implement it into the script. – Joseph Dec 8 '14 at 14:09

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.