Merge pull request #261 from OpenDroneMap/xyz-ptcloud

Xyz ptcloud

Former-commit-id: 70675e7b6c

opendm/types.py

@@ -123,7 +123,6 @@ class ODM_GeoRef(object):
             log.ODM_ERROR('Unknown pole format %s' % _pole)
             return
 
-
     def convert_to_las(self, _file, pdalXML):
 
         if not self.epsg:
@@ -139,7 +138,7 @@ class ODM_GeoRef(object):
                    'xml': pdalXML}
 
         # call txt2las
-        #system.run('{bin}/txt2las -i {f_in} -o {f_out} -skip 30 -parse xyzRGBssss ' \
+        # system.run('{bin}/txt2las -i {f_in} -o {f_out} -skip 30 -parse xyzRGBssss ' \
         #           '-set_scale 0.01 0.01 0.01 -set_offset {east} {north} 0 '  \
         #           '-translate_xyz 0 -epsg {epsg}'.format(**kwargs))
         #           
@@ -166,12 +165,11 @@ class ODM_GeoRef(object):
         pipelineXml += '  </Writer>'
         pipelineXml += '</Pipeline>'
 
-        
         with open(pdalXML, 'w') as f:
             f.write(pipelineXml)
 
         # call pdal 
-        system.run('{bin}/pdal pipeline -i {xml} --readers.ply.filename={f_in} ' \
+        system.run('{bin}/pdal pipeline -i {xml} --readers.ply.filename={f_in} '
                    '--writers.las.filename={f_out}'.format(**kwargs))
 
     def utm_to_latlon(self, _file, _photo, idx):
@@ -185,8 +183,8 @@ class ODM_GeoRef(object):
                    'y': gcp.y + self.utm_north_offset,
                    'z': gcp.z }
 
-        latlon = system.run_and_return('echo {x} {y} | cs2cs +proj=utm ' \
+        latlon = system.run_and_return('echo {x} {y} | cs2cs +proj=utm '
-            '+datum={datum} +ellps={datum} +zone={zone} +units=m +to '   \
+            '+datum={datum} +ellps={datum} +zone={zone} +units=m +to '
             '+proj=latlon +ellps={datum}'.format(**kwargs)).split()
 
         # Example: 83d18'16.285"W
@@ -292,6 +290,7 @@ class ODM_GeoRef(object):
         # update EPSG
         self.epsg = self.calculate_EPSG(self.utm_zone, self.utm_pole)
 
+
 class ODM_Tree(object):
     def __init__(self, root_path):
         ### root path to the project
@@ -336,14 +335,6 @@ class ODM_Tree(object):
             self.odm_texturing, 'odm_textured_model.obj')
         self.odm_textured_model_mtl = io.join_paths(
             self.odm_texturing, 'odm_textured_model.mtl')
-        self.odm_textured_model_txt_geo = io.join_paths(
-            self.odm_texturing, 'odm_textured_model_geo.txt')
-        self.odm_textured_model_ply_geo = io.join_paths(
-            self.odm_texturing, 'odm_textured_model_geo.ply')
-        self.odm_textured_model_obj_geo = io.join_paths(
-            self.odm_texturing, 'odm_textured_model_geo.obj')
-        self.odm_textured_model_mtl_geo = io.join_paths(
-            self.odm_texturing, 'odm_textured_model_geo.mtl')
         self.odm_texuring_log = io.join_paths(
             self.odm_texturing, 'odm_texturing_log.txt')
 
@@ -358,6 +349,16 @@ class ODM_Tree(object):
             self.odm_georeferencing, 'odm_georeferencing_utm_log.txt')
         self.odm_georeferencing_log = io.join_paths(
             self.odm_georeferencing, 'odm_georeferencing_log.txt')
+        self.odm_georeferencing_model_txt_geo = io.join_paths(
+            self.odm_georeferencing, 'odm_georeferencing_model_geo.txt')
+        self.odm_georeferencing_model_ply_geo = io.join_paths(
+            self.odm_georeferencing, 'odm_georeferenced_model.ply')
+        self.odm_georeferencing_model_obj_geo = io.join_paths(
+            self.odm_georeferencing, 'odm_textured_model_geo.obj')
+        self.odm_georeferencing_model_mtl_geo = io.join_paths(
+            self.odm_georeferencing, 'odm_textured_model_geo.mtl')
+        self.odm_georeferencing_xyz_file = io.join_paths(
+            self.odm_georeferencing, 'odm_georeferenced_model.csv')
         self.odm_georeferencing_pdal = io.join_paths(
             self.odm_georeferencing, 'pipeline.xml')
 
