GCP with OpenSfM, refactoring

opendm/gcp.py

@@ -1,12 +1,15 @@
 import glob
 import os
 from opendm import log
+from opendm import location
+from pyproj import CRS
 
 class GCPFile:
     def __init__(self, gcp_path):
         self.gcp_path = gcp_path
         self.entries = []
-        self.srs = ""
+        self.raw_srs = ""
+        self.srs = None
         self.read()
     
     def read(self):
@@ -16,7 +19,8 @@ class GCPFile:
     
             lines = map(str.strip, contents.split('\n'))
             if lines:
-                self.srs = lines[0] # SRS
+                self.raw_srs = lines[0] # SRS
+                self.srs = location.parse_srs_header(self.raw_srs)
 
                 for line in lines[1:]:
                     if line != "" and line[0] != "#":
@@ -47,6 +51,18 @@ class GCPFile:
     def exists(self):
         return 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:
+            point = list(self.entries_dict())[0]
+            longlat = CRS.from_epsg("4326")
+            lat, lon = location.transform(self.srs, longlat, point['x'], point['y'])
+            utm_zone, hemisphere = location.get_utm_zone_and_hemisphere_from(lon, lat)
+            return "WGS84 UTM %s%s" % (utm_zone, hemisphere)
+
     def make_filtered_copy(self, gcp_file_output, images_dir, min_images=3):
         """
         Creates a new GCP file from an existing GCP file includes
@@ -62,7 +78,7 @@ class GCPFile:
 
         files = map(os.path.basename, glob.glob(os.path.join(images_dir, "*")))
 
-        output = [self.srs]
+        output = [self.raw_srs]
         files_found = 0
         
         for entry in self.entries_dict():

opendm/location.py

@@ -1,6 +1,6 @@
 import math
 from opendm import log
-from pyproj import Proj
+from pyproj import Proj, Transformer, CRS
 
 def extract_utm_coords(photos, images_path, output_coords_file):
     """
