How to script shapefile editing with open source tools?

Question

I would really like not to disturb anybody and find my answer by a Google search, but I don't really know where to start.

I have a .shp vector layer with several points on streets, each having a list of attributes (about 40K). I want to create a line between certain points.

For example: I at the corner of a street, one point has an attribute "East", the one facing it an attribute "West". I want to streamline a process where (for example) every point with the attribute "West" from a particular street (included in the attribute) is connected to the corresponding "East" point.

In short, how to script modifications, based on attributes, on a .shp?

The output has to be a shapefile.

Could you please guide me to some relevant documentation please? I don't have ArcGIS.

Thanks a lot!

Answer 1

You could do that fairly easily in python using pyshp: http://code.google.com/p/pyshp/

The way I understood your question you are working with a point shapefile and turning it into a line shapefile.

The code would go something like this:

import shapefile

# Read your points shapefile

r = shapefile.Reader("streetPoints.shp")

shapes = r.shapes()

records = r.records()

# Create the output shapefile

w = shapefile.Writer(shapefile.POLYLINE)

# Loop through the attributes and create line segments.
# There's probably better logic for this but it's a start

segments = []

# Get the field names from the dbf field descriptors
fieldnames = [n[0] for n in r.fields]

st = fieldnames.index("STREET")
d = fieldnames.index("END")
opp = {"WEST":"EAST", "NORTH":"SOUTH"}

# Build the connected segments based on corresponding points
for rec1 in records:
    street = rec[st]
    e = opp[rec[d]]
    for rec2 in records:
        if street in rec2 and e in rec2:
            p1 = shapes[records.index(rec1)]
            p2 = shapes[records.index(rec2)]
            connect = [p1, p2, street]
            if connect not in segments:
                segments.append(connect)

w.field("STREET")
for seg in segments:
    # add the start and end points
    w.poly(parts=[ [ [[seg[0]], [seg[1]] ] ]
    # add the street name as the attribute
    w.record(seg[2])

# Produce a new shapefile
w.save("streetLines")

This example is overly simplistic and could be made more efficient with a little work but that's the basic idea. If this is a process you plan to repeat often it would definitely be worth some additional work and maybe even creating an executable with PyInstaller if you need to give it to other people.

Answer 2

This site has a lot of tutorials and examples on how to use GDAL/OGR. Another place to look at more complicated python examples is to look directly at the GDAL code repository.

Answer 3

I'm not sure how easy this will be for your example but I've used GeoScript for this sort of thing in the past. It gives you easy spatial operations in JavaScript, Python, Scala, and Groovy.