@@ -368,8 +369,3 @@ class ODM_Tree(object):
         self.odm_orthophoto_log = io.join_paths(self.odm_orthophoto, 'odm_orthophoto_log.txt')
         self.odm_orthophoto_tif_log = io.join_paths(self.odm_orthophoto, 'gdal_translate_log.txt')
 
-
-
-
-
-        

scripts/odm_georeferencing.py

@@ -1,4 +1,5 @@
 import ecto
+import csv
 
 from opendm import io
 from opendm import log
@@ -6,6 +7,7 @@ from opendm import types
 from opendm import system
 from opendm import context
 
+
 class ODMGeoreferencingCell(ecto.Cell):
     def declare_params(self, params):
         params.declare("gcp_file", 'path to the file containing the ground control '
@@ -38,8 +40,9 @@ class ODMGeoreferencingCell(ecto.Cell):
         if not self.params.use_gcp and \
            not io.file_exists(tree.odm_georeferencing_coords):
             
-            log.ODM_WARNING('Warning: No coordinates file. ' \
+            log.ODM_WARNING('Warning: No coordinates file. '
-                'Generating coordinates file in: %s' % tree.odm_georeferencing_coords)
+                            'Generating coordinates file in: %s'
+                            % tree.odm_georeferencing_coords)
             try:
                 # odm_georeference definitions
                 kwargs = {
@@ -51,13 +54,13 @@ class ODMGeoreferencingCell(ecto.Cell):
                 }
 
                 # run UTM extraction binary
-                system.run('{bin}/odm_extract_utm -imagesPath {imgs}/ '      \
+                system.run('{bin}/odm_extract_utm -imagesPath {imgs}/ '
-                    '-imageListFile {imgs_list} -outputCoordFile {coords} '  \
+                           '-imageListFile {imgs_list} -outputCoordFile {coords} '
                            '-logFile {log}'.format(**kwargs))
 
             except Exception, e:
-                log.ODM_ERROR('Could not generate GCP file from images metadata.'  \
+                log.ODM_ERROR('Could not generate GCP file from images metadata.'
-                    'Consider rerunning with argument --odm_georeferencing-useGcp' \
+                              'Consider rerunning with argument --odm_georeferencing-useGcp'
                               ' and provide a proper GCP file')
                 log.ODM_ERROR(e)
                 return ecto.QUIT
@@ -69,8 +72,8 @@ class ODMGeoreferencingCell(ecto.Cell):
                      (args['rerun_from'] is not None and
                       'odm_georeferencing' in args['rerun_from'])
 
-        if not io.file_exists(tree.odm_textured_model_obj_geo) or \
+        if not io.file_exists(tree.odm_georeferencing_model_obj_geo) or \
-           not io.file_exists(tree.odm_textured_model_ply_geo) or rerun_cell:
+           not io.file_exists(tree.odm_georeferencing_model_ply_geo) or rerun_cell:
 
             # odm_georeference definitions
             kwargs = {
@@ -82,9 +85,9 @@ class ODMGeoreferencingCell(ecto.Cell):
                 'pc': tree.pmvs_model,
                 'log': tree.odm_georeferencing_log,
                 'coords': tree.odm_georeferencing_coords,
-                'pc_geo': tree.odm_textured_model_ply_geo,
+                'pc_geo': tree.odm_georeferencing_model_ply_geo,
-                'geo_sys': tree.odm_textured_model_txt_geo,
+                'geo_sys': tree.odm_georeferencing_model_txt_geo,
-                'model_geo': tree.odm_textured_model_obj_geo,
+                'model_geo': tree.odm_georeferencing_model_obj_geo,
                 'size': self.params.img_size,
                 'gcp': io.join_paths(tree.root_path, self.params.gcp_file),
 
@@ -93,22 +96,17 @@ class ODMGeoreferencingCell(ecto.Cell):
             if self.params.use_gcp and \
                io.file_exists(tree.odm_georeferencing_coords):
 
-                system.run('{bin}/odm_georef -bundleFile {bundle} -inputCoordFile {coords} ' \
+                system.run('{bin}/odm_georef -bundleFile {bundle} -inputCoordFile {coords} '
-                    '-bundleResizedTo {size} -inputFile {model} -outputFile {model_geo} '    \
+                           '-bundleResizedTo {size} -inputFile {model} -outputFile {model_geo} '
-                    '-inputPointCloudFile {pc} -outputPointCloudFile {pc_geo} '              \
+                           '-inputPointCloudFile {pc} -outputPointCloudFile {pc_geo} '
-                    '-logFile {log} -georefFileOutputPath {geo_sys} -gcpFile {gcp} '         \
+                           '-logFile {log} -georefFileOutputPath {geo_sys} -gcpFile {gcp} '
                            '-outputCoordFile {coords}'.format(**kwargs))
             else:
-                system.run('{bin}/odm_georef -bundleFile {bundle} -inputCoordFile {coords} ' \
+                system.run('{bin}/odm_georef -bundleFile {bundle} -inputCoordFile {coords} '
-                    '-inputFile {model} -outputFile {model_geo} '                            \
+                           '-inputFile {model} -outputFile {model_geo} '
-                    '-inputPointCloudFile {pc} -outputPointCloudFile {pc_geo} '              \
+                           '-inputPointCloudFile {pc} -outputPointCloudFile {pc_geo} '
                            '-logFile {log} -georefFileOutputPath {geo_sys}'.format(**kwargs))
 
-        else:
-            log.ODM_WARNING('Found a valid georeferenced model in: %s' \
-                % tree.odm_textured_model_ply_geo)
-
-
             # update images metadata
             geo_ref = types.ODM_GeoRef()
             geo_ref.parse_coordinate_system(tree.odm_georeferencing_coords)
@@ -117,153 +115,29 @@ class ODMGeoreferencingCell(ecto.Cell):
                 geo_ref.utm_to_latlon(tree.odm_georeferencing_latlon, photo, idx)
 
             # convert ply model to LAS reference system
-        geo_ref.convert_to_las(tree.odm_textured_model_ply_geo, tree.odm_georeferencing_pdal)
+            geo_ref.convert_to_las(tree.odm_georeferencing_model_ply_geo,
+                                   tree.odm_georeferencing_pdal)
 
+            # XYZ point cloud output
+            log.ODM_INFO("Creating geo-referenced CSV file (XYZ format, can be used with GRASS to create DEM)")
+            with open(tree.odm_georeferencing_xyz_file, "wb") as csvfile:
+                csvfile_writer = csv.writer(csvfile, delimiter=",")
+                reachedpoints = False
+                with open(tree.odm_georeferencing_model_ply_geo) as f:
+                    for lineNumber, line in enumerate(f):
+                        if reachedpoints:
+                            tokens = line.split(" ")
+                            csv_line = [float(tokens[0])+geo_ref.utm_east_offset,
+                                        float(tokens[1])+geo_ref.utm_north_offset,
+                                        tokens[2]]
+                            csvfile_writer.writerow(csv_line)
+                        if line.startswith("end_header"):
+                            reachedpoints = True
+            csvfile.close()
+
+        else:
+            log.ODM_WARNING('Found a valid georeferenced model in: %s'
+                            % tree.odm_textured_model_ply_geo)
 
         log.ODM_INFO('Running OMD Georeferencing Cell - Finished')
         return ecto.OK if args['end_with'] != 'odm_georeferencing' else ecto.QUIT
-
-
-def odm_georeferencing():
-    """Run odm_georeferencing"""
-    print "\n  - running georeferencing - " + now()
-
-    os.chdir(jobOptions["jobDir"])
-    try:
-        os.mkdir(jobOptions["jobDir"] + "/odm_georeferencing")
-    except:
-        pass
-
-    if not args.odm_georeferencing_useGcp:
-        run("\"" + BIN_PATH + "/odm_extract_utm\" -imagesPath " + jobOptions["srcDir"] + "/ -imageListFile " \
-            + jobOptions["jobDir"] + "/pmvs/list.rd.txt -outputCoordFile " + jobOptions["jobDir"]            \
-            + "/odm_georeferencing/coordFile.txt")
-        
-        run("\"" + BIN_PATH + "/odm_georef\" -bundleFile " + jobOptions["jobDir"]                            \
-            + "/pmvs/bundle.rd.out -inputCoordFile " + jobOptions["jobDir"]                                  \
-            + "/odm_georeferencing/coordFile.txt -inputFile " + jobOptions["jobDir"]                         \
-            + "-results/odm_texturing/odm_textured_model.obj -outputFile " + jobOptions["jobDir"]            \
-            + "-results/odm_texturing/odm_textured_model_geo.obj -inputPointCloudFile "                      \
-            + jobOptions["jobDir"] + "-results/option-0000.ply -outputPointCloudFile " + jobOptions["jobDir"] \
-            + "-results/option-0000_georef.ply -logFile " + jobOptions["jobDir"]                             \
-            + "/odm_georeferencing/odm_georeferencing_log.txt -georefFileOutputPath " + jobOptions["jobDir"] \
-            + "-results/odm_texturing/odm_textured_model_geo_georef_system.txt")
-        
-    elif os.path.isfile(jobOptions["srcDir"] + "/" + args.odm_georeferencing_gcpFile):
-        run("\"" + BIN_PATH + "/odm_georef\" -bundleFile " + jobOptions["jobDir"]                             \
-            + "/pmvs/bundle.rd.out -gcpFile " + jobOptions["srcDir"] + "/" + args.odm_georeferencing_gcpFile  \
-            + " -imagesPath " + jobOptions["srcDir"] + "/ -imagesListPath " + jobOptions["jobDir"]            \
-            + "/pmvs/list.rd.txt -bundleResizedTo " + str(jobOptions["resizeTo"]) + " -inputFile "            \
-            + jobOptions["jobDir"] + "-results/odm_texturing/odm_textured_model.obj -outputFile "             \
-            + jobOptions["jobDir"] + "-results/odm_texturing/odm_textured_model_geo.obj -outputCoordFile "    \
-            + jobOptions["jobDir"] + "/odm_georeferencing/coordFile.txt -inputPointCloudFile "                \
-            + jobOptions["jobDir"] + "-results/option-0000.ply -outputPointCloudFile " + jobOptions["jobDir"] \
-            + "-results/option-0000_georef.ply -logFile " + jobOptions["jobDir"]                              \
-            + "/odm_georeferencing/odm_georeferencing_log.txt -georefFileOutputPath " + jobOptions["jobDir"]  \
-            + "-results/odm_texturing/odm_textured_model_geo_georef_system.txt")
-    else:
-        print "Warning: No GCP file. Consider rerunning with argument --odm_georeferencing-useGcp false --start-with odm_georeferencing"
-        print "Skipping orthophoto"
-        args.end_with = "odm_georeferencing"
-
-    if "csString" not in jobOptions:
-        parse_coordinate_system()
-
-    if "csString" in jobOptions and "utmEastOffset" in jobOptions and "utmNorthOffset" in jobOptions:
-        images = []
-        with open(jobOptions["jobDir"] + "/pmvs/list.rd.txt") as f:
-            images = f.readlines()
-
-        if len(images) > 0:
-            with open(jobOptions["jobDir"] + "/odm_georeferencing/coordFile.txt") as f:
-                for lineNumber, line in enumerate(f):
-                    if lineNumber >= 2 and lineNumber - 2 < len(images):
-                        tokens = line.split(' ')
-                        
-                        if len(tokens) >= 3:
-                            x = float(tokens[0])
-                            y = float(tokens[1])
-                            z = float(tokens[2])
-                            filename = images[lineNumber - 2]
-                            
-                            run("echo " + str(x + jobOptions["utmEastOffset"]) + " "                          \
-                                + str(y + jobOptions["utmNorthOffset"]) + " " + str(z)                        \
-                                + " | cs2cs " + jobOptions["csString"] + " +to +datum=WGS84 +proj=latlong > " \
-                                + jobOptions["jobDir"] + "/odm_georeferencing/latlong.txt")
-                            
-                            with open(jobOptions["jobDir"] + "/odm_georeferencing/latlong.txt") as latlongFile:
-                                latlongLine = latlongFile.readline()
-                                tokens = latlongLine.split()
-                                if len(tokens) >= 2:
-                                    exifGpsInfoWritten = False
-
-                                    lonString = tokens[0]  # Example: 83d18'16.285"W
-                                    latString = tokens[1]  # Example: 41d2'11.789"N
-                                    altString = ""
-                                    if len(tokens) > 2:
- 
-                                        altString = tokens[2] # Example: 0.998
-
-                                    tokens = re.split("[d '\"]+", lonString)
-                                    if len(tokens) >= 4:
-                                        lonDeg = tokens[0]
-                                        lonMin = tokens[1]
-                                        lonSec = tokens[2]
-                                        lonSecFrac = fractions.Fraction(lonSec)
-                                        lonSecNumerator = str(lonSecFrac._numerator)
-                                        lonSecDenominator = str(lonSecFrac._denominator)
-                                        lonRef = tokens[3]
-
-                                        tokens = re.split("[d '\"]+", latString)
-                                        if len(tokens) >= 4:
-                                            latDeg = tokens[0]
-                                            latMin = tokens[1]
-                                            latSec = tokens[2]
-                                            latSecFrac = fractions.Fraction(latSec)
-                                            latSecNumerator = str(latSecFrac._numerator)
-                                            latSecDenominator = str(latSecFrac._denominator)
-                                            latRef = tokens[3]
-
-                                            exivCmd = "exiv2 -q"
-                                            exivCmd += " -M\"set Exif.GPSInfo.GPSLatitude " + latDeg + "/1 "         \
-                                                + latMin + "/1 " + latSecNumerator + "/" + latSecDenominator + "\""
-                                            
-                                            exivCmd += " -M\"set Exif.GPSInfo.GPSLatitudeRef " + latRef + "\""
-                                            
-                                            exivCmd += " -M\"set Exif.GPSInfo.GPSLongitude " + lonDeg + "/1 "         \
-                                                + lonMin + "/1 " + lonSecNumerator + "/" + lonSecDenominator + "\""
-                                            
-                                            exivCmd += " -M\"set Exif.GPSInfo.GPSLongitudeRef " + lonRef + "\""
-
-                                            altNumerator = arcDenominator = 0  # BUG: arcDenominator is never used
-                                            
-                                            if altString:
-                                                altFrac = fractions.Fraction(altString)
-                                                altNumerator = str(altFrac._numerator)
-                                                altDenominator = str(altFrac._denominator)
-                                                exivCmd += " -M\"set Exif.GPSInfo.GPSAltitude " + altNumerator + "/" + altDenominator + "\""
-                                                exivCmd += " -M\"set Exif.GPSInfo.GPSAltitudeRef 0\""
-
-                                            exivCmd += " " + filename
-                                            run(exivCmd)
-                                            exifGpsInfoWritten = True
-
-                                    if not exifGpsInfoWritten:
-                                        print("    Warning: Failed setting EXIF GPS info for " \
-                                              + filename + " based on " + latlongLine)
-
-    if "epsg" in jobOptions and "utmEastOffset" in jobOptions and "utmNorthOffset" in jobOptions:
-        lasCmd = "\"" + BIN_PATH + "/txt2las\" -i " + jobOptions["jobDir"] +                                         \
-            "-results/option-0000_georef.ply -o " + jobOptions["jobDir"]                                             \
-            + "-results/pointcloud_georef.laz -skip 30 -parse xyzRGBssss -set_scale 0.01 0.01 0.01 -set_offset "     \
-            + str(jobOptions["utmEastOffset"]) + " " + str(jobOptions["utmNorthOffset"]) + " 0 -translate_xyz "      \
-            + str(jobOptions["utmEastOffset"]) + " " + str(jobOptions["utmNorthOffset"])                             \
-            + " 0 -epsg " + str(jobOptions["epsg"])
-        
-        print("    Creating geo-referenced LAS file (expecting warning)...")
-        
-        run(lasCmd)
- 
-
-    if args['--end-with'] != "odm_georeferencing":
-        odm_orthophoto()