Python - 在光栅单元上循环，速度极慢

不是最精明的python用户，我一直在试图找到一种解决方案来加速我的代码。我有2个光栅文件，都具有相同的尺寸和范围。一个栅格是一条从激光雷达图像中提取的河流，其高程值范围很窄，河流文件中的所有其他值都是0.所以基本上，除了河流宽度（具有高程数据的像素）之外，99％河流栅格的值为0.另一个栅格文件是从激光雷达图像中提取的，并在整个范围内具有高程数据。Python - 在光栅单元上循环，速度极慢

我的计划是检查没有0值的河流文件中的每个像素，将其加1（如果河水水位上升1米），并将其与相邻像素高程数据（基本上LIDAR标高的像素位置）。如果它更高，给它一个值1（这意味着这个像素被淹没）。我想出了一些代码，但它看起来很可怕，而且非常缓慢（油漆干得比代码运行得更快）。所以我正在寻找一种加快速度的方法。我看矢量化，numpy等，但还不能理解它（可能是因为我已经在工作超过9小时..）

所以，在这里，它的任何建议将非常多赞赏：

river = arcpy.Raster(r'C:\Flood.gdb\Clip_River1') 
lidar = arcpy.Raster(r'C:\Flood.gdb\Clip_Lidar_1') 

arrayRiver = arcpy.RasterToNumPyArray(river,nodata_to_value=0) 
(rHeight, rWidth)=arrayRiver.shape 

arrayLidar = arcpy.RasterToNumPyArray(lidar,nodata_to_value=0) 
(lHeight, lWidth)=arrayLidar.shape 

# extent of the 2 rasters: columns 3822, rows 10129 

flood = 1 
while flood == 1: 
    flood = 0 
    for row in range(0,rHeight-1): 
     for col in range(0,rWidth-1): 
      if arrayRiver.item(row,col) <> 0 and arrayRiver[row,col] <> 1: 
       if arrayLidar[row,col-1] <> 1 and arrayLidar.item(row,col-1) < arrayRiver.item(row,col)+1: 
        arrayLidar[row,col-1] = 1 
        arrayRiver[row,col-1] = arrayRiver.item(row,col) 
        flood = 1 

       [....doing the same concept for each possible neighboring pixel] 

       if arrayLidar[row+1,col+1] <> 1 and arrayLidar.item(row+1,col+1) < arrayRiver.item(row,col)+1: 
        arrayLidar[row+1,col+1] = 1 
        arrayRiver[row+1,col+1] = arrayRiver.item(row,col) 
        flood = 1 

      arrayRiver[row,col] = 1 

newRaster = arcpy.NumPyArrayToRaster(arrayLidar,lowerLeft,cellSize,value_to_nodata=0) 
newRaster.save(r"C:\Flood.gdb\newRastaRiver_small") 

newEmptyRaster = arcpy.NumPyArrayToRaster(arrayEmpty,lowerLeft,cellSize,value_to_nodata=0) 
newEmptyRaster.save(r"C:\Flood.gdb\newEmptyRastaRiver_small")

来源

2015-04-30 derBrain

在做任何事情之前，你应该做一个分析，抛出整个if语句，看看你的代码在本地生效了多长时间，然后你就可以看到主要瓶颈在哪里。而且你没有添加单元格的尺寸，所以它不清楚瓶颈可能来自哪里。 – user1767754

我会问这个gis.stackexchange.com。一般来说，不要循环访问数组。作为一个例子，这里有一个最近的问题，我问了每个元素的所有邻居没有循环：http://stackoverflow.com/questions/29925684/generalize-stacking-of-array-elements-neighbors-into-3 -d-array –

好的，我删除了整个if语句，整个代码（包括栅格到数组语句）在8.2秒内运行。光栅的尺寸是3822列，10129行 – derBrain

这只是一个部分的答案，但应该让你开始。

使用numpy排除（掩码）低值的所有像素。

arrayRiver = arcpy.RasterToNumPyArray(river,nodata_to_value=0) 
# Create a any array makring pixels less than zero. 
mask = arrayRiver > 0 
arrayRiver = arrayRiver[mask] 
arrayLidar = arrayLidar[mask]

你arrayRiver & arrayLidar现在要小很多，和你的代码的其余部分应该会更快。顺便说一句，你的示例代码似乎没有做一个arrayLidar数组。

来源

2015-04-30 22:05:02

另外，欢迎来到Stack Overflow。 –

感谢您的输入Don。我试过你的建议，但看起来 'arrayRiver = arrayRiver [mask]' line会创建一个列表（如果我错了，请纠正我），并且我无法访问原始位置（行，列）它再次与lidarArray [row，col]标高比较。 – derBrain

Python - 在光栅单元上循环，速度极慢

回答

相关问题