diff --git a/opendm/dem/commands.py b/opendm/dem/commands.py index dcd77dc8..2f2a7b10 100644 --- a/opendm/dem/commands.py +++ b/opendm/dem/commands.py @@ -183,6 +183,7 @@ def create_dem(input_point_cloud, dem_type, output_type='max', radiuses=['0.56'] vrt_path = os.path.abspath(os.path.join(outdir, "merged.vrt")) run('gdalbuildvrt "%s" "%s"' % (vrt_path, '" "'.join(map(lambda t: t['filename'], tiles)))) + geotiff_tmp_path = os.path.abspath(os.path.join(outdir, 'merged.tmp.tif')) geotiff_path = os.path.abspath(os.path.join(outdir, 'merged.tif')) # Build GeoTIFF @@ -190,16 +191,25 @@ def create_dem(input_point_cloud, dem_type, output_type='max', radiuses=['0.56'] 'max_memory': get_max_memory(), 'threads': max_workers if max_workers else 'ALL_CPUS', 'vrt': vrt_path, - 'geotiff': geotiff_path + 'geotiff': geotiff_path, + 'geotiff_tmp': geotiff_tmp_path } if gapfill: + # Sometimes, for some reason gdal_fillnodata.py + # behaves strangely when reading data directly from a .VRT + # so we need to convert to GeoTIFF first. + run('gdal_translate ' + '-co NUM_THREADS={threads} ' + '--config GDAL_CACHEMAX {max_memory}% ' + '{vrt} {geotiff_tmp}'.format(**kwargs)) + run('gdal_fillnodata.py ' '-co NUM_THREADS={threads} ' '--config GDAL_CACHEMAX {max_memory}% ' '-b 1 ' '-of GTiff ' - '{vrt} {geotiff}'.format(**kwargs)) + '{geotiff_tmp} {geotiff}'.format(**kwargs)) else: run('gdal_translate ' '-co NUM_THREADS={threads} ' @@ -209,6 +219,7 @@ def create_dem(input_point_cloud, dem_type, output_type='max', radiuses=['0.56'] post_process(geotiff_path, output_path) os.remove(geotiff_path) + if os.path.exists(geotiff_tmp_path): os.remove(geotiff_tmp_path) if os.path.exists(vrt_path): os.remove(vrt_path) for t in tiles: if os.path.exists(t['filename']): os.remove(t['filename']) diff --git a/opendm/mesh.py b/opendm/mesh.py index d17e4470..703ce434 100644 --- a/opendm/mesh.py +++ b/opendm/mesh.py @@ -74,12 +74,22 @@ def dem_to_points(inGeotiff, outPointCloud, verbose=False): return outPointCloud -def dem_to_mesh_gridded(inGeotiff, outPointCloud, maxVertexCount, verbose=False): +def dem_to_mesh_gridded(inGeotiff, outMesh, maxVertexCount, verbose=False): log.ODM_INFO('Creating mesh from DSM: %s' % inGeotiff) + mesh_path, mesh_filename = os.path.split(outMesh) + # mesh_path = path/to + # mesh_filename = odm_mesh.ply + + basename, ext = os.path.splitext(mesh_filename) + # basename = odm_mesh + # ext = .ply + + outMeshDirty = os.path.join(mesh_path, "{}.dirty{}".format(basename, ext)) + kwargs = { 'bin': context.dem2mesh_path, - 'outfile': outPointCloud, + 'outfile': outMeshDirty, 'infile': inGeotiff, 'maxVertexCount': maxVertexCount, 'verbose': '-verbose' if verbose else '' @@ -90,7 +100,25 @@ def dem_to_mesh_gridded(inGeotiff, outPointCloud, maxVertexCount, verbose=False) '-maxVertexCount {maxVertexCount} ' ' {verbose} '.format(**kwargs)) - return outPointCloud + # Cleanup and reduce vertex count if necessary + # (as dem2mesh cannot guarantee that we'll have the target vertex count) + cleanupArgs = { + 'bin': context.odm_modules_path, + 'outfile': outMesh, + 'infile': outMeshDirty, + 'max_vertex': maxVertexCount, + 'verbose': '-verbose' if verbose else '' + } + + system.run('{bin}/odm_cleanmesh -inputFile {infile} ' + '-outputFile {outfile} ' + '-removeIslands ' + '-decimateMesh {max_vertex} {verbose} '.format(**cleanupArgs)) + + # Delete intermediate results + os.remove(outMeshDirty) + + return outMesh def screened_poisson_reconstruction(inPointCloud, outMesh, depth = 8, samples = 1, maxVertexCount=100000, pointWeight=4, threads=context.num_cores, verbose=False):