tileManagerDemo2.py

Author

lwiniwar

Created: 22.03.2016 Copyright: (c) lwiniwar 2016 Licence: opals license

## @package python.demo.tileManagerDemo2
# 
# @brief Demo script for a complex tiling with tileManager
# 
# Takes a set of pointclouds and a set of rasters, creates buffered tiles, smoothes the raster,
# calculated normalized height in the point clouds, cuts the point clouds to net size, merges them.
# 
# 
 
import sys, os
import glob
 
import pprint
 
import opals
from opals import StatFilter, AddInfo, Import
 
from opals.tools import tileManager
from opals.tools.doxygen import brief
from opals.tools.opalsargparse import PkgArgumentParser, Boolean, PathArgument
 
def run_tile(tile, searchRad, tempdir):
    smoothfile = os.path.join(tempdir, tile['tile_name'], 'dtm_smooth.tif')
    # create smoothed file
    sf = StatFilter.StatFilter(inFile=tile['dtm'],
                               outFile=smoothfile,
                               feature="mean",
                               kernelShape=opals.Types.KernelShape.circle,
                               kernelSize=int(searchRad / 1.0))  # pixel size is 1
    sf.run(reset=True)
    # calculate normalized Z based on smoothed file
    add = AddInfo.AddInfo(inFile=tile['pointcloud'],
                          gridFile=smoothfile,
                          attribute='NormalizedZ = Z - r[0]',
                          resampling=opals.Types.ResamplingMethod.bilinear)
    add.run(reset=True)
 
 
@opals.main
def main(options, logger):
    inp = []  #build list of input files (supporting wildcards)
    for i in options.i:
        inp.extend(glob.glob(str(i)))
    if len(inp) == 0:
        raise Exception("No input files found")
    rasters = [] #build list of raster files (supporting wildcards)
    for r in options.r:
        rasters.extend(glob.glob(str(r)))
    if len(rasters) == 0:
        raise Exception("No raster files found")
 
    inp_files = [{'files': inp, # List of Files (relative or full path, no */?)
                  'type': 'odm',  # either odm or raster
                  'name': 'pointcloud',  # any name that can be used as a dict key and a (windows) filename
                  'required': True},
                 {'files': rasters,
                  'type': 'raster',
                  'name': 'dtm',
                  'required': True}
                 ]
    tiling_concept = 'stripsPolygons.shp'  # use existing shapefile
    naming_concept = 'shp.id'  # follow naming of the shapefile
    tempdir = 'temp'
 
    searchRad = options.d
 
 
    # prepare tileManager
    tm = tileManager.tileManager(logger, inp_files, tiling_concept, naming_concept, tempdir, buffer=searchRad)
    # get number of elements - this already cuts the files!
    lix = len(tm)
 
    outfiles = []
    # iterate over the tiles (no more cutting needed here)
    for idx, tile in enumerate(tm):
        print("Now processing tile %s/%s:" % (idx, lix))
        pprint.pprint(tile)
        run_tile(tile, searchRad=searchRad, tempdir=tempdir)
        _, vectorfiles = tm.cutToNetSize(tile, vectorfiles=tile['pointcloud'], overwrite=False)
        outfiles.append(vectorfiles[0])
    # merge the output files into one mosaic
    imp = Import.Import(inFile=outfiles,
                        outFile='mosaic.odm')
    imp.run()
 
if __name__ == '__main__':
    parser = PkgArgumentParser(usage="%(prog)s [options]\n" + brief(__doc__))
    parser.add_argument("-i", type=PathArgument, remarks='mandatory', nargs='*', help="input point cloud files")
    parser.add_argument("-r", type=PathArgument, remarks='mandatory', nargs='*', help="input raster files")
    parser.add_argument("-d", type=float, default=5.0,  help="smoothing radius")
    options = parser.parse_args(sys.argv[1:])
    main(options)