1. Add a topocentric to geocoords converter, allowing to convert every points directly.

2. update odm_georeferencing stage to convert and generate ply, laz, and texture model in UTM coordinate system. Some minor change to use pdal python api instead of command line tool

opendm/georef.py

@@ -2,32 +2,8 @@ import numpy as np
 from numpy import ndarray
 from typing import Tuple
 from pyproj import Proj
-import argparse
 from opensfm.geo import TopocentricConverter
 
-def parse_pdal_args(pdal_args: dict) -> argparse.Namespace:
-    def validate_arg(name: str, data_type: type):
-        if name not in pdal_args:
-            raise ValueError(f"PDAL arguments should contain {name}.")
-        try:
-             data_type(pdal_args[name])
-        except ValueError:
-            raise ValueError(f"PDAL argument {name} should be of type {data_type}.")
-        return data_type(pdal_args[name])
-    return argparse.Namespace(
-        reflat=validate_arg('reflat', float),
-        reflon=validate_arg('reflon', float),
-        refalt=validate_arg('refalt', float),
-        x_offset=validate_arg('x_offset', float),
-        y_offset=validate_arg('y_offset', float),
-        a_srs=validate_arg('a_srs', str)
-    )
-
-def topocentric_to_georef_pdal(ins, outs):
-    args = parse_pdal_args(ins)
-    outs['X'], outs['Y'], outs['Z'] = topocentric_to_georef(args.reflat, args.reflon, args.refalt, args.a_srs, ins['X'], ins['Y'], ins['Z'], args.x_offset, args.y_offset)
-    return True
-
 def topocentric_to_georef(
         reflat: float, 
         reflon: float, 
@@ -46,4 +22,35 @@ def topocentric_to_georef(
     return easting - x_offset, northing - y_offset, alt
 
 
+class TopocentricToProj:
+    def __init__(self, reflat:float, reflon:float, refalt:float, a_srs:str):
+        self.reference = TopocentricConverter(reflat, reflon, refalt)
+        self.projection = Proj(a_srs)
+        
+    def convert_array(self, arr:ndarray, offset_x:float=0, offset_y:float=0):
+        x, y, z = arr['X'], arr['Y'], arr['Z']
+        easting, northing, alt = self.convert_points(x, y, z, offset_x, offset_y)
+        arr['X'] = easting
+        arr['Y'] = northing
+        arr['Z'] = alt
+        return arr
+    
+    def convert_points(self, x:ndarray, y:ndarray, z:ndarray, offset_x:float=0, offset_y:float=0):
+        lat, lon, alt = self.reference.to_lla(x, y, z)
+        easting, northing = self.projection(lon, lat)
+        return easting - offset_x, northing - offset_y, alt
+        
+    def convert_obj(self, input_obj:str, output_obj:str, offset_x:float=0, offset_y:float=0):
+        with open(input_obj, 'r') as fin:
+            with open(output_obj, 'w') as fout:
+                lines = fin.readlines()
+                for line in lines:
+                    if line.startswith("v "):
+                        v = np.fromstring(line.strip()[2:] + " 1",  sep=' ', dtype=float)
+                        vt = self.convert_points(v[0], v[1], v[2], offset_x, offset_y)
+                        fout.write("v " + " ".join(map(str, list(vt))) + '\n')
+                    else:
+                        fout.write(line)
+
+
     

opendm/osfm.py

@@ -631,7 +631,7 @@ class OSFMContext:
     
     def reference(self):
         ds = DataSet(self.opensfm_project_path)
-        return ds.load_reference_lla()
+        return ds.load_reference()
 
     def name(self):
         return os.path.basename(os.path.abspath(self.path("..")))    

opendm/types.py

@@ -374,7 +374,7 @@ class ODM_Tree(object):
         self.odm_25dmeshing_log = os.path.join(self.odm_meshing, 'odm_25dmeshing_log.txt')
 
         # texturing
-        self.odm_textured_model_obj_topo = os.path.join(self.odm_texturing, 'odm_textured_model_topocentric.obj')
+        self.odm_textured_model_obj_topo = 'odm_textured_model_topocentric.obj'
         self.odm_textured_model_obj = 'odm_textured_model_geo.obj'
         self.odm_textured_model_glb = 'odm_textured_model_geo.glb'
 

stages/mvstex.py

@@ -71,7 +71,6 @@ class ODMMvsTexStage(types.ODM_Stage):
 
             odm_textured_model_obj = os.path.join(r['out_dir'], tree.odm_textured_model_obj_topo)
             unaligned_obj = io.related_file_path(odm_textured_model_obj, postfix="_unaligned")
-
             if not io.file_exists(odm_textured_model_obj) or self.rerun():
                 log.ODM_INFO('Writing MVS Textured file in: %s'
                               % odm_textured_model_obj)
@@ -94,7 +93,7 @@ class ODMMvsTexStage(types.ODM_Stage):
                 # mvstex definitions
                 kwargs = {
                     'bin': context.mvstex_path,
-                    'out_dir': os.path.join(r['out_dir'], "odm_textured_model_geo"),
+                    'out_dir': os.path.join(r['out_dir'], tree.odm_textured_model_obj_topo.split('.obj')[0]),
                     'model': r['model'],
                     'dataTerm': 'gmi',
                     'outlierRemovalType': 'gauss_clamping',

stages/odm_georeferencing.py

@@ -8,7 +8,6 @@ import json
 import zipfile
 import math
 from collections import OrderedDict
-from numpy import rec
 from pyproj import CRS
 import pdal
 
@@ -25,6 +24,7 @@ from opendm.osfm import OSFMContext
 from opendm.boundary import as_polygon, export_to_bounds_files
 from opendm.align import compute_alignment_matrix, transform_point_cloud, transform_obj
 from opendm.utils import np_to_json
+from opendm.georef import TopocentricToProj
 
 class ODMGeoreferencingStage(types.ODM_Stage):
     def process(self, args, outputs):
@@ -120,79 +120,50 @@ class ODMGeoreferencingStage(types.ODM_Stage):
             # prepare pipeline stage for topocentric to georeferenced conversion
             octx = OSFMContext(tree.opensfm)
             reference = octx.reference()
-            pdalargs = {
+            converter = TopocentricToProj(reference.lat, reference.lon, reference.alt, reconstruction.georef.proj4())
-                'reflat': reference.lat,
-                'reflon': reference.lon,
-                'refalt': reference.alt,
-                'a_srs': reconstruction.georef.proj4(),
-                'x_offset': reconstruction.georef.utm_east_offset,
-                'y_offset': reconstruction.georef.utm_north_offset
-            }
-            topo_to_georef_def = {
-                "type": "filters.python",
-                "module": "opendm.georef",
-                "function": "topocentric_to_georef_pdal",
-                "pdalargs": json.dumps(pdalargs)
-            }
         
         if not io.file_exists(tree.filtered_point_cloud) or self.rerun():
             log.ODM_INFO("Georeferecing filtered point cloud")
             if reconstruction.is_georeferenced():
-                ply_georeferencing_pipeline = [
+                pipeline = pdal.Reader.ply(tree.filtered_point_cloud_topo).pipeline()
-                    {
+                pipeline.execute()
-                        "type": "readers.ply",
+                arr = pipeline.arrays[0]
-                        "filename": tree.filtered_point_cloud_topo
+                arr = converter.convert_array(
-                    },
+                    arr,
-                    topo_to_georef_def,
+                    reconstruction.georef.utm_east_offset,
-                    {
+                    reconstruction.georef.utm_north_offset
-                        "type": "writers.ply",
+                )
-                        "filename": tree.filtered_point_cloud
+                pipeline = pdal.Writer.ply(tree.filtered_point_cloud).pipeline(arr)
-                    }
-                ]
-                pipeline = pdal.Pipeline(json.dumps(ply_georeferencing_pipeline))
                 pipeline.execute()
             else:
                 shutil.copy(tree.filtered_point_cloud_topo, tree.filtered_point_cloud)
 
-        if not io.file_exists(tree.odm_textured_model_geo_obj) or self.rerun():
+        if not io.file_exists(os.path.join(tree.odm_texturing, tree.odm_textured_model_obj)) or self.rerun():
             log.ODM_INFO("Georeferecing textured model")
+            obj_in = os.path.join(tree.odm_texturing, tree.odm_textured_model_obj_topo)
+            obj_out = os.path.join(tree.odm_texturing, tree.odm_textured_model_obj)
             if reconstruction.is_georeferenced():
-                obj_georeferencing_pipeline = [
+                converter.convert_obj(
-                    {
+                    obj_in, 
-                        "type": "readers.obj",
+                    obj_out, 
-                        "filename": tree.odm_textured_model_obj
+                    reconstruction.georef.utm_east_offset, 
-                    },
+                    reconstruction.georef.utm_north_offset
-                    topo_to_georef_def,
+                )
-                    {
-                        "type": "writers.obj",
-                        "filename": tree.odm_textured_model_geo_obj
-                    }
-                ]
-                pipeline = pdal.Pipeline(json.dumps(obj_georeferencing_pipeline))
-                pipeline.execute()
             else:
-                shutil.copy(tree.odm_textured_model_obj, tree.odm_textured_model_geo_obj)
+                shutil.copy(obj_in, obj_out)
         
         if not io.file_exists(tree.odm_georeferencing_model_laz) or self.rerun():
-            las_georeferencing_pipeline = [
+            pipeline = pdal.Pipeline()
-                {
+            pipeline |= pdal.Reader.ply(tree.filtered_point_cloud)
-                    "type": "readers.las",
+            pipeline |= pdal.Filter.ferry(dimensions="views => UserData")
-                    "filename": tree.filtered_point_cloud
-                },
-                {
-                    "type": "filters.ferry",
-                    "dimensions": "views => UserData"
-                }
-            ]
             
             if reconstruction.is_georeferenced():
                 log.ODM_INFO("Georeferencing point cloud")
 
                 utmoffset = reconstruction.georef.utm_offset()
-                las_georeferencing_pipeline.append({
+                pipeline |= pdal.Filter.transformation(
-                    "type": "filters.transformation",
+                    matrix=f"1 0 0 {utmoffset[0]} 0 1 0 {utmoffset[1]} 0 0 1 0 0 0 0 1"
-                    "matrix": f"1 0 0 {utmoffset[0]} 0 1 0 {utmoffset[1]} 0 0 1 0 0 0 0 1"
+                )
-                })
 
                 # Establish appropriate las scale for export
                 las_scale = 0.001
@@ -216,7 +187,6 @@ class ODMGeoreferencingStage(types.ODM_Stage):
                     log.ODM_INFO("No point_cloud_stats.json found. Using default las scale: %s" % las_scale)
 
                 las_writer_def = {
-                    "type": "writers.las",
                     "filename": tree.odm_georeferencing_model_laz,
                     "a_srs": reconstruction.georef.proj4(),
                     "offset_x": utmoffset[0],
@@ -230,9 +200,6 @@ class ODMGeoreferencingStage(types.ODM_Stage):
                 if reconstruction.has_gcp() and io.file_exists(gcp_geojson_zip_export_file):
                     if os.path.getsize(gcp_geojson_zip_export_file) <= 65535:
                         log.ODM_INFO("Embedding GCP info in point cloud")
-                        params += [
-                            '--writers.las.vlrs="{\\\"filename\\\": \\\"%s\\\", \\\"user_id\\\": \\\"ODM\\\", \\\"record_id\\\": 2, \\\"description\\\": \\\"Ground Control Points (zip)\\\"}"' % gcp_geojson_zip_export_file.replace(os.sep, "/")
-                        ]
                         las_writer_def["vlrs"] = json.dumps(
                             {
                                 "filename": gcp_geojson_zip_export_file.replace(os.sep, "/"),
@@ -245,9 +212,10 @@ class ODMGeoreferencingStage(types.ODM_Stage):
                     else:
                         log.ODM_WARNING("Cannot embed GCP info in point cloud, %s is too large" % gcp_geojson_zip_export_file)
 
-                las_georeferencing_pipeline.append(las_writer_def)
+                pipeline |= pdal.Writer.las(
+                    **las_writer_def
+                )
     
-                pipeline = pdal.Pipeline(json.dumps(las_georeferencing_pipeline))
                 pipeline.execute()
 
                 self.update_progress(50)
@@ -282,11 +250,9 @@ class ODMGeoreferencingStage(types.ODM_Stage):
                     export_to_bounds_files(outputs['boundary'], reconstruction.get_proj_srs(), bounds_json, bounds_gpkg)
             else:
                 log.ODM_INFO("Converting point cloud (non-georeferenced)")
-                las_georeferencing_pipeline.append({
+                pipeline |= pdal.Writer.las(
-                    "type": "writers.las",
+                    tree.odm_georeferencing_model_laz
-                    "filename": tree.odm_georeferencing_model_laz
+                )
-                })
-                pipeline = pdal.Pipeline(json.dumps(las_georeferencing_pipeline))
                 pipeline.execute()