19. مكتبة تحليل الشبكات

تلميح

تحتاج مقتطفات الشيفرة في هذه الصفحة إلى عمليات الاستيراد التالية إذا كنت خارج وحدة تحكم pyqgis:

from qgis.core import (
  QgsVectorLayer,
  QgsPointXY,
)

يمكن استخدام مكتبة تحليل الشبكات من أجل:

  • إنشاء مخطط رياضي (graph) من البيانات الجغرافية (طبقات متجهية من الخطوط المتعددة)

  • تطبيق الطرق الأساسية من نظرية المخططات (حاليًا خوارزمية Dijkstra فقط)

باختصار، يمكن وصف حالة الاستخدام النموذجية كما يلي:

  1. إنشاء مخطط من البيانات الجغرافية (عادةً طبقة متجهية من الخطوط المتعددة)

  2. تشغيل تحليل المخطط

  3. استخدام نتائج التحليل (مثل عرضها مرئيًا)

19.1. بناء المخطط

أول ما عليك فعله --- هو تجهيز بيانات الإدخال، أي تحويل طبقة متجهية إلى مخطط. وستستخدم جميع الإجراءات اللاحقة هذا المخطط لا الطبقة.

يمكننا استخدام أي طبقة متجهية من الخطوط المتعددة مصدرًا. تصبح عُقد الخطوط المتعددة رؤوسًا في المخطط، وتصبح مقاطع الخطوط المتعددة أضلاعًا (edges) فيه. وإذا كان لعدة عُقد الإحداثيات نفسها فإنها تمثّل رأسًا واحدًا في المخطط. لذا فإن خطين يشتركان في عقدة يصبحان متصلين أحدهما بالآخر.

بالإضافة إلى ذلك، يمكن أثناء إنشاء المخطط "تثبيت" ("ربط") أي عدد من النقاط الإضافية بالطبقة المتجهية المدخلة. ولكل نقطة إضافية يُعثر على مطابق --- أقرب رأس في المخطط أو أقرب ضلع فيه. وفي الحالة الأخيرة يُقسَم الضلع ويُضاف رأس جديد.

يمكن استخدام البيانات الوصفية للطبقة المتجهية وطول الضلع خصائصَ للضلع.

يتم التحويل من طبقة متجهية إلى مخطط باستخدام النمط البرمجي Builder. ويُبنى المخطط باستخدام ما يُسمّى الموجِّه (Director). ولا يوجد حاليًا سوى موجِّه واحد: QgsVectorLayerDirector. يحدد الموجِّه الإعدادات الأساسية التي ستُستخدم لبناء مخطط من طبقة متجهية خطية، ويستخدمها الباني (builder) لإنشاء المخطط. وحاليًا، كما في حالة الموجِّه، يوجد بانٍ واحد فقط: QgsGraphBuilder، الذي ينشئ كائنات QgsGraph. وقد ترغب في تطبيق بُناة خاصين بك لبناء مخطط متوافق مع مكتبات مثل BGL أو NetworkX.

لحساب خصائص الأضلاع يُستخدم النمط البرمجي الاستراتيجية (strategy). ولا تتوفر حاليًا سوى الاستراتيجية QgsNetworkDistanceStrategy (التي تأخذ طول المسار في الحسبان) و QgsNetworkSpeedStrategy (التي تراعي السرعة أيضًا). ويمكنك تطبيق استراتيجيتك الخاصة التي تستخدم جميع المعاملات اللازمة.

حان وقت الخوض في العملية.

  1. أولًا، لاستخدام هذه المكتبة يجب استيراد وحدة التحليل:

    from qgis.analysis import *

  2. ثم بعض الأمثلة على إنشاء موجِّه:

    1# Don't use information about road direction from layer attributes,
2# all roads are treated as two-way
3director = QgsVectorLayerDirector(
4    vectorLayer, -1, "", "", "", QgsVectorLayerDirector.DirectionBoth
5)

    1# Use field with index 5 as source of information about road direction.
2# one-way roads with direct direction have attribute value "yes",
3# one-way roads with reverse direction have the value "1", and accordingly
4# bidirectional roads have "no". By default roads are treated as two-way.
5# This scheme can be used with OpenStreetMap data
6director = QgsVectorLayerDirector(
7    vectorLayer, 5, "yes", "1", "no", QgsVectorLayerDirector.DirectionBoth
8)

    لبناء موجِّه، يجب أن نمرّر طبقة متجهية ستُستخدم مصدرًا لبنية المخطط، ومعلومات عن الحركة المسموح بها على كل مقطع طريق (حركة في اتجاه واحد أو في الاتجاهين، اتجاه مباشر أو معاكس). ويبدو الاستدعاء هكذا (تجد مزيدًا من التفاصيل عن المعاملات في qgis.analysis.QgsVectorLayerDirector):

    1director = QgsVectorLayerDirector(
2    vectorLayer,
3    directionFieldId,
4    directDirectionValue,
5    reverseDirectionValue,
6    bothDirectionValue,
7    defaultDirection,
8)

  3. ثم يلزم إنشاء استراتيجية لحساب خصائص الأضلاع

    1# The index of the field that contains information about the edge speed
