From 321c92691cca2540f14b58cf660700be80106ca0 Mon Sep 17 00:00:00 2001
From: Piero Toffanin <pt@masseranolabs.com>
Date: Wed, 24 Mar 2021 17:22:44 +0000
Subject: [PATCH] Added pc2dem module

---
 contrib/pc2dem/README.md | 44 +++++++++++++++++++++++++++++++
 contrib/pc2dem/pc2dem.py | 57 ++++++++++++++++++++++++++++++++++++++++
 2 files changed, 101 insertions(+)
 create mode 100644 contrib/pc2dem/README.md
 create mode 100755 contrib/pc2dem/pc2dem.py

diff --git a/contrib/pc2dem/README.md b/contrib/pc2dem/README.md
new file mode 100644
index 00000000..2bcc8115
--- /dev/null
+++ b/contrib/pc2dem/README.md
@@ -0,0 +1,44 @@
+# Point Cloud To DEM
+
+Convert point clouds (LAS, LAZ, PLY, and any other format compatible with [PDAL](https://pdal.io/stages/readers.html) to GeoTIFF elevation models using ODM's open source algorithms.
+
+![image](https://user-images.githubusercontent.com/1951843/112354653-492a5100-8ca3-11eb-9f21-4dda4cae976f.png)
+
+This tool includes method to perform efficient and scalable gapfill interpolation and is the same method used by ODM's processing pipeline. It is offered here as a standalone module for processing individual point clouds.
+
+## Usage
+
+```
+docker run -ti --rm -v /home/youruser/folder_with_point_cloud:/input --entrypoint /code/contrib/pc2dem/pc2dem.py opendronemap/odm /input/point_cloud.las [flags]
+```
+
+The result (`dsm.tif` or `dtm.tif`) will be stored in the same folder as the input point cloud. See additional `flags` you can pass at the end of the command above:
+
+```
+usage: pc2dem.py [-h] [--type {dsm,dtm}] [--resolution RESOLUTION]
+                 [--gapfill-steps GAPFILL_STEPS]
+                 point_cloud
+
+Generate DEMs from point clouds using ODM's algorithm.
+
+positional arguments:
+  point_cloud           Path to point cloud file (.las, .laz,
+                        .ply)
+
+optional arguments:
+  -h, --help            show this help message and exit
+  --type {dsm,dtm}      Type of DEM. Default: dsm
+  --resolution RESOLUTION
+                        Resolution in m/px. Default: 0.05
+  --gapfill-steps GAPFILL_STEPS
+                        Number of steps used to fill areas with
+                        gaps. Set to 0 to disable gap filling.
+                        Starting with a radius equal to the output
+                        resolution, N different DEMs are generated
+                        with progressively bigger radius using the
+                        inverse distance weighted (IDW) algorithm
+                        and merged together. Remaining gaps are
+                        then merged using nearest neighbor
+                        interpolation. Default: 3
+
+```
diff --git a/contrib/pc2dem/pc2dem.py b/contrib/pc2dem/pc2dem.py
new file mode 100755
index 00000000..0441565b
--- /dev/null
+++ b/contrib/pc2dem/pc2dem.py
@@ -0,0 +1,57 @@
+#!/usr/bin/env python3
+# Author: Piero Toffanin
+# License: AGPLv3
+
+import os
+import sys
+sys.path.insert(0, os.path.join("..", "..", os.path.dirname(__file__)))
+
+import argparse
+import multiprocessing
+from opendm.dem import commands
+
+parser = argparse.ArgumentParser(description='Generate DEMs from point clouds using ODM\'s algorithm.')
+parser.add_argument('point_cloud',
+                type=str,
+                help='Path to point cloud file (.las, .laz, .ply)')
+parser.add_argument('--type',
+                    type=str,
+                    choices=("dsm", "dtm"),
+                    default="dsm",
+                    help="Type of DEM. Default: %(default)s")
+parser.add_argument('--resolution',
+                type=float,
+                default=0.05,
+                help='Resolution in m/px. Default: %(default)s')
+parser.add_argument('--gapfill-steps',
+                    default=3,
+                    type=int,
+                    help='Number of steps used to fill areas with gaps. Set to 0 to disable gap filling. '
+                            'Starting with a radius equal to the output resolution, N different DEMs are generated with '
+                            'progressively bigger radius using the inverse distance weighted (IDW) algorithm '
+                            'and merged together. Remaining gaps are then merged using nearest neighbor interpolation. '
+                            'Default: %(default)s')
+args = parser.parse_args()
+
+if not os.path.exists(args.point_cloud):
+    print("%s does not exist" % args.point_cloud)
+    exit(1)
+
+outdir = os.path.dirname(args.point_cloud)
+
+radius_steps = [args.resolution / 2.0]
+for _ in range(args.gapfill_steps - 1):
+    radius_steps.append(radius_steps[-1] * 2) # 2 is arbitrary, maybe there's a better value?
+
+commands.create_dem(args.point_cloud,
+                    args.type,
+                    output_type='idw' if args.type == 'dtm' else 'max',
+                    radiuses=list(map(str, radius_steps)),
+                    gapfill=args.gapfill_steps > 0,
+                    outdir=outdir,
+                    resolution=args.resolution,
+                    decimation=1,
+                    verbose=True,
+                    max_workers=multiprocessing.cpu_count(),
+                    keep_unfilled_copy=False
+                )
\ No newline at end of file