Python:改进逐 block 代码以读取数百万点的建议

Question

我用 Python 写了一段代码来读取 *.las 文件。 *las 文件是特殊的 ascii 文件，其中每一行都是 x,y,z 点的值。

我的函数读取 N。点数并检查它们是否在具有 points_inside_poly 的多边形内。

我有以下问题:

1. 当我到达文件末尾时，我收到此消息:LASException:“LASReader_GetPointAt”中的 LASError:点下标超出范围 因为点数低于 block 维度。我不知道如何解决这个问题。
2. a = [file_out.write(c[m]) for m in xrange(len(c))] 我使用 a = 以避免视频打印.是否正确？
3. 在 c = [chunk[l] for l in index] 中，我创建了一个新列表 c 因为我不确定替换新 block 是否是明智的解决方案(例如:chunk = [chunk[l] for l in index])。
4. 在语句if else...else 中我使用pass。这是正确的选择吗？

真的很感谢你的帮助。重要的是要从专家那里改进听力建议!!!!

import shapefile
import numpy
import numpy as np
from numpy import nonzero
from liblas import file as lasfile
from shapely.geometry import Polygon
from matplotlib.nxutils import points_inside_poly  


# open shapefile (polygon)
sf = shapefile.Reader(poly)
shapes = sf.shapes()
# extract vertices
verts = np.array(shapes[0].points,float)

# open las file
f = lasfile.File(inFile,None,'r') # open LAS
# read "header"
h = f.header

# create a file where store the points
file_out = lasfile.File(outFile,mode='w',header= h)


chunkSize = 100000
for i in xrange(0,len(f), chunkSize):
    chunk = f[i:i+chunkSize]

    x,y = [],[]

    # extraxt x and y value for each points
    for p in xrange(len(chunk)):
        x.append(chunk[p].x)
        y.append(chunk[p].y)

    # zip all points 
    points = np.array(zip(x,y))
    # create an index where are present the points inside the polygon
    index = nonzero(points_inside_poly(points, verts))[0]

    # if index is not empty do this otherwise "pass"
    if len(index) != 0:
        c = [chunk[l] for l in index] #Is It correct to create a new list or i can replace chunck?
        # save points
        a = [file_out.write(c[m]) for m in xrange(len(c))] #use a = in order to avoid video print. Is it correct?
    else:
        pass #Is It correct to use pass?

f.close()
file_out.close()

@Roland Smith 提出并由 Gianni 更改的代码

f = lasfile.File(inFile,None,'r') # open LAS
h = f.header
# change the software id to libLAS
h.software_id = "Gianni"
file_out = lasfile.File(outFile,mode='w',header= h)
f.close()
sf = shapefile.Reader(poly) #open shpfile
shapes = sf.shapes()
for i in xrange(len(shapes)):
    verts = np.array(shapes[i].points,float)
    inside_points = [p for p in lasfile.File(inFile,None,'r') if pnpoly(p.x, p.y, verts)]
    for p in inside_points:
        file_out.write(p)
f.close()
file_out.close()

我使用了这些解决方案:1) 读取 f = lasfile.File(inFile,None,'r') 并在读取 head 之后因为我需要在 *.las 输出文件中2)关闭文件3) 我使用 inside_points = [p for p in lasfile.File(inFile,None,'r') if pnpoly(p.x, p.y, verts)] 而不是

with lasfile.File(inFile, None, 'r') as f:
...     inside_points = [p for p in f if pnpoly(p.x, p.y, verts)]
...     

因为我总是收到这个错误信息

追溯(最近一次调用):文件“”，第 1 行，位于属性错误:_退出_

Answer 1

关于 (1):

首先，为什么要使用 block ？只需将 lasfile 用作迭代器(如 tutorial 所示)，并一次处理一个点。下面应该使用 pnpoly 将多边形内的所有点写入输出文件在列表推导式而不是 points_inside_poly 中运行。

from liblas import file as lasfile
import numpy as np
from matplotlib.nxutils import pnpoly

with lasfile.File(inFile, None, 'r') as f:
    inside_points = (p for p in f if pnpoly(p.x, p.y, verts))
    with lasfile.File(outFile,mode='w',header= h) as file_out:
        for p in inside_points:
            file_out.write(p)

上面的五行应该替换整个大的 for 循环。让我们一一回顾:

• with lasfile.File(inFile...:使用此构造意味着当 with block 完成时文件将自动关闭。
• 现在是最重要的部分，完成所有工作的生成器表达式(() 之间的部分)。它遍历输入文件(for p in f)。多边形内的每个点(if pnpoly(p.x, p.y, verts))都被添加到生成器中。
• 我们为输出文件使用另一个with block
• 和所有点(for p in inside_points，这是使用了生成器)
• 被写入输出文件(file_out.write(p))

因为此方法仅将多边形内部的点添加到列表中，所以您不会在不需要的点上浪费内存!

如果上面显示的方法不起作用，您应该仅使用 block 。使用 block 时，您应该正确处理异常。例如:

from liblas import LASException

chunkSize = 100000
for i in xrange(0,len(f), chunkSize):
    try:
        chunk = f[i:i+chunkSize]
    except LASException:
        rem = len(f)-i
        chunk = f[i:i+rem]

关于 (2):抱歉，我不明白您要在这里完成什么。 “视频打印”是什么意思？

关于 (3):由于您不再使用原始 block ，因此可以重新使用该名称。意识到在 python 中，“变量”只是一个 nametag .

关于 (4):您没有使用 else，因此请完全忽略它。