오랜만에 LiDAR 데이터 관련해서 공부한 기록을 남기게 되었다..너무나 감사하게도 우아한 형제들의 로보틱스lab 팀의 계약직으로 들어갈 수 있게 되었다. 입사하기 전에 최대한 공부하고 들어가보자! 먼저 오늘 분석에서 사용할 데이터셋은 "Autzen Stadium" LiDAR 데이터이다.
Autzen Stadium은 미국 오리건 유진에 위치한 야외 미식축구 경기장이다. 검색해보니 1967년에 완공되었고 약 54,000명의 인원을 수용할 수 있다.
1️⃣ Python 환경 준비하기
conda 환경 설정
GIS 분석을 위해 conda 환경을 사용하겠습니다. 먼저 새로운 콘다 가상환경을 생성하고 필수 라이브러리만 간단하게 설치합니다.
# 새로운 환경 생성 후 라이브러리 설치
conda create -n lidar -c conda-forge python=3.12 pdal python-pdal tiledb tiledb-py geopandas pandas numpy wget pyarrow
# 주피터 환경 설정
conda install -c conda-forge notebook
# 커널 연결해주기
conda install ipykernel
python -m ipykernel install --user --name lidar --display-name "Python (lidar)"
# 주피터 노트북 접속
jupyter notebook설치한 라이브러리는 다음과 같습니다.
`pdal`, `python-pdal`, `tiledb`, `tiledb-py`, `geopandas`, `pandas`, `numpy`, `wget`
주피터 노트북에 접속했다면 작업을 진행할 폴더, ipynb 파일을 생성해주면 된다.
파이썬 라이브러리, 데이터 불러오기
import tiledb
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import json
import subprocess
import pyarrowlidar 데이터는 pdal github에 존재하는 autzen.laz 파일을 사용하도록 하겠다.
파일 다운로드 방법 1. 직접 다운로드
링크 : https://github.com/PDAL/data/tree/main/workshop
data/workshop at main · PDAL/data
Example point cloud data for PDAL testing and evaluation - PDAL/data
github.com
해당 링크에 접속하여 autzen.laz 파일을 다운받으면 된다. 다운받은 파일은 노트북 파일이 존재하는 폴더에 넣어주면 된다.
파일 다운로드 방법 2. 셀에서 바로 다운로드(wget)
!wget "https://github.com/PDAL/data/blob/main/workshop/autzen.laz?raw=true" -O autzen.laz
다운로드가 완료되면 폴더에 autzen.laz파일이 생성되었을 것이다.
2️⃣ 데이터 확인하기
데이터가 다운로드 되었으니 PDAL을 사용해서 데이터셋 정보를 확인해보도록 하자. 아래 코드를 실행시키고 잠시 기다리면 된다.
!pdal info autzen.laz | jq .더보기
{
"file_size": 56350988,
"filename": "autzen.laz",
"now": "2026-03-20T19:51:54+0900",
"pdal_version": "2.10.0 (git-version: Release)",
"reader": "readers.las",
"stats": {
"bbox": {
"EPSG:4326": {
"bbox": {
"maxx": -123.0619741,
"maxy": 44.0627857,
"maxz": 615.26,
"minx": -123.0755564,
"miny": 44.04972421,
"minz": 406.14
},
"boundary": {
"type": "Polygon",
"coordinates": [
[
[
-123.075000938064562,
44.049724209884076,
406.14
],
[
-123.075556425524823,
44.062491619980655,
406.14
],
[
-123.062526791255408,
44.062785699800301,
615.26
],
[
-123.061974110327611,
44.050018226731467,
615.26
],
[
-123.075000938064562,
44.049724209884076,
406.14
]
]
]
}
},
"native": {
"bbox": {
"maxx": 639003.73,
"maxy": 853537.66,
"maxz": 615.26,
"minx": 635577.79,
"miny": 848882.15,
"minz": 406.14
},
"boundary": {
"type": "Polygon",
"coordinates": [
[
[
635577.790000000037253,
848882.150000000023283,
406.14
],
[
635577.790000000037253,
853537.660000000032596,
406.14
],
[
639003.729999999981374,
853537.660000000032596,
615.26
],
[
639003.729999999981374,
848882.150000000023283,
615.26
],
[
635577.790000000037253,
848882.150000000023283,
406.14
]
]
]
}
}
},
"statistic": [
{
"average": 637294.1783,
"count": 10653336,
"maximum": 639003.73,
"minimum": 635577.79,
"name": "X",
"position": 0,
"stddev": 948.0423248,
"variance": 898784.2496
},
{
"average": 851247.6953,
"count": 10653336,
"maximum": 853537.66,
"minimum": 848882.15,
"name": "Y",
"position": 1,
"stddev": 1296.496683,
"variance": 1680903.65
},
{
"average": 434.1025002,
"count": 10653336,
"maximum": 615.26,
"minimum": 406.14,
"name": "Z",
"position": 2,
"stddev": 24.67860061,
"variance": 609.033328
},
{
"average": 77.14742312,
"count": 10653336,
"maximum": 254,
"minimum": 0,
"name": "Intensity",
"position": 3,
"stddev": 62.62418437,
"variance": 3921.788468
},
{
"average": 1.17801438,
"count": 10653336,
"maximum": 4,
"minimum": 1,
"name": "ReturnNumber",
"position": 4,
"stddev": 0.4653415116,
"variance": 0.2165427224
},
{
"average": 1.358579791,
"count": 10653336,
"maximum": 4,
"minimum": 1,
"name": "NumberOfReturns",
"position": 5,
"stddev": 0.6656062384,
"variance": 0.4430316645
},
{
"average": 0.4989654884,
"count": 10653336,
"maximum": 1,
"minimum": 0,
"name": "ScanDirectionFlag",
"position": 6,
"stddev": 0.4999989533,
"variance": 0.2499989533
},
{
"average": 0,
"count": 10653336,
"maximum": 0,
"minimum": 0,
"name": "EdgeOfFlightLine",
"position": 7,
"stddev": 0,
"variance": 0
},
{
"average": 1.256686262,
"count": 10653336,
"maximum": 2,
"minimum": 1,
"name": "Classification",
"position": 8,
"stddev": 0.4368048111,
"variance": 0.190798443
},
{
"average": -0.812061405,
"count": 10653336,
"maximum": 20,
"minimum": -21,
"name": "ScanAngleRank",
"position": 9,
"stddev": 8.48431566,
"variance": 71.98361221
},
{
"average": 126.4052859,
"count": 10653336,
"maximum": 156,
"minimum": 115,
"name": "UserData",
"position": 10,
"stddev": 3.832798034,
"variance": 14.69034077
},
{
"average": 7329.903705,
"count": 10653336,
"maximum": 7334,
"minimum": 7326,
"name": "PointSourceId",
"position": 11,
"stddev": 2.149008526,
"variance": 4.618237645
},
{
"average": 121.3214254,
"count": 10653336,
"maximum": 255,
"minimum": 35,
"name": "Red",
"position": 12,
"stddev": 45.56261892,
"variance": 2075.952243
},
{
"average": 126.2526972,
"count": 10653336,
"maximum": 255,
"minimum": 49,
"name": "Green",
"position": 13,
"stddev": 36.85449467,
"variance": 1358.253777
},
{
"average": 111.2207554,
"count": 10653336,
"maximum": 255,
"minimum": 49,
"name": "Blue",
"position": 14,
"stddev": 31.95559955,
"variance": 1021.160343
},
{
"average": 247608.4011,
"count": 10653336,
"maximum": 249783.703,
"minimum": 245369.8966,
"name": "GpsTime",
"position": 15,
"stddev": 1176.138454,
"variance": 1383301.664
},
{
"average": 0,
"count": 10653336,
"maximum": 0,
"minimum": 0,
"name": "Synthetic",
"position": 16,
"stddev": 0,
"variance": 0
},
{
"average": 0,
"count": 10653336,
"maximum": 0,
"minimum": 0,
"name": "KeyPoint",
"position": 17,
"stddev": 0,
"variance": 0
},
{
"average": 0,
"count": 10653336,
"maximum": 0,
"minimum": 0,
"name": "Withheld",
"position": 18,
"stddev": 0,
"variance": 0
},
{
"average": 0,
"count": 10653336,
"maximum": 0,
"minimum": 0,
"name": "Overlap",
"position": 19,
"stddev": 0,
"variance": 0
}
]
}
}데이터셋 정보를 확인하면 해당 데이터의 좌표계, bbox좌표, 포인트 수, PDAL 버전, 파일 사이즈 등 기본적인 메타데이터를 확인할 수 있다. 이제 본격적으로 TileDB를 활용해보자.
3️⃣ LAS 데이터 TileDB에 로드하기
먼저 tiledb 경로와 데이터 무결성을 위해 초기화를 진행하고 시작하자.
# 경로 설정해주기
output_array = '~/autzen_tiledb'
# 중복 생성 방지와 무결성을 위해 초기화 해주고 진행
try:
shutil.rmtree(output_array)
except:
pass이제 las 데이터를 PDAL Pipeline을 활용해서 TileDB에 로드해보자.
pipeline = [
{
'type': 'readers.las',
'filename': 'autzen.laz'
},
{
'type': 'filters.stats'
},
{
'type': 'writers.tiledb',
'array_name': f"{output_array}",
'chunk_size': 1000000
}
]
with open("pipeline.json", "w") as f:
json.dump(pipeline, f)
print("pipeline 생성 완료")파이프라인이 생성되었으면 이제 pdal에서 pipeline을 실행시켜보자
# pdal의 pipeline을 실행
!pdal pipeline -i pipeline.json --nostream파이프라인이 실행되었다면 우리가 방금 생성한 TileDB의 배열 구조를 살펴보도록 하겠다. 먼저 이 배열은 3개의 차원(X,Y,Z)과 13개의 속성(intensity, color(RGB)으로 구성되었다. 특히 TileDB의 경우 3차원 데이터셋에서 매우 빠른 데이터 슬라이싱 능력을 제공하기 때문에 3차원 좌표 데이터를 다루기에 적합하다.
또한 데이터셋을 포함하는 3차원 boundary box가 존재하기 때문에 해당 데이터가 속한 구역의 최솟값과 최댓값 범위를 파악하기 용이하다.
%%time
with tiledb.open(output_array) as arr:
print(f"Non-empty domain: {arr.nonempty_domain()}")
print(arr.schema)위 코드를 실행해서 생성한 배열의 스키마 구조를 출력하고 13개의 속성에 어떤 값이 존재하는지 확인해보자.
더보기
Non-empty domain: ((635577.79, 639003.73), (848882.15, 853537.66), (406.14, 615.26))
ArraySchema(
domain=Domain(*[
Dim(name='X', domain=(635576.79, 639004.73), tile=0.0, dtype='float64', filters=FilterList([FloatScaleFilter(factor=0.01,offset=0.0,bytewidth=4), DeltaFilter(reinterp_dtype=int32), BitShuffleFilter(), ZstdFilter(level=7), ])),
Dim(name='Y', domain=(848881.15, 853538.66), tile=0.0, dtype='float64', filters=FilterList([FloatScaleFilter(factor=0.01,offset=0.0,bytewidth=4), DeltaFilter(reinterp_dtype=int32), BitShuffleFilter(), ZstdFilter(level=7), ])),
Dim(name='Z', domain=(405.14, 616.26), tile=0.0, dtype='float64', filters=FilterList([FloatScaleFilter(factor=0.01,offset=0.0,bytewidth=4), DeltaFilter(reinterp_dtype=int32), BitShuffleFilter(), ZstdFilter(level=7), ])),
]),
attrs=[
Attr(name='Intensity', dtype='uint16', var=False, nullable=False, enum_label=None, filters=FilterList([DeltaFilter(reinterp_dtype=None), ZstdFilter(level=5), ])),
Attr(name='Classification', dtype='uint8', var=False, nullable=False, enum_label=None, filters=FilterList([ZstdFilter(level=5), ])),
Attr(name='Synthetic', dtype='uint8', var=False, nullable=False, enum_label=None),
Attr(name='KeyPoint', dtype='uint8', var=False, nullable=False, enum_label=None),
Attr(name='Withheld', dtype='uint8', var=False, nullable=False, enum_label=None),
Attr(name='Overlap', dtype='uint8', var=False, nullable=False, enum_label=None),
Attr(name='ScanAngleRank', dtype='float32', var=False, nullable=False, enum_label=None, filters=FilterList([ZstdFilter(level=5), ])),
Attr(name='UserData', dtype='uint8', var=False, nullable=False, enum_label=None, filters=FilterList([ZstdFilter(level=5), ])),
Attr(name='PointSourceId', dtype='uint16', var=False, nullable=False, enum_label=None, filters=FilterList([ZstdFilter(level=5), ])),
Attr(name='GpsTime', dtype='float64', var=False, nullable=False, enum_label=None, filters=FilterList([DeltaFilter(reinterp_dtype=int64), BitWidthReductionFilter(window=256), ZstdFilter(level=7), ])),
Attr(name='Red', dtype='uint16', var=False, nullable=False, enum_label=None, filters=FilterList([DeltaFilter(reinterp_dtype=None), BitWidthReductionFilter(window=256), ZstdFilter(level=7), ])),
Attr(name='Green', dtype='uint16', var=False, nullable=False, enum_label=None, filters=FilterList([DeltaFilter(reinterp_dtype=None), BitWidthReductionFilter(window=256), ZstdFilter(level=7), ])),
Attr(name='Blue', dtype='uint16', var=False, nullable=False, enum_label=None, filters=FilterList([DeltaFilter(reinterp_dtype=None), BitWidthReductionFilter(window=256), ZstdFilter(level=7), ])),
Attr(name='BitFields', dtype='uint16', var=False, nullable=False, enum_label=None, filters=FilterList([ZstdFilter(level=5), ])),
],
cell_order='hilbert',
tile_order=None,
capacity=100000,
sparse=True,
allows_duplicates=True,
)
CPU times: user 8.38 ms, sys: 12.1 ms, total: 20.5 ms
Wall time: 29.2 ms- 데이터의 실제 범위 : ((635577.79, 639003.73), (848882.15, 853537.66), (406.14, 615.26))
- 13개 속성값 :
| 순번 | 속성 이름 (Attribute) | 데이터 타입 | 설명 |
| 1 | Intensity | uint16 | 반사 강도: 레이저가 물체에 맞고 돌아온 세기. 지표면의 재질(아스팔트, 풀숲 등)에 따라 값이 달라집니다. |
| 2 | Classification | uint8 | 분류: 점의 종류를 숫자로 구분 (예: 2=지면, 5=높은 식생, 6=건물 등). |
| 3 | Synthetic | uint8 | 가상 포인트: 실제 스캔 데이터가 아니라 사후 처리에 의해 생성된 점인지 여부. |
| 4 | KeyPoint | uint8 | 중요 지점: 모델링 등에서 제거하면 안 되는 중요한 구조적 포인트인지 여부. |
| 5 | Withheld | uint8 | 제외 대상: 노이즈나 오류로 판단되어 분석에서 제외해야 할 점인지 여부. |
| 6 | Overlap | uint8 | 중복 지점: 여러 번의 스캔 경로가 겹치는 구간에서 발생한 중복 데이터 여부. |
| 7 | ScanAngleRank | float32 | 스캔 각도: 레이저가 발사될 때 나간 각도 (-90 ~ +90도). |
| 8 | UserData | uint8 | 사용자 정의 데이터: 분석가가 특정 목적으로 할당해둔 임의의 숫자 값. |
| 9 | PointSourceId | uint16 | 소스 ID: 이 점을 쏜 비행 회차나 레이저 장비의 고유 번호. |
| 10 | GpsTime | float64 | 측정 시간: 각 포인트가 위성으로부터 기록된 정확한 시간 정보. |
| 11 | Red | uint16 | 빨강(R): 포인트에 입혀진 RGB 색상 중 Red 값. |
| 12 | Green | uint16 | 초록(G): 포인트에 입혀진 RGB 색상 중 Green 값. |
| 13 | Blue | uint16 | 파랑(B): 포인트에 입혀진 RGB 색상 중 Blue 값. |
| + | BitFields | uint16 | 비트 필드: 리턴 번호(Return Number) 등 여러 플래그 정보를 비트 단위로 합쳐놓은 값. |
4️⃣ TileDB에서 LiDAR 데이터 슬라이싱
다음으로는 TileDB의 배열 구조에서 데이터를 슬라이싱 해보겠다. 앞서 언급했듯이 TileDB는 3차원 배열 구조에서 슬라이싱에 매우 효율적이므로 슬라이싱에 대해 확인해보자. 결과는 pandas dataframe 형태로 저장된다.
%%time
with tiledb.open(output_array) as arr:
df = arr.df[636800:637800, 851000:853000, 406.14:615.26]
df
간단하게 코드를 해석해보면,
- `arr.df[...]` : TileDB 데이터를 Pandas DataFrame으로 바로 변환해서 가져온다.
- `[636800:637800, 851000:853000, 406.14:615.26]` :
- 첫 번째 범위 : X좌표(동서 방향으로 1,000m 구간)
- 두 번째 범위 : Y좌표(남북 방향으로 2,000m 구간)
- 세 번째 범위 : Z좌표(고도 406.14m ~ 615.26m 구간)
- 위 세 가지 범위가 겹쳐지는 3차원 직육면체 공간에 있는 점들만 슬라이싱하였다.
TileDB는 이러한 힐베르트 곡선 인덱싱 기법 덕분에, 해당 영역의 데이터가 하드디스크 어디에 위치해 있는지 알고 있기 때문에 해당 영역만 데이터를 가져올 수 있다. 그렇다면 어떤 작업에 특화되었는지도 알아보면
- `df['Classification'].value_range()`를 통해서 이 영역에 건물이 많은지 나무가 많은지 통계를 낼 수 있다.
- `df[['Red', 'Green', 'Blue']]` 값을 이용해서 이 영역만 따로 3D 시각화를 할 수 있다.
5️⃣ 시각화
당연히 포인트 시각화가 가능하다. 먼저 matplotlib을 사용해 3차원 포인트를 시각화 해보겠다.
matplotlib 시각화
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
# 3차원 포인트 시각화 함수 생성
def plot(pts):
rgbs = np.array(list(zip(pts['Red'], pts['Green'], pts['Blue']))) / 255.0
fig = plt.figure(figsize=(10, 10))
ax = fig.add_subplot(111, projection='3d')
ax.ticklabel_format(useOffset=False)
ax.scatter(pts['X'], pts['Y'], pts['Z'], c=rgbs)
ax.set_xlabel('x', fontsize=20, labelpad=20)
ax.set_ylabel('y', fontsize=20, labelpad=35)
ax.set_zlabel('z', fontsize=20, labelpad=25)
ax.set_title('Autzen', fontsize=20, pad=20)
ax.view_init(60, 96)
ax.tick_params(axis='y', pad=20)
ax.tick_params(axis='z', pad=10)
plt.show()- `rgps = np.array(list(zip(...))) / 255.0` : LiDAR 데이터의 RGB값은 보통 0~255 범위인데, Matplotlib의 scatter 함수는 색상 값을 0.0~1.0 사이로 받는다. 따라서 LiDAR 데이터의 RGB값을 제대로 표현하기 위해 255.0으로 나누어 정규화 작업을 수행한다.
- `ax.view_init(60, 96)` : 카메라의 시점(각도)설정. 상하 60도, 좌우 96도 방향에서 내려다보는 구도
BabylonJS 시각화
인터렉티브한 시각화를 원한다면 BabylonJS를 사용하면 된다. 단, 주의할 점은 conda-forge에는 pybabylonjs가 존재하지 않으므로 `pip install pybabylonjs`를 셀에서 실행해주었다.
from pybabylonjs import Show as show
# 데이터를 딕셔너리 형태로 변환
data = {
'X': df['X'][::10],
'Y': df['Y'][::10],
'Z': df['Z'][::10],
'Red': df['Red'][::10] / 255.0,
'Green': df['Green'][::10] / 255.0,
'Blue': df['Blue'][::10] / 255.0,
}
show.point_cloud(
source="dict",
data=data,
point_size=2,
width=800,
height=600
)위 코드를 통해 babylonjs를 실행할 수 있다. 조금 오래 걸리긴 하지만 잠시 기다리면 인터렉티브한 시각화가 가능하다.
'GIS > 01. GIS TIL' 카테고리의 다른 글
| 파이썬에서 Rasterio 기반 래스터 데이터 다루기 (0) | 2026.02.03 |
|---|---|
| 파이썬에서 GeoPandas 기반 벡터 데이터 다루기 (0) | 2026.02.02 |
| GIS 분석을 위한 환경설정(with conda, uv) (0) | 2026.02.01 |
| PostGIS 공간쿼리 기초 연습 (2) (0) | 2026.01.30 |
| PostGIS 공간쿼리 기초 연습 (1) (0) | 2026.01.29 |