QGIS大栅格数据处理脚本优化及多CPU调用问题咨询

我正在处理大范围高分辨率栅格数据，现有脚本在小栅格文件上运行效果良好，但应用于大栅格文件时耗时极久。处理流程如下：

1. 从DEM栅格文件生成aspect（坡向）和slope（坡度）栅格
2. 将aspect和slope栅格转换为多边形图层
3. 导出重要的aspect与slope组合
4. 将导出图层合并为最终的“NO-GO-ZONES（禁入区）”图层

当前使用的代码如下：

import pandas as pd
from os.path import join, normpath
import time

start = time.time()

path = 'C:/Users/tlind/Dropbox/Documents/Temp/'

# iface.addRasterLayer(path+'riktigk_raster.tif')

processing.run("native:slope", {'INPUT':path+'riktig.tif','Z_FACTOR':1,'OUTPUT':path+'slope.tif'})
# iface.addRasterLayer(path+'slope.tif')
#
processing.run("native:aspect", {'INPUT':path+'riktig.tif','Z_FACTOR':1,'OUTPUT':path+'aspect.tif'})
# iface.addRasterLayer(path+'aspect.tif')

processing.run("gdal:merge", {'INPUT':['C:/Users/tlind/Dropbox/Documents/Temp/aspect.tif','C:/Users/tlind/Dropbox/Documents/Temp/slope.tif'],'PCT':True,'SEPARATE':True,'NODATA_INPUT':None,'NODATA_OUTPUT':None,'OPTIONS':'','EXTRA':'','DATA_TYPE':5,'OUTPUT':path+'merge.tif'})

processing.run("native:pixelstopolygons", {'INPUT_RASTER':path+'merge.tif','RASTER_BAND':1,'FIELD_NAME':'ASPECT','OUTPUT':path+'aspect.shp'})

processing.run("native:pixelstopolygons", {'INPUT_RASTER':path+'merge.tif','RASTER_BAND':2,'FIELD_NAME':'SLOPE','OUTPUT':path+'slope.shp'})

aspect_layer = iface.addVectorLayer(path+'aspect.shp', "", "ogr")
slope_layer = iface.addVectorLayer(path+'slope.shp', "", "ogr")

pv_apsect = aspect_layer.dataProvider()
pv_apsect.addAttributes([QgsField('ID', QVariant.Double)])

aspect_layer.updateFields()

expression = QgsExpression('$id')

context = QgsExpressionContext()
context.appendScopes(QgsExpressionContextUtils.globalProjectLayerScopes(aspect_layer))

with edit(aspect_layer):
    for f in aspect_layer.getFeatures():
        context.setFeature(f)
        f['ID'] = expression.evaluate(context)
        aspect_layer.updateFeature(f)
    
pv_slope = slope_layer.dataProvider()
pv_slope.addAttributes([QgsField('ID', QVariant.Double)])

slope_layer.updateFields()

expression = QgsExpression('$id')

context = QgsExpressionContext()
context.appendScopes(QgsExpressionContextUtils.globalProjectLayerScopes(slope_layer))

with edit(slope_layer):
    for f in slope_layer.getFeatures():
        context.setFeature(f)
        f['ID'] = expression.evaluate(context)
        slope_layer.updateFeature(f)
    
processing.run("native:joinattributestable", {'INPUT':path+'aspect.shp','FIELD':'ID','INPUT_2':path+'slope.shp','FIELD_2':'ID','FIELDS_TO_COPY':[],'METHOD':1,'DISCARD_NONMATCHING':False,'PREFIX':'','OUTPUT':path+'aspect_slope.shp'})
aspect_slope_layer = iface.addVectorLayer(path+'aspect_slope.shp', "", "ogr")

slope_list = [4.2, 4.6, 5.1, 5.7, 6.4, 7, 9.5, 12, 15, 18.2, 22.5, 30]
aspect_list_1 = [[0, 15], [15,30], [30, 45], [45, 60], [60, 75], [75, 90], [90, 105], [105, 120], [120, 135], [135, 150], [150, 165], [165, 180]]
aspect_list_2 = [[180, 195], [195, 210], [210, 225], [225, 240], [240, 255], [255, 270], [270, 285], [285, 300], [300, 315], [315, 330], [330, 345], [345, 360]]
# aspect_list_3 = aspect_list_1+aspect_list_2

for i in range(len(slope_list)):
    aspect_list_1[i].append(slope_list[i])

for i in range(len(slope_list)):
    aspect_list_2[i].append(slope_list[::-1][i])

