diff --git a/opendm/boundary.py b/opendm/boundary.py index bff5a5b9..a969873f 100644 --- a/opendm/boundary.py +++ b/opendm/boundary.py @@ -1,22 +1,67 @@ import fiona +import fiona.crs import os +import io +import json +from opendm import system +from pyproj import CRS +from opendm.location import transformer +from opendm.utils import double_quote -def load_boundary(boundary_json, reproject_to_crs=None): - if not os.path.isfile(boundary_json): - raise IOError("Cannot load boundary file: %s does not exist" % boundary_json) - - with fiona.open(boundary_json, 'r') as src: +def load_boundary(boundary_json, reproject_to_proj4=None): + with fiona.open(io.BytesIO(json.dumps(boundary_json).encode('utf-8')), 'r') as src: if len(src) != 1: - raise IOError("Boundary file must have a single polygon (found: %s)" % len(src)) + raise IOError("Boundary must have a single polygon (found: %s)" % len(src)) geom = src[0]['geometry'] if geom['type'] != 'Polygon': - raise IOError("Boundary file must have a polygon feature (found: %s)" % geom['type']) + raise IOError("Boundary must have a polygon feature (found: %s)" % geom['type']) - coords = geom['coordinates'] + rings = geom['coordinates'] - if reproject_to_crs is not None: - pass + if len(rings) == 0: + raise IOError("Boundary geometry has no rings") - return coords \ No newline at end of file + coords = rings[0] + if len(coords) == 0: + raise IOError("Boundary geometry has no coordinates") + + dimensions = len(coords[0]) + + if reproject_to_proj4 is not None: + t = transformer(CRS.from_proj4(fiona.crs.to_string(src.crs)), + CRS.from_proj4(reproject_to_proj4)) + coords = [t.TransformPoint(*c)[:dimensions] for c in coords] + + return coords + +def as_polygon(boundary): + return "POLYGON((" + ", ".join([" ".join(map(str, c)) for c in boundary]) + "))" + +def export_to_bounds_files(boundary, proj4, bounds_json_file, bounds_gpkg_file): + with open(bounds_json_file, "w") as f: + f.write(json.dumps({ + "type": "FeatureCollection", + "name": "bounds", + "features": [{ + "type": "Feature", + "properties": {}, + "geometry": { + "type": "Polygon", + "coordinates": [boundary] + } + }] + })) + + if os.path.isfile(bounds_gpkg_file): + os.remove(bounds_gpkg_file) + + kwargs = { + 'proj4': proj4, + 'input': double_quote(bounds_json_file), + 'output': double_quote(bounds_gpkg_file) + } + + system.run('ogr2ogr -overwrite -f GPKG -a_srs "{proj4}" {output} {input}'.format(**kwargs)) + diff --git a/opendm/cropper.py b/opendm/cropper.py index 0728d007..a6ba7532 100644 --- a/opendm/cropper.py +++ b/opendm/cropper.py @@ -5,6 +5,7 @@ from opendm.point_cloud import export_summary_json from osgeo import ogr import json, os from opendm.concurrency import get_max_memory +from opendm.utils import double_quote class Cropper: def __init__(self, storage_dir, files_prefix = "crop"): @@ -41,9 +42,9 @@ class Cropper: try: kwargs = { - 'gpkg_path': gpkg_path, - 'geotiffInput': original_geotiff, - 'geotiffOutput': geotiff_path, + 'gpkg_path': double_quote(gpkg_path), + 'geotiffInput': double_quote(original_geotiff), + 'geotiffOutput': double_quote(geotiff_path), 'options': ' '.join(map(lambda k: '-co {}={}'.format(k, gdal_options[k]), gdal_options)), 'warpOptions': ' '.join(warp_options), 'max_memory': get_max_memory() @@ -252,10 +253,13 @@ class Cropper: bounds_gpkg_path = os.path.join(self.storage_dir, '{}.bounds.gpkg'.format(self.files_prefix)) + if os.path.isfile(bounds_gpkg_path): + os.remove(bounds_gpkg_path) + # Convert bounds to GPKG kwargs = { - 'input': bounds_geojson_path, - 'output': bounds_gpkg_path, + 'input': double_quote(bounds_geojson_path), + 'output': double_quote(bounds_gpkg_path), 'proj4': pc_proj4 } diff --git a/opendm/orthophoto.py b/opendm/orthophoto.py index 02f00ade..72f99954 100644 --- a/opendm/orthophoto.py +++ b/opendm/orthophoto.py @@ -70,7 +70,7 @@ def generate_kmz(orthophoto_file, output_file=None, outsize=None): '--config GDAL_CACHEMAX %s%% ' % (orthophoto_file, output_file, bandparam, get_max_memory())) def post_orthophoto_steps(args, bounds_file_path, orthophoto_file, orthophoto_tiles_dir): - if args.crop > 0: + if args.crop > 0 or args.boundary: Cropper.crop(bounds_file_path, orthophoto_file, get_orthophoto_vars(args), keep_original=not args.optimize_disk_space, warp_options=['-dstalpha']) if args.build_overviews and not args.cog: diff --git a/opendm/osfm.py b/opendm/osfm.py index fc4752f7..df729b7d 100644 --- a/opendm/osfm.py +++ b/opendm/osfm.py @@ -510,10 +510,11 @@ def get_submodel_argv(args, submodels_path = None, submodel_name = None): tweaking --crop if necessary (DEM merging makes assumption about the area of DEMs and their euclidean maps that require cropping. If cropping is skipped, this leads to errors.) removing --gcp (the GCP path if specified is always "gcp_list.txt") reading the contents of --cameras + reading the contents of --boundary """ assure_always = ['orthophoto_cutline', 'dem_euclidean_map', 'skip_3dmodel', 'skip_report'] remove_always = ['split', 'split_overlap', 'rerun_from', 'rerun', 'gcp', 'end_with', 'sm_cluster', 'rerun_all', 'pc_csv', 'pc_las', 'pc_ept', 'tiles', 'copy-to', 'cog'] - read_json_always = ['cameras'] + read_json_always = ['cameras', 'boundary'] argv = sys.argv diff --git a/opendm/types.py b/opendm/types.py index 38fccb8d..7369ad1c 100644 --- a/opendm/types.py +++ b/opendm/types.py @@ -69,7 +69,7 @@ class ODM_Reconstruction(object): return self.georef is not None def has_gcp(self): - return self.is_georeferenced() and self.gcp is not None + return self.is_georeferenced() and self.gcp is not None and self.gcp.exists() def georeference_with_gcp(self, gcp_file, output_coords_file, output_gcp_file, output_model_txt_geo, rerun=False): if not io.file_exists(output_coords_file) or not io.file_exists(output_gcp_file) or rerun: diff --git a/stages/dataset.py b/stages/dataset.py index 6c76a9b2..f90c5257 100644 --- a/stages/dataset.py +++ b/stages/dataset.py @@ -9,7 +9,7 @@ from opendm import system from opendm.geo import GeoFile from shutil import copyfile from opendm import progress - +from opendm import boundary def save_images_database(photos, database_file): with open(database_file, 'w') as f: @@ -154,3 +154,10 @@ class ODMLoadDatasetStage(types.ODM_Stage): reconstruction.save_proj_srs(os.path.join(tree.odm_georeferencing, tree.odm_georeferencing_proj)) outputs['reconstruction'] = reconstruction + + # Try to load boundaries + if args.boundary: + if reconstruction.is_georeferenced(): + outputs['boundary'] = boundary.load_boundary(args.boundary, reconstruction.get_proj_srs()) + else: + log.ODM_WARNING("Reconstruction is not georeferenced, but boundary file provided (will ignore boundary file)") diff --git a/stages/odm_georeferencing.py b/stages/odm_georeferencing.py index e5739e42..03e4a84f 100644 --- a/stages/odm_georeferencing.py +++ b/stages/odm_georeferencing.py @@ -17,7 +17,7 @@ from opendm.cropper import Cropper from opendm import point_cloud from opendm.multispectral import get_primary_band_name from opendm.osfm import OSFMContext - +from opendm.boundary import as_polygon, export_to_bounds_files class ODMGeoreferencingStage(types.ODM_Stage): def process(self, args, outputs): @@ -128,6 +128,12 @@ class ODMGeoreferencingStage(types.ODM_Stage): '--writers.las.vlrs="{\\\"filename\\\": \\\"%s\\\", \\\"user_id\\\": \\\"ODM_GCP\\\", \\\"description\\\": \\\"Ground Control Points (GML)\\\"}"' % gcp_gml_export_file.replace(os.sep, "/") ] + if 'boundary' in outputs: + stages.append("crop") + params += [ + '--filters.crop.polygon="%s"' % as_polygon(outputs['boundary']) + ] + system.run(cmd + ' ' + ' '.join(stages) + ' ' + ' '.join(params)) self.update_progress(50) @@ -152,6 +158,14 @@ class ODMGeoreferencingStage(types.ODM_Stage): except: log.ODM_WARNING("Cannot calculate crop bounds! We will skip cropping") args.crop = 0 + + if 'boundary' in outputs and args.crop == 0: + log.ODM_INFO("Using boundary JSON as cropping area") + + bounds_base, _ = os.path.splitext(tree.odm_georeferencing_model_laz) + bounds_json = bounds_base + ".bounds.geojson" + bounds_gpkg = bounds_base + ".bounds.gpkg" + export_to_bounds_files(outputs['boundary'], reconstruction.get_proj_srs(), bounds_json, bounds_gpkg) else: log.ODM_INFO("Converting point cloud (non-georeferenced)") system.run(cmd + ' ' + ' '.join(stages) + ' ' + ' '.join(params))