Study case 2: Monastery St. Anna & Scharfeneck Ruin, AT¶
About the dataset¶
This LiDAR dataset covers an area in the nature park Mannersdorf-Wüste in Eastern Austria and includes the monastery St. Anna and the Scharfeneck Ruin. It has a point density of approximately 9 pts/m$^2$.
1. Import afwizard, set parameters and load your dataset¶
First, we import the afwizard library:
import afwizard as af
Now, we have to set the LAStools directory to the path where we installed LAStools. This may be a different location on your computer, so you will have to provide the correct path:
af.set_lastools_directory("C:/LAStools")
The next steps are to create the directory "filters", where we will later save our developed filter pipelines. After that, we set the directory where the input data is located and where the output data will be written. In case, the directory is not existing, it will be created.
af.set_data_directory("filters", create_dir=True)
af.set_data_directory("data", create_dir=True)
Now we can load our first dataset by file name and additionally provide the spatial reference system. Make sure that the data are copied to your data directory specified just before.
ds = af.DataSet(filename="StA_last.laz", spatial_reference="EPSG:31256")
We remove any existing classification from the dataset.
ds = af.remove_classification(ds)
2. Restricting datasets¶
In the case of this dataset, there is only one segment. Still, we have to load a GeoJSON file containing a pre-prepared segmentation (in this case with just one segment) of the LiDAR tile. In case you want to produce your own segmentation, you can use e.g. QGIS: create a new shapefile layer, draw the segments (the shapes must not overlap and there must not be areas in your dataset not covered by a segment) and export the layer as GeoJSON.
segmentation = af.load_segmentation(
"StA_segment.geojson", spatial_reference="EPSG:31256"
)
Since the segments in the segmentation may cover large areas and contain many points, we can sample our segments using the restrict method of afwizard. This will make the process of creating filter pipelines much faster!
The aim is to select characteristic subareas running the following command. The command will open a GUI showing the segment on a map. If you want to see a draft hillshade version of the DTM, click on "Visualize".
Click the small polygon on the sidebar on the left to draw a segmentation polygon for an area of interest. When you are done, finish by clicking "Finalize". Your restricted dataset is then saved to the variable rds. Save the data by executing the next code cell "(saved = rds.save(...))". If you do not draw a polygon, the entire tile will be used in the next steps.
rds = ds.restrict(segmentation_overlay=segmentation)
saved = rds.save("data/sample_StA.las", overwrite=True)
3. Creating filter pipelines¶
Now we can create a filter pipelines for the restricted dataset we created before. This is done using the interactive pipeline_tuning function. Here, you can select from different outlier filtering and ground filtering algorithms and combine them in a sequential filter pipeline. You can preview the result of different versions of your pipelines and delete filtering steps if it turns out they are not needed or need different settings.
sampledataset1 = af.DataSet(filename="data/sample_StA.las", spatial_reference="EPSG:31256")
The following command will open the GUI for creating a filter pipeline. In the centre, the point cloud of the spatial sample is visualized as hillshade, slope or a combination of both. Visualization settings can be tuned using the fields in the right column.
In the center, a tab (#0) with an image is displayed showing the unclassified data as hillshade.
In the left column, a filter can be added to the pipeline. Click on "Pipeline stages" and then on "+ Add entry" and select the filter algorithm of your choice (here: lasground_new (LASTools)). Finally, you can set the parameters (offset, spike, step). Here, we let them at their defaut values. (0, 0.5, 1). Clicking "Preview" button on top of the right column will trigger the filtering of the sample using the current filter and parameter strings and display the result as a new window in tab "#1" in the centre column (with the preset hillshade at 45° altitude and 0.2 m resolution).
pipeline1 = af.pipeline_tuning(sampledataset1)
af.set_data_directory("filters", create_dir=True)
af.save_filter(pipeline1, "filters/StA_woodland.json")
4. Mapping segmentations to Filter Pipelines¶
Now we map the created filter pipelines to the segmentations.
af.add_filter_library(path="filters", recursive=False)
Using the selection window below, select the pipeline you created before. Finish the interactive step by clicking "Finalize".
pipelines = af.select_pipeline_from_library("filters")
In this stage, you select which pipeline to assign to which segment. As we have only one segment, on the top right, you leave the dropdown menu to "ID". In the dropdown menu below the segment, select the corresponding filter you created. When you are done with assigning, click "Finalize". Then run the next code cell to save the segments together with their assigned filter pipeline.
assigned_segmentation = af.assign_pipeline(
ds, segmentation=segmentation, pipelines=pipelines
)
af.set_data_directory("data", create_dir=True)
assigned_segmentation.save("data/StA_segments_assigned.geojson")
5. Adaptive filtering by running AFwizard¶
Finally, the adaptive filtering has to be calculated. This is done on a command line basis. Therefore, open a new Tab using the "+" symbol next to the top of this window and cklick on "Terminal".
ATTENTION: the directories refer relatively to the directory of the prompt in the terminal!
In the example below, the prompt is executed from .../trail-groundfiltering.
Filesystem is:
.../trail-groundfiltering.../trail-groundfiltering/data.../trail-groundfiltering/filters
Note:
- file
StA_last.lazis in directorydata - file
StA_segments_assigned.geojsonis in directorydata - Directory
outputwill be created indata/output - filters are stored in
filters
afwizard --dataset=data/StA_last.laz --dataset-crs=EPSG:31256 --segmentation=data/StA_segment_assigned.geojson --segmentation-crs=EPSG:31256 --output-dir=data/output --library filters --lastools-dir="C:\LAStools"
!afwizard --dataset=data/StA_last.laz --dataset-crs=EPSG:31256 --segmentation=data/StA_segments_assigned.geojson --segmentation-crs=EPSG:31256 --output-dir=data/output --library filters --lastools-dir="C:/LAStools"