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 know this kind of question has been asked, in part, before - specifically: Seeing QGIS Python Commands and: How is object created as line using createline stored in .tab file? but this is a specific plea for help!

I've lost my (ancient) version of MapInfo, and so am attempting to recreate some analysis in QGIS. I relied heavily on the following MapBasic Command:

set coordsys table "yourtable"
update "yourtable" set obj = createline(x1, y1, x2, y2)

It joined every single point to every single point in the table and retained attribute data, thus giving me the ability to subsequently analyse the file. I loved this command for its' sheer simplicity.

I'm currently attempting to teach myself some Python, but from a very low knowledge base - so, I was wondering if anyone had seen something similar scripted in Python?!?

UPDATE: Thanks to Nathan and gene, comments really helpful!Household Income movers >$104,000 per year

share|improve this question
    
I'm guessing you called this in a loop right? Otherwise all the objects would have the same line geometry. (Note: Been a while since I have done MapBasic.) –  Nathan W Apr 2 '13 at 6:00
    
Nathan, the above statement doesn't have to be called in a loop if the values X1, ..., Y2 refer to column names. In this case each record will get it's own line –  Peter Horsbøll Møller Apr 2 '13 at 6:21
    
@PeterHorsbøllMøller ahh yeah I remember now. Been a while :) –  Nathan W Apr 2 '13 at 6:24
    
Yup, it was that simple ;) –  Andrew Tice Apr 2 '13 at 22:27

2 Answers 2

Functionally this code does the same thing (in QGIS 2.0):

vl = QgsVectorLayer("LineString", "lines", "memory")
vl.startEditing()
dp = vl.dataProvider()
layer = QgsMapLayerRegistry.instance().mapLayersByName("points")[0]
dp.addAttributes(layer.pendingFields())
for feature in layer.getFeatures():
    x1, y1, x2, y2 = (feature["x1"].toFloat()[0],
                      feature["y1"].toFloat()[0],
                      feature["x2"].toFloat()[0],
                      feature["y2"].toFloat()[0])

    f = QgsFeature(layer.pendingFields())
    line = QgsGeometry.fromPolyline([QgsPoint(x1,y1), QgsPoint(x2,y2)])
    f.setGeometry(line)
    f.setAttributes(feature.attributes())
    vl.addFeature(f)
vl.commitChanges()
QgsMapLayerRegistry.instance().addMapLayer(vl)

A few things to remember:

  1. QGIS can only have one geometry type per layer hence the new memory line layer
  2. QGIS doesn't have DSL like MapBasic. You are dealing with raw API calls.

Under the hood MapInfo is taking the SQL and translating it into something like the above (although very different and more complicated) so you don't have to worry about it.

What we can do however is write a wrapper around this in order to make things easier for the future (code pending)

share|improve this answer
    
Thanks Nathan - worked well. Not an issue holding the original points as a separate layer as I'm using a different analytic process for the data they contain.... also I kind of knew that there was something horrendous and unseen going on when executing the SQL command... just didn't want to think about it! –  Andrew Tice Apr 2 '13 at 22:30
    
Well SQL is there for a reason. I have had dreams of building something like that for QGIS but there is a lot of work involved. –  Nathan W Apr 2 '13 at 22:49
    
I have considered writing a small Python DSL around all this stuff to make it easier for people. –  Nathan W Apr 3 '13 at 0:05

You can even use the original coordinates of the point shapefile

 # iterator, iterate over pairs of points in a list, a layer
`def pair_points(iterable):
     iterator = iter(iterable)
     prev = None
     item = iterator.next() 
     for next in iterator:
         yield (item,next)
         prev = item
         item = next

  `# using the iterator 
  `for pt1,pt2 in pair_points(layer.getFeatures()):
       ptA = pt1.geometry().asPoint()
       ptB = pt2.geometry().asPoint()
       print ptA, ptB, 
       .....
       line = QgsGeometry.fromPolyline([QgsPoint(ptA), QgsPoint(ptB)])
       print line.asPolyline()
       .....

Example with one of my points layers

(271927,155249) (272361,153856) line: [(271927,155249), (272361,153856)]
(272361,153856) (272689,152802) line: [(272361,153856), (272689,152802)]
....
share|improve this answer
    
Thanks gene - many different approaches to the same problem! Going to have to get cracking on teaching myself this. Knowing what is going on in a piece of code, I find, is quite a useful way of teaching. –  Andrew Tice Apr 2 '13 at 22:34
    
Although this only gets the next feature in the layer, which may or may not be the one you want to connect to. I really like the pair_points function though. –  Nathan W Apr 2 '13 at 22:52

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.