From 7a8502f58536d360c1c1c01b795393e9eb8eb907 Mon Sep 17 00:00:00 2001 From: Piero Toffanin Date: Thu, 20 Jun 2019 15:39:49 -0400 Subject: [PATCH] GCP with OpenSfM, refactoring --- opendm/gcp.py | 22 +++- opendm/location.py | 47 +++++++- opendm/types.py | 166 ++++++++++++--------------- stages/dataset.py | 29 ++--- stages/odm_georeferencing.py | 9 +- stages/odm_orthophoto.py | 22 ++-- stages/run_opensfm.py | 4 +- test.sh | 6 +- tests/assets/gcp_epsg_valid.txt | 3 + tests/assets/gcp_utm_north_valid.txt | 4 + tests/test_gcp.py | 21 ++++ 11 files changed, 195 insertions(+), 138 deletions(-) create mode 100644 tests/assets/gcp_epsg_valid.txt create mode 100644 tests/assets/gcp_utm_north_valid.txt create mode 100644 tests/test_gcp.py diff --git a/opendm/gcp.py b/opendm/gcp.py index 9a034d55..cd311d17 100644 --- a/opendm/gcp.py +++ b/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(): diff --git a/opendm/location.py b/opendm/location.py index 69ae2b47..5a9d9ce7 100644 --- a/opendm/location.py +++ b/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 \ No newline at end of file diff --git a/opendm/types.py b/opendm/types.py index 8918c143..e9c2e5fb 100644 --- a/opendm/types.py +++ b/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. diff --git a/stages/dataset.py b/stages/dataset.py index daae9021..641ebb8a 100644 --- a/stages/dataset.py +++ b/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 diff --git a/stages/odm_georeferencing.py b/stages/odm_georeferencing.py index f846e161..d27650ef 100644 --- a/stages/odm_georeferencing.py +++ b/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: diff --git a/stages/odm_orthophoto.py b/stages/odm_orthophoto.py index 61f5c655..34bc9527 100644 --- a/stages/odm_orthophoto.py +++ b/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, diff --git a/stages/run_opensfm.py b/stages/run_opensfm.py index 8f3e12de..b5f8d69a 100644 --- a/stages/run_opensfm.py +++ b/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) \ No newline at end of file diff --git a/test.sh b/test.sh index 1ebe2f67..91fe5153 100755 --- a/test.sh +++ b/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 diff --git a/tests/assets/gcp_epsg_valid.txt b/tests/assets/gcp_epsg_valid.txt new file mode 100644 index 00000000..1836a362 --- /dev/null +++ b/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 diff --git a/tests/assets/gcp_utm_north_valid.txt b/tests/assets/gcp_utm_north_valid.txt new file mode 100644 index 00000000..567d3930 --- /dev/null +++ b/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 diff --git a/tests/test_gcp.py b/tests/test_gcp.py new file mode 100644 index 00000000..45d4d53b --- /dev/null +++ b/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() \ No newline at end of file