diff --git a/opendm/clipper.py b/opendm/clipper.py deleted file mode 100644 index 3b21874a..00000000 --- a/opendm/clipper.py +++ /dev/null @@ -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 - diff --git a/opendm/cropper.py b/opendm/cropper.py new file mode 100644 index 00000000..14c108ad --- /dev/null +++ b/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 + diff --git a/scripts/odm_dem.py b/scripts/odm_dem.py index 6b81b094..eb30c07a 100644 --- a/scripts/odm_dem.py +++ b/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): @@ -62,9 +61,7 @@ class ODMDEMCell(ecto.Cell): if (args.dtm and not io.file_exists(dtm_output_filename)) or \ (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) - kwargs['site'] = '-s {}'.format(bounds_buffer_path) + bounds_shapefile_path = os.path.join(tree.odm_georeferencing, 'odm_georeferenced_model.bounds.shp') + 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) diff --git a/scripts/odm_georeferencing.py b/scripts/odm_georeferencing.py index dd327dfc..a7e2cf9b 100644 --- a/scripts/odm_georeferencing.py +++ b/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, ...) diff --git a/scripts/odm_orthophoto.py b/scripts/odm_orthophoto.py index 4174af42..b1470e96 100644 --- a/scripts/odm_orthophoto.py +++ b/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 = {