9

I'm trying to recreate an Origin-Destination plot like this:

enter image description here

I have managed to munge the data into a MSOA to LAD table and can draw a map like this for one of the origin MSOA.

enter image description here

Which once you allow for the (now ridiculous) distances people in the Peak District commute to work is close.

But I quite like the effect that the author has achieved by "splaying" out the lines. Obviously, with flows of 522 and 371, I can't go for a single line per commuter but it would be nice to produce a proportional arc of lines to show the number of people making the trip.

I thought that I would be able to use the Geometry Generator but without a loop construct, I can't seem to make headway.

  • This ESRI tool may be of interest to you, or at least a spring board for code ideas on creating a "wedge" of lines. – Hornbydd Jul 31 '18 at 13:50
  • And when a line symbolizes, let's say 50 (100, 200) commuters, per line? With some python code (or the geometry generator, I'm not sure about) you could rotate the lines (x/50) with a distinct amount. – Stefan Aug 1 '18 at 6:03
5
+50

A great challenge!

This answer primarily uses the Geometry generator and was written in QGIS 3.2. QGIS crashed (without me having saved!) just after I first built the lines and I almost gave up, but the recently used expression list saved the day - another bonus to using the Geometry generator

I started with two point sets, one source and three destinations. The destinations are labelled with the counts:

Initial points

I then generated lines connecting the source point to all destinations using a virtual layer using the following code:

SELECT d.Count_MF, Makeline( s.geometry, d.geometry) 'geometry' 
  FROM Source AS s JOIN Destinations AS d

Connected points

Then I used the following Geometry generator expression to style the lines:

 intersection(
   geom_from_wkt( 
     'MULTILINESTRING ((' ||  $x_at( 0)  || ' ' || $y_at( 0)  || ', ' || 
     array_to_string(
       array_remove_at( string_to_array( regexp_replace(
             geom_to_wkt(nodes_to_points( tapered_buffer(  $geometry ,0, "Count_MF" * 200, floor("Count_MF" / 10)),true)),
             '[\\(\\)]','')),0)
     , ') , ('  ||  $x_at( 0)  || ' ' || $y_at( 0)  || ', ' )
    || '))')
    ,buffer( point_n(  $geometry ,1), $length))

This takes each line and applies the following steps:

  1. Generates a tapered buffer going from zero width at the source to a width scaled by the destination count at the destination end. The buffer point density is also scaled by the destination count attribute.
  2. The buffer polygon's vertices are converted to points (this is possibly superfluous), and then exported to WKT, and the brackets are removed using a regex, before converting into an array
  3. The array is then expanded back out to a WKT string for a multilinestring, inserting in the source point's coordinates plus relevant formatting - this creates a separate line for each of the extracted vertices connected to the source point
  4. The WKT is converted back to an geometry object and finally intersected with the buffer of the source point to clip them back to the circle that the destination point sits on (see the output of a tapered_buffer to understand why this is necessary)

Fans

In writing up the steps, I realise that the conversion to and from an array is unecessary, and all the WKT manipulation can be done with regexes. This expression is below, and if the tapered_array function can be replaced with a different one then this could also be used in QGIS 2.18.

intersection(
   geom_from_wkt(
    'MULTILINESTRING ((' ||  $x_at( 0)  || ' ' || $y_at( 0)  || ', ' ||
  replace(
    regexp_replace(
      regexp_replace(
        geom_to_wkt(tapered_buffer(  $geometry ,0, "Count_MF" * 200, floor("Count_MF" / 10))),
      '^[^,]*,',''),
    ',[^,]*$',''),
  ',',') , ('  ||  $x_at( 0)  || ' ' || $y_at( 0)  || ', ')
  || '))')
,buffer( point_n(  $geometry ,1), $length))
6

Your question made me curious.

This solution works only for QGIS 2.x in the Python Console

Like mentioned in my comment here is my idea to create the arc of lines with Python.

I have two point layer:

i. One holding the capital (id, capital)

ii. One holding the towns (id, town, commuters)

The amount of commuters are "separated into banknotes" and these will be the lines which build up the arc. So 371 commuters are a combination of 3x100, 1x50, 2x10 and 1x1 and in total 7 banknotes. Afterwards the lines are styled by a rule based styling.

Here is the code:

from qgis.gui import *
from qgis.utils import *
from qgis.core import *
from PyQt4 import QtGui, uic
from PyQt4.QtGui import *
from PyQt4.QtCore import *