2attributeId = 1
3# Default speed value
4defaultValue = 50
5# Conversion from speed to metric units ('1' means no conversion)
6toMetricFactor = 1
7strategy = QgsNetworkSpeedStrategy(attributeId, defaultValue, toMetricFactor)

  4. وإبلاغ الموجِّه بهذه الاستراتيجية

    director = QgsVectorLayerDirector(vectorLayer, -1, "", "", "", 3)
director.addStrategy(strategy)

  5. يمكننا الآن استخدام الباني الذي سينشئ المخطط، باستخدام مُنشئ الصنف QgsGraphBuilder.

    # only CRS is set, all other values are defaults
builder = QgsGraphBuilder(vectorLayer.crs())

  6. ويمكننا أيضًا تعريف عدة نقاط ستُستخدم في التحليل. على سبيل المثال:

    startPoint = QgsPointXY(1179720.1871, 5419067.3507)
endPoint = QgsPointXY(1180616.0205, 5419745.7839)

  7. أصبح كل شيء جاهزًا الآن، لذا يمكننا بناء المخطط و"ربط" هذه النقاط به:

    tiedPoints = director.makeGraph(builder, [startPoint, endPoint])

    قد يستغرق بناء المخطط بعض الوقت (بحسب عدد المعالم في الطبقة وحجم الطبقة). والمتغير tiedPoints قائمة بإحداثيات النقاط "المربوطة".

  8. عند انتهاء عملية البناء يمكننا الحصول على المخطط واستخدامه في التحليل:

    graph = builder.graph()

  9. باستخدام الشيفرة التالية يمكننا الحصول على فهارس رؤوس نقاطنا:

    startId = graph.findVertex(tiedPoints[0])
endId = graph.findVertex(tiedPoints[1])

19.2. تحليل المخطط

يُستخدم تحليل الشبكات للإجابة عن سؤالين: أي الرؤوس متصلة، وكيف نجد أقصر مسار. ولحل هاتين المسألتين توفّر مكتبة تحليل الشبكات خوارزمية Dijkstra.

تجد خوارزمية Dijkstra أقصر طريق من أحد رؤوس المخطط إلى جميع الرؤوس الأخرى، وقيم معاملات التحسين. ويمكن تمثيل النتائج على هيئة شجرة أقصر المسارات.

شجرة أقصر المسارات مخطط موجَّه موزون (أو بدقة أكبر شجرة) له الخصائص التالية:

  • رأس واحد فقط ليس له أضلاع واردة — وهو جذر الشجرة

  • جميع الرؤوس الأخرى لها ضلع وارد واحد فقط

  • إذا كان الوصول إلى الرأس B ممكنًا من الرأس A، فإن المسار من A إلى B هو المسار الوحيد المتاح وهو الأمثل (الأقصر) في هذا المخطط

للحصول على شجرة أقصر المسارات استخدم الدالتين shortestTree() و dijkstra() من الصنف QgsGraphAnalyzer. ويُوصى باستخدام الدالة dijkstra() لأنها أسرع وتستخدم الذاكرة بكفاءة أكبر.

تفيد الدالة shortestTree() عندما تريد التجوّل في شجرة أقصر المسارات. وهي تنشئ دائمًا كائن مخطط جديدًا (QgsGraph) وتقبل ثلاثة متغيرات:

  • source --- المخطط المُدخل

  • startVertexIdx --- فهرس النقطة على الشجرة (جذر الشجرة)

  • criterionNum --- رقم خاصية الضلع المراد استخدامها (يبدأ من 0).

tree = QgsGraphAnalyzer.shortestTree(graph, startId, 0)

للدالة dijkstra() الوسائط نفسها، لكنها تُرجع صفًّا (tuple) من مصفوفتين:

  • في المصفوفة الأولى، يحتوي العنصر n على فهرس الضلع الوارد، أو -1 إذا لم تكن هناك أضلاع واردة.

  • في المصفوفة الثانية، يحتوي العنصر n على المسافة من جذر الشجرة إلى الرأس n، أو DOUBLE_MAX إذا تعذّر الوصول إلى الرأس n من الجذر.

