This is not a beginners book on how to learn Python. This is the book you need, being the bookish type, if you use Python every day and would like to learn to finer details of how to use the standard modules effectively. This is a great book to just browse and start reading anywhere. I promise you that you will learn a lot.

I haven’t taken a whole lot of similar classes but enough to know which ones are a failure, and which ones are worth your time. The one for FME Desktop 2011 was one of the best I’ve taken. In fact, I couldn’t help myself and had to rave to the friends at Safe about it. So let me quote from my own email:

First of all, let me share with you that I found the FME Desktop class to be top-notch ! I have taken a few similar classes or workshops in the past and this one stands out. Not only did I get the impression that Robyn really knows her stuff never mind the occasional FATAL ERROR   but the format of the class (30 minutes teaching, 10 minutes hands-on), plus the quality of the manual, made for a real good package that I would recommend to anyone interested in FME. Well done ! —  That said, I think FME is fantastic. I think I mentioned to both of you that I have pulled my hair out trying to wrap my brain around ESRI’s ModelBuilder more than once. FME makes more sense to me, is a lot more intuitive to my thinking and the work flows I deal with …

So, there you have it. If you haven’t taken a look at FME and deal with lots of ETL (extraction-transformation-loading) procedures, then this may be for you. Funny thing is (and I think both instructor and manual shared this with us), they say, that 90% of all FME transformations consist of shapefile >>> shapefile conversion e.g. reprojection, adding attributes… What does that say about ArcGIS user friendliness or annoyance ? I’ve certainly tried to use ModelBuilder more than once, and every time I found that is was much easier and straight forward to just write things out in Python. With FME, I felt that – while maybe not easy to use – it’s very powerful and much more intuitive.

The table has the following headers:

UWI | NAME | MD | TVD | EWOFFSET | NSOFFSET | DIP | AZM | TVDSS | XPATH | YPATH | SURNAME | SURSTAT

and provides a row of data for each survey depth, with UWI(API) repeating for all rows belonging to the same well.

To read this into Python, I did this:

import csv

### Reference to Input CSV File
infile  = open('c:/temp/dirEXPORT.csv', "rb")
reader = csv.reader(infile)

rownum = 0
well = ""
Wells = []
points = []

### Read CSV File, Line by Line
for row in reader:

    if rownum == 0: # First line in file: read headers
        colcount = 0
        for column in row:
            print column,
            if column == "UWI":
                uwi = colcount
            if column == "XPATH":
                xpath = colcount
            if column == "YPATH":
                ypath = colcount
            if column == "MD":
                z = colcount
            colcount = colcount + 1
            rownum=+1

    else:

        print "\n"

        ### Check if row (line) is for the same well as row above
        ### if so, read coordinates from field and assign to xyz as list
        if row[uwi] == well:
            xyz = [row[xpath],row[ypath],row[z]]
            points.append(xyz)
            print row[uwi], xyz
            rownum=+1

        ### If row does not belong to the same well as the the last row,
### then append coordinate list and well name to Wells list, reset
### coordinates list (=[]), and assign the coordinates to the new well
        else:
            Wells.append([well,points])
            print "New Well"
            points = []
            well = row[uwi]
            xyz = [row[xpath],row[ypath],row[z]]
            points.append(xyz)
            print row[uwi], xyz
            rownum=+1
infile.close

Any remaining print statements allow me to watch what’s going on when the program is run. It’s my crude way of debugging code.
Once all the data is loaded into a list, I start some arcpy magic:

import arcpy
from arcpy import env
import os

env.overwriteOutput = True
env.workspace = "C:/temp"

cnt = 0

mypath = env.workspace
outputFile = r"wormtracks.shp"
template = r"c:\temp\template.shp"
fieldLength = 25

sr = arcpy.Describe(template).spatialReference

arcpy.CreateFeatureclass_management(mypath,outputFile,"POLYLINE","","","",sr)
arcpy.AddField_management(outputFile, 'API', "TEXT", "", "", fieldLength) # add API as new table field

rows = arcpy.InsertCursor(outputFile)
array_container = arcpy.Array()

