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

@@ -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)

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]:

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)

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:

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)