Study case 1: Punta Križa, HR¶
About the dataset¶
This LiDAR dataset covers an area at the south coast of Cres Island in Croatia, to the west of Punta Križa. 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="PK_last.laz", spatial_reference="EPSG:25833")
We remove any existing classification from the dataset.
ds = af.remove_classification(ds)
2. Restricting datasets¶
Since the dataset covers both land and water, the goal is to split our dataset in two segments to process these areas separately. To help us, we load a GeoJSON file containing a pre-prepared segmentation of the LiDAR tile. In case you want to produce your own segmentation, you can use e.g. QGIS: create a shapefile, 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 Shapefile as GeoJSON.
segmentation = af.load_segmentation(
"PK_segments.geojson", spatial_reference="EPSG:25833"
)
Since the segments in the segmentation may cover large areas and contain many points, we 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, one for each segment running following command. The command will open a GUI showing the segments 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 the land area. 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(...)). Then repeat the restrict and the save step for the water area as well. For the water area, the filename is changed to "data/sample_uw.las" when saving (see below).
rds = ds.restrict(segmentation_overlay=segmentation)
saved = rds.save("data/sample_stonewall.las", overwrite=True)
The same procedure is now repeated for the submerged terrain:
rds = ds.restrict(segmentation_overlay=segmentation)
saved = rds.save("data/sample_uw.las", overwrite=True)
3. Creating filter pipelines¶
Now we can create targeted filter pipelines for the two datasets 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. We do this step twice, for each of our two segments.
sampledataset1 = af.DataSet(filename="data/sample_stonewall.las", spatial_reference="EPSG:25833")
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). Before finalizing a finished filtering pipeline, fill in the "Pipeline metadata" in the dropdown on the left.
pipeline1 = af.pipeline_tuning(sampledataset1)
af.save_filter(pipeline1, "filters/PK_stonewall.json")
Now we repeat the procedure for the other segment using our sample from the terrestrial area.
sampledataset2 = af.DataSet(filename="sample_uw.las", spatial_reference="EPSG:25833")
pipeline2 = af.pipeline_tuning(sampledataset2)
af.save_filter(pipeline2, "filters/PK_uw.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 two pipelines you created earlier. Important: You need to select all pipelines that you will need to attach to the segment classes in the following step. Select multiple pipelines by pressing "Ctrl" ("Strg" on German keyboards). As 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. On the top right, switch in the dropdown menu to "segment". You will see that there are now two segments, "UW" for "underwater" and "veg" for "vegetation". You can also visualize their outline by clicking on the marker button next to the segment name. In the dropdown menu below each 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
)
assigned_segmentation.save("data/PK_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
PK_last.lazis in directorydata - file
PK_segments_assigned.geojsonis in directorydata - Directory
outputwill be created indata/output - filters are stored in
filters
afwizard --dataset=data/PK_last.laz --dataset-crs=EPSG:25833 --segmentation=data/PK_segments_assigned.geojson --segmentation-crs=EPSG:25833 --output-dir=data/output --library filters --lastools-dir="C:\LAStools"
We can also execute the command line tool by using the !command syntax in the Jupyter notebook. The exclamation mark tells Jupyter that it should pass the command to the shell instead of executing with Python. Let's try it out below.
!afwizard --dataset=data/PK_last.laz --dataset-crs=EPSG:25833 --segmentation=data/PK_segments_assigned.geojson --segmentation-crs=EPSG:25833 --output-dir=data/output --library filters --lastools-dir="C:/LAStools"