ArcPy+多处理:错误:无法保存栅格数据集 [英] arcpy + multiprocessing: error: could not save Raster dataset

查看:10
本文介绍了ArcPy+多处理:错误:无法保存栅格数据集的处理方法,对大家解决问题具有一定的参考价值,需要的朋友们下面随着小编来一起学习吧!

问题描述

我正在使用arcpy对栅格进行一些数字运算,希望使用multiprocessing包来加快速度。基本上,我需要遍历元组列表,使用每个元组进行一些栅格计算,并将一些输出写入文件。我的输入由一个数据栅格(水深测量)、一个定义区域的栅格和一个由两个浮点(水面高程、深度)组成的元组组成。我的过程由一个函数computeplane组成,该函数获取一个元组并运行一系列栅格计算以生成五个栅格(总栅格、滨海栅格、表层栅格、次表层栅格和深栅格栅格),然后为每个栅格调用processtable函数以使用arcpy.sa.ZonalStatisticsAsTable将值写入DBF,使用arcpy.AddField_management添加一些字段,使用arcpy.TableToTable_conversion将DBF转换为CSV,最后使用arcpy.Delete_management删除DBF文件。

基于一些other SO posts,我将我的代码包装在main()中,以便multiprocessing.Pool可以很好地发挥作用。我使用main()创建元组集,使用pool.map进行多处理。我使用tempfile包为DBF文件选择名称以避免名称冲突;CSV文件名保证不会与其他线程冲突。

我已经使用for循环测试了我的代码,它工作得很好,但当我尝试使用pool.map时,我收到

运行时错误:错误010240:无法将栅格数据集保存到 C:UsersMichaelAppDataLocalESRIDesktop10.4SpatialAnalystlesst_ras 使用输出格式网格。

这里发生了什么?这个错误没有出现在代码的非多进程版本中,我也不会写出任何地方的栅格-话又说回来,我不知道arcpy如何处理中间栅格(我肯定不认为它会让它们留在内存中-它们太大了)。我需要告诉arcpy它处理栅格计算的方式才能使多处理工作吗?我已在下面包含了我的python文件。

import arcpy
arcpy.CheckOutExtension("Spatial")
import arcpy.sa
import numpy
import multiprocessing
import tempfile
bathymetry_path = r'C:/GIS workspace/RRE/habitat.gdb/bathymetry_NGVD_meters'
zones_path = r'C:/GIS workspace/RRE/habitat.gdb/markerzones_meters'
table_folder = r'C:/GIS workspace/RRE/zonetables'
bathymetry = arcpy.sa.Raster(bathymetry_path)
zones = arcpy.sa.Raster(zones_path)

def processtable(raster_obj, zone_file, w, z, out_folder, out_file):
  temp_name = "/" + next(tempfile._get_candidate_names()) + ".dbf"
  arcpy.sa.ZonalStatisticsAsTable(zone_file, 'VALUE', raster_obj, out_folder + temp_name, "DATA", "SUM")
  arcpy.AddField_management(out_folder + temp_name, "wse", 'TEXT')
  arcpy.AddField_management(out_folder + temp_name, "z", 'TEXT')
  arcpy.CalculateField_management(out_folder + temp_name, "wse", "'" + str(w) + "'", "PYTHON")
  arcpy.CalculateField_management(out_folder + temp_name, "z", "'" + str(z) + "'", "PYTHON")
  arcpy.TableToTable_conversion(out_folder + temp_name, out_folder, out_file)
  arcpy.Delete_management(out_folder + temp_name)

