From 946b39ad7bba348853e38933bcdd3534df7cbf0b Mon Sep 17 00:00:00 2001
From: Stephen Mather <1174901+smathermather@users.noreply.github.com>
Date: Thu, 8 Aug 2024 22:14:23 -0400
Subject: [PATCH] Las resolution tied to GSD (#1788)

* When possible, ensure las resolution retains detail based on estimated density of data rather than just fixed value
---
 stages/odm_georeferencing.py | 30 +++++++++++++++++++++++++++---
 1 file changed, 27 insertions(+), 3 deletions(-)

diff --git a/stages/odm_georeferencing.py b/stages/odm_georeferencing.py
index a478ee1db..de16c374b 100644
--- a/stages/odm_georeferencing.py
+++ b/stages/odm_georeferencing.py
@@ -6,6 +6,7 @@
 import fiona.crs
 import json
 import zipfile
+import math
 from collections import OrderedDict
 from pyproj import CRS
 
@@ -125,13 +126,35 @@ def process(self, args, outputs):
 
                 stages.append("transformation")
                 utmoffset = reconstruction.georef.utm_offset()
+
+                # Establish appropriate las scale for export
+                las_scale = 0.001
+                filtered_point_cloud_stats = tree.path("odm_filterpoints", "point_cloud_stats.json")
+                # Function that rounds to the nearest 10
+                # and then chooses the one below so our
+                # las scale is sensible
+                def powerr(r):
+                    return pow(10,round(math.log10(r))) / 10
+
+                if os.path.isfile(filtered_point_cloud_stats):
+                    try:
+                        with open(filtered_point_cloud_stats, 'r') as stats:
+                             las_stats = json.load(stats)
+                             spacing = powerr(las_stats['spacing'])
+                             log.ODM_INFO("las scale calculated as the minimum of 1/10 estimated spacing or %s, which ever is less." % las_scale)
+                             las_scale = min(spacing, 0.001)
+                    except Exception as e:
+                        log.ODM_WARNING("Cannot find file point_cloud_stats.json. Using default las scale: %s" % las_scale)
+                else:
+                    log.ODM_INFO("No point_cloud_stats.json found. Using default las scale: %s" % las_scale)
+
                 params += [
                     f'--filters.transformation.matrix="1 0 0 {utmoffset[0]} 0 1 0 {utmoffset[1]} 0 0 1 0 0 0 0 1"',
                     f'--writers.las.offset_x={reconstruction.georef.utm_east_offset}' ,
                     f'--writers.las.offset_y={reconstruction.georef.utm_north_offset}',
-                    '--writers.las.scale_x=0.001',
-                    '--writers.las.scale_y=0.001',
-                    '--writers.las.scale_z=0.001',
+                    f'--writers.las.scale_x={las_scale}',
+                    f'--writers.las.scale_y={las_scale}',
+                    f'--writers.las.scale_z={las_scale}',
                     '--writers.las.offset_z=0',
                     f'--writers.las.a_srs="{reconstruction.georef.proj4()}"' # HOBU this should maybe be WKT
                 ]
@@ -257,3 +280,4 @@ def transform_textured_model(obj):
             os.remove(tree.filtered_point_cloud)
 
 
+