Automatic cropping

opendm/clipper.py

@@ -1,114 +0,0 @@
-from opendm import context
-from opendm import system
-from osgeo import ogr
-import json, os
-
-def run(command):
-    env_paths = [context.superbuild_bin_path]
-    return system.run(command, env_paths)
-
-class Clipper:
-    def __init__(self, storageDirectory, filesPrefix = "clip"):
-        self.storageDirectory = storageDirectory
-        self.filesPrefix = filesPrefix
-
-    def path(self, suffix):
-        return os.path.join(self.storageDirectory, '{}.{}'.format(self.filesPrefix, suffix))
-
-    def create_buffer_geojson(self, pointCloudPath, bufferDistance = 0):
-        # Use PDAL to dump boundary information
-        # then read the information back
-
-        boundary_file_path = self.path('boundary.json')
-
-        run('pdal info --boundary --filters.hexbin.edge_length=1 --filters.hexbin.threshold=0 {0} > {1}'.format(pointCloudPath,  boundary_file_path))
-        
-        pc_geojson_boundary_feature = None
-
-        with open(boundary_file_path, 'r') as f:
-            json_f = json.loads(f.read())
-            pc_geojson_boundary_feature = json_f['boundary']['boundary_json']
-
-        if pc_geojson_boundary_feature is None: raise RuntimeError("Could not determine point cloud boundaries")
-
-        # Write bounds to GeoJSON
-        buffer_geojson_path = self.path('bounds.geojson')
-        with open(buffer_geojson_path, "w") as f:
-            f.write(json.dumps({
-                "type": "FeatureCollection",
-                "features": [{
-                    "type": "Feature",
-                    "geometry": pc_geojson_boundary_feature
-                }]
-            }))
-
-        # Create a convex hull around the boundary
-        # as to encompass the entire area (no holes)    
-        driver = ogr.GetDriverByName('GeoJSON')
-        ds = driver.Open(buffer_geojson_path, 0) # ready-only
-        layer = ds.GetLayer()
-
-        # Collect all Geometry
-        geomcol = ogr.Geometry(ogr.wkbGeometryCollection)
-        for feature in layer:
-            geomcol.AddGeometry(feature.GetGeometryRef())
-
-        # Calculate convex hull
-        convexhull = geomcol.ConvexHull()
-
-        # If buffer distance is specified
-        # Create two buffers, one shrinked by
-        # N + 3 and then that buffer expanded by 3
-        # so that we get smooth corners. \m/
-        BUFFER_SMOOTH_DISTANCE = 3
-
-        if bufferDistance > 0:
-            convexhull = convexhull.Buffer(-(bufferDistance + BUFFER_SMOOTH_DISTANCE))
-            convexhull = convexhull.Buffer(BUFFER_SMOOTH_DISTANCE)
-
-        # Save to a new file
-        buffer_geojson_path = self.path('buffer.geojson')
-        if os.path.exists(buffer_geojson_path):
-            driver.DeleteDataSource(buffer_geojson_path)
-
-        out_ds = driver.CreateDataSource(buffer_geojson_path)
-        layer = out_ds.CreateLayer("convexhull", geom_type=ogr.wkbPolygon)
-
-        feature_def = layer.GetLayerDefn()
-        feature = ogr.Feature(feature_def)
-        feature.SetGeometry(convexhull)
-        layer.CreateFeature(feature)
-        feature = None
-
-        # Save and close data sources
-        out_ds = ds = None
-
-        return buffer_geojson_path
-
-
-    def create_buffer_shapefile(self, pointCloudPath, bufferDistance = 0):
-        buffer_geojson_path = self.create_buffer_geojson(pointCloudPath, bufferDistance)
-
-        summary_file_path = os.path.join(self.storageDirectory, '{}.summary.json'.format(self.filesPrefix))
-        run('pdal info --summary {0} > {1}'.format(pointCloudPath, summary_file_path))
-        
-        pc_proj4 = None
-        with open(summary_file_path, 'r') as f:
-            json_f = json.loads(f.read())
-            pc_proj4 = json_f['summary']['srs']['proj4']
-
-        if pc_proj4 is None: raise RuntimeError("Could not determine point cloud proj4 declaration")
-
-        bounds_shapefile_path = os.path.join(self.storageDirectory, '{}.buffer.shp'.format(self.filesPrefix))
-
-        # Convert bounds to Shapefile
-        kwargs = {
-            'input': buffer_geojson_path,
-            'output': bounds_shapefile_path,
-            'proj4': pc_proj4
-        }
-
-        run('ogr2ogr -overwrite -a_srs "{proj4}" {output} {input}'.format(**kwargs))
-
-        return bounds_shapefile_path
-

