如何使用liblas获取存储在las point的"数据"字段中的内容？

Question

我正在使用多脉冲激光雷达数据，该数据收集飞行路径内多条直线上的点。我正在尝试确定las文件中各个行的名称和数量。我在python中使用liblas模块。

我找到了这个文档，它解释了存储在las文件中的不同字段。它在页面的最底部提到了一个数据字段（get_data和set_data）。

标头中的"点数据格式"one_answers"点数据记录长度"为此"数据"字段留出了空间。我的标题说我为数据字段预留了28个字节，数据字段中存储了28个值。第19个值（至少在来自两个不同传感器的两个数据集中）指的是行号。我在单脉冲数据中有一个值，在多脉冲数据中是4。

我想知道这个领域中存储的内容是否有一个标准，或者它是否是专有的。

此外，作为获取每条扫描线名称的一种方法，我编写了以下代码：

import liblas
from liblas import file as lasfile
# Get parameters
las_file = r"E:Testing0101.las"
f = lasfile.File(las_file, mode='r')
line_list = []
counter = 0
for p in f:
    line_num = p.data[18]
    if line_num not in line_list:
        line_list.append(line_num)
    counter += 1
print line_list

结果出现以下错误：

Traceback (most recent call last):
  File "D:ToolsPython_Scriptspoint_info.py", line 46, in <module>
    line_num = p.data[18]
  File "C:Python27ArcGIS10.1libsite-packagesliblaspoint.py", line 560, in get_data
    length = self.header.data_record_length
  File "C:Python27ArcGIS10.1libsite-packagesliblaspoint.py", line 546, in get_header
    return header.Header(handle=core.las.LASPoint_GetHeader(self.handle))
WindowsError: [Error -529697949] Windows Error 0xE06D7363

有人知道更多关于las点/表头中存储的行号的信息吗？有人能解释这个错误吗？在我得到错误之前，它似乎分配了将近2gb的内存。我在winxp上，所以我猜这是一个内存错误，但我不明白为什么访问这个"数据"字段会占用内存。非常感谢您的帮助。

Answer 1

我不会假装自己是这方面的专家，但我对GIS数据很感兴趣，所以这引起了我的兴趣。我在Fedora 19系统上安装了liblas及其依赖项，并使用了liblas附带的示例数据文件。

使用你的代码，我遇到了同样的问题，看着我所有的记忆都被消耗殆尽。我不知道为什么会发生这种情况——也许是不需要的引用阻碍了垃圾收集器的工作，正如我们所希望的那样。这可能是可以修复的，但我不会尝试

我确实注意到了liblas模块的一些有趣的功能，并决定尝试一下。我相信你能得到你想要的数据。

打开文件后，请查看标题中的XML描述。

h = f.get_header()
print(h.get_xml())

这很难查看（可以随意使用xml.dom.minidom或lxml.etree），但在我的示例文件中，它显示了点数据的字节布局（我的也有28个字节）。在我的例子中，偏移量18是分配给点源ID的一个短（2字节）。您应该能够使用p.data[18:19]、p.get_data()[18:19]、p.point_source_id或p.get_point_source_id()检索到它。不幸的是，data引用占用了内存，p.point_source_id有一个bug（提交给开发人员的bug修复拉取请求）。如果我们将您的代码更改为使用最后一种访问方法，那么一切似乎都很正常。所以，在你的for循环中试试这个：

for p in f:
    line_num = p.get_point_source_id()
    if line_num not in line_list:
        line_list.append(line_num)
    counter += 1

注意

counter == h.get_count()

如果您只想要一组唯一的点源ID值。。。

line_set = set(p.get_point_source_id() for p in f)

希望您的数据值也可用作p.get_point_source_id()。请在评论中告诉我它是如何为您工作的。干杯