kopia lustrzana https://github.com/OpenDroneMap/ODM
Point cloud, orthophoto boundaries working
Piero Toffanin
rodzic
8791b74b73
commit
5dc2e224ce
|
|
@ -1,22 +1,67 @@
|
|||
import fiona
|
||||
import fiona.crs
|
||||
import os
|
||||
import io
|
||||
import json
|
||||
from opendm import system
|
||||
from pyproj import CRS
|
||||
from opendm.location import transformer
|
||||
from opendm.utils import double_quote
|
||||
|
||||
def load_boundary(boundary_json, reproject_to_crs=None):
|
||||
if not os.path.isfile(boundary_json):
|
||||
raise IOError("Cannot load boundary file: %s does not exist" % boundary_json)
|
||||
|
||||
with fiona.open(boundary_json, 'r') as src:
|
||||
def load_boundary(boundary_json, reproject_to_proj4=None):
|
||||
with fiona.open(io.BytesIO(json.dumps(boundary_json).encode('utf-8')), 'r') as src:
|
||||
if len(src) != 1:
|
||||
raise IOError("Boundary file must have a single polygon (found: %s)" % len(src))
|
||||
raise IOError("Boundary must have a single polygon (found: %s)" % len(src))
|
||||
|
||||
geom = src[0]['geometry']
|
||||
|
||||
if geom['type'] != 'Polygon':
|
||||
raise IOError("Boundary file must have a polygon feature (found: %s)" % geom['type'])
|
||||
raise IOError("Boundary must have a polygon feature (found: %s)" % geom['type'])
|
||||
|
||||
coords = geom['coordinates']
|
||||
rings = geom['coordinates']
|
||||
|
||||
if reproject_to_crs is not None:
|
||||
pass
|
||||
if len(rings) == 0:
|
||||
raise IOError("Boundary geometry has no rings")
|
||||
|
||||
return coords
|
||||
coords = rings[0]
|
||||
if len(coords) == 0:
|
||||
raise IOError("Boundary geometry has no coordinates")
|
||||
|
||||
dimensions = len(coords[0])
|
||||
|
||||
if reproject_to_proj4 is not None:
|
||||
t = transformer(CRS.from_proj4(fiona.crs.to_string(src.crs)),
|
||||
CRS.from_proj4(reproject_to_proj4))
|
||||
coords = [t.TransformPoint(*c)[:dimensions] for c in coords]
|
||||
|
||||
return coords
|
||||
|
||||
def as_polygon(boundary):
|
||||
return "POLYGON((" + ", ".join([" ".join(map(str, c)) for c in boundary]) + "))"
|
||||
|
||||
def export_to_bounds_files(boundary, proj4, bounds_json_file, bounds_gpkg_file):
|
||||
with open(bounds_json_file, "w") as f:
|
||||
f.write(json.dumps({
|
||||
"type": "FeatureCollection",
|
||||
"name": "bounds",
|
||||
"features": [{
|
||||
"type": "Feature",
|
||||
"properties": {},
|
||||
"geometry": {
|
||||
"type": "Polygon",
|
||||
"coordinates": [boundary]
|
||||
}
|
||||
}]
|
||||
}))
|
||||
|
||||
if os.path.isfile(bounds_gpkg_file):
|
||||
os.remove(bounds_gpkg_file)
|
||||
|
||||
kwargs = {
|
||||
'proj4': proj4,
|
||||
'input': double_quote(bounds_json_file),
|
||||
'output': double_quote(bounds_gpkg_file)
|
||||
}
|
||||
|
||||
system.run('ogr2ogr -overwrite -f GPKG -a_srs "{proj4}" {output} {input}'.format(**kwargs))
|
||||
|
||||
|
|
|
|||
|
|
@ -5,6 +5,7 @@ from opendm.point_cloud import export_summary_json
|
|||
from osgeo import ogr
|
||||
import json, os
|
||||
from opendm.concurrency import get_max_memory
|
||||
from opendm.utils import double_quote
|
||||
|
||||
class Cropper:
|
||||
def __init__(self, storage_dir, files_prefix = "crop"):
|
||||
|
|
@ -41,9 +42,9 @@ class Cropper:
|
|||
|
||||
try:
|
||||
kwargs = {
|
||||
'gpkg_path': gpkg_path,
|
||||
'geotiffInput': original_geotiff,
|
||||
'geotiffOutput': geotiff_path,
|
||||
'gpkg_path': double_quote(gpkg_path),
|
||||
'geotiffInput': double_quote(original_geotiff),
|
||||
'geotiffOutput': double_quote(geotiff_path),
|
||||
'options': ' '.join(map(lambda k: '-co {}={}'.format(k, gdal_options[k]), gdal_options)),
|
||||
'warpOptions': ' '.join(warp_options),
|
||||
'max_memory': get_max_memory()
|
||||
|
|
@ -252,10 +253,13 @@ class Cropper:
|
|||
|
||||
bounds_gpkg_path = os.path.join(self.storage_dir, '{}.bounds.gpkg'.format(self.files_prefix))
|
||||
|
||||
if os.path.isfile(bounds_gpkg_path):
|
||||
os.remove(bounds_gpkg_path)
|
||||
|
||||
# Convert bounds to GPKG
|
||||
kwargs = {
|
||||
'input': bounds_geojson_path,
|
||||
'output': bounds_gpkg_path,
|
||||
'input': double_quote(bounds_geojson_path),
|
||||
'output': double_quote(bounds_gpkg_path),
|
||||
'proj4': pc_proj4
|
||||
}
|
||||
|
||||
|
|
|
|||
|
|
@ -70,7 +70,7 @@ def generate_kmz(orthophoto_file, output_file=None, outsize=None):
|
|||
'--config GDAL_CACHEMAX %s%% ' % (orthophoto_file, output_file, bandparam, get_max_memory()))
|
||||
|
||||
def post_orthophoto_steps(args, bounds_file_path, orthophoto_file, orthophoto_tiles_dir):
|
||||
if args.crop > 0:
|
||||
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'])
|
||||
|
||||
if args.build_overviews and not args.cog:
|
||||
|
|
|
|||
|
|
@ -510,10 +510,11 @@ def get_submodel_argv(args, submodels_path = None, submodel_name = None):
|
|||
tweaking --crop if necessary (DEM merging makes assumption about the area of DEMs and their euclidean maps that require cropping. If cropping is skipped, this leads to errors.)
|
||||
removing --gcp (the GCP path if specified is always "gcp_list.txt")
|
||||
reading the contents of --cameras
|
||||
reading the contents of --boundary
|
||||
"""
|
||||
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', 'cog']
|
||||
read_json_always = ['cameras']
|
||||
read_json_always = ['cameras', 'boundary']
|
||||
|
||||
argv = sys.argv
|
||||
|
||||
|
|
|
|||
|
|
@ -69,7 +69,7 @@ class ODM_Reconstruction(object):
|
|||
return self.georef is not None
|
||||
|
||||
def has_gcp(self):
|
||||
return self.is_georeferenced() and self.gcp is not None
|
||||
return self.is_georeferenced() and self.gcp is not None and self.gcp.exists()
|
||||
|
||||
def georeference_with_gcp(self, gcp_file, output_coords_file, output_gcp_file, output_model_txt_geo, rerun=False):
|
||||
if not io.file_exists(output_coords_file) or not io.file_exists(output_gcp_file) or rerun:
|
||||
|
|
|
|||
|
|
@ -9,7 +9,7 @@ from opendm import system
|
|||
from opendm.geo import GeoFile
|
||||
from shutil import copyfile
|
||||
from opendm import progress
|
||||
|
||||
from opendm import boundary
|
||||
|
||||
def save_images_database(photos, database_file):
|
||||
with open(database_file, 'w') as f:
|
||||
|
|
@ -154,3 +154,10 @@ class ODMLoadDatasetStage(types.ODM_Stage):
|
|||
|
||||
reconstruction.save_proj_srs(os.path.join(tree.odm_georeferencing, tree.odm_georeferencing_proj))
|
||||
outputs['reconstruction'] = reconstruction
|
||||
|
||||
# Try to load boundaries
|
||||
if args.boundary:
|
||||
if reconstruction.is_georeferenced():
|
||||
outputs['boundary'] = boundary.load_boundary(args.boundary, reconstruction.get_proj_srs())
|
||||
else:
|
||||
log.ODM_WARNING("Reconstruction is not georeferenced, but boundary file provided (will ignore boundary file)")
|
||||
|
|
|
|||
|
|
@ -17,7 +17,7 @@ from opendm.cropper import Cropper
|
|||
from opendm import point_cloud
|
||||
from opendm.multispectral import get_primary_band_name
|
||||
from opendm.osfm import OSFMContext
|
||||
|
||||
from opendm.boundary import as_polygon, export_to_bounds_files
|
||||
|
||||
class ODMGeoreferencingStage(types.ODM_Stage):
|
||||
def process(self, args, outputs):
|
||||
|
|
@ -128,6 +128,12 @@ class ODMGeoreferencingStage(types.ODM_Stage):
|
|||
'--writers.las.vlrs="{\\\"filename\\\": \\\"%s\\\", \\\"user_id\\\": \\\"ODM_GCP\\\", \\\"description\\\": \\\"Ground Control Points (GML)\\\"}"' % gcp_gml_export_file.replace(os.sep, "/")
|
||||
]
|
||||
|
||||
if 'boundary' in outputs:
|
||||
stages.append("crop")
|
||||
params += [
|
||||
'--filters.crop.polygon="%s"' % as_polygon(outputs['boundary'])
|
||||
]
|
||||
|
||||
system.run(cmd + ' ' + ' '.join(stages) + ' ' + ' '.join(params))
|
||||
|
||||
self.update_progress(50)
|
||||
|
|
@ -152,6 +158,14 @@ class ODMGeoreferencingStage(types.ODM_Stage):
|
|||
except:
|
||||
log.ODM_WARNING("Cannot calculate crop bounds! We will skip cropping")
|
||||
args.crop = 0
|
||||
|
||||
if 'boundary' in outputs and args.crop == 0:
|
||||
log.ODM_INFO("Using boundary JSON as cropping area")
|
||||
|
||||
bounds_base, _ = os.path.splitext(tree.odm_georeferencing_model_laz)
|
||||
bounds_json = bounds_base + ".bounds.geojson"
|
||||
bounds_gpkg = bounds_base + ".bounds.gpkg"
|
||||
export_to_bounds_files(outputs['boundary'], reconstruction.get_proj_srs(), bounds_json, bounds_gpkg)
|
||||
else:
|
||||
log.ODM_INFO("Converting point cloud (non-georeferenced)")
|
||||
system.run(cmd + ' ' + ' '.join(stages) + ' ' + ' '.join(params))
|
||||
|
|
|
|||
Ładowanie…
Reference in New Issue