From c679d400c820c4cfd38d61168667284f158f9fc0 Mon Sep 17 00:00:00 2001 From: Merten Fermont Date: Thu, 12 Oct 2023 22:06:53 +0200 Subject: [PATCH 1/3] Tiler zoom level is calculated from GSD Instead of hardcoding a value, calculate the maximum zoomlevel in which there is still an increase in detail. By using the configured orthophoto resolution or GSD. The higher the latitude, the higher the resolution will be of the tile. Resulting in a chance of generating useless tiles, as there is no compensation for this. At the moment it'll use the worst-case resolution from the equator. Zoom level calulation from: https://wiki.openstreetmap.org/wiki/Zoom_levels --- opendm/orthophoto.py | 4 ++-- opendm/tiles/tiler.py | 21 +++++++++++++++------ stages/odm_dem.py | 2 +- stages/odm_orthophoto.py | 12 ++++++------ stages/splitmerge.py | 4 ++-- 5 files changed, 26 insertions(+), 17 deletions(-) diff --git a/opendm/orthophoto.py b/opendm/orthophoto.py index 79865c718..319a504cb 100644 --- a/opendm/orthophoto.py +++ b/opendm/orthophoto.py @@ -85,7 +85,7 @@ def generate_kmz(orthophoto_file, output_file=None, outsize=None): system.run('gdal_translate -of KMLSUPEROVERLAY -co FORMAT=PNG "%s" "%s" %s ' '--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): +def post_orthophoto_steps(args, bounds_file_path, orthophoto_file, orthophoto_tiles_dir, resolution): 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']) @@ -99,7 +99,7 @@ def post_orthophoto_steps(args, bounds_file_path, orthophoto_file, orthophoto_ti generate_kmz(orthophoto_file) if args.tiles: - generate_orthophoto_tiles(orthophoto_file, orthophoto_tiles_dir, args.max_concurrency) + generate_orthophoto_tiles(orthophoto_file, orthophoto_tiles_dir, args.max_concurrency, resolution) if args.cog: convert_to_cogeo(orthophoto_file, max_workers=args.max_concurrency, compression=args.orthophoto_compression) diff --git a/opendm/tiles/tiler.py b/opendm/tiles/tiler.py index 6b36f2098..f091e7146 100644 --- a/opendm/tiles/tiler.py +++ b/opendm/tiles/tiler.py @@ -1,16 +1,25 @@ import os import sys +import math from opendm import log from opendm import system from opendm import io -def generate_tiles(geotiff, output_dir, max_concurrency): +def generate_tiles(geotiff, output_dir, max_concurrency, resolution): + circumference_earth_cm = 2*math.pi*637_813_700 + px_per_tile = 256 + resolution_equator_cm = circumference_earth_cm/px_per_tile + zoom = math.ceil(math.log(resolution_equator_cm/resolution, 2)) + + min_zoom = 5 # 4.89 km/px + max_zoom = min(zoom, 30) # No deeper zoom than 30 (0.0146 cm/px) + gdal2tiles = os.path.join(os.path.dirname(__file__), "gdal2tiles.py") - system.run('%s "%s" --processes %s -z 5-21 -n -w none "%s" "%s"' % (sys.executable, gdal2tiles, max_concurrency, geotiff, output_dir)) + system.run('%s "%s" --processes %s -z %s-%s -n -w none "%s" "%s"' % (sys.executable, gdal2tiles, max_concurrency, min_zoom, max_zoom, geotiff, output_dir)) -def generate_orthophoto_tiles(geotiff, output_dir, max_concurrency): +def generate_orthophoto_tiles(geotiff, output_dir, max_concurrency, resolution): try: - generate_tiles(geotiff, output_dir, max_concurrency) + generate_tiles(geotiff, output_dir, max_concurrency, resolution) except Exception as e: log.ODM_WARNING("Cannot generate orthophoto tiles: %s" % str(e)) @@ -37,10 +46,10 @@ def generate_colored_hillshade(geotiff): log.ODM_WARNING("Cannot generate colored hillshade: %s" % str(e)) return (None, None, None) -def generate_dem_tiles(geotiff, output_dir, max_concurrency): +def generate_dem_tiles(geotiff, output_dir, max_concurrency, resolution): try: colored_dem, hillshade_dem, colored_hillshade_dem = generate_colored_hillshade(geotiff) - generate_tiles(colored_hillshade_dem, output_dir, max_concurrency) + generate_tiles(colored_hillshade_dem, output_dir, max_concurrency, resolution) # Cleanup for f in [colored_dem, hillshade_dem, colored_hillshade_dem]: diff --git a/stages/odm_dem.py b/stages/odm_dem.py index d9029eea0..ca242152b 100755 --- a/stages/odm_dem.py +++ b/stages/odm_dem.py @@ -126,7 +126,7 @@ def process(self, args, outputs): pseudogeo.add_pseudo_georeferencing(dem_geotiff_path) if args.tiles: - generate_dem_tiles(dem_geotiff_path, tree.path("%s_tiles" % product), args.max_concurrency) + generate_dem_tiles(dem_geotiff_path, tree.path("%s_tiles" % product), args.max_concurrency, resolution) if args.cog: convert_to_cogeo(dem_geotiff_path, max_workers=args.max_concurrency) diff --git a/stages/odm_orthophoto.py b/stages/odm_orthophoto.py index 44a691de4..cf362bb30 100644 --- a/stages/odm_orthophoto.py +++ b/stages/odm_orthophoto.py @@ -28,10 +28,10 @@ def process(self, args, outputs): if not io.file_exists(tree.odm_orthophoto_tif) or self.rerun(): - resolution = 1.0 / (gsd.cap_resolution(args.orthophoto_resolution, tree.opensfm_reconstruction, - ignore_gsd=args.ignore_gsd, - ignore_resolution=(not reconstruction.is_georeferenced()) and args.ignore_gsd, - has_gcp=reconstruction.has_gcp()) / 100.0) + resolution = gsd.cap_resolution(args.orthophoto_resolution, tree.opensfm_reconstruction, + ignore_gsd=args.ignore_gsd, + ignore_resolution=(not reconstruction.is_georeferenced()) and args.ignore_gsd, + has_gcp=reconstruction.has_gcp()) # odm_orthophoto definitions kwargs = { @@ -39,7 +39,7 @@ def process(self, args, outputs): 'log': tree.odm_orthophoto_log, 'ortho': tree.odm_orthophoto_render, 'corners': tree.odm_orthophoto_corners, - 'res': resolution, + 'res': 1.0 / (resolution/100.0), 'bands': '', 'depth_idx': '', 'inpaint': '', @@ -127,7 +127,7 @@ def process(self, args, outputs): os.path.join(tree.odm_orthophoto, "odm_orthophoto_cut.tif"), blend_distance=20, only_max_coords_feature=True) - orthophoto.post_orthophoto_steps(args, bounds_file_path, tree.odm_orthophoto_tif, tree.orthophoto_tiles) + orthophoto.post_orthophoto_steps(args, bounds_file_path, tree.odm_orthophoto_tif, tree.orthophoto_tiles, resolution) # Generate feathered orthophoto also if args.orthophoto_cutline: diff --git a/stages/splitmerge.py b/stages/splitmerge.py index 1dca676bf..7bdc1a6c8 100644 --- a/stages/splitmerge.py +++ b/stages/splitmerge.py @@ -266,7 +266,7 @@ def process(self, args, outputs): orthophoto_vars = orthophoto.get_orthophoto_vars(args) orthophoto.merge(all_orthos_and_ortho_cuts, tree.odm_orthophoto_tif, orthophoto_vars) - orthophoto.post_orthophoto_steps(args, merged_bounds_file, tree.odm_orthophoto_tif, tree.orthophoto_tiles) + orthophoto.post_orthophoto_steps(args, merged_bounds_file, tree.odm_orthophoto_tif, tree.orthophoto_tiles, args.orthophoto_resolution) elif len(all_orthos_and_ortho_cuts) == 1: # Simply copy log.ODM_WARNING("A single orthophoto/cutline pair was found between all submodels.") @@ -306,7 +306,7 @@ def merge_dems(dem_filename, human_name): log.ODM_INFO("Created %s" % dem_file) if args.tiles: - generate_dem_tiles(dem_file, tree.path("%s_tiles" % human_name.lower()), args.max_concurrency) + generate_dem_tiles(dem_file, tree.path("%s_tiles" % human_name.lower()), args.max_concurrency, args.orthophoto_resolution) if args.cog: convert_to_cogeo(dem_file, max_workers=args.max_concurrency) From a89803c2eb211354bd1b6124d42a5565d5ac3f65 Mon Sep 17 00:00:00 2001 From: Merten Fermont Date: Sun, 15 Oct 2023 23:26:56 +0200 Subject: [PATCH 2/3] Use dem instead of orthophoto resolution for generating DEM tiles --- stages/splitmerge.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/stages/splitmerge.py b/stages/splitmerge.py index 7bdc1a6c8..fef7a2321 100644 --- a/stages/splitmerge.py +++ b/stages/splitmerge.py @@ -306,7 +306,7 @@ def merge_dems(dem_filename, human_name): log.ODM_INFO("Created %s" % dem_file) if args.tiles: - generate_dem_tiles(dem_file, tree.path("%s_tiles" % human_name.lower()), args.max_concurrency, args.orthophoto_resolution) + generate_dem_tiles(dem_file, tree.path("%s_tiles" % human_name.lower()), args.max_concurrency, args.dem_resolution) if args.cog: convert_to_cogeo(dem_file, max_workers=args.max_concurrency) From f70e55c9eba6cebdffa4aea7b5a5c14c591c8b7c Mon Sep 17 00:00:00 2001 From: Merten Fermont Date: Mon, 16 Oct 2023 09:59:48 +0200 Subject: [PATCH 3/3] Limit maximum tiler zoom level to 23 --- opendm/tiles/tiler.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/opendm/tiles/tiler.py b/opendm/tiles/tiler.py index f091e7146..349c1ed5d 100644 --- a/opendm/tiles/tiler.py +++ b/opendm/tiles/tiler.py @@ -12,7 +12,7 @@ def generate_tiles(geotiff, output_dir, max_concurrency, resolution): zoom = math.ceil(math.log(resolution_equator_cm/resolution, 2)) min_zoom = 5 # 4.89 km/px - max_zoom = min(zoom, 30) # No deeper zoom than 30 (0.0146 cm/px) + max_zoom = min(zoom, 23) # No deeper zoom than 23 (1.86 cm/px at equator) gdal2tiles = os.path.join(os.path.dirname(__file__), "gdal2tiles.py") system.run('%s "%s" --processes %s -z %s-%s -n -w none "%s" "%s"' % (sys.executable, gdal2tiles, max_concurrency, min_zoom, max_zoom, geotiff, output_dir))