PCD¶

Using region of interest (ROI, e.g. plot boundary), cropping each of them from large digital orthomosaic (DOM), digital surface model (DSM), and point cloud of whole field, without using any GIS or point cloud processing software.

荷花池

包和数据准备。¶

导入 easyidp 包并使用示例的最简单方法是：

[1]:

import easyidp as idp

lotus = idp.data.Lotus()

如果您是第一次运行，它将自动从 Google Drive 下载 3.3GB 的数据集，请参考 Data 了解更多详情。

从 Shapefile 读取 ROI¶

以下代码将加载并打开 Lotus 数据集中的地块边界 shapefile，shp 文件看起来像这样（红色多边形）：

荷花池

[2]:

roi = idp.ROI(lotus.shp, name_field='plot_id')

[shp][proj] Use projection [WGS 84] for loaded shapefile [plots.shp]

Read shapefile [plots.shp]: 100%|██████████| 112/112 [00:00<00:00, 2347.62it/s]

[3]:

roi

[3]:

<easyidp.ROI> with 112 items
[0]     N1W1
array([[139.54052962,  35.73475194],
       [139.54055106,  35.73475596],
       [139.54055592,  35.73473843],
       [139.54053438,  35.73473446],
       [139.54052962,  35.73475194]])
[1]     N1W2
array([[139.54053488,  35.73473289],
       [139.54055632,  35.73473691],
       [139.54056118,  35.73471937],
       [139.54053963,  35.73471541],
       [139.54053488,  35.73473289]])
...
[110]   S4E6
array([[139.54090456,  35.73453742],
       [139.540926  ,  35.73454144],
       [139.54093086,  35.7345239 ],
       [139.54090932,  35.73451994],
       [139.54090456,  35.73453742]])
[111]   S4E7
array([[139.54090986,  35.73451856],
       [139.54093129,  35.73452258],
       [139.54093616,  35.73450504],
       [139.54091461,  35.73450107],
       [139.54090986,  35.73451856]])

另见

有关加载 shapefile 时参数的更多详细信息，请参考从 shapefile 加载 ROI，或 python API easyidp.ROI，和 easyidp.shp.read_shp。

读取和裁剪 geotiff (DOM/DSM)¶

首先，通过以下方式打开 DOM geotiff 文件：

[4]:

lotus.metashape.dom

[4]:

PosixPath('/Users/hwang/Library/Application Support/easyidp.data/2017_tanashi_lotus/170531.Lotus.outputs/170531.Lotus_dom.tif')

[5]:

dom = idp.GeoTiff(lotus.metashape.dom)

另见

GeoTiff 是 easyidp 定义的类，包含几个必需的信息。请查看 python API easyidp.GeoTiff 了解更多信息

然而，在这种情况下，ROI 和 GeoTiff 不共享相同的地理坐标：

[6]:

print("ROI.CRS: ", roi.crs.name, "\nDOM.CRS: ", dom.crs.name)

ROI.CRS:  WGS 84
DOM.CRS:  WGS 84 / UTM zone 54N

因此需要将 ROI 转换为与 GeoTiff 相同的 CRS，更多详细信息请参考从 shapefile 加载 ROI

[7]:

roi.change_crs(dom.crs)

然后使用以下函数从整个区域 GeoTIff 中裁剪每个 ROI（图）：

[8]:

dom_parts = roi.crop(dom)

Crop roi from geotiff [170531.Lotus_dom.tif]: 100%|██████████| 112/112 [00:03<00:00, 30.07it/s]

输出的 dom_parts 是一个字典，使用地块标签作为键，裁剪的 imarray 作为值：

[9]:

import matplotlib.pyplot as plt

[10]:

plt.imshow(dom_parts['N1W1'])

[10]:

<matplotlib.image.AxesImage at 0x7fbf288812e0>

../_images/jupyter_crop_outputs_21_1.png

如果您想保存裁剪的 GeoTiff，请在裁剪时传递 save_folder 参数

>>> dom_parts = roi.crop(dom, save_folder=r"expected\save\folder")

它会将所有裁剪的部分保存到带有地理偏移的 GeoTiff 文件中（您可以几乎完美地将裁剪的 DOM 与原始 DOM 重叠）

未来工作

目前不能仅将单个输出 dom_parts["N1W1"] 保存为具有正确地理位置的标准 GeoTiff 文件，而无需之前的 save_folder 批量保存，但未来将支持通过 dom_part["N1W1"].save("path\to\save\N1W1.tif") 直接保存此类文件

裁剪 DSM 的步骤与 DOM 相同，这里忽略。

读取和裁剪点云¶

点云也使用与 GeoTiff 相同的过程

另见

PointCloud 是 easyidp 定义的类，包含几个必需的信息。请查看 python API easyidp.PointCloud 了解更多信息

[11]:

pcd = idp.PointCloud(lotus.metashape.pcd)

检查点云值：

[12]:

pcd

[12]:

                   x            y        z  r    g    b        nx      ny      nz
       0  368016.142  3955509.827   94.017  129  132  147  nodata  nodata  nodata
       1  368016.142  3955509.827   94.017  129  132  147  nodata  nodata  nodata
       2  368016.142  3955509.827   94.017  129  132  147  nodata  nodata  nodata
     ...     ...          ...      ...      ...  ...  ...     ...     ...     ...
10158190  368056.169  3955481.725   97.425  174  163  149  nodata  nodata  nodata
10158191  368056.51   3955481.904   97.44   166  153  139  nodata  nodata  nodata
10158192  368056.538  3955481.486   97.447  97   112  72   nodata  nodata  nodata

并裁剪：

注意

请确保 ROI 和点云之间的 CRS 相同。

[13]:

pcd_parts = roi.crop(pcd)

Crop roi from point cloud [170531.Lotus.laz]: 100%|██████████| 112/112 [00:30<00:00,  3.61it/s]

输出的 pcd_parts 是一个字典，使用地块标签作为键，裁剪的点云作为值：

[14]:

pcd_parts["N1W1"]

[14]:

                x            y        z  r       g       b           nx      ny      nz
    0  368017.889  3955510.812   97.214  nodata  nodata  nodata  nodata  nodata  nodata
    1  368017.888  3955510.8     97.218  nodata  nodata  nodata  nodata  nodata  nodata
    2  368017.897  3955510.964   97.22   nodata  nodata  nodata  nodata  nodata  nodata
  ...     ...          ...      ...      ...     ...     ...        ...     ...     ...
35745  368020.003  3955509.546   97.558  nodata  nodata  nodata  nodata  nodata  nodata
35746  368020.046  3955509.665   97.436  nodata  nodata  nodata  nodata  nodata  nodata
35747  368019.914  3955509.974   97.444  nodata  nodata  nodata  nodata  nodata  nodata

同样，您可以传递 save_folder 参数来保存裁剪的点云

>>> pcd_parts = roi.crop(pcd, save_folder=r"expected\save\folder")

或者您可以通过以下方式仅保存一个点云：

>>> pcd_part["N1W1"].save(r"path\to\save\N1W1.ply")