Generating LiDAR flight line boundaries with PDAL's tindex

A common task after an airborne LiDAR data collection is to determine the boundaries of the LiDAR LAS file strips, to visualize where they cover approximately. The screenshot below shows an example of how these LiDAR strips look like, colored by elevation. 

Example of LiDAR data strips

The boundaries can be generated in various formats using PDAL's tindex command. Using this command is simple, as shown in the example steps below:

1. Open up the OSGeo4W Command Prompt. Change the directory to the folder containing the *.laz files.
   
   

   

   A folder containing LiDAR *.las file strips

   
   
2. At the prompt, type in the command:
   
   C:\> pdal tindex boundary.sqlite *.laz -f SQLite --lyr_name "bnd"
   
   where boundary.sqlite is the output file
   *.laz specifies the input files using a wildcard
   -f SQLite is the output format
   --lyr_name "bnd" specifies the layer name to append to the output name
3. Press RETURN.
   
   The boundaries are generated.
   
   
4. Optional. After the process is completed, drag and drop the output file e.g. boundary.sqlite onto QGIS.
   
   The LiDAR strips boundary_bnd_polygon are displayed.
   
   
   The output boundaries have the following attributes as shown below.
   

Android Studio Emulator: How I resolved a Google Maps SupportMapFragment GLThread exception

While attempting to open up a Google Maps SupportMapFragment in the Android Studio Emulator, the Activity crashed with the following error messages:

[ 08-07 00:50:44.423  3266: 3523 W/         ]
                                                                     Unrecognized GLES max version string in extensions: ANDROID_EMU_CHECKSUM_HELPER_v1 ANDROID_EMU_dma_v1 
08-07 00:50:44.429 3266-3523/com.dom925.volcamr.nz E/EGL_emulation: rcCreateContext returned 0
08-07 00:50:44.429 3266-3523/com.dom925.volcamr.nz E/EGL_emulation: tid 3523: eglCreateContext(1434): error 0x3003 (EGL_BAD_ALLOC)
                                                                    
                                                                    --------- beginning of crash
08-07 00:50:44.429 3266-3523/com.dom925.volcamr.nz E/AndroidRuntime: FATAL EXCEPTION: GLThread 276
                                                                     Process: com.dom925.volcamr.nz, PID: 3266
                                                                     java.lang.RuntimeException: createContext failed: 12291
                                                                         at com.google.maps.api.android.lib6.gmm6.vector.az.a(:com.google.android.gms.DynamiteModulesB:29)
                                                                         at com.google.maps.api.android.lib6.gmm6.vector.ba.f(:com.google.android.gms.DynamiteModulesB:158)
                                                                         at com.google.maps.api.android.lib6.gmm6.vector.ba.run(:com.google.android.gms.DynamiteModulesB:11)

The same Android app has no issue running on a real handset.

After some fiddling around with the emulator settings, the problem was resolved by changing the emulator's graphics emulated performance from automatic to either Hardware GLES 2.0 or Software GLES 2.0 as shown below.

PDAL: Colorize a LAS file with multiple GeoTiff images

A LiDAR LAS file may not necessarily share the same bounds as a single raster image file. More often the case, to cover the LAS file more than one raster image file is needed. One way to resolve this is to merge all the raster images into a single mosaic file. Another way may be to keep the raster images as they are and load them in one by one (just-in-time fashion) to colorize the LAS file. This workflow is a little more complicated but it can be done using PDAL's pipeline processing mechanism.

In this post, PDAL will be used to colorize a LAS file using multiple tiles of GeoTIFF image files - four to be exact. You just need to create a pipeline that loads in the LAS file, apply the RGB values from each image file individually and write the colored points into an output LAS file. The following sections illustrate this workflow.

Here's how the source LAS and GeoTIFF files look like.

The source LiDAR LAS file displayed by elevation

The source GeoTIFF raster image files

Create a PDAL pipeline JSON file

1. Using a text editor, type in the JSON syntax to colorize a LAS file using PDAL's colorization filter, as shown below.
   
   
   where
   uncompahgre.laz is the source LAS file,
   tile0.tif, tile1.tif, tile2.tif and tile3.tif are the source image files,
   and color.laz is the output LAS file name
   
2. Save the text into a file e.g. colorLas.json.
   

A sample pipeline JSON is show in the listing below.

{
    "pipeline": [
        "uncompahgre.laz",
        {
            "type": "filters.colorization",
            "raster": "tile0.tif"
        },
        {
            "type": "filters.colorization",
            "raster": "tile1.tif"
        },
        {
            "type": "filters.colorization",
            "raster": "tile2.tif"
        },
        {
            "type": "filters.colorization",
            "raster": "tile3.tif"
        },
        {
            "type": "writers.las",
            "compression": "true",
            "minor_version": "2",
            "dataformat_id": "3",
            "filename":"color.laz"
        }
    ]
}

Run the colorization process

1. Open up a OSGeo4W Command Prompt.
2. At the prompt, type in the pdal pipeline command:
   
   c:\> pdal pipeline colorLas.json
   where colorLas.json is the pipeline JSON file created in the previous section.
3. Run the command.
   
   Processing messages appear.
   
   
   The output colorized LAS file color.laz is generated.
4. Optional. Display the resultant colored LAS file in a LAS Viewer.
   
   

Using Orfeo ToolBox to combine grayscale TIFF images into RGB and CIR composites

For working with geo-referenced TIFF images, Orfeo ToolBox (OTB) provides a useful convenience function otbcli_concatenateImages to combine separate grayscale images into a composite file. Examples of separate bands of grayscale images are shown below.

For instance, the separate grayscale bands representing the red, green and blue channels can be combined into a single RGB composite image. To do this using OTB, the following steps can be done.

1. Open up a OSGeo4W Command Prompt.
2. At the prompt, type in the otbcli_concatenateImages command with options:
   
   C:\> otbcli_concatenateImages -il band_red.tif band_green.tif band_blue.tif -out rgb.tif uint16
   
   where
   band_red.tif, band_green.tif, and band_blue.tif are the input files in the correct order,
   rgb.tif is the output RGB composite file
   uint16 is the output file datatype. By default it is float32 but that is too much. 
3. Run the command.
   
   Processing messages appear. The bands are merged into the output file rgb.tif.
   
4. Optional. Display the resultant file in a suitable viewer.
   

Similarly, to combine the separate grayscale bands into a composite Colored Infra-Red CIR file, the following steps can be done.

1. Open up a OSGeo4W Command Prompt.
2. Type in the otbcli_concatenateImages command.
   
   C:\> otbcli_concatenateImages -il band_nir.tif band_red.tif band_green.tif -out cir.tif uint16
   Where
   band_nir.tif, band_red.tif, and band_green.tif are the ordered input bands
   cir.tif is the output file name
   uint16 is the output file data type
3. Run the command.
   
   Processing messages appear. The bands are merged into the output file cir.tif.
   
4. Optional. Display the resultant file in a suitable viewer.
   

Use PDAL to generate a DEM from a LiDAR LAS file

PDAL can be used to generate a GeoTIFF digital elevation model (DEM) from a LiDAR point cloud LAS file via the GDAL writer driver. The resultant elevation GeoTIFF file may not be displayed as what is expected in some software such as Global Mapper because the PDAL generated file can have one or more bands of data including an alpha channel, depending on the PDAL options used; refer to https://www.pdal.io/stages/writers.gdal.html for the types of output bands that can be written to the DEM GeoTIFF file. For normal display, it might be necessary to separate the bands into individual files - in this post, the OTB (Orfeo Tool Box) otbcli_splitImage executable is used.

The following steps show how to generate a GeoTIFF DEM and split the resultant file into separate bands.

Generate the DEM

1. Open up a text editor. Type in the following JSON text. Save the text into a file e.g. pdal_dtm.json
   
   
   
   Note: the JSON pipeline text will do the following:
   (a) Load in a LAS file e.g. autzen.laz
   (b) Filter away non-ground classified point
   (c) Write only the mean Z values to a GeoTIFF file dem.tif with a resolution of 10 meters
   
2. Open up the OSGeo4W Command Prompt.
   
   The Command Prompt appears.
   
3. Type in the PDAL command to process the pipeline. Run the command.
   
   C:\> pdal pipeline pdal_dtm.json
   
   The GeoTIFF file dem.tif is generated.
   

Splitting the bands

Displaying the resultant dem.tif containing only the mean Z band in QGIS will show something resembling a grayscale terrain with masked areas, as shown in the screenshot below.

However, displaying the resultant GeoTIFF in Global Mapper will show a black rectangle, as shown below, due to the presence of an alpha channel. This need to be removed for proper interpretation of the elevation data.

To split the bands, do the following:

1. In the OSGeo4W Command Window, type in the command:
   
   C:\> otbcli_splitimage -in dem.tif -out split.tif
   
   
   
   The input file dem.tif is split into separate files with a numbered suffix starting from 0 e.g. split_0.tif.
   
2. Now, the file can be displayed properly in Global Mapper.
   

Applying a simple X,Y shift or translation to a raster GeoTIFF file using GDAL

Sometimes when comparing a raster ortho-mosaic GeoTIFF file to ground control points (GCPs), the raster file may appear to be slightly shifted relative to the ground control points. To resolve this problem correctly, it may be necessary to regenerate the ortho-mosaic after adjusting the tie points. But sometimes, if the shift is linear or you want to do a fast correction, then the open source software GDAL may be able to solve the problem for you.

To shift a raster GeoTIFF file e.g. input.tif with GDAL, do the following:

1. Open up a Windows Command Prompt.
   
2. Type in the gdal_translate command with the a_ullr option:
   
   C:\> gdal_translate -a_ullr 760726.437 4557390.415 772996.466 4540826.364 input.tif translate.tif
   
   Note: 
   translate.tif is the output filename
   -a_ullr specifies the new upper left X, upper left Y, lower right X, lower right Y coordinates  
   
3. Run the command.
   
   Processing messages appear. The file is shifted to the new coordinates.
   
   
   The file is shifted to the new coordinates.
   

Mask a LAS file using PDAL and QGIS

The Point Data Abstraction Library (PDAL) comes with a useful function to mask or crop out a LiDAR LAS file with one or more polygons. The example screenshot below shows a point cloud overlaid with a red polygon in the upper right corner, which outlines the desired area of the point cloud to be cropped.

The cropping can be done using PDAL's crop filter but it requires the cropping polygon to be specified in the Well Known Text (WKT) string format. This is a bit of pain but can be overcome using a few methods, one of which is shown below using QGIS and the Plain Geometry Editor plugin.

Define the cropping polygon

1. In QGIS, draw a new polygon e.g. mask, as shown below.
   
   
   
   Note: The mask should be created in the same coordinate system as the LAS file. 
2. Click the Plain Geometry Editor icon in the toolbar (red circle above). Click on the clipping polygon.
   
   The Plain Geometry Editor dialog box appears.
   
   
   Note: Install the Plain Geometry Editor plugin if the icon is not available.
3. In the text field, select and copy all the polygon WKT text string into the Windows clipboard.  

Create a PDAL processing pipeline JSON file

1. In a text editor, type in something similar to the example below.
   
2. From the Windows Clipboard, paste the WKT string from the previous section to the "polygon" attribute and surround it with double quote '"' characters.
   
3. Save the JSON file e.g. process.json.

{
  "pipeline":[
    "autzen.laz",
    {
      "type":"filters.crop",
      "polygon": "Polygon ((638500.66904077248182148 853359.34703735215589404, 638869.71793351718224585 853365.15883093187585473, 638881.34152069012634456 853208.24040409235749394, 638677.92874516302254051 853199.52271371346432716, 638878.43562390014994889 852818.85023379256017506, 638733.14078423334285617 852842.09740813483949751, 638611.09311891463585198 852975.76866063219495118, 638416.39803376235067844 853211.14630088687408715, 638500.66904077248182148 853359.34703735215589404))"
    },
    {
      "type":"writers.las",
      "filename":"file-cropped.las"
    }
  ]
}

