Automate Data Preparation
< Previous section Next section >
ImportantPrerequisite
If you have not yet setup LIS Pro 3D for usage with Python (specifically the PySAGA api), please have a look at our Automation Guide and specifically the section Installation of LIS Pro 3D with Python. There, you will find a step-by-step description on how to setup your development environment!
TipRecommendation: Complete the Automation Guide First!
It is recommended to complete the Automation Guide before proceeding with this section, in order to have more background knowledge about what we are doing here.
Folder structure
For this tutorial, we assume the following folder structure:
├── data
| ├── lidar
| | ├── pointclouds
| | | ├── 2684000_1246000.laz
| | | └── ...
| | └── DTMs
| | ├── 2684000_1246000.tif
| | └── ...
| └── vector
| └── Amtliche_Vermessung_-_Datenmodell_ZH_-Standard-_-OGD
| | ├── Bo_BoFlaeche_A.shp
| | └── ...
└── scripts
├── run_script.py
├── project_config.py
├── 01_data_preparation.py
└── 02_data_processing.pyScripts
project_config.py01_data_preparation.py
project_config.py
import os
import sys
# set path to PySAGA and LIS Pro 3D tools
os.environ["SAGA_PATH"] = "C:\Program Files\Laserdata\LIS Pro 3D"
sys.path.insert(0, os.environ["SAGA_PATH"])
# -----------------------------------------------------------------------------
# processing directories
BASE_DIR = r"C:\...\city-modeller"
## these folders should exist
LAZ_DIR = os.path.join(BASE_DIR, "data", "lidar", "pointclouds")
DTM_DIR = os.path.join(BASE_DIR, "data", "lidar", "DTMs")
VECT_DIR = os.path.join(BASE_DIR, "data", "vector")
INPUT_VECT_DIR = os.path.join(
VECT_DIR, "Amtliche_Vermessung_-_Datenmodell_ZH_-Standard-_-OGD"
)
## these folders have to be created
PROCESS_DIR = os.path.join(BASE_DIR, "processing")
SHAPES_DIR = os.path.join(PROCESS_DIR, "shapes")
GML_DIR = os.path.join(PROCESS_DIR, "gml")
JSON_DIR = os.path.join(PROCESS_DIR, "json")
# -----------------------------------------------------------------------------
# -----------------------------------------------------------------------------
# create (non-existing) dirs
os.makedirs(PROCESS_DIR, exist_ok=True)
os.makedirs(SHAPES_DIR, exist_ok=True)
os.makedirs(GML_DIR, exist_ok=True)
os.makedirs(JSON_DIR, exist_ok=True)
# -----------------------------------------------------------------------------
# -----------------------------------------------------------------------------
# global settings: files, tool parameters, ...
INPUT_POLYGON_LAYER = os.path.join(INPUT_VECT_DIR, "Bo_BoFlaeche_A.shp")
BUILDING_FOOTPRINTS_RAW = os.path.join(VECT_DIR, "Building_Footprints.shp")
BUILDING_FOOTPRINTS_USED = os.path.join(
VECT_DIR, "Building_Footprints_Used.shp"
)
AOI = os.path.join(VECT_DIR, "AOI.shp")
TILESHAPE_250 = os.path.join(VECT_DIR, "Tiles_250.shp")
LASVF = os.path.join(LAZ_DIR, "tiles.lasvf")
VRT = os.path.join(DTM_DIR, "tiles.vrt")
TILE_SIZE_PROCESSING = 250.0
PROJ = "+proj=somerc +lat_0=46.95240555555556 +lon_0=7.439583333333333 +k_0=1 +x_0=2600000 +y_0=1200000 +ellps=bessel +towgs84=674.374,15.056,405.346,0,0,0,0 +units=m +no_defs"
# parameters needed for the processing script
OVERLAP = 5
## params clustering
RADIUS_CLUSTER = 2.0
MIN_POINTS_CLUSTER = 15
MIN_CLUSTERSIZE = 200
## params segmentation
TILESIZE_SEGMENTATION = 10
OVERLAP_SEGMENTATION = 5
SEGMENT_MINSIZE = 20
ROBUST_DIST_THRES = 0.150000
RANSAC_ITER = 50
RANSAC_SEED = 2
REGION_RADIUS = 0.750000
REGION_MINSIZE = 20
REGION_MAXANGLE = 10.000000
REGION_MAXOFFSET = 0.500000
## params filtering
FILTER_MIN_SEGSIZE = 50
FILTER_MAX_SLOPE = 85
## params assign roofsegment to buildings
PERCENTAGE_INSIDE = 10
## params building modeller
MODEL_MAX_SLOPE = 85
MIN_POINTS_ROOF = 20
MAX_EDGELENGTH = 1.500000
MAX_LINEDIST = 0.300000
SNAP_DIST = 0.750000
ITERATIONS_ROOF_GAPS = 5
MIN_BUILDING_HEIGHT = 0.100000
MAX_AREA_RECTIFY = 2
## params export
EPSG = 2056
EXTERNALID = "OBJID"Change the BASE_DIR to the working directory you are using on your own machine!
01_data_preparation.py
- Create LAS/LAZ index
- Create a point cloud catalogue (lasvf) for the available point clouds
- Create a raster catalogue (VRT) for the available DTMs
- Create a shape of the extent of the entire catalogue
- Set projection (Spatial Reference) for the shape
- Create an AOI from the overview
- Select and save relevant building footprints for the AOI
- Create a tiling scheme for processing with 250m tilesize
import os
import project_config as cfg
from PySAGA import saga_api, lis_cmd
from PySAGA.tools import *
def prepare_data():
# get list of point cloud files (laz-files)
laz_list_file = os.path.join(cfg.LAZ_DIR, "lasvf_input_file_list.txt")
laz_files = lis_cmd.Create_File_List(
cfg.LAZ_DIR, ".las;.laz", laz_list_file
)
if not lis_io.Create_LASLAZ_Index(
FILES="",
INPUT_FILE_LIST=laz_files,
INDEX_ALL=True,
TILE_SIZE=0.000000,
THRESHOLD=1000,
MIN_POINTS=100000,
MAX_INTERVALS=-20,
):
return lis_cmd.Error_Return("Failed to execute tool")
# create virtual point cloud for seamless access
if not lis_virtual.Create_Virtual_LASLAZ_Dataset(
FILES="",
INPUT_FILE_LIST=laz_files,
FILENAME=cfg.LASVF,
METHOD_PATHS="relative",
COLOR_DEPTH="16 bit",
):
return lis_cmd.Error_Return("Failed to execute tool")
# get list of DTM files (tif-files)
dem_list_file = os.path.join(cfg.DTM_DIR, "vrt_input_file_list.txt")
dtm_files = lis_cmd.Create_File_List(
cfg.DTM_DIR, ".tif;.TIFF", dem_list_file
)
# create virtual raster table for seamless access
if not io_gdal.Create_Virtual_Raster_VRT(
FILES="",
FILE_LIST=dtm_files,
VRT_NAME=cfg.VRT,
RESAMPLING="nearest",
RESOLUTION="highest",
):
return lis_cmd.Error_Return("Failed to execute tool")
# load polygons and select building footprints
input_polygons = saga_api.SG_Get_Data_Manager().Add_Shapes(
cfg.INPUT_POLYGON_LAYER
)
if not input_polygons:
return lis_cmd.Error_Return("Failed to load input polygons")
buildings_raw = (
saga_api.SG_Get_Data_Manager().Add_Shapes()
) # output data object
if not shapes_tools.Select_by_Attributes_Numerical_Expression(
SHAPES=input_polygons,
COPY=buildings_raw,
FIELD="ART",
EXPRESSION="a < 8",
INVERSE=False,
USE_NODATA=False,
METHOD="new selection",
POSTJOB="copy",
):
return lis_cmd.Error_Return("Failed to execute tool")
# add building enumeration
if not table_tools.Field_Enumeration(
INPUT=buildings_raw,
FIELD="<not set>",
ENUM="<not set>",
NAME="proc_id", # enumeration field name
ORDER="ascending",
):
return lis_cmd.Error_Return("Failed to execute tool")
# generate polygon shapes of physical tiling of DTM files
catalogue = saga_api.SG_Get_Data_Manager().Add_Shapes()
if not io_gdal.Create_Raster_Catalogue_from_Virtual_Raster_VRT(
CATALOGUE=catalogue, VRT_FILE=cfg.VRT
):
return lis_cmd.Error_Return("Failed to execute tool")
catalogue.Get_Projection().Create(saga_api.CSG_String(cfg.PROJ))
# get the total extent of all DTM files
total_extent = saga_api.SG_Get_Data_Manager().Add_Shapes()
if not shapes_tools.Get_Shapes_Extents(
SHAPES=catalogue, EXTENTS=total_extent, OUTPUT="all shapes"
):
return lis_cmd.Error_Return("Failed to execute tool")
total_extent.Get_Projection().Create(saga_api.CSG_String(cfg.PROJ))
total_extent.Save(cfg.AOI)
# select all buildings within the total extent
if not shapes_tools.Select_by_Location(
SHAPES=buildings_raw,
LOCATIONS=total_extent,
CONDITION="are completely within",
METHOD="new selection",
):
return lis_cmd.Error_Return("Failed to execute tool")
# copy selected buildings to new layer
buildings_used = saga_api.SG_Get_Data_Manager().Add_Shapes()
if not shapes_tools.Copy_Selection_to_New_Shapes_Layer(
INPUT=buildings_raw, OUTPUT=buildings_used
):
return lis_cmd.Error_Return("Failed to execute tool")
buildings_used.Get_Projection().Create(saga_api.CSG_String(cfg.PROJ))
buildings_used.Save(cfg.BUILDING_FOOTPRINTS_USED)
# create spatial chunks for processing
processing_units = saga_api.SG_Get_Data_Manager().Add_Shapes()
if not lis_tools.Create_Spatial_Processing_Units(
SHP_EXTENT=total_extent,
SHP_UNITS=processing_units,
TILE_SIZE=250.000000,
):
return lis_cmd.Error_Return("Failed to execute tool")
processing_units.Get_Projection().Create(saga_api.CSG_String(cfg.PROJ))
processing_units.Save(cfg.TILESHAPE_250)
# free memory ressources
saga_api.SG_Get_Data_Manager().Delete()
return True
if __name__ == "__main__":
# -----------------------------------------------------------------------------
# start the logging
starttime = lis_cmd.Start_Logging(False)
# -----------------------------------------------------------------------------
# call the processing function
result = prepare_data()
if not result:
lis_cmd.Error_Exit("Data preparation failed.", starttime)
# -----------------------------------------------------------------------------
# stop logging
lis_cmd.Stop_Logging(starttime)