@@ -54,6 +54,9 @@ 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):
+    transformer = Transformer.from_crs(from_srs, to_srs)
+    return transformer.transform(x, y)
 
 def get_utm_zone_and_hemisphere_from(lon, lat):
     """
@@ -84,3 +87,45 @@ def convert_to_utm(lon, lat, alt, utm_zone, hemisphere):
     x,y = p(lon, lat)
     return [x, y, alt]
 
+def parse_srs_header(header):
+    """
+    Parse a header coming from GCP or coordinate file
+    :param header (str) line
+    :return Proj object
+    """
+    log.ODM_DEBUG('Parsing SRS header: %s' % header)
+    header = header.strip()
+    ref = header.split(' ')
+    try:
+        if ref[0] == 'WGS84' and ref[1] == 'UTM':
+            datum = ref[0]
+            utm_pole = (ref[2][len(ref[2]) - 1]).upper()
+            utm_zone = int(ref[2][:len(ref[2]) - 1])
+            
+            proj_args = {
+                'zone': utm_zone, 
+                'datum': datum
+            }
+
+            proj4 = '+proj=utm +zone={zone} +datum={datum} +no_defs=True'
+            if utm_pole == 'S':
+                proj4 += ' +south=True'
+
+            srs = CRS.from_proj4(proj4.format(**proj_args))
+        elif '+proj' in header:
+            srs = CRS.from_proj4(header.strip('\''))
+        elif header.lower().startswith("epsg:"):
+            srs = CRS.from_epsg(header.lower()[5:])
+        else:
+            log.ODM_ERROR('Could not parse coordinates. Bad SRS supplied: %s' % header)
+    except RuntimeError as e:
+        log.ODM_ERROR('Uh oh! There seems to be a problem with your coordinates/GCP file.\n\n'
+                            'The line: %s\n\n'
+                            'Is not valid. Projections that are valid include:\n'
+                            ' - EPSG:*****\n'
+                            ' - WGS84 UTM **(N|S)\n'
+                            ' - Any valid proj4 string (for example, +proj=utm +zone=32 +north +ellps=WGS84 +datum=WGS84 +units=m +no_defs)\n\n'
+                            'Modify your input and try again.' % header)
+        raise RuntimeError(e)
+    
+    return srs

opendm/types.py

@@ -5,7 +5,9 @@ import os
 from fractions import Fraction
 from opensfm.exif import sensor_string
 from opendm import get_image_size
-from pyproj import Proj
+from opendm import location
+from opendm.gcp import GCPFile
+from pyproj import CRS
 
 import log
 import io
@@ -87,121 +89,103 @@ class ODM_Photo:
 
 
 class ODM_Reconstruction(object):
-    """docstring for ODMReconstruction"""
-
-    def __init__(self, photos, projstring = None, coords_file = None):
-        self.photos = photos    # list of ODM_Photos
-        self.projection = None  # Projection system the whole project will be in
+    def __init__(self, photos):
+        self.photos = photos
         self.georef = None
-        if projstring:
-            self.projection = self.set_projection(projstring)
-            self.georef = ODM_GeoRef(self.projection)
-        else:
-            self.projection = self.parse_coordinate_system(coords_file)
-            if self.projection:
-                self.georef = ODM_GeoRef(self.projection)
 
-    def parse_coordinate_system(self, _file):
-        """Write attributes to jobOptions from coord file"""
+    def is_georeferenced(self):
+        return self.georef is not None
+
+    def georeference_with_gcp(self, gcp_file, output_coords_file, reload_coords=False):
+        if not io.file_exists(output_coords_file) or reload_coords:
+            gcp = GCPFile(gcp_file)
+            if gcp.exists():
+                # Create coords file
+                with open(output_coords_file, 'w') as f:
+                    coords_header = gcp.wgs84_utm_zone()
+                    f.write(coords_header + "\n")
+                    log.ODM_DEBUG("Generated coords file from GCP: %s" % coords_header)
+            else:
+                log.ODM_WARNING("GCP file does not exist: %s" % gcp_file)
+                return
+        else:
+            log.ODM_INFO("Coordinates file already exist: %s" % output_coords_file)
+        
+        self.georef = ODM_GeoRef.FromCoordsFile(output_coords_file)
+        return self.georef
+
+    def georeference_with_gps(self, images_path, output_coords_file, reload_coords=False):
+        try:
+            if not io.file_exists(output_coords_file) or reload_coords:
+                location.extract_utm_coords(photos, tree.dataset_raw, output_coords_file)
+            else:
+                log.ODM_INFO("Coordinates file already exist: %s" % output_coords_file)
+            
+            self.georef = ODM_GeoRef.FromCoordsFile(output_coords_file)
+        except:
+            log.ODM_WARNING('Could not generate coordinates file. An orthophoto will not be generated.')
+
+        return self.georef
+
+    def save_proj_srs(self, file):
+        # Save proj to file for future use (unless this 
+        # dataset is not georeferenced)
+        if self.is_georeferenced():
+            with open(file, 'w') as f:
+                f.write(self.georef.proj4())
+
+class ODM_GeoRef(object):
+    @staticmethod
+    def FromProj(projstring):
+        return ODM_GeoRef(CRS.from_proj4(projstring))
+
+    @staticmethod
+    def FromCoordsFile(coords_file):
         # check for coordinate file existence
-        if not io.file_exists(_file):
-            log.ODM_WARNING('Could not find file %s' % _file)
+        if not io.file_exists(coords_file):
+            log.ODM_WARNING('Could not find file %s' % coords_file)
             return
 
-        with open(_file) as f:
+        srs = None
+
+        with open(coords_file) as f:
             # extract reference system and utm zone from first line.
             # We will assume the following format:
             # 'WGS84 UTM 17N' or 'WGS84 UTM 17N \n'
             line = f.readline().rstrip()
-            log.ODM_DEBUG('Line: %s' % line)
-            ref = line.split(' ')
-            # match_wgs_utm = re.search('WGS84 UTM (\d{1,2})(N|S)', line, re.I)
-            try:
-                if ref[0] == 'WGS84' and ref[1] == 'UTM':  # match_wgs_utm:
-                    datum = ref[0]
-                    utm_pole = (ref[2][len(ref[2]) - 1]).upper()
-                    utm_zone = int(ref[2][:len(ref[2]) - 1])
+            srs = location.parse_srs_header(line)
 
-                    proj_args = {
-                        'proj': "utm", 
-                        'zone': utm_zone, 
-                        'datum': datum,
-                        'no_defs': True
-                    }
-                    if utm_pole == 'S':
-                        proj_args['south'] = True
+        return ODM_GeoRef(srs)
 
-                    return Proj(**proj_args)
-                elif '+proj' in line:
-                    return Proj(line.strip('\''))
-                elif 'epsg' in line.lower():
-                    return Proj(init=line)
-                else:
-                    log.ODM_ERROR('Could not parse coordinates. Bad CRS supplied: %s' % line)
-            except RuntimeError as e:
-                log.ODM_ERROR('Uh oh! There seems to be a problem with your GCP file.\n\n'
-                                    'The line: %s\n\n'
-                                    'Is not valid. Projections that are valid include:\n'
-                                    ' - EPSG:*****\n'
-                                    ' - WGS84 UTM **(N|S)\n'
-                                    ' - Any valid proj4 string (for example, +proj=utm +zone=32 +north +ellps=WGS84 +datum=WGS84 +units=m +no_defs)\n\n'
-                                    'Modify your GCP file and try again.' % line)
-                raise RuntimeError(e)
-
-
-    def set_projection(self, projstring):
-        try:
-            return Proj(projstring)
-        except RuntimeError:
-            log.ODM_EXCEPTION('Could not set projection. Please use a proj4 string')
-
-
-class ODM_GeoRef(object):
-    """docstring for ODMUtmZone"""
-
-    def __init__(self, projection):
-        self.projection = projection
-        self.epsg = None
+    def __init__(self, srs):
+        self.srs = srs
         self.utm_east_offset = 0
         self.utm_north_offset = 0
         self.transform = []
-        self.gcps = []
 
-    def coord_to_fractions(self, coord, refs):
-        deg_dec = abs(float(coord))
-        deg = int(deg_dec)
-        minute_dec = (deg_dec - deg) * 60
-        minute = int(minute_dec)
+    def proj4(self):
+        return self.srs.to_proj4()
+    
+    def valid_utm_offsets(self):
+        return self.utm_east_offset and self.utm_north_offset
 
-        sec_dec = (minute_dec - minute) * 60
-        sec_dec = round(sec_dec, 3)
-        sec_denominator = 1000
-        sec_numerator = int(sec_dec * sec_denominator)
-        if float(coord) >= 0:
-            latRef = refs[0]
-        else:
-            latRef = refs[1]
-
-        output = str(deg) + '/1 ' + str(minute) + '/1 ' + str(sec_numerator) + '/' + str(sec_denominator)
-        return output, latRef
-
-    def extract_offsets(self, _file):
-        if not io.file_exists(_file):
-            log.ODM_ERROR('Could not find file %s' % _file)
+    def extract_offsets(self, geo_sys_file):
+        if not io.file_exists(geo_sys_file):
+            log.ODM_ERROR('Could not find file %s' % geo_sys_file)
             return
 
-        with open(_file) as f:
+        with open(geo_sys_file) as f:
             offsets = f.readlines()[1].split(' ')
             self.utm_east_offset = float(offsets[0])
             self.utm_north_offset = float(offsets[1])
 
-    def parse_transformation_matrix(self, _file):
-        if not io.file_exists(_file):
-            log.ODM_ERROR('Could not find file %s' % _file)
+    def parse_transformation_matrix(self, matrix_file):
+        if not io.file_exists(matrix_file):
+            log.ODM_ERROR('Could not find file %s' % matrix_file)
             return
 
         # Create a nested list for the transformation matrix
-        with open(_file) as f:
+        with open(matrix_file) as f:
             for line in f:
                 # Handle matrix formats that either
                 # have leading or trailing brakets or just plain numbers.

stages/dataset.py

@@ -6,7 +6,6 @@ from opendm import io
 from opendm import types
 from opendm import log
 from opendm import system
-from opendm import location
 from shutil import copyfile
 from opendm import progress
 
@@ -100,27 +99,13 @@ class ODMLoadDatasetStage(types.ODM_Stage):
 
         log.ODM_INFO('Found %s usable images' % len(photos))
 
-        # append photos to cell output
-        if not self.params.get('proj'):
-            if tree.odm_georeferencing_gcp:
-                outputs['reconstruction'] = types.ODM_Reconstruction(photos, coords_file=tree.odm_georeferencing_gcp)
-            else:
-                # Generate UTM from images
-                try:
-                    if not io.file_exists(tree.odm_georeferencing_coords) or self.rerun():
-                        location.extract_utm_coords(photos, tree.dataset_raw, tree.odm_georeferencing_coords)
-                    else:
-                        log.ODM_INFO("Coordinates file already exist: %s" % tree.odm_georeferencing_coords)
-                except:
-                    log.ODM_WARNING('Could not generate coordinates file. '
-                                    'Ignore if there is a GCP file')
+        # Create reconstruction object
+        reconstruction = types.ODM_Reconstruction(photos)
 
-                outputs['reconstruction'] = types.ODM_Reconstruction(photos, coords_file=tree.odm_georeferencing_coords)
+        if tree.odm_georeferencing_gcp:
+            reconstruction.georeference_with_gcp(tree.odm_georeferencing_gcp, tree.odm_georeferencing_coords)
         else:
-            outputs['reconstruction'] = types.ODM_Reconstruction(photos, projstring=self.params.get('proj'))
+            reconstruction.georeference_with_gps(tree.dataset_raw, tree.odm_georeferencing_coords, reload_coords=self.rerun())
 
-        # Save proj to file for future use (unless this 
-        # dataset is not georeferenced)
-        if outputs['reconstruction'].projection:
-            with open(io.join_paths(tree.odm_georeferencing, tree.odm_georeferencing_proj), 'w') as f:
-                f.write(outputs['reconstruction'].projection.srs)
+        reconstruction.save_proj_srs(io.join_paths(tree.odm_georeferencing, tree.odm_georeferencing_proj))
+        outputs['reconstruction'] = reconstruction

stages/odm_georeferencing.py

@@ -20,7 +20,6 @@ class ODMGeoreferencingStage(types.ODM_Stage):
         doPointCloudGeo = True
         transformPointCloud = True
         verbose = '-verbose' if self.params.get('verbose') else ''
-        geo_ref = reconstruction.georef
 
         runs = [{
             'georeferencing_dir': tree.odm_georeferencing,
@@ -73,8 +72,8 @@ class ODMGeoreferencingStage(types.ODM_Stage):
                 if transformPointCloud:
                     kwargs['pc_params'] = '-inputPointCloudFile {input_pc_file} -outputPointCloudFile {output_pc_file}'.format(**kwargs)
 
-                    if geo_ref and geo_ref.projection and geo_ref.projection.srs:
-                        kwargs['pc_params'] += ' -outputPointCloudSrs %s' % pipes.quote(geo_ref.projection.srs)
+                    if reconstruction.is_georeferenced():
+                        kwargs['pc_params'] += ' -outputPointCloudSrs %s' % pipes.quote(reconstruction.georef.proj4())
                     else:
                         log.ODM_WARNING('NO SRS: The output point cloud will not have a SRS.')
                 else:
@@ -99,9 +98,7 @@ class ODMGeoreferencingStage(types.ODM_Stage):
                     doPointCloudGeo = False # skip the rest of the georeferencing
 
                 if doPointCloudGeo:
-                    # update images metadata
-                    geo_ref.extract_offsets(odm_georeferencing_model_txt_geo_file)
-                    reconstruction.georef = geo_ref
+                    reconstruction.georef.extract_offsets(odm_georeferencing_model_txt_geo_file)
 
                     # XYZ point cloud output
                     if args.pc_csv:

stages/odm_orthophoto.py

@@ -38,19 +38,17 @@ class ODMOrthoPhotoStage(types.ODM_Stage):
 
             # Check if the georef object is initialized
             # (during a --rerun this might not be)
-            # TODO: we should move this to a more central
-            # location (perhaps during the dataset initialization)
-            if georef and not georef.utm_east_offset:
+            # TODO: this should be moved to a more central location?
+            if reconstruction.is_georeferenced() and not reconstruction.georef.valid_utm_offsets():
                 georeferencing_dir = tree.odm_georeferencing if args.use_3dmesh and not args.skip_3dmodel else tree.odm_25dgeoreferencing
                 odm_georeferencing_model_txt_geo_file = os.path.join(georeferencing_dir, tree.odm_georeferencing_model_txt_geo)
 
                 if io.file_exists(odm_georeferencing_model_txt_geo_file):
-                    georef.extract_offsets(odm_georeferencing_model_txt_geo_file)
+                    reconstruction.georef.extract_offsets(odm_georeferencing_model_txt_geo_file)
                 else:
                     log.ODM_WARNING('Cannot read UTM offset from {}. An orthophoto will not be generated.'.format(odm_georeferencing_model_txt_geo_file))
 
-
-            if georef:
+            if reconstruction.is_georeferenced():
                 if args.use_3dmesh:
                     kwargs['model_geo'] = os.path.join(tree.odm_texturing, tree.odm_georeferencing_model_obj_geo)
                 else:
@@ -69,7 +67,7 @@ class ODMOrthoPhotoStage(types.ODM_Stage):
             # Create georeferenced GeoTiff
             geotiffcreated = False
 
-            if georef and georef.projection and georef.utm_east_offset and georef.utm_north_offset:
+            if reconstruction.is_georeferenced() and reconstruction.georef.valid_utm_offsets():
                 ulx = uly = lrx = lry = 0.0
                 with open(tree.odm_orthophoto_corners) as f:
                     for lineNumber, line in enumerate(f):
@@ -77,13 +75,13 @@ class ODMOrthoPhotoStage(types.ODM_Stage):
                             tokens = line.split(' ')
                             if len(tokens) == 4:
                                 ulx = float(tokens[0]) + \
-                                    float(georef.utm_east_offset)
+                                    float(reconstruction.georef.utm_east_offset)
                                 lry = float(tokens[1]) + \
-                                    float(georef.utm_north_offset)
+                                    float(reconstruction.georef.utm_north_offset)
                                 lrx = float(tokens[2]) + \
-                                    float(georef.utm_east_offset)
+                                    float(reconstruction.georef.utm_east_offset)
                                 uly = float(tokens[3]) + \
-                                    float(georef.utm_north_offset)
+                                    float(reconstruction.georef.utm_north_offset)
                 log.ODM_INFO('Creating GeoTIFF')
 
                 orthophoto_vars = orthophoto.get_orthophoto_vars(args)
@@ -94,7 +92,7 @@ class ODMOrthoPhotoStage(types.ODM_Stage):
                     'lrx': lrx,
                     'lry': lry,
                     'vars': ' '.join(['-co %s=%s' % (k, orthophoto_vars[k]) for k in orthophoto_vars]),
-                    'proj': georef.projection.srs,
+                    'proj': reconstruction.georef.proj4(),
                     'png': tree.odm_orthophoto_file,
                     'tiff': tree.odm_orthophoto_tif,
                     'log': tree.odm_orthophoto_tif_log,

stages/run_opensfm.py

@@ -84,7 +84,7 @@ class ODMOpenSfMStage(types.ODM_Stage):
 
         self.update_progress(90)
 
-        if reconstruction.georef and (not io.file_exists(tree.opensfm_transformation) or self.rerun()):
-            octx.run('export_geocoords --transformation --proj \'%s\'' % reconstruction.georef.projection.srs)
+        if reconstruction.is_georeferenced() and (not io.file_exists(tree.opensfm_transformation) or self.rerun()):
+            octx.run('export_geocoords --transformation --proj \'%s\'' % reconstruction.georef.proj4())
         else:
             log.ODM_WARNING("Will skip exporting %s" % tree.opensfm_transformation)

test.sh

@@ -1 +1,5 @@
-python -m unittest discover tests "test_*.py"
+if [ ! -z "$1" ]; then
+	python -m unittest discover tests "test_$1.py"
+else
+	python -m unittest discover tests "test_*.py"
+fi

tests/assets/gcp_epsg_valid.txt

@@ -0,0 +1,3 @@
+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_utm_north_valid.txt

@@ -0,0 +1,4 @@
+WGS84 UTM 16N
+609925.8180680283 4950688.771811823 171.662982 1331 350 DJI_0002.JPG
+609925.8180680283 4950688.771811823 171.662982 2028 87 DJI_0003.JPG
+609925.8180680283 4950688.771811823 171.662982 2101 1181 DJI_0004.JPG

tests/test_gcp.py

@@ -0,0 +1,21 @@
+import time
+import unittest
+import os
+from opendm.gcp import GCPFile
+
+class TestRemote(unittest.TestCase):
+    def setUp(self):
+        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):
+        gcp = GCPFile("tests/assets/gcp_epsg_valid.txt")
+        self.assertTrue(gcp.exists())
+        self.assertEqual(gcp.wgs84_utm_zone(), "WGS84 UTM 16N")
+
+if __name__ == '__main__':
+    unittest.main()