opendm/cropper.py

@@ -0,0 +1,180 @@
+from opendm import context
+from opendm import system
+from opendm import log
+from osgeo import ogr
+import json, os
+
+def run(command):
+    env_paths = [context.superbuild_bin_path]
+    return system.run(command, env_paths)
+
+class Cropper:
+    def __init__(self, storage_dir, files_prefix = "crop"):
+        self.storage_dir = storage_dir
+        self.files_prefix = files_prefix
+
+    def path(self, suffix):
+        """
+        @return a path relative to storage_dir and prefixed with files_prefix
+        """
+        return os.path.join(self.storage_dir, '{}.{}'.format(self.files_prefix, suffix))
+
+    @staticmethod
+    def crop(shapefile_path, geotiff_path, gdal_options):
+        if not os.path.exists(shapefile_path) or not os.path.exists(geotiff_path):
+            log.ODM_WARNING("Either {} or {} does not exist, will skip cropping.".format(shapefile_path, geotiff_path))
+            return geotiff_path
+
+        # Rename original file
+        # path/to/odm_orthophoto.tif --> path/to/odm_orthophoto.original.tif
+        
+        path, filename = os.path.split(geotiff_path)
+        # path = path/to
+        # filename = odm_orthophoto.tif
+
+        basename, ext = os.path.splitext(filename)
+        # basename = odm_orthophoto
+        # ext = .tif
+
+        original_geotiff = os.path.join(path, "{}.original{}".format(basename, ext))
+        os.rename(geotiff_path, original_geotiff)
+
+        try:
+            kwargs = {
+                'shapefile_path': shapefile_path,
+                'geotiffInput': original_geotiff,
+                'geotiffOutput': geotiff_path,
+                'options': ' '.join(map(lambda k: '-co {}={}'.format(k, gdal_options[k]), gdal_options))
+            }
+
+            run('gdalwarp -cutline {shapefile_path} '
+                '-crop_to_cutline '
+                '{options} '
+                '{geotiffInput} '
+                '{geotiffOutput} '.format(**kwargs))
+
+        except Exception as e:
+            log.ODM_WARNING('Something went wrong while cropping: {}'.format(e.message))
+            
+            # Revert rename
+            os.rename(original_geotiff, geotiff_path)
+
+        return geotiff_path
+
+    def create_bounds_geojson(self, pointcloud_path, buffer_distance = 0):
+        """
+        Compute a buffered polygon around the data extents (not just a bounding box)
+        of the given point cloud.
+
+        @return filename to GeoJSON containing the polygon
+        """
+        if not os.path.exists(pointcloud_path):
+            log.ODM_WARNING('Point cloud does not exist, cannot generate shapefile bounds {}'.format(pointcloud_path))
+            return ''
+
+        # Use PDAL to dump boundary information
+        # then read the information back
+
+        boundary_file_path = self.path('boundary.json')
+
+        run('pdal info --boundary --filters.hexbin.edge_length=1 --filters.hexbin.threshold=0 {0} > {1}'.format(pointcloud_path,  boundary_file_path))
+        
+        pc_geojson_boundary_feature = None
+
+        with open(boundary_file_path, 'r') as f:
+            json_f = json.loads(f.read())
+            pc_geojson_boundary_feature = json_f['boundary']['boundary_json']
+
+        if pc_geojson_boundary_feature is None: raise RuntimeError("Could not determine point cloud boundaries")
+
+        # Write bounds to GeoJSON
+        bounds_geojson_path = self.path('bounds.geojson')
+        with open(bounds_geojson_path, "w") as f:
+            f.write(json.dumps({
+                "type": "FeatureCollection",
+                "features": [{
+                    "type": "Feature",
+                    "geometry": pc_geojson_boundary_feature
+                }]
+            }))
+
+        # Create a convex hull around the boundary
+        # as to encompass the entire area (no holes)    
+        driver = ogr.GetDriverByName('GeoJSON')
+        ds = driver.Open(bounds_geojson_path, 0) # ready-only
+        layer = ds.GetLayer()
+
+        # Collect all Geometry
+        geomcol = ogr.Geometry(ogr.wkbGeometryCollection)
+        for feature in layer:
+            geomcol.AddGeometry(feature.GetGeometryRef())
+
+        # Calculate convex hull
+        convexhull = geomcol.ConvexHull()
+
+        # If buffer distance is specified
+        # Create two buffers, one shrinked by
+        # N + 3 and then that buffer expanded by 3
+        # so that we get smooth corners. \m/
+        BUFFER_SMOOTH_DISTANCE = 3
+
+        if buffer_distance > 0:
+            convexhull = convexhull.Buffer(-(buffer_distance + BUFFER_SMOOTH_DISTANCE))
+            convexhull = convexhull.Buffer(BUFFER_SMOOTH_DISTANCE)
+
+        # Save to a new file
+        bounds_geojson_path = self.path('bounds.geojson')
+        if os.path.exists(bounds_geojson_path):
+            driver.DeleteDataSource(bounds_geojson_path)
+
+        out_ds = driver.CreateDataSource(bounds_geojson_path)
+        layer = out_ds.CreateLayer("convexhull", geom_type=ogr.wkbPolygon)
+
+        feature_def = layer.GetLayerDefn()
+        feature = ogr.Feature(feature_def)
+        feature.SetGeometry(convexhull)
+        layer.CreateFeature(feature)
+        feature = None
+
+        # Save and close data sources
+        out_ds = ds = None
+
+        return bounds_geojson_path
+
+
+    def create_bounds_shapefile(self, pointcloud_path, buffer_distance = 0):
+        """
+        Compute a buffered polygon around the data extents (not just a bounding box)
+        of the given point cloud.
+        
+        @return filename to Shapefile containing the polygon
+        """
+        if not os.path.exists(pointcloud_path):
+            log.ODM_WARNING('Point cloud does not exist, cannot generate shapefile bounds {}'.format(pointcloud_path))
+            return ''
+
+        bounds_geojson_path = self.create_bounds_geojson(pointcloud_path, buffer_distance)
+
+        summary_file_path = os.path.join(self.storage_dir, '{}.summary.json'.format(self.files_prefix))
+        run('pdal info --summary {0} > {1}'.format(pointcloud_path, summary_file_path))
+        
+        pc_proj4 = None
+        with open(summary_file_path, 'r') as f:
+            json_f = json.loads(f.read())
+            pc_proj4 = json_f['summary']['srs']['proj4']
+
+        if pc_proj4 is None: raise RuntimeError("Could not determine point cloud proj4 declaration")
+
+        bounds_shapefile_path = os.path.join(self.storage_dir, '{}.bounds.shp'.format(self.files_prefix))
+
+        # Convert bounds to Shapefile
+        kwargs = {
+            'input': bounds_geojson_path,
+            'output': bounds_shapefile_path,
+            'proj4': pc_proj4
+        }
+
+        run('ogr2ogr -overwrite -a_srs "{proj4}" {output} {input}'.format(**kwargs))
+
+        return bounds_shapefile_path
+