for lyr in QgsMapLayerRegistry.instance().mapLayers().values():
    if lyr.name() == "capital":
        capital_layer = lyr

for lyr in QgsMapLayerRegistry.instance().mapLayers().values():
    if lyr.name() == "town":
        town_layer = lyr

    # creating the memory layer
d_lyr = QgsVectorLayer('LineString', 'distance', 'memory')
QgsMapLayerRegistry.instance().addMapLayer(d_lyr)
prov = d_lyr.dataProvider()
prov.addAttributes( [ QgsField("id", QVariant.Int), QgsField("banknote",QVariant.Int)])

    # function to create the banknotes
def banknoteOutput(number):
    number_list = []
    number_list.append(number)
    banknote_count = []
    temp_list = []
    banknote_list = []
    for n in number_list:
        total_sum = 0
        total = int(n/100)
        total_sum = total_sum + total
        if total > 0:
            banknote_count.append([total, 100])
        n = n-(total*100)
        total = int(n/50)
        total_sum = total_sum + total
        if total > 0:
            banknote_count.append([total, 50])
        n = n-(total*50)
        total = int(n/10)
        total_sum = total_sum + total
        if total > 0:
            banknote_count.append([total, 10])
        n = n-(total*10)
        total = int(n/5)
        total_sum = total_sum + total
        if total > 0:
            banknote_count.append([total, 5])
        n = n-(total*5)
        total = int(n/1)
        total_sum = total_sum + total
        if total > 0:
            banknote_count.append([total, 1])
        for i in banknote_count:
            temp_list.append(i*i[0])
        banknote_list = [item for sublist in temp_list for item in sublist][1::2]
        return banknote_list

        # creating lines with the amount of banknotes
for capital in capital_layer.getFeatures():
    for town in town_layer.getFeatures():
        commuter_splitting = banknoteOutput(town['commuters'])
        for i,banknote in enumerate(commuter_splitting):
            angle = 2
            distance = QgsDistanceArea()
            distance.measureLine(capital.geometry().asPoint(), town.geometry().asPoint())
            vect = QgsFeature()
            vect.setGeometry(QgsGeometry.fromPolyline([capital.geometry().asPoint(), town.geometry().asPoint()]))
            vect.geometry().rotate(0+(i*angle), capital.geometry().asPoint())
            vect.setAttributes([int(town["id"]), int(banknote)])
            prov.addFeatures([vect])

d_lyr.updateExtents()
d_lyr.triggerRepaint()
d_lyr.updateFields()

The result could look like this:

enter image description here

UPDATE: distinction male/female

Results in 4 Memory Layer.

from qgis.gui import *
from qgis.utils import *
from qgis.core import *
from PyQt4 import QtGui, uic
from PyQt4.QtGui import *
from PyQt4.QtCore import *

for lyr in QgsMapLayerRegistry.instance().mapLayers().values():
    if lyr.name() == "capital":
        capital_layer = lyr

for lyr in QgsMapLayerRegistry.instance().mapLayers().values():
    if lyr.name() == "town":
        town_layer = lyr

    # function to create the banknotes
def banknoteOutput(number):
    number_list = []
    number_list.append(number)
    banknote_count = []
    temp_list = []
    banknote_list = []
    for n in number_list:
        total_sum = 0
        total = int(n/100)
        total_sum = total_sum + total
        if total > 0:
            banknote_count.append([total, 100])
        n = n-(total*100)
        total = int(n/50)
        total_sum = total_sum + total
        if total > 0:
            banknote_count.append([total, 50])
        n = n-(total*50)
        total = int(n/10)
        total_sum = total_sum + total
        if total > 0:
            banknote_count.append([total, 10])
        n = n-(total*10)
        total = int(n/5)
        total_sum = total_sum + total
        if total > 0:
            banknote_count.append([total, 5])
        n = n-(total*5)
        total = int(n/1)
        total_sum = total_sum + total
        if total > 0:
            banknote_count.append([total, 1])
        for i in banknote_count:
            temp_list.append(i*i[0])
        banknote_list = [item for sublist in temp_list for item in sublist][1::2]
        return banknote_list

    # creating the male memory layer
cmt_male = QgsVectorLayer('LineString', 'Commuters_Male', 'memory')
QgsMapLayerRegistry.instance().addMapLayer(cmt_male)
prov_male = cmt_male.dataProvider()
prov_male.addAttributes( [ QgsField("id", QVariant.Int), QgsField("banknote",QVariant.Int)])

    # creating the male polygon memory layer
