How to make your first QGIS custom form on Windows

For a long time I wanted to give a go at creating a QGIS plugin. Nothing in particular, just a ‘Hello World’ plugin or similar so I could understand how the whole QGIS/python thing hangs together. Seasoned readers of my blog may have noticed by now that I mainly (only?) work on Windows OS. And let me tell you, trying to find exactly what is needed in terms of python libraries for being able to create a simple form in QGIS proved to be a big challenge. .

But I have figured it out and those are my notes in case you also want to give it a try. They don’t go as far as creating a plugin. They merely create a custom form you can use for editing a layer. But I consider this a HUGE step for me nevertheless!

I was running QGIS 1.8 on a Windows 7 64-bit laptop. I also had ArcGIS 10 installed and for Python development I use PyScripter. This posed an interesting problem of how to use PyScripter for both ArcGIS and QGIS but I had solved this earlier

So next, I apparently had to use PyQt4 for being able to create forms within QGIS. But DO NOT use the version for Python 2.7 from this link as -as I found out at my expense- its not compatible with QGIS. The safest way is to use the OSGeo4W Installer. Select the Advanced install and select to download and install the packages as shown in the next screenshot. Not entirely sure that ALL these are required but they work for me so….

py_packagesNote the last item on the list (SIP) which the version I have is the previous from the current one. I *think* you need this version otherwise you get runtime errors when trying to load the PyQt4 libraries i.e. “RuntimeError: the sip module implements API vXX to vXX but the PyQt4.QtCore module requires API vXX

Yes- it was a very trial and error job…

Moving on, you can then open PyScripter for Python 2.7 and check if you can load the PyQt4 libraries by running on the python prompt:

>>> from PyQt4.QtCore import *
>>> from PyQt4.QtGui import *

If all goes well you should not get any errors. To be on the safe side and ensure you have the same setup in the QGIS python environment (it may sound an unecessary step but trust me, its not) run the above command on the QGIS python window as well.

And so, finally, you are ready to make your custom form. But for this I simply followed Nathan’s example from his QGIS blog

Next step would be to create my Hello World plugin. On a next post….

 
Advertisements

Identifying internal polygons using Arcpy

This is just a small python script that will identify any polygons that are within other polygons in a shapefile. It has come handy for me for a project I was working on where I had to identify the “internal” polygons. The main trick here is to compare each polygon with the outer boundary of every other

# -*- coding: utf-8 -*-
#!/usr/bin/env python
import arcpy
#-------------------------------------------------------------------------------
# Name:        flag_int_polygons
# Purpose:     Identifies interior polygons
#
#-------------------------------------------------------------------------------
def main():

    output_folder = u"C:\temp"
    inFC=u"C:\temp\polys.shp"

    # Check if the compare shapefile already exists and if so, recreate it
    if arcpy.Exists(output_folder+"/tmpCompare.shp")==True:
        arcpy.Delete_management(output_folder+"/tmpCompare.shp")
    arcpy.CopyFeatures_management(inFC, output_folder+"/tmpCompare.shp")

    try:
        arcpy.AddField_management(output_folder+"/tmpCompare.shp", "LABEL", "TEXT", "", "", 2, "", "", "REQUIRED")
    except:
        pass

     # Get shape column
    shapeName = arcpy.Describe(inFC).shapeFieldName
    shapeName1 = arcpy.Describe(output_folder+"/tmpCompare.shp").shapeFieldName

     # Loop through input features
    rows=arcpy.SearchCursor(inFC)

    # Start loop through input fc
    r=1
    for row in rows:
        partnum = 0
        parcel=row.getValue(shapeName)
        exBoundary=getPolyBoundary(parcel)
        fid=row.getValue("FID")
        rows1=arcpy.UpdateCursor(output_folder+"/tmpCompare.shp")
        for row1 in rows1:
            fid1=row1.getValue("FID")
            parcel1=row1.getValue(shapeName1)
            if fid1 <> fid:
                if parcel1.within(exBoundary)==True:
                    print "Feature " + str(fid1) + " is within feature " + str(fid)
            rows1.updateRow(row1)

        del row1
        del rows1
    del row
    del rows