(tree, cost) = QgsGraphAnalyzer.dijkstra(graph, startId, 0)

فيما يلي شيفرة بسيطة جدًا لعرض شجرة أقصر المسارات باستخدام المخطط المُنشأ بالدالة shortestTree() أو الدالة dijkstra() (حدّد طبقة خطية في اللوحة Layers واستبدل الإحداثيات بإحداثياتك).

تحذير

استخدم هذه الشيفرة مثالًا فقط، فهي تنشئ كثيرًا من كائنات QgsRubberBand وقد تكون بطيئة مع مجموعات البيانات الكبيرة.

 1from qgis.core import *
 2from qgis.gui import *
 3from qgis.analysis import *
 4from qgis.PyQt.QtCore import *
 5from qgis.PyQt.QtGui import *
 6
 7vectorLayer = QgsVectorLayer(
 8    "testdata/network.gpkg|layername=network_lines", "lines"
 9)
10director = QgsVectorLayerDirector(
11    vectorLayer, -1, "", "", "", QgsVectorLayerDirector.DirectionBoth
12)
13strategy = QgsNetworkDistanceStrategy()
14director.addStrategy(strategy)
15builder = QgsGraphBuilder(vectorLayer.crs())
16
17pStart = QgsPointXY(1179661.925139, 5419188.074362)
18tiedPoint = director.makeGraph(builder, [pStart])
19pStart = tiedPoint[0]
20
21graph = builder.graph()
22
23idStart = graph.findVertex(pStart)
24
25tree = QgsGraphAnalyzer.shortestTree(graph, idStart, 0)
26
27i = 0
28while i < tree.edgeCount():
29    rb = QgsRubberBand(iface.mapCanvas())
30    rb.setColor(Qt.red)
31    rb.addPoint(tree.vertex(tree.edge(i).fromVertex()).point())
32    rb.addPoint(tree.vertex(tree.edge(i).toVertex()).point())
33    i = i + 1

19.2.1. إيجاد أقصر المسارات

لإيجاد المسار الأمثل بين نقطتين يُتّبع النهج التالي. تُربط كلتا النقطتين (نقطة البداية A ونقطة النهاية B) بالمخطط عند بنائه. ثم نبني باستخدام الدالة shortestTree() أو dijkstra() شجرة أقصر المسارات بحيث يكون جذرها نقطة البداية A. ونجد في الشجرة نفسها نقطة النهاية B، ثم نبدأ التجوّل في الشجرة من النقطة B إلى النقطة A. ويمكن كتابة الخوارزمية كاملة كما يلي:

1assign T = B
2while T != B
3    add point T to path
4    get incoming edge for point T
5    look for point TT, that is start point of this edge
6    assign T = TT
7add point A to path

عند هذه المرحلة يكون لدينا المسار على هيئة قائمة معكوسة من الرؤوس (تُسرد الرؤوس بترتيب معكوس من نقطة النهاية إلى نقطة البداية) التي ستُزار أثناء السير على هذا المسار.

فيما يلي شيفرة نموذجية لوحدة تحكم Python في QGIS (قد تحتاج إلى تحميل طبقة خطية وتحديدها في جدول المحتويات واستبدال الإحداثيات في الشيفرة بإحداثياتك) تستخدم الدالة shortestTree() أو dijkstra():

 1from qgis.core import *
 2from qgis.gui import *
 3from qgis.analysis import *
 4
 5from qgis.PyQt.QtCore import *
 6from qgis.PyQt.QtGui import *
 7
 8vectorLayer = QgsVectorLayer(
 9    "testdata/network.gpkg|layername=network_lines", "lines"
10)
11director = QgsVectorLayerDirector(
12    vectorLayer, -1, "", "", "", QgsVectorLayerDirector.DirectionBoth
13)
14strategy = QgsNetworkDistanceStrategy()
15director.addStrategy(strategy)
16
17builder = QgsGraphBuilder(vectorLayer.sourceCrs())
18
19startPoint = QgsPointXY(1179661.925139, 5419188.074362)
20endPoint = QgsPointXY(1180942.970617, 5420040.097560)
21
22tiedPoints = director.makeGraph(builder, [startPoint, endPoint])
23tStart, tStop = tiedPoints
24
25graph = builder.graph()
26idxStart = graph.findVertex(tStart)
27
28tree = QgsGraphAnalyzer.shortestTree(graph, idxStart, 0)
29
30idxStart = tree.findVertex(tStart)
31idxEnd = tree.findVertex(tStop)
32
33if idxEnd == -1:
34    raise Exception("No route!")
35
36# Add last point
37route = [tree.vertex(idxEnd).point()]
38
39# Iterate the graph
40while idxEnd != idxStart:
41    edgeIds = tree.vertex(idxEnd).incomingEdges()
42    if len(edgeIds) == 0:
43        break
44    edge = tree.edge(edgeIds[0])
45    route.insert(0, tree.vertex(edge.fromVertex()).point())
46    idxEnd = edge.fromVertex()
47
48# Display
49rb = QgsRubberBand(iface.mapCanvas())
50rb.setColor(Qt.green)
51
52# This may require coordinate transformation if project's CRS
53# is different from layer's CRS
54for p in route:
55    rb.addPoint(p)

