From 73bf8ca2a65731d7ad7980e0ee45e2dd2d1c42e9 Mon Sep 17 00:00:00 2001 From: Piero Toffanin Date: Mon, 29 Apr 2019 14:01:55 -0400 Subject: [PATCH] Merge DEM step finished --- opendm/dem/merge.py | 32 ++++++++++++++++---------------- opendm/dem/utils.py | 9 +++++++++ scripts/odm_dem.py | 10 ++-------- scripts/splitmerge.py | 30 +++++++++++++++++++++++++++--- 4 files changed, 54 insertions(+), 27 deletions(-) create mode 100644 opendm/dem/utils.py diff --git a/opendm/dem/merge.py b/opendm/dem/merge.py index fa3ba656..44741930 100644 --- a/opendm/dem/merge.py +++ b/opendm/dem/merge.py @@ -9,7 +9,7 @@ import os from rasterio import logging -def euclidean_merge_dems(input_dems, output_dem): +def euclidean_merge_dems(input_dems, output_dem, creation_options={}): """ Based on https://github.com/mapbox/rio-merge-rgba and ideas from Anna Petrasova @@ -22,13 +22,7 @@ def euclidean_merge_dems(input_dems, output_dem): inputs = [] bounds=None precision=7 - creation_options={} - log_ = logging.getLogger() - debug_enabled = log_.isEnabledFor(logging.DEBUG) - if debug_enabled: - logging.disable(logging.DEBUG) - existing_dems = [] for dem in input_dems: if not io.file_exists(dem): @@ -40,6 +34,12 @@ def euclidean_merge_dems(input_dems, output_dem): log.ODM_WARNING("No input DEMs, skipping euclidean merge.") return + # Suppress DEBUG logging for rasterio operations + log_ = logging.getLogger() + debug_enabled = log_.isEnabledFor(logging.DEBUG) + if debug_enabled: + logging.disable(logging.DEBUG) + with rasterio.open(existing_dems[0]) as first: src_nodata = 0 nodatavals = first.get_nodatavals() @@ -110,10 +110,10 @@ def euclidean_merge_dems(input_dems, output_dem): profile["transform"] = output_transform profile["height"] = output_height profile["width"] = output_width - profile["tiled"] = True - profile["blockxsize"] = 512 - profile["blockysize"] = 512 - profile["compress"] = 'lzw' + profile["tiled"] = creation_options.get('TILED', 'YES') == 'YES' + profile["blockxsize"] = creation_options.get('BLOCKXSIZE', 512) + profile["blockysize"] = creation_options.get('BLOCKYSIZE', 512) + profile["compress"] = creation_options.get('COMPRESS', 'LZW') profile["nodata"] = src_nodata # Creation opts @@ -133,14 +133,10 @@ def euclidean_merge_dems(input_dems, output_dem): dst_count = first.count dst_shape = (dst_count, dst_rows, dst_cols) log.ODM_DEBUG("Temp shape: %r", dst_shape) - # dstarr = np.full(dst_shape, src_nodata, dtype=dtype) dstarr = np.zeros(dst_shape, dtype=dtype) distsum = np.zeros(dst_shape, dtype=dtype) - # Read up srcs until - # a. everything is data; i.e. no nodata - # b. no sources left for src_d, src_e in sources: # The full_cover behavior is problematic here as it includes # extra pixels along the bottom right when the sources are @@ -187,5 +183,9 @@ def euclidean_merge_dems(input_dems, output_dem): for _, euclidean_geotiff in inputs: if io.file_exists(euclidean_geotiff): os.remove(euclidean_geotiff) - + + # Restore logging + if debug_enabled: + logging.disable(logging.NOTSET) + return output_dem \ No newline at end of file diff --git a/opendm/dem/utils.py b/opendm/dem/utils.py new file mode 100644 index 00000000..45118f42 --- /dev/null +++ b/opendm/dem/utils.py @@ -0,0 +1,9 @@ + +def get_dem_vars(args): + return { + 'TILED': 'YES', + 'COMPRESS': 'LZW', + 'BLOCKXSIZE': 512, + 'BLOCKYSIZE': 512, + 'NUM_THREADS': args.max_concurrency, + } diff --git a/scripts/odm_dem.py b/scripts/odm_dem.py index a311ac21..41d94470 100644 --- a/scripts/odm_dem.py +++ b/scripts/odm_dem.py @@ -7,7 +7,7 @@ from opendm import system from opendm import context from opendm import types from opendm import gsd -from opendm.dem import commands +from opendm.dem import commands, utils from opendm.cropper import Cropper class ODMDEMStage(types.ODM_Stage): @@ -80,13 +80,7 @@ class ODMDEMStage(types.ODM_Stage): if args.crop > 0: bounds_file_path = os.path.join(tree.odm_georeferencing, 'odm_georeferenced_model.bounds.gpkg') if os.path.exists(bounds_file_path): - Cropper.crop(bounds_file_path, os.path.join(odm_dem_root, "{}.tif".format(product)), { - 'TILED': 'YES', - 'COMPRESS': 'LZW', - 'BLOCKXSIZE': 512, - 'BLOCKYSIZE': 512, - 'NUM_THREADS': self.params.get('max_concurrency') - }) + Cropper.crop(bounds_file_path, os.path.join(odm_dem_root, "{}.tif".format(product)), utils.get_dem_vars(args)) else: log.ODM_WARNING('Found existing outputs in: %s' % odm_dem_root) else: diff --git a/scripts/splitmerge.py b/scripts/splitmerge.py index 0c56131f..68c8725b 100644 --- a/scripts/splitmerge.py +++ b/scripts/splitmerge.py @@ -6,7 +6,8 @@ from opendm import types from opendm import io from opendm import system from opendm import orthophoto -from opendm.dem import pdal +from opendm.dem import pdal, utils +from opendm.dem.merge import euclidean_merge_dems from opensfm.large import metadataset from opendm.cropper import Cropper from opendm.concurrency import get_max_memory, get_max_memory_mb @@ -237,10 +238,33 @@ class ODMMergeStage(types.ODM_Stage): else: log.ODM_WARNING("Found merged orthophoto in %s" % tree.odm_orthophoto_tif) - # Merge DEM + # Merge DEMs - # TODO: crop DEM if necessary + def merge_dems(dem_filename, human_name): + dem_file = os.path.join(tree.path("odm_dem", dem_filename)) + if not io.file_exists(dem_file) or self.rerun(): + all_dems = get_submodel_paths(tree.submodels_path, "odm_dem", dem_filename) + log.ODM_INFO("Merging %ss" % human_name) + + # Merge + dem_vars = utils.get_dem_vars(args) + euclidean_merge_dems(all_dems, dem_file, dem_vars) + if io.file_exists(dem_file): + # Crop + if args.crop > 0: + Cropper.crop(merged_bounds_file, dem_file, dem_vars) + log.ODM_INFO("Created %s" % dem_file) + else: + log.ODM_WARNING("Cannot merge %s, %s was not created" % (human_name, dem_file)) + else: + log.ODM_WARNING("Found merged %s in %s" % (human_name, dsm_file)) + + if args.dsm: + merge_dems("dsm.tif", "DSM") + + if args.dtm: + merge_dems("dtm.tif", "DTM") # Stop the pipeline short! We're done. self.next_stage = None