#-------------------------------------------------------------------------------
# Name:        getPolyBoundary
# Purpose:     Returns the external polygon boundary
#
#-------------------------------------------------------------------------------
def getPolyBoundary(geom):
    partno=0
    pntObj = arcpy.Point()
    arrayObj = arcpy.Array()

    for part in geom:
        pntCounter=0;
        # Loop through points in part
        for pnt in geom.getPart(partno):
            if pnt:
                pntObj.X=pnt.X
                pntObj.Y=pnt.Y
                arrayObj.add(pntObj)
            else:
                break
        pntCounter=pntCounter+1
    bound=arcpy.Polygon(arrayObj)
    #print str(bound)
    return bound
if __name__ == '__main__':
    main()

Configure PyScripter to use with QGIS (and still use arcpy) on Windows

As many of you may already know, PyScripter is an excellent IDE for Python. For those of you who haven’t tried it… well, you don’t know what you are missing! I was using PyScripter quite happily with ArcGIS 10 and arcpy without ever having to tweak anything in PyScripter. All was well until I decided to explore PyQGIS, QGIS’ scripting language. I could not, for the life of me, get this initial statement to work in PyScripter although it was working fine in the QGIS Python console:

import qgis.core

I was getting the “usual” error:

Traceback (most recent call last):
File “<interactive input>”, line 1, in <module>
ImportError: No module named qgis.core

Now, arcpy is using Python 2.6 whereas QGIS (I was using the latest 1.8 Lisboa version) is using 2.7. One would have thought that since PyScripter allows you to work with different versions of Python this wouldn’t be such a big deal. Yet it was.  Below I am outlining the steps I took to solve it.

One thing to note here is that I installed QGIS 1.8 using the standalone windows installer. Things may have been different if I used the OSeo4W installer which apparently sets the environment for you….but I can’t be sure.

One of the problems with the standalone installer is that it installs Python 2.7 but does not create any registry entries- which is what PyScripter mainly uses. So the first thing to do is to create the relevant registry entries for 2.7. Assuming you have installed QGIS in the default folder (C:\Program Files (x86)\Quantum GIS Lisboa for a 64bit machine), you can use this registry settings file (edit the file accordingly if using 32-bit Windows).

Next, copy the python27.dll from C:\Program Files (x86)\Quantum GIS Lisboa\bin to c:\Windows\SysWOW64 (or \system32 for 32-bit machines)