19.2.2. مناطق الإتاحة

منطقة الإتاحة للرأس A هي المجموعة الجزئية من رؤوس المخطط التي يمكن الوصول إليها من الرأس A، بحيث لا تتجاوز تكلفة المسارات من A إلى هذه الرؤوس قيمة معينة.

يمكن توضيح ذلك أكثر بالمثال التالي: "توجد محطة إطفاء. ما أجزاء المدينة التي يمكن لسيارة الإطفاء الوصول إليها خلال 5 دقائق؟ 10 دقائق؟ 15 دقيقة؟". الإجابات عن هذه الأسئلة هي مناطق الإتاحة لمحطة الإطفاء.

لإيجاد مناطق الإتاحة يمكننا استخدام الدالة dijkstra() من الصنف QgsGraphAnalyzer. ويكفي مقارنة عناصر مصفوفة التكلفة بقيمة محددة مسبقًا. فإذا كانت cost[i] أصغر من القيمة المحددة أو تساويها، فإن الرأس i يقع داخل منطقة الإتاحة، وإلا فهو خارجها.

المسألة الأصعب هي الحصول على حدود منطقة الإتاحة. فالحد الأدنى هو مجموعة الرؤوس التي لا يزال الوصول إليها ممكنًا، والحد الأعلى هو مجموعة الرؤوس التي لا يمكن الوصول إليها. والأمر في الواقع بسيط: إنه حد الإتاحة المبني على أضلاع شجرة أقصر المسارات التي يمكن الوصول إلى رأسها المصدر ولا يمكن الوصول إلى رأسها الهدف.

إليك مثالًا:

 1director = QgsVectorLayerDirector(
 2  vectorLayer, -1, "", "", "", QgsVectorLayerDirector.DirectionBoth
 3)
 4strategy = QgsNetworkDistanceStrategy()
 5director.addStrategy(strategy)
 6builder = QgsGraphBuilder(vectorLayer.crs())
 7
 8
 9pStart = QgsPointXY(1179661.925139, 5419188.074362)
10delta = iface.mapCanvas().getCoordinateTransform().mapUnitsPerPixel() * 1
11
12rb = QgsRubberBand(iface.mapCanvas())
13rb.setColor(Qt.green)
14rb.addPoint(QgsPointXY(pStart.x() - delta, pStart.y() - delta))
15rb.addPoint(QgsPointXY(pStart.x() + delta, pStart.y() - delta))
16rb.addPoint(QgsPointXY(pStart.x() + delta, pStart.y() + delta))
17rb.addPoint(QgsPointXY(pStart.x() - delta, pStart.y() + delta))
18
19tiedPoints = director.makeGraph(builder, [pStart])
20graph = builder.graph()
21tStart = tiedPoints[0]
22
23idStart = graph.findVertex(tStart)
24
25(tree, cost) = QgsGraphAnalyzer.dijkstra(graph, idStart, 0)
26
27upperBound = []
28r = 1500.0
29i = 0
30tree.reverse()
31
32while i < len(cost):
33    if cost[i] > r and tree[i] != -1:
34        outVertexId = graph.edge(tree [i]).toVertex()
35        if cost[outVertexId] < r:
36            upperBound.append(i)
37    i = i + 1
38
39for i in upperBound:
40    centerPoint = graph.vertex(i).point()
41    rb = QgsRubberBand(iface.mapCanvas())
42    rb.setColor(Qt.red)
43    rb.addPoint(QgsPointXY(centerPoint.x() - delta, centerPoint.y() - delta))
44    rb.addPoint(QgsPointXY(centerPoint.x() + delta, centerPoint.y() - delta))
45    rb.addPoint(QgsPointXY(centerPoint.x() + delta, centerPoint.y() + delta))
46    rb.addPoint(QgsPointXY(centerPoint.x() - delta, centerPoint.y() + delta))