Merge pull request #1883 from pierotofy/cogfix

Fix COGs

opendm/orthophoto.py

@@ -13,6 +13,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 opendm.utils import add_raster_meta_tags
 from osgeo import gdal
 from osgeo import ogr
 
@@ -166,7 +167,7 @@ def generate_tfw(orthophoto_file):
         log.ODM_WARNING("Cannot create .tfw for %s: %s" % (orthophoto_file, str(e)))
 
 
-def post_orthophoto_steps(args, bounds_file_path, orthophoto_file, orthophoto_tiles_dir, resolution):
+def post_orthophoto_steps(args, bounds_file_path, orthophoto_file, orthophoto_tiles_dir, resolution, reconstruction, tree, embed_gcp_meta=False):
     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'])
 
@@ -179,6 +180,8 @@ def post_orthophoto_steps(args, bounds_file_path, orthophoto_file, orthophoto_ti
     if args.orthophoto_kmz:
         generate_kmz(orthophoto_file)
 
+    add_raster_meta_tags(orthophoto_file, reconstruction, tree, embed_gcp_meta=embed_gcp_meta)
+
     if args.tiles:
         generate_orthophoto_tiles(orthophoto_file, orthophoto_tiles_dir, args.max_concurrency, resolution)
 

opendm/utils.py

@@ -1,8 +1,12 @@
 import os, shutil
 import numpy as np
 import json
+import rasterio
+from osgeo import gdal
+from datetime import datetime
+
 from opendm import log
-from opendm.photo import find_largest_photo_dims
+from opendm.photo import find_largest_photo_dims, find_mean_utc_time
 from osgeo import gdal
 from opendm.arghelpers import double_quote
 
@@ -113,4 +117,42 @@ def np_to_json(arr):
     return json.dumps(arr, cls=NumpyEncoder)
 
 def np_from_json(json_dump):
-    return np.asarray(json.loads(json_dump))
+    return np.asarray(json.loads(json_dump))
+
+def add_raster_meta_tags(raster, reconstruction, tree, embed_gcp_meta=True):
+    try:
+        if os.path.isfile(raster):
+            mean_capture_time = find_mean_utc_time(reconstruction.photos)
+            mean_capture_dt = None
+            if mean_capture_time is not None:
+                mean_capture_dt = datetime.fromtimestamp(mean_capture_time).strftime('%Y:%m:%d %H:%M:%S') + '+00:00'
+
+            log.ODM_INFO("Adding TIFFTAGs to {}".format(raster))
+            with rasterio.open(raster, 'r+') as rst:
+                if mean_capture_dt is not None:
+                    rst.update_tags(TIFFTAG_DATETIME=mean_capture_dt)
+                rst.update_tags(TIFFTAG_SOFTWARE='ODM {}'.format(log.odm_version()))
+
+            if embed_gcp_meta:
+                # Embed GCP info in 2D results via
+                # XML metadata fields
+                gcp_gml_export_file = tree.path("odm_georeferencing", "ground_control_points.gml")
+
+                if reconstruction.has_gcp() and os.path.isfile(gcp_gml_export_file):
+                    gcp_xml = ""
+
+                    with open(gcp_gml_export_file) as f:
+                        gcp_xml = f.read()
+
+                    ds = gdal.Open(raster)
+                    if ds is not None:
+                        if ds.GetMetadata('xml:GROUND_CONTROL_POINTS') is None or self.rerun():
+                            ds.SetMetadata(gcp_xml, 'xml:GROUND_CONTROL_POINTS')
+                            ds = None
+                            log.ODM_INFO("Wrote xml:GROUND_CONTROL_POINTS metadata to %s" % raster)
+                        else:
+                            log.ODM_WARNING("Already embedded ground control point information")
+                    else:
+                        log.ODM_WARNING("Cannot open %s for writing, skipping GCP embedding" % raster)
+    except Exception as e:
+        log.ODM_WARNING("Cannot write raster meta tags to %s: %s" % (raster, str(e)))

stages/odm_dem.py

@@ -12,6 +12,8 @@ from opendm.cropper import Cropper
 from opendm import pseudogeo
 from opendm.tiles.tiler import generate_dem_tiles
 from opendm.cogeo import convert_to_cogeo
+from opendm.utils import add_raster_meta_tags
+
 
 class ODMDEMStage(types.ODM_Stage):
     def process(self, args, outputs):
@@ -86,6 +88,8 @@ class ODMDEMStage(types.ODM_Stage):
 
                     if pseudo_georeference:
                         pseudogeo.add_pseudo_georeferencing(dem_geotiff_path)
+                    
+                    add_raster_meta_tags(dem_geotiff_path, reconstruction, tree, embed_gcp_meta=not outputs['large'])
 
                     if args.tiles:
                         generate_dem_tiles(dem_geotiff_path, tree.path("%s_tiles" % product), args.max_concurrency, resolution)

stages/odm_orthophoto.py

@@ -132,7 +132,8 @@ class ODMOrthoPhotoStage(types.ODM_Stage):
                     else:
                         log.ODM_INFO("Not a submodel run, skipping mask raster generation")
 
-                orthophoto.post_orthophoto_steps(args, bounds_file_path, tree.odm_orthophoto_tif, tree.orthophoto_tiles, resolution)
+                orthophoto.post_orthophoto_steps(args, bounds_file_path, tree.odm_orthophoto_tif, tree.orthophoto_tiles, resolution, 
+                    reconstruction, tree, not outputs["large"])
 
                 # Generate feathered orthophoto also
                 if args.orthophoto_cutline and submodel_run:

stages/odm_postprocess.py

@@ -1,12 +1,8 @@
 import os
-import rasterio
 
-from datetime import datetime
-from osgeo import gdal
 from opendm import io
 from opendm import log
 from opendm import types
-from opendm import photo
 from opendm.utils import copy_paths, get_processing_results_paths
 from opendm.ogctiles import build_3dtiles
 
@@ -17,54 +13,6 @@ class ODMPostProcess(types.ODM_Stage):
 
         log.ODM_INFO("Post Processing")
 
-        rasters = [tree.odm_orthophoto_tif,
-                    tree.path("odm_dem", "dsm.tif"),
-                    tree.path("odm_dem", "dtm.tif")]
-
-        mean_capture_time = photo.find_mean_utc_time(reconstruction.photos)
-        mean_capture_dt = None
-        if mean_capture_time is not None:
-            mean_capture_dt = datetime.fromtimestamp(mean_capture_time).strftime('%Y:%m:%d %H:%M:%S') + '+00:00'
-
-        # Add TIFF tags
-        for product in rasters:
-            if os.path.isfile(product):
-                log.ODM_INFO("Adding TIFFTAGs to {}".format(product))
-                with rasterio.open(product, 'r+') as rst:
-                    if mean_capture_dt is not None:
-                        rst.update_tags(TIFFTAG_DATETIME=mean_capture_dt)
-                    rst.update_tags(TIFFTAG_SOFTWARE='ODM {}'.format(log.odm_version()))
-
-        # GCP info
-        if not outputs['large']:
-            # TODO: support for split-merge?
-
-            # Embed GCP info in 2D results via
-            # XML metadata fields
-            gcp_gml_export_file = tree.path("odm_georeferencing", "ground_control_points.gml")
-
-            if reconstruction.has_gcp() and io.file_exists(gcp_gml_export_file):
-                skip_embed_gcp = False
-                gcp_xml = ""
-
-                with open(gcp_gml_export_file) as f:
-                    gcp_xml = f.read()
-
-                for product in rasters:
-                    if os.path.isfile(product):
-                        ds = gdal.Open(product)
-                        if ds is not None:
-                            if ds.GetMetadata('xml:GROUND_CONTROL_POINTS') is None or self.rerun():
-                                ds.SetMetadata(gcp_xml, 'xml:GROUND_CONTROL_POINTS')
-                                ds = None
-                                log.ODM_INFO("Wrote xml:GROUND_CONTROL_POINTS metadata to %s" % product)
-                            else:
-                                skip_embed_gcp = True
-                                log.ODM_WARNING("Already embedded ground control point information")
-                                break
-                        else:
-                            log.ODM_WARNING("Cannot open %s for writing, skipping GCP embedding" % product)
-
         if getattr(args, '3d_tiles'):
             build_3dtiles(args, tree, reconstruction, self.rerun())
 

stages/splitmerge.py

@@ -13,7 +13,7 @@ from opendm.cropper import Cropper
 from opendm.remote import LocalRemoteExecutor
 from opendm.shots import merge_geojson_shots, merge_cameras
 from opendm import point_cloud
-from opendm.utils import double_quote
+from opendm.utils import double_quote, add_raster_meta_tags
 from opendm.tiles.tiler import generate_dem_tiles
 from opendm.cogeo import convert_to_cogeo
 from opendm import multispectral
@@ -263,7 +263,8 @@ 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, args.orthophoto_resolution)
+                        orthophoto.post_orthophoto_steps(args, merged_bounds_file, tree.odm_orthophoto_tif, tree.orthophoto_tiles, args.orthophoto_resolution,
+                            reconstruction, tree, False)
                     elif len(all_orthos_and_ortho_cuts) == 1:
                         # Simply copy
                         log.ODM_WARNING("A single orthophoto/cutline pair was found between all submodels.")
@@ -305,6 +306,8 @@ class ODMMergeStage(types.ODM_Stage):
                         if args.tiles:
                             generate_dem_tiles(dem_file, tree.path("%s_tiles" % human_name.lower()), args.max_concurrency, args.dem_resolution)
                         
+                        add_raster_meta_tags(dem_file, reconstruction, tree, embed_gcp_meta=False)
+
                         if args.cog:
                             convert_to_cogeo(dem_file, max_workers=args.max_concurrency)
                     else: