Add support for generating cloud optimized GeoTIFFs (#1292)

* Add support for generating cloud optimized GeoTIFFs

* Generate COGs in DEM stage, not in create_dem

* Bump version

* Add cogeo.py

VERSION

@@ -1 +1 @@
-2.5.2
+2.5.3

opendm/cogeo.py

@@ -0,0 +1,63 @@
+import os
+import shutil
+from opendm import system
+from opendm.concurrency import get_max_memory
+from opendm import io
+from opendm import log
+
+def convert_to_cogeo(src_path, blocksize=256, max_workers=1):
+    """
+    Guarantee that the .tif passed as an argument is a Cloud Optimized GeoTIFF (cogeo)
+    The file is destructively converted into a cogeo.
+    If the file cannot be converted, the function does not change the file
+    :param src_path: path to GeoTIFF
+    :return: True on success
+    """
+
+    if not os.path.isfile(src_path):
+        logger.warning("Cannot convert to cogeo: %s (file does not exist)" % src_path)
+        return False
+
+    log.ODM_INFO("Optimizing %s as Cloud Optimized GeoTIFF" % src_path)
+
+
+    tmpfile = io.related_file_path(src_path, postfix='_cogeo')
+    swapfile = io.related_file_path(src_path, postfix='_cogeo_swap')
+
+    kwargs = {
+        'threads': max_workers if max_workers else 'ALL_CPUS',
+        'blocksize': blocksize,
+        'max_memory': get_max_memory(),
+        'src_path': src_path,
+        'tmpfile': tmpfile,
+    }
+
+    try:
+        system.run("gdal_translate "
+                "-of COG "
+                "-co NUM_THREADS={threads} "
+                "-co BLOCKSIZE={blocksize} "
+                "-co COMPRESS=deflate "
+                "-co BIGTIFF=IF_SAFER "
+                "-co RESAMPLING=NEAREST "
+                "--config GDAL_CACHEMAX {max_memory}% "
+                "--config GDAL_NUM_THREADS {threads} "
+                "\"{src_path}\" \"{tmpfile}\" ".format(**kwargs))
+    except Exception as e:
+        log.ODM_WARNING("Cannot create Cloud Optimized GeoTIFF: %s" % str(e))
+
+    if os.path.isfile(tmpfile):
+        shutil.move(src_path, swapfile) # Move to swap location
+
+        try:
+            shutil.move(tmpfile, src_path)
+        except IOError as e:
+            log.ODM_WARNING("Cannot move %s to %s: %s" % (tmpfile, src_path, str(e)))
+            shutil.move(swapfile, src_path) # Attempt to restore
+
+        if os.path.isfile(swapfile):
+            os.remove(swapfile)
+
+        return True
+    else:
+        return False

opendm/config.py

@@ -581,6 +581,13 @@ def config(argv=None, parser=None):
                         default=False,
                         help='Build orthophoto overviews for faster display in programs such as QGIS. Default: %(default)s')
 
+    parser.add_argument('--cog',
+                        action=StoreTrue,
+                        nargs=0,
+                        default=False,
+                        help='Create Cloud-Optimized GeoTIFFs instead of normal GeoTIFFs. Default: %(default)s')
+
+
     parser.add_argument('--verbose', '-v',
                         action=StoreTrue,
                         nargs=0,

opendm/dem/commands.py

@@ -270,7 +270,7 @@ def process_tile(q):
         if os.path.exists(cleanup_file): os.remove(cleanup_file)
     for t in tiles:
         if os.path.exists(t['filename']): os.remove(t['filename'])
-    
+
     log.ODM_INFO('Completed %s in %s' % (output_file, datetime.now() - start))
 
 

opendm/orthophoto.py

@@ -12,6 +12,7 @@
 from rasterio.mask import mask
 from opendm import io
 from opendm.tiles.tiler import generate_orthophoto_tiles
+from opendm.cogeo import convert_to_cogeo
 from osgeo import gdal
 
 
@@ -72,7 +73,7 @@ def post_orthophoto_steps(args, bounds_file_path, orthophoto_file, orthophoto_ti
     if args.crop > 0:
         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:
+    if args.build_overviews and not args.cog:
         build_overviews(orthophoto_file)
 
     if args.orthophoto_png:
@@ -84,6 +85,8 @@ def post_orthophoto_steps(args, bounds_file_path, orthophoto_file, orthophoto_ti
     if args.tiles:
         generate_orthophoto_tiles(orthophoto_file, orthophoto_tiles_dir, args.max_concurrency)
 
+    if args.cog:
+        convert_to_cogeo(orthophoto_file, max_workers=args.max_concurrency)
 
 def compute_mask_raster(input_raster, vector_mask, output_raster, blend_distance=20, only_max_coords_feature=False):
     if not os.path.exists(input_raster):

opendm/osfm.py

@@ -451,7 +451,7 @@ def get_submodel_argv(args, submodels_path = None, submodel_name = None):
         reading the contents of --cameras
     """
     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']
+    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']
 
     argv = sys.argv

stages/odm_dem.py

@@ -11,6 +11,7 @@
 from opendm.cropper import Cropper
 from opendm import pseudogeo
 from opendm.tiles.tiler import generate_dem_tiles
+from opendm.cogeo import convert_to_cogeo
 
 class ODMDEMStage(types.ODM_Stage):
     def process(self, args, outputs):
@@ -126,6 +127,9 @@ def process(self, args, outputs):
                     if args.tiles:
                         generate_dem_tiles(dem_geotiff_path, tree.path("%s_tiles" % product), args.max_concurrency)
 
+                    if args.cog:
+                        convert_to_cogeo(dem_geotiff_path, max_workers=args.max_concurrency)
+
                     progress += 30
                     self.update_progress(progress)
             else:

stages/splitmerge.py

@@ -19,6 +19,7 @@
 from opendm import point_cloud
 from opendm.utils import double_quote
 from opendm.tiles.tiler import generate_dem_tiles
+from opendm.cogeo import convert_to_cogeo
 
 class ODMSplitStage(types.ODM_Stage):
     def process(self, args, outputs):
@@ -337,6 +338,9 @@ def merge_dems(dem_filename, human_name):
 
                         if args.tiles:
                             generate_dem_tiles(dem_file, tree.path("%s_tiles" % human_name.lower()), args.max_concurrency)
+
+                        if args.cog:
+                            convert_to_cogeo(dem_file, max_workers=args.max_concurrency)
                     else:
                         log.ODM_WARNING("Cannot merge %s, %s was not created" % (human_name, dem_file))