scripts/odm_dem.py

@@ -6,7 +6,6 @@ from opendm import log
 from opendm import system
 from opendm import context
 from opendm import types
-from opendm.clipper import Clipper
 
 
 class ODMDEMCell(ecto.Cell):
@@ -63,8 +62,6 @@ class ODMDEMCell(ecto.Cell):
                 (args.dsm and not io.file_exists(dsm_output_filename)) or \
                 rerun_cell:
                  
-                clipper = Clipper(odm_dem_root, 'odm_georeferenced_model')
-
                 # Process with lidar2dems
                 terrain_params_map = {
                     'flatnonforest': (1, 3), 
@@ -83,8 +80,9 @@ class ODMDEMCell(ecto.Cell):
                 }
 
                 if args.crop > 0:
-                    bounds_buffer_path = clipper.create_buffer_shapefile(tree.odm_georeferencing_model_las, args.crop)
+                    bounds_shapefile_path = os.path.join(tree.odm_georeferencing, 'odm_georeferenced_model.bounds.shp')
-                    kwargs['site'] = '-s {}'.format(bounds_buffer_path)
+                    if os.path.exists(bounds_shapefile_path):
+                        kwargs['site'] = '-s {}'.format(bounds_shapefile_path)
 
                 l2d_params = '--slope {slope} --cellsize {cellsize} ' \
                              '{verbose} ' \
