Wednesday, June 28, 2017

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.
      "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))"

  • 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.

Wednesday, June 21, 2017

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

    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.

Monday, June 12, 2017

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

    -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.

Monday, June 5, 2017

Extract GPS tags from photographs into a CSV file using Exiftool

There is a nice command line utility Exiftool from 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

    -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.

Monday, May 29, 2017

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

    -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

    -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

Monday, May 22, 2017

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;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
    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.

Sunday, May 14, 2017

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.

Monday, May 1, 2017

Create labels from multi-columns in QGIS

It is quite simple to create labels from one or more database columns in QGIS; just that you have to type in the correct syntax as shown below.

  1. Run QGIS. Display some layers, e.g. stations.

  2. In the Layer Panel, mouse right click on the layer to be labelled, e.g. stations.

    A pop up menu appears

  3. Choose Properties.

    The Layer Properties dialog box appears.

  4. Choose Labels. In the combo box, select Show labels for this layer.

  5.  In the Label with field, click the Expression icon.

    The Expression dialog box appears.

  6. In the Expression tab, type in the syntax to form the multi-column label e.g.
    name || ' ' || code.

    Note: name and code are two columns from the layer database table, and || is the string concatenation operator and ' ' is a space character.
  7. Click OK. And Click OK to apply the changes.

    The labels appear.
  8. If multi-line labels are desired, then the newline character '\n' needs to be concatenated with the database fields in the Expression dialog box, as shown below.
    name || '\n' || code

    And the resultant map look like this.

Thursday, March 23, 2017

LiDAR data coordinate readout in Cloud Compare

It is useful to obtain a readout of LiDAR data points in Cloud Compare, perhaps to determine the height value, or the intensity value of a LiDAR point. This post illustrates how to read out coordinate values in Cloud Compare.
  1. In Cloud Compare, load in a LAS file, e.g. Serpent Mound Model.las.

    The Global Shift/scale dialog box appears.

    Note: Real LiDAR data have large coordinate values requiring double data types, but most 3-D software like Cloud Compare works with coordinates in float data types internally. So in order to preserve the LiDAR coordinates in Cloud Compare, it may be necessary to apply some coordinate offset and/or scaling so that the LiDAR data points can fit into float data types. So when you load a LAS file in Cloud Compare, it will calculate and recommend an offset/scaling values before importing, as shown in the screen shot below.
  2. Click Yes.

    The LAS file is loaded; its data points are globally shifted and scaled.
  3. In the toolbar, click the Point picking icon.

    A point picking toolbar appears on the top right corner.
  4. By default, the first command (Select one point and displays its info) on the left is active. If not, click on it.
  5. Click on one LAS point.

    A balloon appears showing the X, Y, Z and a scalar value.

    Note: the xl, yl, zl on the left are local float data type coordinates. The xg, yg, zg on the right are the actual LiDAR data point vertices. The Point Source ID is the current scalar value.
If you want to get the read out to show another scalar field, you have to set the desired field before entering the Point Pick command, as described below.
  1. In the Tree pane, select the point cloud node, e.g. Serpent Mound Model.

    The Properties appear below the Tree pane.
  2. In the Scalar Fields | Active field, choose the desired field, e.g. Intensity.

    The point colors change to intensity values.
  3. Now, use the Point Pick command to readout a point.

    The readout balloon now shows the selected scalar field - Intensity.

Monday, March 13, 2017

Batch sample or convert videos to animated GIF files WebApp

I wrote a WebApp VidzGif for batch sampling or converting video files into animated GIF images that works with just a Web browser without uploading the videos to a server for processing.

To use this WebApp, do the following:

  1. Open up the WebApp at with a Web browser.

  2. Optional. Click Settings. Set the desired parameters, e.g. sample interval, duration, etc.
  3. Click the Browse button and select one or more video files. Or drag and drop one or more video files into the dashed box.

    The videos are sampled and rendered into GIFs.
  4. To save out the GIF, right click on the image and choose Save image as.

Monday, March 6, 2017

Interactive cross-section in Cloud Compare

After loading LiDAR point clouds in Cloud Compare, it helps sometimes to cut a cross section or profile to get a feel for the land. There is a couple of Segment commands but the one I use for this purpose is the Cross Section Tool with the graphic interactors.

The following steps show how to use the Cross Section Tool.
  1. Load a point cloud file, e.g. Serpent Mound Model LAS Data.las, into Cloud Compare.
  2. Select the cloud in the DB Tree pane.

  3. Select Tools | Segmentation | Cross Section.

    The graphic interactors (the large 3D arrows pointing to the X, Y, Z axes) and Cross Section box appear
  4. Drag a cone of an interactor to thin the cross section box.

  5. Drag the torus of an interactor to rotate the cross section box.

  6. Adjust the interactors and camera view until the desired orientation is achieved.

    Note: if you wish to segment out the cross section into another cloud, click the Tick icon on the Cross Section tool box. Otherwise, click the X icon to exit the command when done.
Related Posts with Thumbnails