การปรับปรุงประสิทธิภาพเมทริกซ์ระยะทาง PyQGIS

Oct 09 2020

ฉันได้เขียนเครื่องมือประมวลผลแบบกำหนดเองซึ่งฉันใช้ภายในโมเดลกราฟิกขนาดใหญ่ใน QGIS ฉันกำลังมองหาวิธีเร่งความเร็ว

ขั้นตอนแรกของเครื่องมือประมวลผลนี้คือการเชื่อมต่อทุกจุดของเลเยอร์อินพุตหนึ่งกับทุก ๆ จุดของเลเยอร์อินพุตอื่นและวัดความยาวบรรทัด ใช้เวลาไม่นานเพราะโดยปกติทั้งสองชั้นจะมีมากกว่า 10,000 คะแนน แต่ฉันสนใจเฉพาะคู่จุดในระยะห่างที่กำหนดซึ่งกันและกัน ตอนนี้ฉันกำลังจะลบบรรทัดในภายหลังซึ่งยาวเกินไป ฉันคิดว่ามันค่อนข้างไม่ได้ผลดังนั้นจึงอยากถามว่ามีวิธีง่ายๆในการเร่งขั้นตอนนี้หรือไม่โดยใช้ดัชนีเชิงพื้นที่หรืออะไรทำนองนั้น ฉันสามารถสร้างได้โดยใช้index = QgsSpatialIndex(layer.getFeatures())แต่ฉันจะเชื่อมต่อเฉพาะจุดที่ฉันสนใจได้อย่างไร

นี่เป็นส่วนน้อยที่สุดของรหัสปัจจุบันของฉัน:

Memorylayer_VL = QgsVectorLayer(str('LineString?crs=25832')), "tmp_lines", "memory")
Memorylayer_PR = Memorylayer_VL.dataProvider()
maxdistance = 150

# iterate over stop features and create lines
for stop_feat in stop_layer.getFeatures():
    point1 = QgsPoint(stop_feat.geometry().asPoint())
    for source_feat in possibility_layer.getFeatures():
        point2 = QgsPoint(source_feat.geometry().asPoint())
        temp_feat = QgsFeature(fields)
        temp_feat.setGeometry(QgsGeometry.fromPolyline([point1, point2]))
        temp_feat["line_length"] = temp_feat.geometry().length()    
        # Only add lines shorter than maxdistance and longer than 1m (Already simplified like this thanks to @BERAs suggestion in comments) 
        # However, couldn't using a spatial index even have more effect?
        if (temp_feat.geometry().length() < maxdistance) and (tmp_line_feat.geometry().length() > 1):       
            Memorylayer_PR.addFeature(temp_feat)
            Memorylayer_VL.updateFields()
Memorylayer_VL.commitChanges()

# Doing more processing stuff on the just created lines layer....

คำตอบ

3 J.Monticolo Oct 09 2020 at 12:50

ด้านล่างรหัสสำหรับเลือกรายการโซลูชันทั้งหมดระหว่าง 0 (ยกเว้น) และmaxdistance(ยกเว้นด้วย) POSSIBILITYคุณสมบัติทั้งหมดถูกแทรกไว้ในดัชนีเชิงพื้นที่ ( QgsSpatialIndex) ตัวกรองซึ่งต้องขอบคุณดัชนีเชิงพื้นที่โดยมีการใช้maxdistanceสี่เหลี่ยมรอบคุณสมบัติของเลเยอร์STOPเพื่อเพิ่มความเร็วในกระบวนการและเลือก (ตัดกัน) เฉพาะPOSSIBILITYคุณสมบัติในสี่เหลี่ยมจัตุรัส

หากเลเยอร์STOPและPOSSIBILITYเหมือนกันโค้ดด้านล่างจะซ้ำกันและไม่มีประสิทธิภาพมากนักเนื่องจากซ้ำกัน

qgis_project = QgsProject.instance()
layers = {lyr.name(): lyr for lyr in qgis_project.mapLayers().values()}
lyr_from = layers["STOP"]
lyr_to = layers["POSSIBILITY"]
maxdistance = 150

spatial_idx_to = QgsSpatialIndex()
pr_from = lyr_from.dataProvider()
pr_to = lyr_to.dataProvider()
feats_from = pr_from.getFeatures()
feats_to = pr_to.getFeatures()

# insert destination features to the spatial index
spatial_idx_to.addFeatures(feats_to)

result_lines = []
# iterate over origin features
for ft_from in feats_from:
    # find all the nearest feature thanks to the spatial index
    geom_from = ft_from.geometry().asPoint()
    x_from = geom_from.x()
    y_from = geom_from.y()
    # create a rectangle at maxdistance of the from point
    rect = QgsRectangle(
        QgsPointXY(x_from - maxdistance, y_from - maxdistance),
        QgsPointXY(x_from + maxdistance, y_from + maxdistance)
    )
    # get all points in that intersects the rectangle
    distance_ids = spatial_idx_to.intersects(rect)
    for fid in distance_ids:
        geom_to = lyr_to.getFeature(fid).geometry().asPoint()
        # if the distance between origin and destination < maxdistance
        if 0 < geom_from.distance(geom_to) < maxdistance:
            line_feat = QgsFeature()
            # create a PolylineXY from origin and destination PointXY
            line_feat.setGeometry(QgsGeometry.fromPolylineXY([geom_from, geom_to]))
            # append the line feature to the result list
            result_lines.append(line_feat)

# create the memory layer
memlyr = QgsVectorLayer("LineString?crs=epsg:25832", "tmp_lines", "memory")
memlyr_pr = memlyr.dataProvider()
# add all resulting lines features
memlyr_pr.addFeatures(result_lines)
# finally add the line memory layer to the project
qgis_project.addMapLayer(memlyr)