This tutorial is aims at students in geosciences who want to learn how to open and manipulate LiDAR data, and eventually integrate the data to other pipelines. The tutorial is based on the lidR package (https://cran.r-project.org/web/packages/lidR/index.html), which was created by Jean-Romain Roussel, and also contributed to by David Auty, Florien de Boissieu, Andrew Sanchez Meador Jean-Francois Bourdon and Gatziolis Demetrios.The package is distributed with a GPL-3 license.

When working in R with packages and libraries, the first thing to do is to make sure that you have installed the package on your machine. In the present case, we are installing the lidR package, that is a set of tools and algorithms to manipulate and process lidar (ALS and TLS) data.

このチュートリアルは、LiDARデータを開いて操作し、最終的には他のパイプラインにデータを統合する方法を学びたいと考えている地球科学の学生を対象としています。このチュートリアルは lidR パッケージ(https://cran.r-project.org/web/packages/lidR/index.html), which was created by Jean-Romain Roussel, and also contributed to by David Auty, Florien de Boissieu, Andrew Sanchez Meador Jean-Francois Bourdon and Gatziolis Demetrios.

パッケージやライブラリを使ってRで作業する場合、まず最初にやることは、そのパッケージがマシンにインストールされているかどうかを確認することです。今回はlidRパッケージをインストールしていますが、これはlidar(ALSとTLS)データを操作して処理するためのツールとアルゴリズムのセットです。

# The '#' sign is for the html document creation, otherwise it generates an error using the knit package.
# '#'記号はhtml文書作成のためのもので、そうでない場合はknitパッケージを使用してエラーを発生させます。

#install.packages("lidR")
#install.packages("raster") # necessary for lidR
#install.packages("sp")    # same as above

You are then all set up to start working with R on your project. For the students in my class, you are provided in the learning environment with a lidar dataset in .las or .laz format. For those outside the class, you can download for free any small subset of a lidar dataset to get started. これで、あなたのプロジェクトでRを使って作業を始める準備が整いました。私のクラスの学生には、.lasまたは.laz形式のライダデータセットが学習環境に用意されています。クラス外の学生には、ライダデータセットの任意の小さなサブセットを無料でダウンロードして、作業を始めることができます。

library(raster) # necessary for lidR
library(sp)    # necessary for lidR 
library(lidR)  # import the lidR library to your project
library(ggplot2)

setwd("C:/Users/kaiki/Desktop/DESKTOPFILES/LIDAR processing with lidR") # you then want to make sure that you are in the right working space, where your data is located, etc., so that you are all nice and cozy. Please note that this is just a setup on my computer, yours will look different. In R-studio, you can also choose your workspace in the files/more/ option.

las <- readLAS("points.las") # read the pointcloud
summary(las) # this traditional function provides information on the "metadata" of the las file. Another shorter version is print(las).
## class        : LAS (v1.0 format 1)
## memory       : 60 Mb 
## extent       : 320739.8, 321035.2, 5466478, 5466750 (xmin, xmax, ymin, ymax)
## coord. ref.  : NAD83 / UTM zone 11N 
## area         : 80373.01 m²
## points       : 786.4 thousand points
## density      : 9.78 points/m²
## File signature:           LASF 
## File source ID:           0 
## Global encoding:
##  - GPS Time Type: GPS Week Time 
##  - Synthetic Return Numbers: no 
##  - Well Know Text: CRS is GeoTIFF 
##  - Aggregate Model: false 
## Project ID - GUID:        00000000-0000-0000-0000-000000000000 
## Version:                  1.0
## System identifier:        LAStools (c) by rapidlasso GmbH 
## Generating software:      TerraScan + OT 
## File creation d/y:        364/2020
## header size:              227 
## Offset to point data:     1489 
## Num. var. length record:  3 
## Point data format:        1 
## Point data record length: 28 
## Num. of point records:    786433 
## Num. of points by return: 600793 151682 31124 2834 0 
## Scale factor X Y Z:       0.01 0.01 0.01 
## Offset X Y Z:             0 0 0 
## min X Y Z:                320739.8 5466478 746.32 
## max X Y Z:                321035.2 5466750 971.96 
## Variable length records: 
##    Variable length record 1 of 3 
##        Description: GeoTIFF GeoKeyDirectoryTag 
##        Tags:
##           Key 1024 value 1 
##           Key 1025 value 1 
##           Key 1026 value 0 
##           Key 2049 value 31 
##           Key 2054 value 9102 
##           Key 3072 value 26911 
##           Key 3076 value 9001 
##           Key 4096 value 5703 
##           Key 4097 value 37 
##           Key 4098 value 5103 
##           Key 4099 value 9001 
##    Variable length record 2 of 3 
##        Description: GeoTIFF GeoAsciiParamsTag 
##        data:                 NAD83 / UTM zone 11N + VERT_CS|NAD83|NAVD88 height| 
##    Variable length record 3 of 3 
##        Description: OGR variant of OpenGIS WKT SRS 
##        WKT OGC COORDINATE SYSTEM:  [...] (truncated)
# Above, is provided the basic reading of a all file, but you may want to filter your dataset and/or only import some variables out of the las file, and for this, there is the Filter option. You can find all the different filters by running:

readLAS(filter = "-help")
# You can generate a 3D view of the pointcloud using the plot command (will not be generated on this page) 

# plotコマンドを用いて点群の3Dビューを生成することができます (このページでは生成されません)

plot(las)

# you can also create transects between two points with the ggplot2 library:

# ggplot2 ライブラリを使って 2 点間の横断線を作成することもできます。


pt1 <- c(320740, 5466650)
pt2 <- c(321035, 5466650)
transect <- clip_transect(las, pt1, pt2, width = 2, xz = TRUE) 

# this create a transect with all the points in a 2 m radius from the file las, from point 1 to point 2.

# lasファイルから半径2mの範囲内にあるすべての点を、点1から点2までの区間を作成します.

ggplot(transect@data, aes(X,Z, color = Z)) +
geom_point(size = 0.5)  +
coord_equal()  +
theme_minimal() +
scale_color_gradientn(colours = height.colors(50))

pt1 <- c(320740, 5466478)
pt2 <- c(321035, 5466750)
transect <- clip_transect(las, pt1, pt2, width = 2, xz = TRUE) 

ggplot(transect@data, aes(X,Z, color = Z)) +
geom_point(size = 0.2)  +
coord_equal()  +
theme_minimal() +
scale_color_gradientn(colours = height.colors(100))

# in this case, we changed the size of the points and the length of the gradient color scale, as well as the position of point 1 and point 

# 2今回は、点の大きさとグラデーションの色目盛の長さ、点1と点2の位置を変更しました。
# In the present case, we reused point 1 and point 2, but this time, the transects are constructed along a line that is not 2 units, but 10 (we will see more trees and the slope across the transect)

# 今回は点1と点2を再利用していますが、今回は2単位ではなく10単位の線に沿って横断線を構成しています(より多くの樹木と横断線を挟んだ斜面を見ることになります)

transect <- clip_transect(las, pt1, pt2, width = 10, xz = TRUE) 

ggplot(transect@data, aes(X,Z, color = Z)) +
geom_point(size = 0.2)  +
coord_equal()  +
theme_minimal() +
scale_color_gradientn(colours = height.colors(10))

# You can check the change in the width of the "transect" object by plotting it using plot(transect)

# plot(transect) を用いて, "transect" オブジェクトの幅の変化をプロットすることで確認することができます.

Geomorphologists often pester against the vegetation for it is covering rocks, sediments and the morphology (even if we secretely love forest and the vegetation) we are trying to examine. For non-submarine environments, LiDAR technology has brought the possibility to separate the vegetation from the ground return (when it is visible) even under the canopy, allowing the geomorphologist to retrieve precious data. But to extract the topography from all the points and noise, one must use some algorithms. Happily, the lidR package provides such ready-algorithm to classify the pointcloud from ground and not-ground returns.

地形学者は、しばしばそれが岩石、堆積物と我々が調査しようとしている形態学（我々は密かに森林と植生を愛していても）を覆っているため、植生に対して苦言を呈します。海底以外の環境では、LiDAR技術により、樹冠の下でも植生を地上から切り離すことができるようになりました（目に見えている場合）。しかし、すべての点とノイズから地形を抽出するためには、いくつかのアルゴリズムを使用しなければなりません。幸いなことに、lidRパッケージは、地上と非地上のリターンから点群を分類するためのアルゴリズムを提供しています。

las <- classify_ground(las, algorithm = pmf(ws = 5, th = 3))
## Original dataset already contains 194957 ground points. These points were reclassified as 'unclassified' before performing a new ground classification.
# this function was created by Jean-Romain Roussel (Laval University, Quebec, CA) to combine both the transect location and plotting the data.

# この関数はJean-Romain Roussel (Laval University, Quebec, CA)によって作成されました。

plot_crossection <- function(las,
                             p1 = c(min(las@data$X), mean(las@data$Y)),
                             p2 = c(max(las@data$X), mean(las@data$Y)),
                             width = 4, colour_by = NULL)
{
  colour_by <- enquo(colour_by)
  data_clip <- clip_transect(las, p1, p2, width)
  p <- ggplot(data_clip@data, aes(X,Z)) + geom_point(size = 0.5) + coord_equal() + theme_minimal()

  if (!is.null(colour_by))
    p <- p + aes(color = !!colour_by) + labs(color = "")

  return(p)
}

plot_crossection(las, pt1 , pt2, colour_by = factor(Classification))

In order to subset the ground from the vegetation, we can also simply use the filter_ground function, that is also part of the lidR package, and then plot it, or make a cross section out of it.

植生から地面をサブセットするために、lidRパッケージの一部であるfilter_ground関数を使って、地面をプロットしたり、断面を作成したりすることができます。

Ground <- filter_ground(las)
ground_transect <- clip_transect(Ground, pt1, pt2, width = 4, xz = TRUE)
ggplot(ground_transect@data, aes(X,Z, color = Z)) + 
  geom_point(size = 0.5) + 
  coord_equal() + 
  theme_minimal() +
  scale_color_gradientn(colours = height.colors(50))

ggplot(transect@data, aes(X,Z, color = Z)) + 
  geom_point(size = 0.5) + 
  coord_equal() + 
  theme_minimal() +
  scale_color_gradientn(colours = height.colors(50))

From the two transects above, you can then see the version with or without vegetation.

You should be provided by your lecturer or teaching assistant or by myself with a basic dataset, so that you can try to reproduce this set of data by yourself.

上の2つのトランセクトから、植生がある場合とない場合のバージョンを見ることができます。

このデータセットを自分で再現できるように、講師やティーチングアシスタント、あるいは私が基本的なデータセットを用意しておく必要があります。