위성사진 + Elevation + TIN -> RGB PointCloud

작성한 다른 글에서 사용하려고 테스트 데이터를 만드는 과정을 해결하여 정리해 봅니다.

사용하는 데이터셋은 비슷합니다만 GDAL만으로 처리하기에는 실력 부족으로 막혔던 부분을 데이터베이스에 raster데이터를 적재하고 각 포인트의 값을 추출하는 방식으로 해결하였습니다.

QGIS GeoReferencing을 통한 정사영상 좌표 부여

PointCloud의 자료 형식은 많이 있지만 제일 단순하다고 판단된 [ Text : X Y Z R G B ] 형식으로 원천 데이터를 생성하였고 CloudCompare등을 통해 LAS, PLY등으로 변환하여 사용 할 수 있었습니다.

 처음 사용한 자료는 국토정보플랫폼에서 제공하는 공개(DEM) 자료를 사용하였습니다 하지만 이자료는 구역별 높이값을 제공하고 있어서 생성된 결과물에 계단 처럼 층이 저 있었습니다.

정사영상과 계단처럼 표현된 Z 값 PointCloud

 그래서 표고데이터 및 등고선 데이터를 이용하여 PostGIS에서 제공하는 함수를 통해 raster 각 픽셀에 대한 Band 데이터에 접근할 수 있었고, 정사영상 픽셀위치에서 DEM으로 생성한 TIN에서 Z 값을 보간하는 방식을 사용하였습니다.

 PostGIS Intersection상으로는 2D, 3D 결과에 오류가 있어 세 점을 이용한 평면 방정식을 사용하여 Z값을 취하였습니다.

보간된 Z값을 갖는 PointCloud

 수집된 데이터에 대한 특성 파악이 중요하다는것을 알게된 작업이었던것 같습니다. 작업에 사용한 스크립트 입니다. 결과물은 Potree + PotreeConvertor를 이용하여 전처리하였고 CesiumJS를 이용하였습니다.

1 정사 영상 Load : raster2pgsql (EPSG:5186) GeoReferrencing
 
2 국토정보플랫폼 표고, 등고선 Elevation(EPSG:5181) -> TIN : st_delaunaytriangles
   (SHP to PgSQL)
 
create table tmp_rater_el_tin  as
 select (st_dump(st_delaunaytriangles(st_collect(st_setsrid(st_makepoint(st_x(geom), st_y(geom), val), st_srid(geom)))))).geom as geom
     from (
        select (st_pixelaspoints(rast, 1)).*
         from el
         ) a
         ;
 
3.정사영상 raster pixcel geometry + RGB 추출 -> Table
create table tmp_rater_rgb_pt as
select row_number () over() id, geom, r,g,b
 from
 (
  select geom    
             ,st_value(rast, 1, geom ) r
             ,st_value(rast, 2, geom ) g
             ,st_value(rast, 3, geom ) b
           from (
             select rast, (st_pixelaspoints(rast, 1)).*
             from "_test_raster" tr
             ) pt
     ) src
 
4 Mesh + Point Intersect를 통한 Z값 보정 RGB point 자료 생성

* 필요 시 Reprojection(표고 : 5181, 정사영상:5186 )
create table tmp_el_tin as
select
       id
      ,ST_Transform(geom, 5186) geom
from tmp_el_pt
;
            
drop table tmp_raster             ;

--fine Elevetion Z in rgb raster
create table tmp_raster as

 select
         st_x(rgb.geom) x
        ,st_y(rgb.geom) y
/*        ,st_z(st_intersection(el.geom, rgb.geom)) z */
        ,st_z(fnfindzonmeshaspoint2d(el.geom, rgb.geom)) z
        ,rgb.r
        ,rgb.g
        ,rgb.b
        --,st_intersection(rgb.geom, el.geom) ins */
  from  tmp_rater_rgb_pt  rgb
      , tmp_el_tin el
  where st_intersects(el.geom,rgb.geom) = true
 
5 Table Export (Text Blank delimiter)
  x y z r g b
-해더 및 포멧, 인코딩 오류 시 CloudCompare에서 las 변환 처리)

6. PointCloud Viewer : Potree + PotreeConvertor

세점 평면방정식을 이용한 Z 값 보간에 대해서는 보완하여 다음 글에서 소개해 드리겠습니다.