Note: 

• The pipeline JSON file stores the processes to be done in sequence in an array bracketed by the '[' and ']' characters. 
• autzen.laz is the input LAS file for this example
• filters.crop is the process to apply using the "polygon" attribute.
• file-cropped.las is the output LAS file. 

Run the cropping process

1. Open up the OSGeo4W Shell.
   
2. At the prompt, type in the pdal pipeline command:
   
   C:\> pdal pipeline process.json -v 4
   Processing messages appear. The file is cropped.
   
   
3. Optional. Using your preferred LAS Viewer, open up the resultant cropped LAS file.
   
   The cropped file showing only the cropped area is displayed.
   

Combine separate gray scale TIFF images into a single RGB TIFF image

Some multi spectral cameras capture images as individual gray scale TIFF images - each file containing data for a single light wavelength e.g. red, green, blue, or infra-red. Example files are the gray scale TIFF files shown below, each file representing a single color band.

Combining these individual color bands into a single RGB color composite can be done using the free and open source ImageMagick software, as shown in the steps below.

1. Open up a Command Prompt.
2. At the prompt, type in the ImageMagick convert command:
   
   C:\> magick convert -verbose band_red.tif -channel R band_green.tif -channel G band_blue.tif -channel B -combine -channel RGB -alpha off -colorspace sRGB out_rgb.tif
   
   Note:
   band_*.tif are the input TIFF files. Each input file must be followed with the correct channel option.
   out_rgb.tif is the output RGB TIFF file.
   -verbose is used to print out processing messages.
   -combine -channel RGB are for combining the bands.
   -alpha off is used to disable the creation of the alpha channel in the output RGB file.
   -colorspace sRGB is used to set the output color space.
3. Run the command.
   
   Processing messages appear. The output file out_rgb.tif is created.
   
4. Open up the resultant file in an image editor, e.g. GIMP.
   
   The file is displayed in color and the the bands are combined into the image's red, green and blue channels.
   

Using PDAL to classify isolated LiDAR points as noise

LiDAR data often contains noise and it is necessary to identify and/or remove them. An example of a LAS file containing noise in the form of low isolated points beneath the ground is shown in the screen shot below.

Isolated points can be easily identified by using statistical filtering methods, which the PDAL open source software has.

To filter out these points using PDAL, perform the following steps.

1. Open up the OSGeo4W Shell.
   
   The OSGeo4W Shell prompt appears.
2. In the prompt, type in the command:
   
   C:\> pdal translate -i in_noisy.las -o out_filtered.las outlier --filters.outlier.method="statistical" --filters.outlier.mean_k=8 --filters.outlier.multiplier=3.0 -v 4
   
   Note:
   -i in_noisy.las is the input LAS file
   -o out_filtered.las specifies the output LAS file
   outlier tells PDAL to apply the outlier filter
   --filters.outlier.**** options specify the various outlier parameters
   -v 4 indicates the processing messages verbosity level
   
   
3. After running the command, the isolated points are classified as Class 7 - Low Noise points in the output LAS file.
   
   The point cloud colored by classification.
   
   
   The resultant LAS file colored by elevation and with the class 7 - Low noise points turned off.
   

Extract GPS tags from photographs into a CSV file using Exiftool

There is a nice command line utility Exiftool from http://www.sno.phy.queensu.ca/~phil/exiftool/ that can be used to quickly extract out the GPS positions and other tags from photographs and other images into a comma separated values (CSV) file.

The following example uses the Windows executable version of the utility to illustrated the extraction steps.

1. Open up a Windows Command Prompt. In the prompt, type in the command.
   
   C:\> "exiftool(-k).exe" -n -gpslongitude -gpslatitude -gpstimestamp -csv D:\MyDocuments\temp\somephotos
   
   Note:
   -n means to print out as numbers only
   -csv prints out csv including the file path and name
   -D points to the folder directory of the photographs to extract
2. Press RETURN
   
   The extraction values are displayed to the screen.
   
   
3. To output to a file, use the > character to redirect the standard output to a file, e.g. type in the command:
   
   C:\> "exiftool(-k).exe" -n -gpslongitude -gpslatitude -gpstimestamp -csv D:\MyDocuments\temp\somephotos > outgps.csv
   
   The output file outgps.csv is created.
   
4. Display the resultant file in a spreadsheet. Or plot the locations on a map.
   
   

Simple LiDAR ground points classification and segmentation using PDAL

PDAL (Point Data Abstraction Library) comes with a couple of options to segment point clouds by classifying LiDAR ground points (an example unclassified point cloud is shown below) - Simple Morphological Filter (SMRF) or Progressive Morphological Filter (PMF).

I have found the SMRF method to be fast and produce reasonable results while the PMF method seems to take a much longer time to do the job. The steps to run ground classification on a LAS file are describe below.

1. In Windows, open up the OSGeo4W Shell.
   
   The OSGeo4W Shell is displayed.
2. In the OSGeo4W prompt, type in and run the command:
   
   C:\> pdal translate -i unclassified.las -o ground.las smrf -v 4
   
   Notes:
   -i unclassified.las is the input file
   -o ground.las specifies the output file
   smrf is the option to apply the Simple Morphological Filter
   -v 4 is the processing messages verbosity level
   
   Processing messages appear.
   
3. Display the ground classified LAS file in a viewer.
   
   
   
   
4. To use the Progressive Morphological Filter to perform the ground classification, type in the following command:
   
   C:\> pdal translate -i unclassified.las -o ground.las pmf -v 4
   
   Notes:
   -i unclassified.las is the input file
   -o ground.las specifies the output file
   pmf is the option to apply the Progressive Morphological Filter
   -v 4 is the processing messages verbosity level

Merging multiple comma separated values CSV files

If there are many comma separated value CSV files with headers and you want to merge them into a single CSV file, it can be a pain having to do it by hand. Fortunately, there are ways to automate the task. One method which I like is to use the tail command from Unix. For Windows, a tail utility can be downloaded from http://tailforwin32.sourceforge.net/;There are others which a Google search can reveal.

Below are screenshots of a few CSV files to merge.

1. Using a text editor, create a script or batch file. In the editor, type in the tail command to create a new merged file with a header.
   
   tail -n +1 -q green.csv > merge.csv
   
   Note: -n +1 means to start from the first line.
   A single > means to output to a new file
   
2. Type in the commands to append lines from subsequent files without headers to the output file.
   
   tail -n +2 -q orange.csv >> merge.csv
   tail -n +2 -q transfer.csv >> merge.csv#...etc...
   Note: -n +2 means to start from the 2nd line
   >> means to append to the output file
   
3. Run the script or batch file in a Command Prompt.
   
   The CSV files are merged.
   

Resolving Python Module 'numpy' has no 'xxxx' member error message in Visual Studio Code

While writing some Python code in Visual Studio Code that calls some Numpy classes, Pylint error messages appear complaining about non-existing Numpy members, e.g. the log2 method, and red wavy lines underline the offending code, as shown in the screenshot below.


One way to clear this up is to get Pylint to white-list Numpy using Visual Studio Code's Settings.

1. In Visual Studio Code, select File | Preferences | Settings.
   
   The settings.json file is displayed in the editor.
   
   
2. Click the Workspace Settings tab.
   
   The settings appear in the editor.
   
   
3. In the editor, type in the following and save the file.
   
   { "python.linting.pylintArgs" : [ "--extension-pkg-whitelist=numpy" ] }
   
   
   
   
4. Close and reopen the folder in Visual Studio Code.
   
   The missing member messages and the red wavy underlines no longer appear.