for aspect_interval in aspect_list_1:
        start = aspect_interval[0] 
        end = aspect_interval[1]
        slope_loop = aspect_interval[2]
        aspect_slope_layer.selectByExpression('"ASPECT">'+str(start)+' and "ASPECT"<='+str(end)+' and "SLOPE">='+str(slope_loop))
        QgsVectorFileWriter.writeAsVectorFormat(aspect_slope_layer, str(path)+'aspect_slope_'+str(start)+'-'+str(end)+'.shp', "UTF-8", aspect_slope_layer.crs(), "ESRI Shapefile", onlySelected=True)
        # iface.addVectorLayer(str(path)+'aspect_slope_'+str(start)+'-'+str(end)+'.shp', "", "ogr")

for aspect_interval in aspect_list_2:
        start = aspect_interval[0] 
        end = aspect_interval[1]
        slope_loop = aspect_interval[2]
        aspect_slope_layer.selectByExpression('"ASPECT">'+str(start)+' and "ASPECT"<='+str(end)+' and "SLOPE">='+str(slope_loop))
        QgsVectorFileWriter.writeAsVectorFormat(aspect_slope_layer, str(path)+'aspect_slope_'+str(start)+'-'+str(end)+'.shp', "UTF-8", aspect_slope_layer.crs(), "ESRI Shapefile", onlySelected=True)
        # iface.addVectorLayer(str(path)+'aspect_slope_'+str(start)+'-'+str(end)+'.shp', "", "ogr")
    
processing.run("native:mergevectorlayers", {'LAYERS':['C:/Users/tlind/Dropbox/Documents/Temp/aspect_slope_0-15.shp','C:/Users/tlind/Dropbox/Documents/Temp/aspect_slope_105-120.shp','C:/Users/tlind/Dropbox/Documents/Temp/aspect_slope_120-135.shp','C:/Users/tlind/Dropbox/Documents/Temp/aspect_slope_135-150.shp','C:/Users/tlind/Dropbox/Documents/Temp/aspect_slope_15-30.shp','C:/Users/tlind/Dropbox/Documents/Temp/aspect_slope_150-165.shp','C:/Users/tlind/Dropbox/Documents/Temp/aspect_slope_165-180.shp','C:/Users/tlind/Dropbox/Documents/Temp/aspect_slope_180-195.shp','C:/Users/tlind/Dropbox/Documents/Temp/aspect_slope_195-210.shp','C:/Users/tlind/Dropbox/Documents/Temp/aspect_slope_210-225.shp','C:/Users/tlind/Dropbox/Documents/Temp/aspect_slope_225-240.shp','C:/Users/tlind/Dropbox/Documents/Temp/aspect_slope_240-255.shp','C:/Users/tlind/Dropbox/Documents/Temp/aspect_slope_255-270.shp','C:/Users/tlind/Dropbox/Documents/Temp/aspect_slope_270-285.shp','C:/Users/tlind/Dropbox/Documents/Temp/aspect_slope_285-300.shp','C:/Users/tlind/Dropbox/Documents/Temp/aspect_slope_30-45.shp','C:/Users/tlind/Dropbox/Documents/Temp/aspect_slope_300-315.shp','C:/Users/tlind/Dropbox/Documents/Temp/aspect_slope_315-330.shp','C:/Users/tlind/Dropbox/Documents/Temp/aspect_slope_330-345.shp','C:/Users/tlind/Dropbox/Documents/Temp/aspect_slope_345-360.shp','C:/Users/tlind/Dropbox/Documents/Temp/aspect_slope_45-60.shp','C:/Users/tlind/Dropbox/Documents/Temp/aspect_slope_60-75.shp','C:/Users/tlind/Dropbox/Documents/Temp/aspect_slope_75-90.shp','C:/Users/tlind/Dropbox/Documents/Temp/aspect_slope_90-105.shp'],'CRS':None,'OUTPUT':str(path)+'aspect_slope_final.shp'})
iface.addVectorLayer(str(path)+'aspect_slope_final.shp', "", "ogr")

end = time.time()
print("Elapsed time:", (end-start)/60, "minutes.")

目前生成多边形图层前的步骤速度尚可，但添加ID字段的操作耗时较长。想咨询三个问题：

• 能否在使用native:pixelstopolygons工具进行栅格转矢量时直接添加多个字段，而非后续手动添加？
• 在末尾的导出循环中，能否一次性选中所有符合条件的aspect与slope组合并仅导出一次，以提升处理效率？
• 能否让QGIS调用更多CPU核心来加速处理？目前观察到它仅使用单核心。

1. 栅格转矢量时直接生成多字段图层 #

无需分别转换aspect和slope再通过ID连接，改用gdal:polygonize工具可一次性处理多波段栅格，将每个波段的值作为单独字段输出到同一个多边形图层，省去添加ID、表连接的耗时步骤。

替换原有的native:pixelstopolygons及后续ID添加、表连接代码：

# 直接处理多波段栅格，生成包含ASPECT和SLOPE字段的矢量图层
processing.run("gdal:polygonize", {
    'INPUT': path+'merge.tif',
    'BAND': [1, 2],
    'FIELD': ['ASPECT', 'SLOPE'],
    'OUTPUT': path+'aspect_slope_direct.shp'
})
aspect_slope_layer = iface.addVectorLayer(path+'aspect_slope_direct.shp', "", "ogr")