cmt_male_polygon = QgsVectorLayer('Polygon', 'Commuters_Male_Poly', 'memory')
QgsMapLayerRegistry.instance().addMapLayer(cmt_male_polygon)
prov_cmt_male_polygon = cmt_male_polygon.dataProvider()
prov_cmt_male_polygon.addAttributes( [ QgsField("id", QVariant.Int), QgsField("banknote",QVariant.Int)])

    # creating lines with the amount of banknotes
for capital in capital_layer.getFeatures():
    for town in town_layer.getFeatures():
        commuter_splitting = banknoteOutput(town['cmt_male'])
        points = []
        for i,banknote in enumerate(reversed(commuter_splitting)):
            angle = 2
            distance = QgsDistanceArea()
            distance.measureLine(capital.geometry().asPoint(), town.geometry().asPoint())
            vect = QgsFeature()
            vect.setGeometry(QgsGeometry.fromPolyline([capital.geometry().asPoint(), town.geometry().asPoint()]))
            vect.geometry().rotate(0+(i*angle), capital.geometry().asPoint())
            vect.setAttributes([int(town["id"]), int(banknote)])
            points.append(vect.geometry().asPolyline()[1])
            prov_male.addFeatures([vect])
        polygon = QgsFeature()
        points.insert(0,capital.geometry().asPoint())
        points.insert(len(points),capital.geometry().asPoint())
        polygon.setGeometry(QgsGeometry.fromPolygon([points]))
        polygon.setAttributes([1, 2])
        prov_cmt_male_polygon.addFeatures([polygon])

cmt_male.updateExtents()
cmt_male.triggerRepaint()
cmt_male.updateFields()
cmt_male_polygon.updateExtents()
cmt_male_polygon.triggerRepaint()
cmt_male_polygon.updateFields()

    # creating the female memory layer
cmt_female = QgsVectorLayer('LineString', 'Commuters_Female', 'memory')
QgsMapLayerRegistry.instance().addMapLayer(cmt_female)
prov_female = cmt_female.dataProvider()
prov_female.addAttributes( [ QgsField("id", QVariant.Int), QgsField("banknote",QVariant.Int)])

    # creating the female polygon memory layer
cmt_female_polygon = QgsVectorLayer('Polygon', 'Commuters_Female_Poly', 'memory')
QgsMapLayerRegistry.instance().addMapLayer(cmt_female_polygon)
prov_cmt_female_polygon = cmt_female_polygon.dataProvider()
prov_cmt_female_polygon.addAttributes( [ QgsField("id", QVariant.Int), QgsField("banknote",QVariant.Int)])

    # creating lines with the amount of banknotes
for capital in capital_layer.getFeatures():
    for town in town_layer.getFeatures():
        commuter_splitting = banknoteOutput(town['cmt_female'])
        points = []
        for i,banknote in enumerate(commuter_splitting):
            angle = 2
            distance = QgsDistanceArea()
            distance.measureLine(capital.geometry().asPoint(), town.geometry().asPoint())
            vect = QgsFeature()
            vect.setGeometry(QgsGeometry.fromPolyline([capital.geometry().asPoint(), town.geometry().asPoint()]))
            vect.geometry().rotate(-angle-(i*angle), capital.geometry().asPoint())
            vect.setAttributes([int(town["id"]), int(banknote)])
            points.append(vect.geometry().asPolyline()[1])
            prov_female.addFeatures([vect])
        polygon = QgsFeature()
        points.insert(0,capital.geometry().asPoint())
        points.insert(len(points),capital.geometry().asPoint())
        polygon.setGeometry(QgsGeometry.fromPolygon([points]))
        polygon.setAttributes([1, 2])
        prov_cmt_female_polygon.addFeatures([polygon])

cmt_female.updateExtents()
cmt_female.triggerRepaint()
cmt_female.updateFields()
cmt_female_polygon.updateExtents()
cmt_female_polygon.triggerRepaint()
cmt_female_polygon.updateFields()

The result could look like this:enter image description here

One thing that is not ideal from a cartographic point of view:

The size of an arc of line can be irritating at first glance, in the way that a bigger arc could represent more commuters. An arc can be bigger with less commuters (289 commuters / 11 banknotes) than another with more commuters (311 commuters/ 5 banknotes).

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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