def computeplane(wsedepth):
  wse = wsedepth[0]
  depth = wsedepth[1]
  total = bathymetry < depth
  littoral = total & ((wse - bathymetry) < 2)
  surface = total & ~(littoral) & ((total + wse - depth) < (total + 2))
  profundal = total & ((total + wse - depth) > (total + 5))
  subsurface = total & ~(profundal | surface | littoral)
  # zonal statistics table names
  total_name = 'total_w' + str(wse) + '_z' + str(depth) + '.csv'
  littoral_name = 'littoral_w' + str(wse) + '_z' + str(depth) + '.csv'
  surface_name = 'surface_w' + str(wse) + '_z' + str(depth) + '.csv'
  subsurface_name = 'subsurface_w' + str(wse) + '_z' + str(depth) + '.csv'
  profundal_name = 'profundal_w' + str(wse) + '_z' + str(depth) + '.csv'
  # compute zonal statistics
  processtable(total, zones, wse, depth, table_folder, total_name)
  processtable(littoral, zones, wse, depth, table_folder, littoral_name)
  processtable(surface, zones, wse, depth, table_folder, surface_name)
  processtable(profundal, zones, wse, depth, table_folder, profundal_name)
  processtable(subsurface, zones, wse, depth, table_folder, subsurface_name)

def main():
  watersurface = numpy.arange(-15.8, 2.7, 0.1)
  # take small subset of the tuples: watersurface[33:34]  
  wsedepths = [(watersurface[x], watersurface[y]) for x in range(watersurface.size)[33:34] for y in range(watersurface[0:x+1].size)]
  pool = multiprocessing.Pool()
  pool.map(computeplane, wsedepths)
  pool.close()
  pool.join()

if __name__ == '__main__':
  main()

更新

更多的调查表明,这与其说是multiprocessing的问题,不如说是ArcGIS处理栅格的方式的问题。栅格代数结果被写入到默认工作空间中的文件中;在我的例子中,我没有指定文件夹,因此arcpy将栅格写入到某个AppData文件夹。ArcGIS根据您的代数表达式使用诸如Lessth、Lessth_1等基本名称。因为我没有指定工作区,所以所有multiprocessing线程都写入此文件夹。虽然单个arcpy进程可以跟踪名称,但多个进程都试图写入相同的栅格名称并撞上其他进程的锁。

我尝试在computeplane开头创建一个随机工作区(文件gdb),然后在结尾将其删除,但arcpy通常不会及时释放其锁定,并且进程在DELETE语句上崩溃。所以我不确定如何继续。

推荐答案

我也花了很大力气来制作一个可以工作的脚本。

由于您的数据与我习惯的数据略有不同,我无法测试答案,也不知道它是否有效,但可以说一些东西来打开您的思维,使您能够进行多进程处理,并提出一个新的脚本供您测试。

我看到有些东西可以并行化,有些东西不能并行化,或者从多进程中几乎得不到任何好处。

使用多进程和arcpy的好做法是不使用地理数据库,尝试写入实际的‘.tif’文件或‘in_Memory/Tempfile’,然后在计算后删除,在要并行化的函数内部。不要一开始就为空间而烦恼,减少要测试的数据,当它奏效时,通过删除临时文件来改善它,在中间放一些有用的打印。

回到问题上,脚本中有一些东西是不能使用的,比如传递一个arcpy.Raster对象,而不是传递给栅格的路径,并在函数中读取它。栅格的路径是一个字符串,因此它是可拾取的,这意味着它们可以轻松地传递到一个新进程,必须特别注意这一点。

另一件事是声明一个类似于‘TABLE_FOLDER’的变量,并期望您的所有7核内核都返回到第一个内核中的原始进程,以询问它的路径。 多进程之间不容易共享变量,对于原始代码,您必须创建一个独立运行的函数。主模块向空闲的CPU发送的内容是创建一个新的Python进程,其中几乎只有函数声明、传递给它的参数和必要的导入模块,所有这些都准备好执行,它不能回头询问。

对我起作用的还是在一个非常切分的函数(不是我创建的)中使用apply_async,该函数使用进程作为键,apply_async的结果作为值来创建字典,然后将其放入try/except中,以便如果一个进程发生错误,它不会停止主进程的执行。 大概是这样的:

def computeplane_multi(processList):
    pool = Pool(processes=4, maxtasksperchild=10)
    jobs= {}
    for item in processList:
        jobs[item[0]] = pool.apply_async(computeplane, [x for x in item])
    for item,result in jobs.items():
        try:
            result = result.get()
        except Exception as e:
            print(e)
    pool.close()
    pool.join()