2. 一次性筛选并导出所有符合条件的要素 #

将所有条件合并为一个表达式，一次性选中所有目标要素后导出一次，避免多次IO操作和后续图层合并步骤。

实现代码：

# 构建所有筛选条件的表达式
condition_list = []

# 处理aspect_list_1的条件
for start, end, slope_val in aspect_list_1:
    condition_list.append(f'(ASPECT > {start} AND ASPECT <= {end} AND SLOPE >= {slope_val})')

# 处理aspect_list_2的条件
for start, end, slope_val in aspect_list_2:
    condition_list.append(f'(ASPECT > {start} AND ASPECT <= {end} AND SLOPE >= {slope_val})')

# 用OR连接所有条件
full_expression = ' OR '.join(condition_list)

# 一次性选中并导出
aspect_slope_layer.selectByExpression(full_expression)
QgsVectorFileWriter.writeAsVectorFormat(
    aspect_slope_layer, 
    path+'NO-GO-ZONES.shp', 
    "UTF-8", 
    aspect_slope_layer.crs(), 
    "ESRI Shapefile", 
    onlySelected=True
)

3. 启用多核心加速处理 #

QGIS部分工具支持多线程，可通过以下方式启用：

• GDAL工具添加多线程参数：在gdal:merge中添加'-multi'，在gdal:polygonize中添加'-co NUM_THREADS=ALL_CPUS'，示例：

  # 多线程合并栅格
  processing.run("gdal:merge", {
      'INPUT': [path+'aspect.tif', path+'slope.tif'],
      'PCT': True,
      'SEPARATE': True,
      'DATA_TYPE': 5,
      'OUTPUT': path+'merge.tif',
      'EXTRA': '-multi'
  })
  
  # 多线程栅格转矢量
  processing.run("gdal:polygonize", {
      'INPUT': path+'merge.tif',
      'BAND': [1, 2],
      'FIELD': ['ASPECT', 'SLOPE'],
      'OUTPUT': path+'aspect_slope_direct.shp',
      'EXTRA': '-co NUM_THREADS=ALL_CPUS'
  })
  

• QGIS全局设置：在选项 -> 处理 -> 并行处理中勾选启用并行处理并设置最大线程数（仅部分工具支持）。
• 超大规模栅格拆分处理：将栅格切割为小块，用Pythonmultiprocessing模块并行处理后合并结果。

import pandas as pd
from os.path import join, normpath
import time

start = time.time()

path = 'C:/Users/tlind/Dropbox/Documents/Temp/'

# 生成坡度和坡向栅格
processing.run("native:slope", {'INPUT':path+'riktig.tif','Z_FACTOR':1,'OUTPUT':path+'slope.tif'})
processing.run("native:aspect", {'INPUT':path+'riktig.tif','Z_FACTOR':1,'OUTPUT':path+'aspect.tif'})

# 多线程合并栅格
processing.run("gdal:merge", {
    'INPUT': [path+'aspect.tif', path+'slope.tif'],
    'PCT': True,
    'SEPARATE': True,
    'DATA_TYPE': 5,
    'OUTPUT': path+'merge.tif',
    'EXTRA': '-multi'
})

# 多线程栅格转矢量，直接生成包含ASPECT和SLOPE的图层
processing.run("gdal:polygonize", {
    'INPUT': path+'merge.tif',
    'BAND': [1, 2],
    'FIELD': ['ASPECT', 'SLOPE'],
    'OUTPUT': path+'aspect_slope_direct.shp',
    'EXTRA': '-co NUM_THREADS=ALL_CPUS'
})

aspect_slope_layer = iface.addVectorLayer(path+'aspect_slope_direct.shp', "", "ogr")

# 定义条件列表
slope_list = [4.2, 4.6, 5.1, 5.7, 6.4, 7, 9.5, 12, 15, 18.2, 22.5, 30]
aspect_list_1 = [[0, 15], [15,30], [30, 45], [45, 60], [60, 75], [75, 90], [90, 105], [105, 120], [120, 135], [135, 150], [150, 165], [165, 180]]
aspect_list_2 = [[180, 195], [195, 210], [210, 225], [225, 240], [240, 255], [255, 270], [270, 285], [285, 300], [300, 315], [315, 330], [330, 345], [345, 360]]

# 为aspect列表添加坡度阈值
for i in range(len(slope_list)):
    aspect_list_1[i].append(slope_list[i])
for i in range(len(slope_list)):
    aspect_list_2[i].append(slope_list[::-1][i])

# 构建完整筛选表达式
condition_list = []
for start, end, slope_val in aspect_list_1:
    condition_list.append(f'(ASPECT > {start