kopia lustrzana https://github.com/OpenDroneMap/ODM
Piero Toffanin
rodzic
2f941b58f4
commit
369d538509
|
|
@ -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
|
||||
|
||||
|
|
@ -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
|
||||
|
||||
|
|
@ -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)
|
||||
|
|
|
|||
|
|
@ -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, ...)
|
||||
|
|
|
|||
|
|
@ -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 = {
|
||||
|
|
|
|||
Ładowanie…
Reference in New Issue