EPSG code support, GCP files reprojection to UTM

opendm/config.py

@@ -95,17 +95,13 @@ def config():
                        choices=processopts,
                        help=('Can be one of:' + ' | '.join(processopts)))
 
-    parser.add_argument('--video',
+    # parser.add_argument('--video',
-                        metavar='<string>',
+    #                     metavar='<string>',
-                        help='Path to the video file to process')
+    #                     help='Path to the video file to process')
 
-    parser.add_argument('--slam-config',
+    # parser.add_argument('--slam-config',
-                        metavar='<string>',
+    #                     metavar='<string>',
-                        help='Path to config file for orb-slam')
+    #                     help='Path to config file for orb-slam')
-
-    parser.add_argument('--proj',
-                        metavar='<PROJ4 string>',
-                        help='Projection used to transform the model into geographic coordinates')
 
     parser.add_argument('--min-num-features',
                         metavar='<integer>',

opendm/gcp.py

@@ -6,6 +6,9 @@ from pyproj import CRS
 
 class GCPFile:
     def __init__(self, gcp_path):
+        if not gcp_path:
+            raise RuntimeError("gcp_path must be specified")
+
         self.gcp_path = gcp_path
         self.entries = []
         self.raw_srs = ""
@@ -30,39 +33,64 @@ class GCPFile:
                         else:
                             log.ODM_WARNING("Malformed GCP line: %s" % line)
 
-    def entries_dict(self):
+    def iter_entries(self, allowed_filenames=None):
         for entry in self.entries:
+            pe = self.parse_entry(entry)
+
+            if allowed_filenames is None or pe.filename in allowed_filenames:
+                yield pe
+    
+    def parse_entry(self, entry):
+        if entry:
             parts = entry.split()
             x, y, z, px, py, filename = parts[:6]
             extras = " ".join(parts[6:])
-            yield {
+            return GCPEntry(float(x), float(y), float(z), int(px), int(py), filename, extras)
-                'x': x,
-                'y': y,
-                'z': z,
-                'px': px,
-                'py': py,
-                'filename': filename,
-                'extras': extras
-            }
 
-    def entry_dict_to_s(self, entry):
+    def get_entry(self, n):
-        return "{x} {y} {z} {px} {py} {filename} {extras}".format(**entry).rstrip()
+        if n < self.entries_count():
+            return self.parse_entry(self.entries[n])
+
+    def entries_count(self):
+        return len(self.entries)
    
     def exists(self):
-        return self.gcp_path and os.path.exists(self.gcp_path)
+        return bool(self.gcp_path and os.path.exists(self.gcp_path))
 
     def wgs84_utm_zone(self):
         """
         Finds the UTM zone where the first point of the GCP falls into
         :return utm zone string valid for a coordinates header
         """
-        if len(self.entries) > 0:
+        if self.entries_count() > 0:
-            point = list(self.entries_dict())[0]
+            entry = self.get_entry(0)
             longlat = CRS.from_epsg("4326")
-            lat, lon = location.transform(self.srs, longlat, point['x'], point['y'])
+            lon, lat = location.transform2(self.srs, longlat, entry.x, entry.y)
             utm_zone, hemisphere = location.get_utm_zone_and_hemisphere_from(lon, lat)
             return "WGS84 UTM %s%s" % (utm_zone, hemisphere)
 
+    def create_utm_copy(self, gcp_file_output, filenames=None):
+        """
+        Creates a new GCP file from an existing GCP file
+        by optionally including only filenames and reprojecting each point to 
+        a UTM CRS
+        """
+        if os.path.exists(gcp_file_output):
+            os.remove(gcp_file_output)
+
+        output = [self.wgs84_utm_zone()]
+        target_srs = location.parse_srs_header(output[0])
+        transformer = location.transformer(self.srs, target_srs)
+
+        for entry in self.iter_entries(filenames):
+            entry.x, entry.y, entry.z = transformer.TransformPoint(entry.x, entry.y, entry.z)
+            output.append(str(entry))
+
+        with open(gcp_file_output, 'w') as f:
+            f.write('\n'.join(output) + '\n')
+
+        return gcp_file_output
+
     def make_filtered_copy(self, gcp_file_output, images_dir, min_images=3):
         """
         Creates a new GCP file from an existing GCP file includes
@@ -81,9 +109,9 @@ class GCPFile:
         output = [self.raw_srs]
         files_found = 0
         
-        for entry in self.entries_dict():
+        for entry in self.iter_entries():
-            if entry['filename'] in files:
+            if entry.filename in files:
-                output.append(self.entry_dict_to_s(entry))
+                output.append(str(entry))
                 files_found += 1
 
         if files_found >= min_images:
@@ -91,3 +119,19 @@ class GCPFile:
                 f.write('\n'.join(output) + '\n')
 
             return gcp_file_output
+
+class GCPEntry:
+    def __init__(self, x, y, z, px, py, filename, extras=""):
+        self.x = x
+        self.y = y
+        self.z = z
+        self.px = px
+        self.py = py
+        self.filename = filename
+        self.extras = extras
+    
+    def __str__(self):
+        return "{} {} {} {} {} {} {}".format(self.x, self.y, self.z, 
+                                             self.px, self.py, 
+                                             self.filename, 
+                                             self.extras).rstrip()

opendm/location.py

@@ -1,6 +1,7 @@
 import math
 from opendm import log
 from pyproj import Proj, Transformer, CRS
+from osgeo import osr
 
 def extract_utm_coords(photos, images_path, output_coords_file):
     """
@@ -54,9 +55,31 @@ def extract_utm_coords(photos, images_path, output_coords_file):
         for coord in coords:
             f.write("%s %s %s\n" % (coord[0] - dx, coord[1] - dy, coord[2]))
     
-def transform(from_srs, to_srs, x, y):
+def transform2(from_srs, to_srs, x, y):
-    transformer = Transformer.from_crs(from_srs, to_srs)
+    return transformer(from_srs, to_srs).TransformPoint(x, y, 0)[:2]
-    return transformer.transform(x, y)
+
+def transform3(from_srs, to_srs, x, y, z):
+    return transformer(from_srs, to_srs).TransformPoint(x, y, z)
+
+def proj_srs_convert(srs):
+    """
+    Convert a Proj SRS object to osr SRS object
+    """
+    res = osr.SpatialReference()
+    epsg = srs.to_epsg()
+
+    if epsg:
+        res.ImportFromEPSG(epsg)
+    else:
+        proj4 = srs.to_proj4()
+        res.ImportFromProj4(proj4)
+    
+    return res
+
+def transformer(from_srs, to_srs):
+    src = proj_srs_convert(from_srs)
+    tgt = proj_srs_convert(to_srs)
+    return osr.CoordinateTransformation(src, tgt)
     
 def get_utm_zone_and_hemisphere_from(lon, lat):
     """
@@ -107,7 +130,7 @@ def parse_srs_header(header):
                 'datum': datum
             }
 
-            proj4 = '+proj=utm +zone={zone} +datum={datum} +no_defs=True'
+            proj4 = '+proj=utm +zone={zone} +datum={datum} +units=m +no_defs=True'
             if utm_pole == 'S':
                 proj4 += ' +south=True'
 

opendm/osfm.py

@@ -58,7 +58,7 @@ class OSFMContext:
             raise Exception("Reconstruction could not be generated")
 
 
-    def setup(self, args, images_path, photos, gcp_path=None, append_config = [], rerun=False):
+    def setup(self, args, images_path, photos, gcp_path=None, append_config = [], rerun=False): #georeferenced=True, 
         """
         Setup a OpenSfM project
         """
@@ -117,6 +117,10 @@ class OSFMContext:
                 config.append("align_method: orientation_prior")
                 config.append("align_orientation_prior: vertical")
             
+            # if not georeferenced:
+            #     config.append("")
+
+
             if args.use_hybrid_bundle_adjustment:
                 log.ODM_DEBUG("Enabling hybrid bundle adjustment")
                 config.append("bundle_interval: 100")          # Bundle after adding 'bundle_interval' cameras

stages/dataset.py

@@ -103,7 +103,7 @@ class ODMLoadDatasetStage(types.ODM_Stage):
         reconstruction = types.ODM_Reconstruction(photos)
 
         if tree.odm_georeferencing_gcp:
-            reconstruction.georeference_with_gcp(tree.odm_georeferencing_gcp, tree.odm_georeferencing_coords)
+            reconstruction.georeference_with_gcp(tree.odm_georeferencing_gcp, tree.odm_georeferencing_coords, reload_coords=self.rerun())
         else:
             reconstruction.georeference_with_gps(tree.dataset_raw, tree.odm_georeferencing_coords, reload_coords=self.rerun())
 

tests/assets/.gitignore

@@ -0,0 +1 @@
+output/

tests/assets/gcp_epsg_valid.txt

@@ -1,3 +0,0 @@
-EPSG:4326
-44.701151842605441 -85.61322729531112 150.5 1331 350 DJI_0002.JPG
-44.701023803183617 -85.612848642726888 151.5 1331 350 DJI_0003.JPG

tests/assets/gcp_latlon_south.txt

@@ -0,0 +1,3 @@
+EPSG:4326
+105.03506 -5.23270 20 1331 350 DJI_0002.JPG
+

tests/assets/gcp_latlon_valid.txt

@@ -0,0 +1,3 @@
++proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs
+-85.61322729531112 44.701151842605441 150.5 1331 350 DJI_0002.JPG
+-85.612848642726888 44.701023803183617 151.5 1331 350 DJI_0003.JPG

tests/assets/gcp_michigan_feet_valid.txt

@@ -0,0 +1,2 @@
+EPSG:2251
+26607594.75468351 -26813.882366498266 563.1987598425 1331 350 DJI_0002.JPG

tests/test_gcp.py

@@ -1,21 +1,57 @@
 import time
 import unittest
 import os
+import shutil
+
 from opendm.gcp import GCPFile
 
-class TestRemote(unittest.TestCase):
+
+class TestGcp(unittest.TestCase):
     def setUp(self):
-        pass
+        if os.path.exists("tests/assets/output"):
+            shutil.rmtree("tests/assets/output")
+        os.makedirs("tests/assets/output")
+        # pass
 
     def test_utm_north(self):
         gcp = GCPFile("tests/assets/gcp_utm_north_valid.txt")
         self.assertTrue(gcp.exists())
         self.assertEqual(gcp.wgs84_utm_zone(), "WGS84 UTM 16N")
     
-    def test_epsg(self):
+    def test_latlon_south(self):
-        gcp = GCPFile("tests/assets/gcp_epsg_valid.txt")
+        gcp = GCPFile("tests/assets/gcp_latlon_south.txt")
+        self.assertTrue(gcp.exists())
+        self.assertEqual(gcp.wgs84_utm_zone(), "WGS84 UTM 48S")
+    
+    def test_latlon(self):
+        gcp = GCPFile("tests/assets/gcp_latlon_valid.txt")
         self.assertTrue(gcp.exists())
         self.assertEqual(gcp.wgs84_utm_zone(), "WGS84 UTM 16N")
 
+    def test_utm_conversion(self):
+        gcp = GCPFile("tests/assets/gcp_latlon_valid.txt")
+        copy = GCPFile(gcp.create_utm_copy("tests/assets/output/gcp_utm.txt"))
+        self.assertTrue(copy.exists())
+        self.assertEqual(copy.raw_srs, "WGS84 UTM 16N")
+        self.assertEqual(copy.get_entry(0).x, 609865.707705)
+        self.assertEqual(copy.get_entry(0).y, 4950688.36182)
+    
+    def test_utm_conversion_feet(self):
+        gcp = GCPFile("tests/assets/gcp_michigan_feet_valid.txt")
+        copy = GCPFile(gcp.create_utm_copy("tests/assets/output/gcp_utm_z.txt"))
+        self.assertTrue(copy.exists())
+        self.assertEqual(copy.raw_srs, "WGS84 UTM 16N")
+        self.assertEqual(round(copy.get_entry(0).x, 3), 609925.818)
+        self.assertEqual(round(copy.get_entry(0).y, 3), 4950688.772)
+        self.assertEqual(round(copy.get_entry(0).z, 3), 171.663)
+
+    def test_filtered_copy(self):
+        gcp = GCPFile('tests/assets/gcp_latlon_valid.txt')
+        self.assertTrue(gcp.exists())
+        self.assertEqual(gcp.entries_count(), 2)
+        copy = GCPFile(gcp.make_filtered_copy('tests/assets/output/filtered_copy.txt', 'tests/assets/images', min_images=1))
+        self.assertTrue(copy.exists())
+        self.assertEqual(copy.entries_count(), 1)
+
 if __name__ == '__main__':
     unittest.main()