From c679d400c820c4cfd38d61168667284f158f9fc0 Mon Sep 17 00:00:00 2001 From: Merten Fermont Date: Thu, 12 Oct 2023 22:06:53 +0200 Subject: [PATCH] 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 | 23 ++++++++++++++++------- stages/odm_dem.py | 2 +- stages/odm_orthophoto.py | 12 ++++++------ stages/splitmerge.py | 4 ++-- 5 files changed, 27 insertions(+), 18 deletions(-) diff --git a/opendm/orthophoto.py b/opendm/orthophoto.py index 79865c71..319a504c 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 6b36f209..f091e714 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): - 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)) +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)) -def generate_orthophoto_tiles(geotiff, output_dir, max_concurrency): + 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 %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, 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 d9029eea..ca242152 100755 --- a/stages/odm_dem.py +++ b/stages/odm_dem.py @@ -126,7 +126,7 @@ class ODMDEMStage(types.ODM_Stage): 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 44a691de..cf362bb3 100644 --- a/stages/odm_orthophoto.py +++ b/stages/odm_orthophoto.py @@ -28,10 +28,10 @@ class ODMOrthoPhotoStage(types.ODM_Stage): 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 @@ class ODMOrthoPhotoStage(types.ODM_Stage): '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 @@ class ODMOrthoPhotoStage(types.ODM_Stage): 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 1dca676b..7bdc1a6c 100644 --- a/stages/splitmerge.py +++ b/stages/splitmerge.py @@ -266,7 +266,7 @@ class ODMMergeStage(types.ODM_Stage): 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 @@ class ODMMergeStage(types.ODM_Stage): 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)