@@ -103,7 +101,8 @@ class ODMDEMCell(ecto.Cell):
                         args.dem_initial_distance, tree.odm_georeferencing), env_paths)
                 else:
                     log.ODM_INFO("Will skip classification, only DSM is needed")
-                    copyfile(tree.odm_georeferencing_model_las, os.path.join(odm_dem_root, 'bounds-0_l2d_s{slope}c{cellsize}.las'.format(**kwargs)))
+                    l2d_classified_pattern = 'odm_georeferenced_model.bounds-0_l2d_s{slope}c{cellsize}.las' if args.crop > 0 else 'l2d_s{slope}c{cellsize}.las'
+                    copyfile(tree.odm_georeferencing_model_las, os.path.join(odm_dem_root, l2d_classified_pattern.format(**kwargs)))
 
                 products = []
                 if args.dsm: products.append('dsm') 
@@ -136,9 +135,11 @@ class ODMDEMCell(ecto.Cell):
 
                     # Rename final output
                     if product == 'dsm':
-                        os.rename(os.path.join(odm_dem_root, 'bounds-0_dsm.idw.tif'), dsm_output_filename)
+                        dsm_pattern = 'odm_georeferenced_model.bounds-0_dsm.idw.tif' if args.crop > 0 else 'dsm.idw.tif'
+                        os.rename(os.path.join(odm_dem_root, dsm_pattern), dsm_output_filename)
                     elif product == 'dtm':
-                        os.rename(os.path.join(odm_dem_root, 'bounds-0_dtm.idw.tif'), dtm_output_filename)
+                        dtm_pattern = 'odm_georeferenced_model.bounds-0_dsm.idw.tif' if args.crop > 0 else 'dtm.idw.tif'
+                        os.rename(os.path.join(odm_dem_root, dtm_pattern), dtm_output_filename)
 
             else:
                 log.ODM_WARNING('Found existing outputs in: %s' % odm_dem_root)

scripts/odm_georeferencing.py

@@ -7,6 +7,7 @@ from opendm import log
 from opendm import types
 from opendm import system
 from opendm import context
+from opendm.cropper import Cropper
 
 
 class ODMGeoreferencingCell(ecto.Cell):
@@ -188,6 +189,11 @@ class ODMGeoreferencingCell(ecto.Cell):
                                     reachedpoints = True
                     csvfile.close()
 
+                    if args.crop > 0:
+                        log.ODM_INFO("Calculating cropping area and generating bounds shapefile from point cloud")
+                        cropper = Cropper(tree.odm_georeferencing, 'odm_georeferenced_model')
+                        cropper.create_bounds_shapefile(tree.odm_georeferencing_model_las, args.crop)
+
                     # Do not execute a second time, since
                     # We might be doing georeferencing for 
                     # multiple models (3D, 2.5D, ...)

scripts/odm_orthophoto.py

@@ -5,6 +5,7 @@ from opendm import log
 from opendm import system
 from opendm import context
 from opendm import types
+from opendm.cropper import Cropper
 
 
 class ODMOrthoPhotoCell(ecto.Cell):
@@ -128,6 +129,18 @@ class ODMOrthoPhotoCell(ecto.Cell):
                                '-a_srs \"EPSG:{epsg}\" '
                                '{png} {tiff} > {log}'.format(**kwargs))
 
+                    if args.crop > 0:
+                        shapefile_path = os.path.join(tree.odm_georeferencing, 'odm_georeferenced_model.bounds.shp')
+                        Cropper.crop(shapefile_path, tree.odm_orthophoto_tif, {
+                                'TILED': 'NO' if self.params.no_tiled else 'YES',
+                                'COMPRESS': self.params.compress,
+                                'PREDICTOR': '2' if self.params.compress in ['LZW', 'DEFLATE'] else '1',
+                                'BIGTIFF': self.params.bigtiff,
+                                'BLOCKXSIZE': 512,
+                                'BLOCKYSIZE': 512,
+                                'NUM_THREADS': 'ALL_CPUS'
+                            })
+
                     if self.params.build_overviews:
                         log.ODM_DEBUG("Building Overviews")
                         kwargs = {