for well in Wells[1:]: # Ignore first item in Wells list. Poor programming where an empty list is created.
    well_api = well[0] # First item in list is the API

    wellList = well[1] # 2nd item in list is a list of XYZ points

    for pt in wellList:
        point_object = arcpy.Point()
        point_object.X = float(pt[0])
        point_object.Y = float(pt[1])

        print pt[0]
        print pt[1]
        array_container.append(point_object) # I couldn't figure out the difference between array.append and array.add
        print point_object
        del point_object

    row = rows.newRow()
    row.shape = arcpy.Polyline(array_container) # Do not use arcpy.Multipoint which created a really messy polyline w/points ouf of sequence.
    row.API = well_api
    rows.insertRow(row)
    array_container.removeAll()

del row #unlock row
del rows #unlock table

del arcpy
del outputFile

Everywhere on the web, you find people developing new GIS tools for/in Python. But I’ve never really come across a book that ties all this together. All the standard Open Source GIS books and manuals (paper and online) have chapters about Python.

All the major packages have Python support, e.g. Mapserver (Python Mapscript) or the Quantum GIS Python Cookbook for QGIS. There a numerous other projects and packages.

So when I saw this new title: Python Geospatial Development, I got excited. It covers everything from basic GIS concepts to sources of geospatial/GIS data and imagery. It covers MySQL, PostGIS, and Spatialite, mapnik, and finally even setting up a Geodjango application.

This is probably not the kind of book to read from cover to cover if you’ve already worked with some of these projects, but it covers a wide range of available tools and is packed with online resources and links to even more tools and resources. I will try to remember to elaborate once I’ve read a little more.

.

I briefly looked at Dracones, because it’s written in Python. But I quickly got the sense that setup was less than straight forward. It sound like some of the adventures I had with getting any useful results out of (Geo)Django. Ka-Map appears to be pretty popular. But I didn’t like the look of their standard template, and wanted something to work out of the box.

So I turned to pmapper. It comes in a download-ready MS4W package that installs super easily – no problems at all if you’ve got a functioning version of MS4W version installed. I found it helpful to first get the sample project up and running – using the demodata provided. This works with a mapfile and some shapefiles of European countries and cities. Depending where you install (copy ‘n paste) your pmapper directory, you might have to adjust some of the Apache directives in httpd.conf. But if you’ve been through this for the MS4W installation, you won’t have any troubles.

Once you start using your own data, it might be easiest to work with the existing directory structure and slowly start changing/recycling the demodata mapfile. Anyway, now I have a great looking template with a somewhat dynamic legend, pan and continuous zoom, and more. The categories & groups references in config/config_default.xml confused me initially because groups actually refers to GROUP names in the MAPFILE. But with some trial and error, you get all your layers displayed the way you want.

 There  are still a few bugs, though, to be ironed out.

1. The feature search still refuses to work properly.
2. I think I may have made a mistake with my reference map, because the box highlighting my zoom area only adjusts in the north-south directions, but nevre east-west, and coincidentally, my panning now only works north-south but not east-west.

Since pmapper uses PHP and Javascript neither of which I understand well, the challenge now lies in adding some Python Mapscript capabilities so that I can extend the framework myself. I’m thinking about a little drop down input box where you select SHP data that’s added dynamically… I will get to that when I have time.

To find out where you are:

import os
WorkingDirectory = os.getcwd()

To change directories:

import os
os.chdir(“NewWorkingDirectory”)

Really useful is the glob module if you want to read the contents of a directory into a list, and maybe iterate over that list…

import glob
for file in glob.glob(“*.*”):
    print file # print all file names in dir

or

import glob
for shp in glob.glob(“*.shp”):
    print shp # print names of all shapefiles in dir

Obviously, there is much more, e.g. os.path. See here (http://effbot.org/librarybook/os.htm) and here (http://docs.python.org/library/os.path.html).

The script certainly is shorter and looks more straight forward than other Python examples I had looked at. But as Tim Schaub observes, Darren’s Windows script (which I adapted for my purposes) is much more light weight.

Python bindings for OGR are included in the installation for FWTools, and provided your PYTHONPATH in Windows is set up properly, all you need to get going is:

import ogr

For other installations, e.g. OSGEO4W, you might need from osgeo import ogr

There is a script here (page 22) for merging shapefiles. But it works by iterating through first files and then each feature to write each one into a new shapefile. This took forever ! I guess what I was looking for was a program or way to simply take one shapefile and tack it to the “end” of the next without crawling through one field at a time and copying it over to the next file.

Then I read about shapefile.py on Joel Lawhead’s geospatialpython.com. This is a really neat module that allows you to access features and features parts/points in shapefile directly.

import shapefile

r = shapefile.Reader(“yourfile.shp“)
features = r.shapes()
print “Number of features: “,len(features)
cnt=1
for point in features[yourchoice].points:
   print “Coordinates for Point “,cnt, point
   cnt=cnt+1

Works like a char. But again, to simply aggregate many shapefile in one SHP, this would take to long. So for now Darren’s little Windows script to batch process the SHP’s with ogr2ogr has been the best way. I haven’t really used Python for such purposes but I guess you could.

Using OGR, you can merge 2 files like this:

ogr2ogr merge.shp file1.shp
ogr2ogr -update -append merged.shp file2.shp -nln merge

You create a new shapefile merged and copy the contents of file1 into the new shapefile. Then you append the contents of file2 to the newly created shapefile and name it merged.

What do you do if you don’t want to do this 25 times ?

On page 22 of this document (Using GDAL / OGR for Data Processing and Analysis) is a Python script for doing this. But I couldn’t get the Python bindings for GDAL to work. The errors were always with importing _gdal which I understand is the gdal dll and references to _gdal are created by SWIG. I don’t care about the script. I can write my own. But I would like to get GDAL(OGR) working from Python.

So after giving up on this concept temporarily, I resorted to Windows scripting, as modified from this page:

as above, use: ogr2ogr merge.shp file1.shp to create a shapefile merge containing the data of file1

then:

for %f in (*.shp) do (
ogr2ogr -update -append merge.shp %f  -f “esri shapefile” -nln merge )
That worked.

So I’ve cobbled together a simple HTML template with a simple CSS (for the design of a few div’s) and map controls for zoom/pan and adding data layers. Essentially this was based on any number of Mapserver Tutorials, e.g. Mitchell’s Web Mapping book or Beginning Mapserver. Initially, I was having some problems trying to reference the correct directories for Apache and mapserv. For example, I found that the href=”…” for the css <link> didn’t work as a relative link. It had to be an absolute link to the location of the CSS file even if that was located in the same directory as the HTML template.

Also, I remembered the hard way that while you can spread your <form> over multiple <div>’s you can’t do this and expect it to work:

<div1>
   <form>
</div1>
<div2>
   </form>
</div2>                 YIKES !!!

As you can see, right now, I have 3 zoom controls, and (for trial purposes) 4 data layers – counties and interstates to help you zoom to your area of interest, and roads and geology for more detail – that can be turned on and off. As soon as I get around to it, I will add capability to hide certain layers at large scales and add labeling for layers at zoom scales where it makes sense. Then, I would like add a tool that allows the user to add project specific data points, e.g. boreholes stored in SHP files, by selecting from a list of project numbers. In an effort to stay clear of ArcGIS, I will look into using PostGIS for this.

It took me a long time re-create some of the color/hatch patterns for the geology map which I have been typically been producing through an ArcGIS style file. This was good practice for a first Mapserver project, but there may be a better/faster/more accurate way to re-produce this map effects. I remember reading somewhere     that creating maps in Quantum GIS and then saving your styles in mapfile format is one way to go. Here is someone trying to use ArcGIS styles in QGIS.

In the process of wanting to make this an all Open Source project, I had to resort to OGR once again to combine public roads data from 25 counties into one shapefile for Mapserver. Im pretty sure there is a way to use some kind of tiling for this but I haven’t gotten that far. Anyway… new lessons learned with OGR int he next post.