How to filter an elevation from point cloud LIDAR las file and save it as 4th band of RGB band raster tif file

[amuthalingeswaranbose]

I am trying to filter an elevation from a point cloud LIDAR LAS file and save it as the 4th band of an RGB band raster TIFF file. The RGB band raster file has 3 bands, and I want to add the elevation band to it. I am using the following code:

import gdal

def filter_elevation(las_file, output_raster):
"""
Filters an elevation from a point cloud LIDAR LAS file and saves it as the 4th band of an RGB band raster TIFF file.

Args:
las_file: The input LAS file.
output_raster: The output RGB raster file.
"""

Read the LAS file.
las = gdal.Open(las_file)

Get the number of points in the LAS file.
n_points = las.GetRasterCount()

Create a new raster with 4 bands.
output_raster = gdal.GetDriver("GTiff").Create(output_raster, las.RasterXSize, las.RasterYSize, 4)

Extract the elevation values from the LAS file.
elevation_values = las.GetRasterBand(1).ReadAsArray()

Set the elevation values to the 4th band of the output raster.
output_raster.SetRasterBand(3, elevation_values)

Set the NoData value to -9999.
output_raster.SetNoDataValue(-9999)

Write the output raster.
output_raster.Flush()

However, this code is not working. I am getting the following error:

RuntimeError: Unable to create band 4.

I would appreciate any help you can provide.

Steps to reproduce

1. Create a point cloud LIDAR LAS file with 1 band.
2. Run the following code:

import gdal

def filter_elevation(las_file, output_raster):
"""
Filters an elevation from a point cloud LIDAR LAS file and saves it as the 4th band of an RGB band raster TIFF file.

Args:
las_file: The input LAS file.
output_raster: The output RGB raster file.
"""

Read the LAS file.
las = gdal.Open(las_file)

Get the number of points in the LAS file.
n_points = las.GetRasterCount()

Create a new raster with 4 bands.
output_raster = gdal.GetDriver("GTiff").Create(output_raster, las.RasterXSize, las.RasterYSize, 4)

Extract the elevation values from the LAS file.
elevation_values = las.GetRasterBand(1).ReadAsArray()

Set the elevation values to the 4th band of the output raster.
output_raster.SetRasterBand(3, elevation_values)

Set the NoData value to -9999.
output_raster.SetNoDataValue(-9999)

Write the output raster.
output_raster.Flush()

Save the code as a Python file.
Run the Python file.
Expected behavior
The code should successfully filter the elevation from the LAS file and save it as the 4th band of the RGB band raster TIFF file.

Actual behavior
The code throws the following error:

RuntimeError: Unable to create band 4.

Versions

• GDAL version: 3.2.1
• Operating system: MacOS 10.15.7

Additional information

I have attached the LAS file and the RGB band raster TIFF file.
I hope this helps!