Finally, create this batch file to start pyscripter with (note you may change the PROGRA~2 to PROGRA~1 for 32-bit machines.

set OSGEO4W_ROOT=C:\PROGRA~2\QUANTU~1
call%OSGEO4W_ROOT%“\bin\o4w_env.bat
call%OSGEO4W_ROOT%“\apps\grass\grass-6.4.2\etc\env.bat
set GDAL_DRIVER_PATH=%OSGEO4W_ROOT%\bin\gdalplugins\1.9
path %PATH%;%OSGEO4W_ROOT%\apps\qgis\bin
path %PATH%;%OSGEO4W_ROOT%\apps\grass\grass-6.4.2\lib
path %PATH%;”%OSGEO4W_ROOT%\apps\Python27\Scripts\”
set PYTHONPATH=%PYTHONPATH%;%OSGEO4W_ROOT%\apps\qgis\python;
set PYTHONPATH=%PYTHONPATH%;%OSGEO4W_ROOT%\apps\Python27\Lib\site-packages
set QGISPATH=%OSGEO4W_ROOT%\apps\qgis
start “PyScripter” /B “C:\Program Files (x86)\PyScripter\PyScripter.exe” –python27

Happy coding!

Convert a Polyline shapefile to PolylineM using arcpy

I am sure I am re-inventing the wheel here but I couldn’t find any ready-made scripts to do this. Plus, I thought, it was high time I get my feet wet with Python. So- behold of my first python script(!). Its quite simple- it basically takes an input linear shapefile and converts it to a new shapefile with M values. If a field containing the measured line length is given, it will use it to calibrate the M values – otherwise it will use the map (geometry) length so the M value of the last vertex of each line will equal the length of the line.

This is the first of a bunch of ‘useful’ Avenue scripts that I used quite often in the past and converted to Python. Will eventually convert the others and learn more about Python and arcpy in the process.

The script is listed below for your viewing pleasure. It is heavily commented so it should be self explanatory. Or, you can simply download it from here

import arcpy, math, os

# Input feature class
#
inFeatures="C:/mapdata/osm_greek_data/projected/road_sample.shp"

# Fieldname containing the measured length for each line
# if null, then the map length will be used
mlField="meas_len"

# Output feature class settings
out_path = "C:/mapdata/osm_greek_data/projected"
out_name = "road_sampleM.shp"
geometry_type = "POLYLINE"
template = inFeatures
spatial_reference=inFeatures
has_m = "ENABLED"
has_z = "DISABLED"

# Create the new shapefile
arcpy.CreateFeatureclass_management(out_path, out_name, geometry_type, template, has_m, has_z, spatial_reference)

#
# Create a list of fields to use when copying field values from input to output
#
fields = arcpy.ListFields(inFeatures)
fieldList = []
for field in fields:
  fieldList.append(field.name)

# Loop through input features
rows=arcpy.SearchCursor(inFeatures)

# Output features
outFeatures=out_path+"/"+out_name

# Open insertcursor
#
outRows = arcpy.InsertCursor(outFeatures)

# Create point and array objects
#
pntObj = arcpy.Point()
arrayObj = arcpy.Array()

shapeName=arcpy.Describe(inFeatures).shapeFieldName #Get shape field name
featType=arcpy.Describe(inFeatures).shapeType       # Get feature type
isMeasured=arcpy.Describe(inFeatures).hasM  # Check if its a polylineM feature already

# Record counter
i=1

# Loop through input shapefile
for row in rows:
    print "Processing row %d " % i
    feat=row.getValue(shapeName)
    noparts=feat.partCount
    partno = 0
    currentM=0
    # Loop through geometry parts
    for part in feat:
        pntCounter=0;
        segmentDistance=0
        prevX=0
        prevY=0
        # Loop through points in part
        for pnt in feat.getPart(partno):
          if pntCounter == 0: # First point
            prevX=pnt.X
            prevY=pnt.Y
            pntObj.X=pnt.X
            pntObj.Y=pnt.Y
            pntObj.M=0
          else:
            curX=pnt.X
            curY=pnt.Y
            # Calculate segment length
            segmentDistance = math.sqrt(math.pow((curX-prevX),2) + math.pow((curY-prevY),2))
            # Use the measured distance
            if mlField <> "":
                segmentDistance=segmentDistance/feat.length*row.getValue(mlField)
            # Add previous M to segment length
            currentM=currentM + segmentDistance
            pntObj.X=pnt.X
            pntObj.Y=pnt.Y
            pntObj.M=currentM
            prevX=pnt.X
            prevY=pnt.Y
          # Add measured points to array
          arrayObj.add(pntObj)
          pntCounter=pntCounter+1
    # Create new row
    #
    outFeat=outRows.newRow()

    # Copy attributes to new row
    #
    for fieldName in fieldList:
      if fieldName not in  ("FID", "Shape"): # ignore FID and Shape fields
        outFeat.setValue(fieldName, row.getValue(fieldName))

    # Assign array of PointMs to new feature
    outFeat.shape=arrayObj

    # Insert feature
    outRows.insertRow(outFeat)

    # Clear array of points
    #
    arrayObj.removeAll()
    i=i+1

# Remove cursor from memory
del outRows
del rows