3

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

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

我有以下问题：

  1. 当我到达文件末尾时，我收到此消息：LASException: LASError in "LASReader_GetPointAt": point subscript out of range因为点数低于块维度。我不知道如何解决这个问题。
  2. a = [file_out.write(c[m]) for m in xrange(len(c))]我使用a =是为了避免视频打印。这是对的吗？
  3. 在c = [chunk[l] for l in index]我创建一个新列表中c，因为我不确定替换一个新块是明智的解决方案（例如：）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')并在读取头之后，因为我需要在 *.las 输出文件中 2）关闭文件 3）我使用inside_points = [p for p in lasfile.File(inFile,None,'r') if pnpoly(px, py, verts)]而不是

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

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

回溯（最近一次调用）：文件“”，第 1 行，在 AttributeError 中：_退出_

4

1 回答 1

3

关于（1）：

首先，为什么要使用块？只需将 lasfile 用作迭代器（如教程中所示），并一次处理一个点。以下应该通过使用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块完成时自动关闭。
  • 现在是好的部分，完成所有工作的生成器表达式（之间的部分()）。它遍历输入文件 ( for p in f)。多边形 ( ) 内的每个点都if pnpoly(p.x, p.y, verts)被添加到生成器中。
  • with我们为输出文件使用另一个块
  • 和所有的点（for p in inside_points，这是使用发电机）
  • 被写入输出文件 ( file_out.write(p))

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

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

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）：由于您不再使用原始块，您可以重新使用名称。意识到在 python 中，“变量”只是一个nametag。

关于（4）：您没有使用else，所以完全不使用它。

于 2012-10-07T14:31:27.737 回答