回到您的代码,我对您的尝试(并告诉我)做了一些修改:

import arcpy
import numpy
import multiprocessing
import tempfile

bathymetry_path = r'C:/GIS workspace/RRE/habitat.gdb/bathymetry_NGVD_meters'
zones_path = r'C:/GIS workspace/RRE/habitat.gdb/markerzones_meters'
table_folder = r'C:/GIS workspace/RRE/zonetables'

def computeplane(bathymetry_path,zones_path,wsedepth):
    def processtable(raster_obj, zones_path, w, z, out_folder, out_file):
        zone_file = arcpy.sa.Raster(zones_path)
        temp_name = "/" + next(tempfile._get_candidate_names()) + ".dbf"
        arcpy.sa.ZonalStatisticsAsTable(zone_file, 'VALUE', raster_obj, out_folder + temp_name, "DATA", "SUM")
        arcpy.AddField_management(out_folder + temp_name, "wse", 'TEXT')
        arcpy.AddField_management(out_folder + temp_name, "z", 'TEXT')
        arcpy.CalculateField_management(out_folder + temp_name, "wse", "'" + str(w) + "'", "PYTHON")
        arcpy.CalculateField_management(out_folder + temp_name, "z", "'" + str(z) + "'", "PYTHON")
        arcpy.TableToTable_conversion(out_folder + temp_name, out_folder, out_file)
        arcpy.Delete_management(out_folder + temp_name)
    bathymetry = arcpy.sa.Raster(bathymetry_path)
    wse = wsedepth[0]
    depth = wsedepth[1]
    total = bathymetry < depth
    littoral = total & ((wse - bathymetry) < 2)
    surface = total & ~(littoral) & ((total + wse - depth) < (total + 2))
    profundal = total & ((total + wse - depth) > (total + 5))
    subsurface = total & ~(profundal | surface | littoral)
    # zonal statistics table names
    total_name = 'total_w' + str(wse) + '_z' + str(depth) + '.csv'
    littoral_name = 'littoral_w' + str(wse) + '_z' + str(depth) + '.csv'
    surface_name = 'surface_w' + str(wse) + '_z' + str(depth) + '.csv'
    subsurface_name = 'subsurface_w' + str(wse) + '_z' + str(depth) + '.csv'
    profundal_name = 'profundal_w' + str(wse) + '_z' + str(depth) + '.csv'
    # compute zonal statistics
    processtable(total, zones_path, wse, depth, table_folder, total_name)
    processtable(littoral, zones_path, wse, depth, table_folder, littoral_name)
    processtable(surface, zones_path, wse, depth, table_folder, surface_name)
    processtable(profundal, zones_path, wse, depth, table_folder, profundal_name)
    processtable(subsurface, zones_path, wse, depth, table_folder, subsurface_name)
    print('point processed : {},{}'.format(wsedepth[0],wsedepth[1]))


def computeplane_multi(processList):
    pool = Pool(processes=4, maxtasksperchild=10)
    jobs= {}
    for item in processList:
        jobs[item[0]] = pool.apply_async(computeplane, [x for x in item])
    for item,result in jobs.items():
        try:
            result = result.get()
        except Exception as e:
            print(e)
    pool.close()
    pool.join()


def main():
    watersurface = numpy.arange(-15.8, 2.7, 0.1)
    # take small subset of the tuples: watersurface[33:34]  
    wsedepths = [(watersurface[x], watersurface[y]) for x in range(watersurface.size)[33:34] for y in range(watersurface[0:x+1].size)]
    processList = []
    for i in wsedepths:
        processList.append((bathymetry_path,zones_path,i))
    computeplane_multi(processList)

if __name__ == '__main__':
    main()

这篇关于ArcPy+多处理:错误:无法保存栅格数据集的文章就介绍到这了,希望我们推荐的答案对大家有所帮助,也希望大家多多支持IT屋!

查看全文
登录 关闭
扫码关注1秒登录
发送“验证码”获取 | 15天全站免登陆