Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

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!

share|improve this question
add comment

3 Answers 3

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.

If I've misunderstood your question or you're having trouble with the solution feel free to contact me and I can help further.

Cheers, Joel

twitter: spatialpython

web: geospatialpython.com

share|improve this answer
add comment

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.

share|improve this answer
add comment

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.

share|improve this answer
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.