In this tutorial you will learn how to write logfiles and parallelize your processing pipeline for more efficiency.
Introduction
In the previous tutorials we have learned, how we can automate processing tasks with python. In this tutorial you will learn how to do this in a more efficient way with the ability to get logging information in order to locate occuring problems in larger processing.
Step 1: Why Logging
In the moment, we are running our scripts directly within Visual Studio Code. This has the advantage, that we have everything together in one GUI. In the TERMINAL we can monitor our processing job and get both errors and successful executions printed.
However, if you are processing thousands of tiles and your processing job takes multiple days, it is not feasable to sit in front of the computer and wait for errors occuring. Therefore, it is advisable to run your scripts with logging and outside of VS Code.
Step 2: Create a RunScript
In your scripts folder, create a new python script called “run_script.py” and add the following code to the script!
from subprocess import PIPE, call
import os, sys, argparse
logdir = '../log'
# -----------------------------------------------------------------------------
parser = argparse.ArgumentParser(description='run a script with automated logging.')
parser.add_argument('script', type=str, help='the script(s) to execute, e.g. my_script.py or my_script0.py;my_script1.py;my_script2.py')
args = parser.parse_args()
if not os.path.exists(logdir):
os.makedirs(logdir)
arguments = args.script.split(';')
# -----------------------------------------------------------------------------
for argument in arguments:
scriptname, fileext = os.path.splitext(argument)
if fileext == '.py':
command = sys.executable + ' -u ' + argument + ' > ' + logdir + os.sep + argument + '.std.log 2> ' + logdir + os.sep + argument + '.err.log'
print('running: ' + command)
call(command, universal_newlines=True, shell=True)
else:
print('unknown file extension for argument: ' + argument)Step 3: Run your scripts by a RunScript
Go to your Miniforge Prompt and ensure that your conda environment is activated (indicated by (myenv) on the left side)
Run both 01_data_preparation.py and 02_data_processing.py in one go by typing python run_script.py 01_data_preparation.py;02_data_processing.py.
Tip
Start typing the filenames and then use tab on your keyboard in order auto-complete filenames!
In the Miniforge Prompt is indicated which of the selected tools is currently running.
Note
You get no detailed prints here, because all of the processing information is piped into logfiles
If the whole processing pipeline has finished, you see a new line in the miniforge_prompt and a blinking cursor again.
Step 4: Inspect Logfiles
Tip
(For longer processing tasks) If you want to have a look what is currently being processed, you can simply open the logfiles in a text editor and check its status.
You will find a new folder called log in your project folder
If you open this folder, you will find logfiles for each python script that has been executed. You can find std.log and err.log files.
In the std.log files you find the complete log of your script execution. In this case you can see that something went wrong…
In the corresponding err.log only the errors are shown. Thus, you can have a target-orientated overview of all errors. (In this case, there was a problem with the LIS Pro 3D license)
Note
Usually, the err.log files should be empty and have 0 KB. But sometimes there are also errors written to here if a tile has simply no points to process.
Logging Conclusion
You have seen, that you can run your scripts independent from VSCode in the Miniforge Prompt terminal. Using a runscript is very helpful in running scripts on huge projects, running multiple scripts in one go and in providing still a clearly arranged overview of occuring errors.
Step 5: Parallelization introduction
Besides structural improvements of your scripting environment, you have still the possibility of accelerating your processing by parallelization.
This is especially recommended, if you repeat your workflow on multiple tiles (as we do in the the 02_data_processing.py script)
A normal for-loop processes each tile one after another. But if we restructure our code, this will allow running multiple tiles in parallel. The number of parallel threads is limited by your computer system or the maximum number of threads you define as a user.
Warning
The unparallelized processing is usually done in order to save computer memory by chunking down a huge dataset into smaller tiles.
Keep in mind that if you use e.g. 4 threads for parallelized processing, you will use four times as much memory as with the single core processing!
Thus, you have to take care that the amount of data that is processed in parallel still fits into the working memory of your computer (RAM)!
Step 6: Start a parallelized script
Create a new script called 02_data_processing_parallel_tiles.py!
Copy the code of the previous tutorial into it!
import os
import project_config as cfg
from PySAGA import saga_api, lis_cmd
from PySAGA.tools import lis_virtual
from PySAGA.tools import lis_filtering
from PySAGA.tools import lis_classification
from PySAGA.tools import lis_segmentation
from PySAGA.tools import lis_analysis
from PySAGA.tools import lis_arithmetic
from PySAGA.tools import lis_tools
# -----------------------------------------------------------------------------
# the processing pipeline in a single function; the function is called for each
# tile to be processed
# note: this would not be absolutely necessary here, but makes it easier to
# parallelize a script later on
def Task(TaskNum=0):
# -----------------------------------------------------------------------------
# load the shapefile with the tiling scheme used for processing
tileshape = saga_api.SG_Get_Data_Manager().Add_Shapes(cfg.TILESHAPE_100)
if not tileshape.is_Valid():
return lis_cmd.Error_Return('Loading tileshape ' + cfg.TILESHAPE_100 + ' failed')
# -----------------------------------------------------------------------------
# get tile to process
tile = tileshape.Get_Shape(TaskNum)
xmin = tile.Get_Extent().Get_XMin()
ymin = tile.Get_Extent().Get_YMin()
xmax = tile.Get_Extent().Get_XMax()
ymax = tile.Get_Extent().Get_YMax()
# tiles are referenced by their lower-left corner coordinate
lis_cmd.Message_Print('\nProcessing tile {:.0f}_{:.0f} \n'.format(xmin, ymin))
tilename = str(int(xmin)) + '_' + str(int(ymin))
# -----------------------------------------------------------------------------
# start processing pipeline
# -----------------------------------------------------------------------------
# query point cloud (with overlap) from virtual point cloud dataset
pc_out_list = []
if not lis_virtual.Get_Subset_from_Virtual_LASLAZ(PC_OUT=pc_out_list, FILENAME=cfg.LASVF, COPY_ATTR=True,
CONSTRAIN_QUERY=False, FILE_COMPRESSION=True, ERROR_IF_EMPTY=True,
THIN_OUT=True, SPACING=0.01, AOI_XRANGE=str(xmin)+';'+str(xmax),
AOI_YRANGE=str(ymin)+';'+str(ymax), AOI_ADD_OVERLAP=True,
OVERLAP=cfg.TILE_OVERLAP, SKIP_EMPTY_AOIS=True,
METHOD_CLIP='lower and left boundary is inside'):
return lis_cmd.Error_Return('Failed to execute tool on tile ' + tilename)
if len(pc_out_list) < 1:
return lis_cmd.Error_Return('Tile ' + tilename + ' contains no point cloud data, skipping!')
# add the point cloud datasets to the data manager for easier memory management
pc_buff = saga_api.SG_Get_Data_Manager().Add(pc_out_list[0]).asPointCloud()
if pc_buff.Get_Count() < 5:
return lis_cmd.Error_Return('Tile ' + tilename + ' has less than 5 points, skipping ...')
# -----------------------------------------------------------------------------
# filter noise points
if not lis_filtering.Remove_Isolated_Points(PC_IN=pc_buff, COPY_ATTR=True, SEARCH_RADIUS=cfg.SEARCH_RADIUS_OUTLIER,
MAX_POINTS=2, METHOD='label points', ATTR_NAME='noise', CLASSID=7,
CLASSID_UNDEF=1):
return lis_cmd.Error_Return('Failed to execute tool on tile ' + tilename)
# -----------------------------------------------------------------------------
# ground classification
if not lis_classification.Ground_Classification(PC_IN=pc_buff, FIELD_CLASS='noise', COPY_ATTR=True,
INIT_CELLSIZE=cfg.GRD_INITIAL_CELLSIZE, FINAL_CELLSIZE=cfg.GRD_FINAL_CELLSIZE,
FILTER_SEEDS=False, TERRAIN_ANGLE=88, MIN_EDGE_LENGTH=cfg.GRD_FINAL_CELLSIZE,
MAX_ANGLE=cfg.GRD_MAX_ANGLE, MAX_DISTANCE=cfg.GRD_MAX_DISTANCE,
Z_TOLERANCE=0.05, ATTACH_DZ=True, EXCLUDE_CLASSES='7'):
return lis_cmd.Error_Return('Failed to execute tool on tile ' + tilename)
# -----------------------------------------------------------------------------
# segmentation of planar objects for building roof detection
if not lis_segmentation.Segmentation_by_Plane_Growing(PC_IN=pc_buff, FILTER_FIELD='grdclass', COPY_ATTR=True, LOW_RANGE='0;1',
HIGH_RANGE='3;64', METHOD_NN='radius', RADIUS_NN=cfg.SEARCH_RADIUS,
MIN_K=3, ROBUST=True, ROBUST_THRES=0.1, ROBUST_PERCENTAGE=50.0,
METHOD_ROBUST_SEARCH_POINT='search point must be included in best fit plane',
GROWING_RADIUS=cfg.SEARCH_RADIUS / 2.0, GROWING_PLANE_DIST=0.1,
SEG_ITERATIONS=1, REGION_RADIUS=cfg.SEARCH_RADIUS, REGION_MINSIZE=20,
REGION_MAXSIZE=100000, REGION_MAXANGLE=10, REGION_MAXOFFSET=0.2,
REGION_RECALC_NORMAL=100000):
return lis_cmd.Error_Return('Failed to execute tool on tile ' + tilename)
# -----------------------------------------------------------------------------
# building and vegetation classification
if not lis_classification.Enhanced_Point_Cloud_Classification(PC_IN=pc_buff, FIELD_SEGMENTID='segmentid', FIELD_CLASSID_GROUND='grdclass',
FIELD_DZ='dz', COPY_ATTR=True, EXCLUDE_CLASSES='7', METHOD_BUILDING='point cloud features',
MIN_DZ_BUILDING=2.0, GROWING_RADIUS_BUILDING=cfg.SEARCH_RADIUS,
FEATURE_RADIUS_BUILDING=cfg.SEARCH_RADIUS, MIN_NEIGHBOR_RATIO_BUILDING=60.0,
SLOPE_ADAPTIVE_RATIO=False, MIN_POINTS_BUILDING=200, RADIUS_VEGETATION=cfg.SEARCH_RADIUS * 2.0,
MIN_POINTS_VEGETATION=5, MIN_DZ_LOW_VEGETATION=0, MIN_DZ_MEDIUM_VEGETATION=0.3,
MIN_DZ_HIGH_VEGETATION=2, OUTPUT_CALIBRATION=False):
return lis_cmd.Error_Return('Failed to execute tool on tile ' + tilename)
# -----------------------------------------------------------------------------
# clean building facades
if not lis_classification.Clean_Building_Facades(PC_IN=pc_buff, FIELD_FILTERATTR='classid', COPY_ATTR=True, ATTR_SUFFIX='facades',
RADIUS_NN=cfg.SEARCH_RADIUS, BUILDING_CLASS='6;6', LOW_ATTR_RANGE='3;5',
HIGH_ATTR_RANGE='7;64', TARGET_CLASS=1):
return lis_cmd.Error_Return('Failed to execute tool on tile ' + tilename)
# -----------------------------------------------------------------------------
# clean building roofs
if not lis_classification.Clean_Building_Roofs(PC_IN=pc_buff, FIELD_FILTERATTR='classid_facades', COPY_ATTR=True, ATTR_SUFFIX='roofs',
DIMENSION='3D search', RADIUS_NN=2.5, USE_INNER_RADIUS=False, THRES_MAJORITY=25,
TARGET_CLASS=1, FILTER_FACADE=False):
return lis_cmd.Error_Return('Failed to execute tool on tile ' + tilename)
# -----------------------------------------------------------------------------
# power line classification
if not lis_classification.Power_Line_Classification(PC_IN=pc_buff, FIELD_CLASS='classid_facades_roofs', FIELD_THRES='dz',
COPY_ATTR=True, ATTR_SUFFIX='pl', EXCLUDE_CLASSES='7', THRES_RANGE='3;1000',
SEARCH_RADIUS=cfg.SEARCH_RADIUS_PL, CLASSIFY_POLES=False):
return lis_cmd.Error_Return('Failed to execute tool on tile ' + tilename)
# -----------------------------------------------------------------------------
# clustering of power line points for filtering
if not lis_segmentation.Densitybased_Spatial_Clustering(PC_IN=pc_buff, FIELD_ATTR_1='X', FIELD_ATTR_2='Y', FIELD_ATTR_3='Z',
FIELD_CONSTRAINT='classid_pl', COPY_ATTR=True, ATTR_SUFFIX='pl',
CHOICE_DIM='3-dimensional', SEARCH_RADIUS=cfg.SEARCH_RADIUS * 2, MIN_POINTS=1,
RANGE_CONSTRAINT='14;14', TOL_CONSTRAINT='-0.5;0.5'):
return lis_cmd.Error_Return('Failed to execute tool on tile ' + tilename)
# -----------------------------------------------------------------------------
# calculate power line cluster lengths for filtering
if not lis_analysis.Segment_Features(PC_IN=pc_buff, FIELD_SEGMENTID='cluster_pl', COPY_ATTR=True, LENGTH_2D=True):
return lis_cmd.Error_Return('Failed to execute tool on tile ' + tilename)
# -----------------------------------------------------------------------------
# filter power line clusters by minimum length
if not lis_filtering.Attribute_Filter(PC_IN=pc_buff, FIELD_ATTR_1='length2d', COPY_ATTR=True, METHOD='binary flag',
TYPE_1='lower limit (min)', MIN_1=cfg.MIN_WIRE_LENGTH, ATTR_NAME='pl_filtered'):
return lis_cmd.Error_Return('Failed to execute tool on tile ' + tilename)
# -----------------------------------------------------------------------------
# relabel power line classification
if not lis_arithmetic.Point_Cloud_Attribute_Calculator(PC_IN=pc_buff,
FORMULA='ifelse(eq([classid_pl],14),ifelse(eq([pl_filtered],1),14,1),[classid_pl])',
NAME='classid_pl_filtered', TYPE='1 byte unsigned integer',
INCLUDE_NODATA=False):
return lis_cmd.Error_Return('Failed to execute tool on tile ' + tilename)
# -----------------------------------------------------------------------------
# pylon classification
if not lis_classification.Pylon_Classification(PC_IN=pc_buff, FIELD_CLASS='classid_pl_filtered', FIELD_DZ='dz',
COPY_ATTR=True, ATTR_SUFFIX='pylon', MIN_PYLON_HEIGHT=cfg.MIN_PYLON_HEIGHT,
MAX_PYLON_DIAMETER=25.0, RADIUS_NN=3.0, MAX_DIFFERENCE=1.0, MIN_POINTS=200,
SEED_CLASS=14, EXCLUDE_CLASSES='7', OUTPUT_CALIBRATION=False):
return lis_cmd.Error_Return('Failed to execute tool on tile ' + tilename)
# -----------------------------------------------------------------------------
# remove tile overlap
pc_out = saga_api.SG_Get_Data_Manager().Add_PointCloud()
if not lis_tools.Clip_Point_Cloud_with_Extent(PC_IN=pc_buff, PC_OUT=pc_out,
X_EXTENT=str(xmin)+';'+str(xmax), Y_EXTENT=str(ymin)+';'+str(ymax),
METHOD_CLIP='lower and left boundary is inside'):
return lis_cmd.Error_Return('Failed to execute tool on tile ' + tilename)
# -----------------------------------------------------------------------------
# save classified tile (without overlap)
pc_out.Save(cfg.SPCCLASSDIR + os.sep + 'tile_' + tilename + '.sg-pts-z')
# -----------------------------------------------------------------------------
# free memory resources and return success
saga_api.SG_Get_Data_Manager().Delete()
return True
# -----------------------------------------------------------------------------
# -----------------------------------------------------------------------------
# -----------------------------------------------------------------------------
# the main function
#
if __name__ == '__main__':
# -----------------------------------------------------------------------------
# start the logging
starttime = lis_cmd.Start_Logging(False)
# -----------------------------------------------------------------------------
# load the shapefile with the tiling scheme used for processing
tileshape = saga_api.SG_Get_Data_Manager().Add_Shapes(cfg.TILESHAPE_100)
if not tileshape.is_Valid():
lis_cmd.Error_Exit('Loading tileshape ' + cfg.TILESHAPE_100 + ' failed', starttime, True)
lis_cmd.Message_Print('Number of tiles to process: {:d} \n\n'.format(tileshape.Get_Count()))
# -----------------------------------------------------------------------------
# iterate over the tiles and process each tile by calling the task function
# with the processing pipeline; have a look at the "parallel tiles" version of
# this script on how to process each tile on a separate CPU core in parallel
for iTile in range(tileshape.Get_Count()):
result = Task(iTile)
if not result:
lis_cmd.Error_Exit('Processing task failed', starttime)
# -----------------------------------------------------------------------------
# processing completed, now ...
# create virtual point cloud dataset of classified tiles for subsequent processing steps (here LAS/LAZ export in script #3)
inputfiles = lis_cmd.Create_File_List(cfg.SPCCLASSDIR, '.sg-pts;.sg-pts-z', cfg.SPCCLASSDIR + os.sep + 'spcvf_input_file_list.txt')
if not lis_virtual.Create_Virtual_Point_Cloud_Dataset(INPUT_FILE_LIST=inputfiles, FILENAME=cfg.SPCVF_CLASS,
METHOD_PATHS=1, USE_HEADER=True):
lis_cmd.Error_Exit('Failed to execute tool', starttime, True)
# -----------------------------------------------------------------------------
# free memory and stop logging
saga_api.SG_Get_Data_Manager().Delete()
lis_cmd.Stop_Logging(starttime)Step 7.3: Add decorator to your function
Add the processing function wrapper decoration to your python function
Note
the ‘processing function wrapper’ decoration takes care of tool initialization and logging of the parallelized processing
@lis_cmd.proc_function_wrapper
def Task(TaskNum=0):Step 7.5: Final python function
The final python function will look like this:
import os
from PySAGA import saga_api, lis_cmd
import project_config as cfg
from PySAGA.tools import lis_virtual
from PySAGA.tools import lis_filtering
from PySAGA.tools import lis_classification
from PySAGA.tools import lis_segmentation
from PySAGA.tools import lis_analysis
from PySAGA.tools import lis_arithmetic
from PySAGA.tools import lis_tools
# -----------------------------------------------------------------------------
# the processing pipeline in a single function, will be executed per tile (each on a single CPU core)
# note: the 'processing function wrapper' decoration takes care of tool initialization and logging
@lis_cmd.proc_function_wrapper
def Task(TaskNum=0):
# -----------------------------------------------------------------------------
# load the shapefile with the tiling scheme used for processing
tileshape = saga_api.SG_Get_Data_Manager().Add_Shapes(cfg.TILESHAPE_100)
if not tileshape.is_Valid():
return lis_cmd.Error_Return('Loading tileshape ' + cfg.TILESHAPE_100 + ' failed')
# -----------------------------------------------------------------------------
# get tile to process
tile = tileshape.Get_Shape(TaskNum)
xmin = tile.Get_Extent().Get_XMin()
ymin = tile.Get_Extent().Get_YMin()
xmax = tile.Get_Extent().Get_XMax()
ymax = tile.Get_Extent().Get_YMax()
# tiles are referenced by their lower-left corner coordinate
lis_cmd.Message_Print('\nProcessing tile {:.0f}_{:.0f} \n'.format(xmin, ymin))
tilename = str(int(xmin)) + '_' + str(int(ymin))
# -----------------------------------------------------------------------------
# start processing pipeline
# -----------------------------------------------------------------------------
# query point cloud (with overlap) from virtual point cloud dataset
pc_out_list = []
if not lis_virtual.Get_Subset_from_Virtual_LASLAZ(PC_OUT=pc_out_list, FILENAME=cfg.LASVF, COPY_ATTR=True,
CONSTRAIN_QUERY=False, FILE_COMPRESSION=True, ERROR_IF_EMPTY=True,
THIN_OUT=True, SPACING=0.01, AOI_XRANGE=str(xmin)+';'+str(xmax),
AOI_YRANGE=str(ymin)+';'+str(ymax), AOI_ADD_OVERLAP=True,
OVERLAP=cfg.TILE_OVERLAP, SKIP_EMPTY_AOIS=True,
METHOD_CLIP='lower and left boundary is inside'):
return lis_cmd.Error_Return('Failed to execute tool on tile ' + tilename)
if len(pc_out_list) < 1:
return lis_cmd.Error_Return('Tile ' + tilename + ' contains no point cloud data, skipping!')
# add the point cloud datasets to the data manager for easier memory management
pc_buff = saga_api.SG_Get_Data_Manager().Add(pc_out_list[0]).asPointCloud()
if pc_buff.Get_Count() < 5:
return lis_cmd.Error_Return('Tile ' + tilename + ' has less than 5 points, skipping ...')
# -----------------------------------------------------------------------------
# filter noise points
if not lis_filtering.Remove_Isolated_Points(PC_IN=pc_buff, COPY_ATTR=True, SEARCH_RADIUS=cfg.SEARCH_RADIUS_OUTLIER,
MAX_POINTS=2, METHOD='label points', ATTR_NAME='noise', CLASSID=7,
CLASSID_UNDEF=1):
return lis_cmd.Error_Return('Failed to execute tool on tile ' + tilename)
# -----------------------------------------------------------------------------
# ground classification
if not lis_classification.Ground_Classification(PC_IN=pc_buff, FIELD_CLASS='noise', COPY_ATTR=True,
INIT_CELLSIZE=cfg.GRD_INITIAL_CELLSIZE, FINAL_CELLSIZE=cfg.GRD_FINAL_CELLSIZE,
FILTER_SEEDS=False, TERRAIN_ANGLE=88, MIN_EDGE_LENGTH=cfg.GRD_FINAL_CELLSIZE,
MAX_ANGLE=cfg.GRD_MAX_ANGLE, MAX_DISTANCE=cfg.GRD_MAX_DISTANCE,
Z_TOLERANCE=0.05, ATTACH_DZ=True, EXCLUDE_CLASSES='7'):
return lis_cmd.Error_Return('Failed to execute tool on tile ' + tilename)
# -----------------------------------------------------------------------------
# segmentation of planar objects for building roof detection
if not lis_segmentation.Segmentation_by_Plane_Growing(PC_IN=pc_buff, FILTER_FIELD='grdclass', COPY_ATTR=True, LOW_RANGE='0;1',
HIGH_RANGE='3;64', METHOD_NN='radius', RADIUS_NN=cfg.SEARCH_RADIUS,
MIN_K=3, ROBUST=True, ROBUST_THRES=0.1, ROBUST_PERCENTAGE=50.0,
METHOD_ROBUST_SEARCH_POINT='search point must be included in best fit plane',
GROWING_RADIUS=cfg.SEARCH_RADIUS / 2.0, GROWING_PLANE_DIST=0.1,
SEG_ITERATIONS=1, REGION_RADIUS=cfg.SEARCH_RADIUS, REGION_MINSIZE=20,
REGION_MAXSIZE=100000, REGION_MAXANGLE=10, REGION_MAXOFFSET=0.2,
REGION_RECALC_NORMAL=100000):
return lis_cmd.Error_Return('Failed to execute tool on tile ' + tilename)
# -----------------------------------------------------------------------------
# building and vegetation classification
if not lis_classification.Enhanced_Point_Cloud_Classification(PC_IN=pc_buff, FIELD_SEGMENTID='segmentid', FIELD_CLASSID_GROUND='grdclass',
FIELD_DZ='dz', COPY_ATTR=True, EXCLUDE_CLASSES='7', METHOD_BUILDING='point cloud features',
MIN_DZ_BUILDING=2.0, GROWING_RADIUS_BUILDING=cfg.SEARCH_RADIUS,
FEATURE_RADIUS_BUILDING=cfg.SEARCH_RADIUS, MIN_NEIGHBOR_RATIO_BUILDING=60.0,
SLOPE_ADAPTIVE_RATIO=False, MIN_POINTS_BUILDING=200, RADIUS_VEGETATION=cfg.SEARCH_RADIUS * 2.0,
MIN_POINTS_VEGETATION=5, MIN_DZ_LOW_VEGETATION=0, MIN_DZ_MEDIUM_VEGETATION=0.3,
MIN_DZ_HIGH_VEGETATION=2, OUTPUT_CALIBRATION=False):
return lis_cmd.Error_Return('Failed to execute tool on tile ' + tilename)
# -----------------------------------------------------------------------------
# clean building facades
if not lis_classification.Clean_Building_Facades(PC_IN=pc_buff, FIELD_FILTERATTR='classid', COPY_ATTR=True, ATTR_SUFFIX='facades',
RADIUS_NN=cfg.SEARCH_RADIUS, BUILDING_CLASS='6;6', LOW_ATTR_RANGE='3;5',
HIGH_ATTR_RANGE='7;64', TARGET_CLASS=1):
return lis_cmd.Error_Return('Failed to execute tool on tile ' + tilename)
# -----------------------------------------------------------------------------
# clean building roofs
if not lis_classification.Clean_Building_Roofs(PC_IN=pc_buff, FIELD_FILTERATTR='classid_facades', COPY_ATTR=True, ATTR_SUFFIX='roofs',
DIMENSION='3D search', RADIUS_NN=2.5, USE_INNER_RADIUS=False, THRES_MAJORITY=25,
TARGET_CLASS=1, FILTER_FACADE=False):
return lis_cmd.Error_Return('Failed to execute tool on tile ' + tilename)
# -----------------------------------------------------------------------------
# power line classification
if not lis_classification.Power_Line_Classification(PC_IN=pc_buff, FIELD_CLASS='classid_facades_roofs', FIELD_THRES='dz',
COPY_ATTR=True, ATTR_SUFFIX='pl', EXCLUDE_CLASSES='7', THRES_RANGE='3;1000',
SEARCH_RADIUS=cfg.SEARCH_RADIUS_PL, CLASSIFY_POLES=False):
return lis_cmd.Error_Return('Failed to execute tool on tile ' + tilename)
# -----------------------------------------------------------------------------
# clustering of power line points for filtering
if not lis_segmentation.Densitybased_Spatial_Clustering(PC_IN=pc_buff, FIELD_ATTR_1='X', FIELD_ATTR_2='Y', FIELD_ATTR_3='Z',
FIELD_CONSTRAINT='classid_pl', COPY_ATTR=True, ATTR_SUFFIX='pl',
CHOICE_DIM='3-dimensional', SEARCH_RADIUS=cfg.SEARCH_RADIUS * 2, MIN_POINTS=1,
RANGE_CONSTRAINT='14;14', TOL_CONSTRAINT='-0.5;0.5'):
return lis_cmd.Error_Return('Failed to execute tool on tile ' + tilename)
# -----------------------------------------------------------------------------
# calculate power line cluster lengths for filtering
if not lis_analysis.Segment_Features(PC_IN=pc_buff, FIELD_SEGMENTID='cluster_pl', COPY_ATTR=True, LENGTH_2D=True):
return lis_cmd.Error_Return('Failed to execute tool on tile ' + tilename)
# -----------------------------------------------------------------------------
# filter power line clusters by minimum length
if not lis_filtering.Attribute_Filter(PC_IN=pc_buff, FIELD_ATTR_1='length2d', COPY_ATTR=True, METHOD='binary flag',
TYPE_1='lower limit (min)', MIN_1=cfg.MIN_WIRE_LENGTH, ATTR_NAME='pl_filtered'):
return lis_cmd.Error_Return('Failed to execute tool on tile ' + tilename)
# -----------------------------------------------------------------------------
# relabel power line classification
if not lis_arithmetic.Point_Cloud_Attribute_Calculator(PC_IN=pc_buff,
FORMULA='ifelse(eq([classid_pl],14),ifelse(eq([pl_filtered],1),14,1),[classid_pl])',
NAME='classid_pl_filtered', TYPE='1 byte unsigned integer',
INCLUDE_NODATA=False):
return lis_cmd.Error_Return('Failed to execute tool on tile ' + tilename)
# -----------------------------------------------------------------------------
# pylon classification
if not lis_classification.Pylon_Classification(PC_IN=pc_buff, FIELD_CLASS='classid_pl_filtered', FIELD_DZ='dz',
COPY_ATTR=True, ATTR_SUFFIX='pylon', MIN_PYLON_HEIGHT=cfg.MIN_PYLON_HEIGHT,
MAX_PYLON_DIAMETER=25.0, RADIUS_NN=3.0, MAX_DIFFERENCE=1.0, MIN_POINTS=200,
SEED_CLASS=14, EXCLUDE_CLASSES='7', OUTPUT_CALIBRATION=False):
return lis_cmd.Error_Return('Failed to execute tool on tile ' + tilename)
# -----------------------------------------------------------------------------
# remove tile overlap
pc_out = saga_api.SG_Get_Data_Manager().Add_PointCloud()
if not lis_tools.Clip_Point_Cloud_with_Extent(PC_IN=pc_buff, PC_OUT=pc_out,
X_EXTENT=str(xmin)+';'+str(xmax), Y_EXTENT=str(ymin)+';'+str(ymax),
METHOD_CLIP='lower and left boundary is inside'):
return lis_cmd.Error_Return('Failed to execute tool on tile ' + tilename)
# -----------------------------------------------------------------------------
# save classified tile (without overlap)
pc_out.Save(cfg.SPCCLASSDIR + os.sep + 'tile_' + tilename + '.sg-pts-z')
# -----------------------------------------------------------------------------
# free memory resources and return success
saga_api.SG_Get_Data_Manager().Delete()
return TrueStep 8: Add parallelization code
In the Main function, replace this code:
# -----------------------------------------------------------------------------
# iterate over the tiles and process each tile by calling the task function
# with the processing pipeline; have a look at the "parallel tiles" version of
# this script on how to process each tile on a separate CPU core in parallel
for iTile in range(tileshape.Get_Count()):
result = Task(iTile)
if not result:
lis_cmd.Error_Exit('Processing task failed', starttime)with this:
num_threads = 32
if saga_api.SG_OMP_Get_Max_Num_Threads() < num_threads:
num_threads = saga_api.SG_OMP_Get_Max_Num_Threads()
lis_cmd.Run_Parallel_Shapes(Task, tiles, num_threads)
Note
This will replace the normal for-loop by a parallelized processing of the given 10 jobs.
Warning
Set the num_threads with respect to your computer memory. If e.g. the processing per tile will need 1 GB RAM, the 32 parallel processes will need 32 GB RAM in total!
The Main function will look like this:
# -----------------------------------------------------------------------------
# -----------------------------------------------------------------------------
# -----------------------------------------------------------------------------
# the main function
#
if __name__ == '__main__':
# -----------------------------------------------------------------------------
# start the logging
starttime = lis_cmd.Start_Logging(False)
# -----------------------------------------------------------------------------
# load the shapefile with the tiling scheme used for processing
tiles = saga_api.SG_Get_Data_Manager().Add_Shapes(cfg.TILESHAPE_PROCESSING)
if not tiles.is_Valid():
lis_cmd.Error_Exit('Loading tileshape ' + cfg.TILESHAPE_PROCESSING + ' failed', starttime, True)
# -----------------------------------------------------------------------------
# process the Task function in parallel (each tile on a CPU core) with the maximum number of available CPU cores
lis_cmd.Run_Parallel_Shapes(Task, tiles, saga_api.SG_OMP_Get_Max_Num_Threads())
# -----------------------------------------------------------------------------
# parallelized processing completed, now ...
# create virtual point cloud dataset of classified tiles for subsequent processing steps (here LAS/LAZ export in script #3)
inputfiles = lis_cmd.Create_File_List(cfg.SPCCLASSDIR, '.sg-pts;.sg-pts-z', cfg.SPCCLASSDIR + os.sep + 'spcvf_input_file_list.txt')
if not lis_virtual.Create_Virtual_Point_Cloud_Dataset(INPUT_FILE_LIST=inputfiles, FILENAME=cfg.SPCVF_CLASS,
METHOD_PATHS=1, USE_HEADER=True):
lis_cmd.Error_Exit('Failed to execute tool', starttime, True)
# -----------------------------------------------------------------------------
# free memory and stop logging
saga_api.SG_Get_Data_Manager().Delete()
lis_cmd.Stop_Logging(starttime)
Tip
Use the same Main Function code for all your future parallelized python scripts. Only edit the def Task() function.
Step 9: Final Code
The final code will look like this:
import os
import project_config as cfg
from PySAGA import saga_api, lis_cmd
from PySAGA.tools import lis_virtual
from PySAGA.tools import lis_filtering
from PySAGA.tools import lis_classification
from PySAGA.tools import lis_segmentation
from PySAGA.tools import lis_analysis
from PySAGA.tools import lis_arithmetic
from PySAGA.tools import lis_tools
@lis_cmd.proc_function_wrapper
def Task(TaskNum=0):
# -----------------------------------------------------------------------------
# load the shapefile with the tiling scheme used for processing
tileshape = saga_api.SG_Get_Data_Manager().Add_Shapes(cfg.TILESHAPE_100)
if not tileshape.is_Valid():
return lis_cmd.Error_Return('Loading tileshape ' + cfg.TILESHAPE_100 + ' failed')
# -----------------------------------------------------------------------------
# get tile to process
tile = tileshape.Get_Shape(TaskNum)
xmin = tile.Get_Extent().Get_XMin()
ymin = tile.Get_Extent().Get_YMin()
xmax = tile.Get_Extent().Get_XMax()
ymax = tile.Get_Extent().Get_YMax()
# tiles are referenced by their lower-left corner coordinate
lis_cmd.Message_Print('\nProcessing tile {:.0f}_{:.0f} \n'.format(xmin, ymin))
tilename = str(int(xmin)) + '_' + str(int(ymin))
# -----------------------------------------------------------------------------
# start processing pipeline
# -----------------------------------------------------------------------------
# query point cloud (with overlap) from virtual point cloud dataset
pc_out_list = []
if not lis_virtual.Get_Subset_from_Virtual_LASLAZ(PC_OUT=pc_out_list, FILENAME=cfg.LASVF, COPY_ATTR=True,
CONSTRAIN_QUERY=False, FILE_COMPRESSION=True, ERROR_IF_EMPTY=True,
THIN_OUT=True, SPACING=0.01, AOI_XRANGE=str(xmin)+';'+str(xmax),
AOI_YRANGE=str(ymin)+';'+str(ymax), AOI_ADD_OVERLAP=True,
OVERLAP=cfg.TILE_OVERLAP, SKIP_EMPTY_AOIS=True,
METHOD_CLIP='lower and left boundary is inside'):
return lis_cmd.Error_Return('Failed to execute tool on tile ' + tilename)
if len(pc_out_list) < 1:
return lis_cmd.Error_Return('Tile ' + tilename + ' contains no point cloud data, skipping!')
# add the point cloud datasets to the data manager for easier memory management
pc_buff = saga_api.SG_Get_Data_Manager().Add(pc_out_list[0]).asPointCloud()
if pc_buff.Get_Count() < 5:
return lis_cmd.Error_Return('Tile ' + tilename + ' has less than 5 points, skipping ...')
# -----------------------------------------------------------------------------
# filter noise points
if not lis_filtering.Remove_Isolated_Points(PC_IN=pc_buff, COPY_ATTR=True, SEARCH_RADIUS=cfg.SEARCH_RADIUS_OUTLIER,
MAX_POINTS=2, METHOD='label points', ATTR_NAME='noise', CLASSID=7,
CLASSID_UNDEF=1):
return lis_cmd.Error_Return('Failed to execute tool on tile ' + tilename)
# -----------------------------------------------------------------------------
# ground classification
if not lis_classification.Ground_Classification(PC_IN=pc_buff, FIELD_CLASS='noise', COPY_ATTR=True,
INIT_CELLSIZE=cfg.GRD_INITIAL_CELLSIZE, FINAL_CELLSIZE=cfg.GRD_FINAL_CELLSIZE,
FILTER_SEEDS=False, TERRAIN_ANGLE=88, MIN_EDGE_LENGTH=cfg.GRD_FINAL_CELLSIZE,
MAX_ANGLE=cfg.GRD_MAX_ANGLE, MAX_DISTANCE=cfg.GRD_MAX_DISTANCE,
Z_TOLERANCE=0.05, ATTACH_DZ=True, EXCLUDE_CLASSES='7'):
return lis_cmd.Error_Return('Failed to execute tool on tile ' + tilename)
# -----------------------------------------------------------------------------
# segmentation of planar objects for building roof detection
if not lis_segmentation.Segmentation_by_Plane_Growing(PC_IN=pc_buff, FILTER_FIELD='grdclass', COPY_ATTR=True, LOW_RANGE='0;1',
HIGH_RANGE='3;64', METHOD_NN='radius', RADIUS_NN=cfg.SEARCH_RADIUS,
MIN_K=3, ROBUST=True, ROBUST_THRES=0.1, ROBUST_PERCENTAGE=50.0,
METHOD_ROBUST_SEARCH_POINT='search point must be included in best fit plane',
GROWING_RADIUS=cfg.SEARCH_RADIUS / 2.0, GROWING_PLANE_DIST=0.1,
SEG_ITERATIONS=1, REGION_RADIUS=cfg.SEARCH_RADIUS, REGION_MINSIZE=20,
REGION_MAXSIZE=100000, REGION_MAXANGLE=10, REGION_MAXOFFSET=0.2,
REGION_RECALC_NORMAL=100000):
return lis_cmd.Error_Return('Failed to execute tool on tile ' + tilename)
# -----------------------------------------------------------------------------
# building and vegetation classification
if not lis_classification.Enhanced_Point_Cloud_Classification(PC_IN=pc_buff, FIELD_SEGMENTID='segmentid', FIELD_CLASSID_GROUND='grdclass',
FIELD_DZ='dz', COPY_ATTR=True, EXCLUDE_CLASSES='7', METHOD_BUILDING='point cloud features',
MIN_DZ_BUILDING=2.0, GROWING_RADIUS_BUILDING=cfg.SEARCH_RADIUS,
FEATURE_RADIUS_BUILDING=cfg.SEARCH_RADIUS, MIN_NEIGHBOR_RATIO_BUILDING=60.0,
SLOPE_ADAPTIVE_RATIO=False, MIN_POINTS_BUILDING=200, RADIUS_VEGETATION=cfg.SEARCH_RADIUS * 2.0,
MIN_POINTS_VEGETATION=5, MIN_DZ_LOW_VEGETATION=0, MIN_DZ_MEDIUM_VEGETATION=0.3,
MIN_DZ_HIGH_VEGETATION=2, OUTPUT_CALIBRATION=False):
return lis_cmd.Error_Return('Failed to execute tool on tile ' + tilename)
# -----------------------------------------------------------------------------
# clean building facades
if not lis_classification.Clean_Building_Facades(PC_IN=pc_buff, FIELD_FILTERATTR='classid', COPY_ATTR=True, ATTR_SUFFIX='facades',
RADIUS_NN=cfg.SEARCH_RADIUS, BUILDING_CLASS='6;6', LOW_ATTR_RANGE='3;5',
HIGH_ATTR_RANGE='7;64', TARGET_CLASS=1):
return lis_cmd.Error_Return('Failed to execute tool on tile ' + tilename)
# -----------------------------------------------------------------------------
# clean building roofs
if not lis_classification.Clean_Building_Roofs(PC_IN=pc_buff, FIELD_FILTERATTR='classid_facades', COPY_ATTR=True, ATTR_SUFFIX='roofs',
DIMENSION='3D search', RADIUS_NN=2.5, USE_INNER_RADIUS=False, THRES_MAJORITY=25,
TARGET_CLASS=1, FILTER_FACADE=False):
return lis_cmd.Error_Return('Failed to execute tool on tile ' + tilename)
# -----------------------------------------------------------------------------
# power line classification
if not lis_classification.Power_Line_Classification(PC_IN=pc_buff, FIELD_CLASS='classid_facades_roofs', FIELD_THRES='dz',
COPY_ATTR=True, ATTR_SUFFIX='pl', EXCLUDE_CLASSES='7', THRES_RANGE='3;1000',
SEARCH_RADIUS=cfg.SEARCH_RADIUS_PL, CLASSIFY_POLES=False):
return lis_cmd.Error_Return('Failed to execute tool on tile ' + tilename)
# -----------------------------------------------------------------------------
# clustering of power line points for filtering
if not lis_segmentation.Densitybased_Spatial_Clustering(PC_IN=pc_buff, FIELD_ATTR_1='X', FIELD_ATTR_2='Y', FIELD_ATTR_3='Z',
FIELD_CONSTRAINT='classid_pl', COPY_ATTR=True, ATTR_SUFFIX='pl',
CHOICE_DIM='3-dimensional', SEARCH_RADIUS=cfg.SEARCH_RADIUS * 2, MIN_POINTS=1,
RANGE_CONSTRAINT='14;14', TOL_CONSTRAINT='-0.5;0.5'):
return lis_cmd.Error_Return('Failed to execute tool on tile ' + tilename)
# -----------------------------------------------------------------------------
# calculate power line cluster lengths for filtering
if not lis_analysis.Segment_Features(PC_IN=pc_buff, FIELD_SEGMENTID='cluster_pl', COPY_ATTR=True, LENGTH_2D=True):
return lis_cmd.Error_Return('Failed to execute tool on tile ' + tilename)
# -----------------------------------------------------------------------------
# filter power line clusters by minimum length
if not lis_filtering.Attribute_Filter(PC_IN=pc_buff, FIELD_ATTR_1='length2d', COPY_ATTR=True, METHOD='binary flag',
TYPE_1='lower limit (min)', MIN_1=cfg.MIN_WIRE_LENGTH, ATTR_NAME='pl_filtered'):
return lis_cmd.Error_Return('Failed to execute tool on tile ' + tilename)
# -----------------------------------------------------------------------------
# relabel power line classification
if not lis_arithmetic.Point_Cloud_Attribute_Calculator(PC_IN=pc_buff,
FORMULA='ifelse(eq([classid_pl],14),ifelse(eq([pl_filtered],1),14,1),[classid_pl])',
NAME='classid_pl_filtered', TYPE='1 byte unsigned integer',
INCLUDE_NODATA=False):
return lis_cmd.Error_Return('Failed to execute tool on tile ' + tilename)
# -----------------------------------------------------------------------------
# pylon classification
if not lis_classification.Pylon_Classification(PC_IN=pc_buff, FIELD_CLASS='classid_pl_filtered', FIELD_DZ='dz',
COPY_ATTR=True, ATTR_SUFFIX='pylon', MIN_PYLON_HEIGHT=cfg.MIN_PYLON_HEIGHT,
MAX_PYLON_DIAMETER=25.0, RADIUS_NN=3.0, MAX_DIFFERENCE=1.0, MIN_POINTS=200,
SEED_CLASS=14, EXCLUDE_CLASSES='7', OUTPUT_CALIBRATION=False):
return lis_cmd.Error_Return('Failed to execute tool on tile ' + tilename)
# -----------------------------------------------------------------------------
# remove tile overlap
pc_out = saga_api.SG_Get_Data_Manager().Add_PointCloud()
if not lis_tools.Clip_Point_Cloud_with_Extent(PC_IN=pc_buff, PC_OUT=pc_out,
X_EXTENT=str(xmin)+';'+str(xmax), Y_EXTENT=str(ymin)+';'+str(ymax),
METHOD_CLIP='lower and left boundary is inside'):
return lis_cmd.Error_Return('Failed to execute tool on tile ' + tilename)
# -----------------------------------------------------------------------------
# save classified tile (without overlap)
pc_out.Save(cfg.SPCCLASSDIR + os.sep + 'tile_' + tilename + '.sg-pts-z')
# -----------------------------------------------------------------------------
# free memory resources and return success
saga_api.SG_Get_Data_Manager().Delete()
return True
# -----------------------------------------------------------------------------
# -----------------------------------------------------------------------------
# -----------------------------------------------------------------------------
# the main function
#
if __name__ == '__main__':
# -----------------------------------------------------------------------------
# start the logging
starttime = lis_cmd.Start_Logging(False)
# -----------------------------------------------------------------------------
# load the shapefile with the tiling scheme used for processing
tiles = saga_api.SG_Get_Data_Manager().Add_Shapes(cfg.TILESHAPE_100)
if not tiles.is_Valid():
lis_cmd.Error_Exit('Loading tileshape ' + cfg.TILESHAPE_100 + ' failed', starttime, True)
# -----------------------------------------------------------------------------
# process the Task function in parallel (each tile on a CPU core) with the maximum number of available CPU cores
lis_cmd.Run_Parallel_Shapes(Task, tiles, saga_api.SG_OMP_Get_Max_Num_Threads())
# -----------------------------------------------------------------------------
# parallelized processing completed, now ...
# create virtual point cloud dataset of classified tiles for subsequent processing steps (here LAS/LAZ export in script #3)
inputfiles = lis_cmd.Create_File_List(cfg.SPCCLASSDIR, '.sg-pts;.sg-pts-z', cfg.SPCCLASSDIR + os.sep + 'spcvf_input_file_list.txt')
if not lis_virtual.Create_Virtual_Point_Cloud_Dataset(INPUT_FILE_LIST=inputfiles, FILENAME=cfg.SPCVF_CLASS,
METHOD_PATHS=1, USE_HEADER=True):
lis_cmd.Error_Exit('Failed to execute tool', starttime, True)
# -----------------------------------------------------------------------------
# free memory and stop logging
saga_api.SG_Get_Data_Manager().Delete()
lis_cmd.Stop_Logging(starttime)Step 11: Run the script using Runscript
Run the tool in the Miniforge Prompt using the command python run_script.py 02_data_processing_parallel_tiles.py.
If you look into the final logfiles and compare the total processing time of both 02_data_processing.py and 02_data_processing_parallel_tiles.py, you will recognize, that the parallelized version of the same processing pipeline was much faster. This becomes more and more pronounced, the more tiles you have to process.
Conclusion
We are able to design processing pipelines with LIS Pro 3D and implement a python script for its automation. Furthermore we are able to parallelize the whole processing job using python.
If you have finalized all tutorials, you will be able to read, interprete and modify scripts implemented by others.
Note
If you get the whole project folder including, data and scripts folders from a colleague, you will be able to run it by only editing the filepaths defined in project_config.py
Additional Script - Data Delivery
Additionally to this tutorial, you could e.g. add a specified script for data delivery, merging the tiled data into one single block again using a respective tileshape and exporting it into a *.laz file again.
Create a new python script 03_data_delivery.py and add the following code:
import os
import project_config as cfg
from PySAGA import saga_api, lis_cmd
from PySAGA.tools import lis_virtual
from PySAGA.tools import lis_io
# -----------------------------------------------------------------------------
# the processing pipeline in a single function
# note: this would not be absolutely necessary here, but makes it easier to
# parallelize a script later on
def Task(TaskNum=0):
# -----------------------------------------------------------------------------
# load the shapefile with the tiling scheme used for data delivery
tileshape = saga_api.SG_Get_Data_Manager().Add_Shapes(cfg.TILESHAPE_500)
if not tileshape.is_Valid():
lis_cmd.Error_Exit('Loading tileshape ' + cfg.TILESHAPE_500 + ' failed', starttime, True)
# -----------------------------------------------------------------------------
# get tile to process
tile = tileshape.Get_Shape(TaskNum)
xmin = tile.Get_Extent().Get_XMin()
ymin = tile.Get_Extent().Get_YMin()
xmax = tile.Get_Extent().Get_XMax()
ymax = tile.Get_Extent().Get_YMax()
# tiles are referenced by their lower-left corner coordinate
lis_cmd.Message_Print('\nProcessing tile {:.0f}_{:.0f} \n'.format(xmin, ymin))
tilename = str(int(xmin)) + '_' + str(int(ymin))
# --------------------------------------------------------------
# query point cloud (without overlap) from virtual point cloud dataset
pc_out_list = []
if not lis_virtual.Get_Subset_from_Virtual_Point_Cloud(PC_OUT=pc_out_list, FILENAME=cfg.SPCVF_CLASS, COPY_ATTR=True,
CONSTRAIN_QUERY=False, ERROR_IF_EMPTY=True, AOI_XRANGE=str(xmin)+';'+str(xmax),
AOI_YRANGE=str(ymin)+';'+str(ymax), AOI_ADD_OVERLAP=False,
METHOD_CLIP='lower and left boundary is inside'):
return lis_cmd.Error_Return('Failed to execute tool on tile ' + tilename)
if len(pc_out_list) < 1:
return lis_cmd.Error_Return('Tile ' + tilename + ' contains no point cloud data, skipping!')
# add the point cloud datasets to the data manager for easier memory management
pc = saga_api.SG_Get_Data_Manager().Add(pc_out_list[0]).asPointCloud()
# -----------------------------------------------------------------------------
# export point cloud as LAZ
if not lis_io.Export_LASLAZ_File(POINTS=pc,
T='gps-time',
r='number of the return',
n='number of returns of given pulse',
i='intensity',
c='classid_pylon',
a='scan angle',
d='direction of scan flag',
e='edge of flight line flag',
u='user data',
p='point source ID',
FILE=cfg.EXPORTDIR + os.sep + 'tile_' + tilename + '.laz',
FILE_FORMAT='LAS 1.4', FORMAT='0', RGB_RANGE='8 bit',
SCALE_X=0.001, SCALE_Y=0.001, SCALE_Z=0.001):
return lis_cmd.Error_Return('Failed to execute tool on tile ' + tilename)
# -----------------------------------------------------------------------------
# free memory and return success
saga_api.SG_Get_Data_Manager().Delete()
return True
# -----------------------------------------------------------------------------
# -----------------------------------------------------------------------------
# -----------------------------------------------------------------------------
# the main function
#
if __name__ == '__main__':
# -----------------------------------------------------------------------------
# start the logging
starttime = lis_cmd.Start_Logging(False)
# -----------------------------------------------------------------------------
# load the shapefile with the tiling scheme used for data delivery
tileshape = saga_api.SG_Get_Data_Manager().Add_Shapes(cfg.TILESHAPE_500)
if not tileshape.is_Valid():
lis_cmd.Error_Exit('Loading tileshape ' + cfg.TILESHAPE_500 + ' failed', starttime, True)
lis_cmd.Message_Print('Number of tiles to process: {:d} \n\n'.format(tileshape.Get_Count()))
# -----------------------------------------------------------------------------
# iterate over the tiles and process each tile by calling the task function
# with the processing pipeline
for iTile in range(tileshape.Get_Count()):
result = Task(iTile)
if not result:
lis_cmd.Error_Exit('Processing task failed', starttime)
# -----------------------------------------------------------------------------
# free memory and stop logging
saga_api.SG_Get_Data_Manager().Delete()
lis_cmd.Stop_Logging(starttime)