Lidar Samples

Lidar is used to create impactful 3D point clouds that represent terrain and surface features. It can be used to generate digital elevation models (DEMs), contour lines, and terrain models. Users can classify, filter, and colorize lidar points based on their return type, elevation or intensity and can even perform feature extraction. This allows features like buildings, vegetation, powerlines, and more to be extracted from a lidar point cloud.

Exporting an Elevation Layer

Note

Before running this sample, configure these variables:

• export_dir: the directory for resulting files to be exported to

A more complex sample script, showing a common workflow using Global Mapper Pro’s lidar tools: loading lidar data, creating an elevation layer from it, and exporting it to an elevation grid file. This sample grabs any loaded lidar layers, classifies the noise and ground points, and creates a class filter to only allow the use of ground points. An elevation grid is then generated, scaled to the size of the layer’s GlobalRect. Similar dimensions are then used to determine the export bounds for the newly created elevation grid. Then the elevation grid is exported to a directory specified by ‘resource_dir.’

# Resolution of the output dataset in meters
pixel_resolution = 1.0 
lidar_layers = []

# Load all LiDAR layers found into a list
for i in range(num_layers):
    info = gm.GetLayerInfo(all_layers[i])
    if info.mEnabled and info.mTypeName == "LIDAR_LAS":
        lidar_layers.append(all_layers[i])

# If there's no lidar layers, no need to run the rest of the script
if len(lidar_layers):
    # Classify noise points
    print("Classifying noise...")
    noise_setup = gm.GM_ClassifyNoiseSetup_t()
    # Only check unclassified points for noise
    noise_setup.mFlags += gm.GM_ClassifyNoise_UnclassifiedOnly
    # Detect any points that resemble noisy lidar returns and reclassify them (best practice to do this first)
    err_LidarClassifyNoise = gm.LidarClassifyNoise(
        lidar_layers, noise_setup, 0
    ) 

    # Classify the ground points
    print("Classifying ground...")
    # Default settings work well for ground classification
    err_LidarClassifyGround = gm.LidarClassifyGround(
        lidar_layers, gm.GM_ClassifyGroundSetup_t(), 0
    )

    # Create the class filter to only use ground points
    print("Applying class filter...")
    # Save the current settings so they can be restored afterward
    err_GetLidarClassFilter, prev_filter = gm.GetLidarClassFilter()
    # New filter object with all classes disabled
    err_LidarClassFilter_Init, class_filter = gm.LidarClassFilter_Init(False)
    # Enable ground points to pass through the filter
    err_LidarClassFilter_SetClassEnabled, mask = gm.LidarClassFilter_SetClassEnabled(
        class_filter, gm.GM_LidarClass_Ground, True
    )
    # Apply new filter: only ground points will be used in lidar operations now
    err_SetLidarClassFilter1 = gm.SetLidarClassFilter(class_filter) 

    # Get the dimensions of the layer in meters in case it isn't already
    bounds = gm.GetLayerInfo(lidar_layers[0]).mGlobalRect
    # Get the total width and height of the data in meters
    err_CalcDistance1, width = gm.CalcDistance(
        bounds.mMaxX, bounds.mMaxY, bounds.mMinX, bounds.mMaxY, 0
    )
    err_CalcDistance2, height = gm.CalcDistance(
        bounds.mMaxX, bounds.mMaxY, bounds.mMaxX, bounds.mMinY, 0
    )
    groundcrs_width = (abs(bounds.mMaxX - bounds.mMinX) / width / pixel_resolution) 
    groundcrs_height = (abs(bounds.mMaxY - bounds.mMinY) / height / pixel_resolution)

    # Create elevation grid layer from ground points at the approximate resolution
    print("Generating elevation grid...")
    grid_setup = gm.GM_GridGenSetup_t()
    grid_setup.mXRes = groundcrs_width
    grid_setup.mYRes = groundcrs_height
    grid_setup.mGridAlg = gm.GM_GridAlg_BinAverage
    # This value is high enough to fill gaps in most lidar data
    grid_setup.mTightnessMult = 32 
    # Define how the gaps are filled
    gap_fill = gm.GM_FillGapsParams_t()
    gap_fill.mGapAlg = gm.GM_GapFillAlg_IDW
    # Needs to be the same as grid_setup.mGridAlg
    gap_fill.mGridAlg = gm.GM_GridAlg_BinAverage
    grid_setup.mGapFillSetup = gap_fill
    err_GenerateElevationGrid, grid_layer, tin_layer = gm.GenerateElevationGrid(lidar_layers, grid_setup)
    # We don't need the TIN layer
    err_CloseLayer = gm.CloseLayer(tin_layer)

    # Reset the lidar class filter to previous settings
    err_SetLidarClassFilter2 = gm.SetLidarClassFilter(prev_filter)

    # Get the dimensions of grid layer in meters to use as the export size
    bounds = gm.GetLayerInfo(grid_layer).mGlobalRect
    err_CalcDistance3, grid_width = gm.CalcDistance(
        bounds.mMaxX, bounds.mMaxY, bounds.mMinX, bounds.mMaxY, 0
    ) 
    err_CalcDistance4, grid_height = gm.CalcDistance(
        bounds.mMaxX, bounds.mMaxY, bounds.mMaxX, bounds.mMinY, 0
    )
    # We are using the pixel_resolution variable for spacing
    # A higher pixel_resolution value means a bigger pixels, and therefore a smaller file
    pixel_width = round(grid_width / pixel_resolution) 
    pixel_height = round(grid_height / pixel_resolution)

    # Export the new elevation data
    print("Exporting elevation grid layer...")
    err_ExportElevation = gm.ExportElevation(
        export_dir + "ExampleDem.dem",
        gm.GM_Export_DEM,
        grid_layer,
        bounds,
        pixel_width,
        pixel_height,
        0,
        gm.GM_ElevUnit_Centimeters,
    )
else:
    print("There are no lidar layers open.")

Pixels to Points

Note

Before running this sample, configure these variables:

• image_folder_path: the file path to the directory that contains the images for a P2P operation

• gc_points_file: the file path to a text file containing ground control points

• export_dir: the absolute path to the desired export directory

This script demonstrates a pixels to points (P2P) process. First a setup object is configured. This involves assigning the image_folder_path and image mask to use for the imported images. Then the export directories and filenames are established for cloud, image, and mesh outputs that will be produced once the PixelsToPoints function is called. Now, the ground control points (GCPs) are read from a text file, extracting the image name, pixel coordinates, and real-world coordinates. It groups image pixel observations by GCP name and stores them in a dictionary. Then it creates control point objects (GM_P2PControlPoint_t) with the associated image pixels and adds them to the GCP list for use in the Pixels to Points process. The GCP list is then converted into an array so that it can be assigned to ‘aSetup.mControlPointList.’ Finally, the PixelsToPoints function is called, passing in the created setup object.

# The espg code for Mercator which is the correct projection for this dataset
ESPG_Code = 4326
err_LoadProjectionFromEPSGCode, projection = gm.LoadProjectionFromEPSGCode(ESPG_Code)
gm.SetProjection(projection)

# Import a sample for Pixels to Points
err_PixelsToPointsInitSetup, aSetup = gm.PixelsToPointsInitSetup()

# Path to image Folder
aSetup.mImageFolder = image_folder_path
# Image file name
aSetup.mImageFilenameMask = "*.JPG"
# Image reduce factor reduces image size, 2 is half size.
aSetup.mImageReduceFactor = 2

# Cloud Layer Name & Path
aSetup.mCloudLayerDesc = "P2P Cloud Layer"
aSetup.mCloudLayerFilename = export_dir + "p2p_test_cloud_output.gmp"
# Image Layer Name & Path
aSetup.mImageLayerDesc = "P2P Image Layer"
aSetup.mImageLayerFilename = export_dir + "p2p_test_image_output.gmp"
# Mesh Layer Name & Path
aSetup.mMeshLayerDesc = "P2P Mesh Layer"
aSetup.mMeshLayerFilename = export_dir + "p2p_test_mesh_output.obj"
# Sets Resampling Method to Nearest Neighbors
aSetup.mImageSamplingMethod = 0

# Get the ground control points from an external file
gcp_file = open(gc_points_file, "r")
gcp_text = gcp_file.readlines()
gcp_file.close()
# Define variable pix_arr null
pix_arr = None
# iterate through P2P_GCP text file
gcp_list = []

for line in gcp_text:
    parts = line.strip().split(",")
    if len(parts) < 7:
        continue

    gcp_name = parts[0]
    image_name = parts[1]
    pixel_x = float(parts[2])
    pixel_y = float(parts[3])
    world_x = float(parts[4])
    world_y = float(parts[5])
    elev = float(parts[6])

    # Create a new control point object for each line
    cp = gm.GM_P2PControlPoint_t()
    cp.mName = gcp_name
    # Set lat/lon/elev
    geo = gm.GM_Point3D_t()
    geo.mX, geo.mY, geo.mElev = world_x, world_y, elev
    cp.mLatLonPos = geo
    # Create the pixel object
    pix = gm.GM_P2PImagePixel_t()
    pix.mImageFilename = image_name
    pix.mPixelPos = gm.GM_Point_t(pixel_x, pixel_y)

    # Create an array and assign the single pixel
    pixel_array = gm.GM_P2PImagePixel_array(1)
    pixel_array[0] = pix
    cp.mImagePixelList = pixel_array
    cp.mImagePixelCount = 1
    gcp_list.append(cp)

# Variable GCP Array is defined as list GCPS
# Converts ground control points from list format to array for
# use with function gm.PixelsToPointsInitSetup()
gcp_array = gm.GM_P2PControlPoint_array(len(gcp_list))
for i in range(len(gcp_list)):
    gcp_array[i] = gcp_list[i]

# Attribute mControlPointList is assigned to gcp_array
# Attribute mControlPointCount is the length of list GCPS
aSetup.mControlPointList = gcp_array
aSetup.mControlPointCount = len(gcp_list)
# aReserved assigned to default, 0
aReserved = 0
# PixelsToPoints function is called
err_PixelstoPoints, info = gm.PixelsToPoints(aSetup, aReserved)

Colorize a Lidar Layer

This sample demonstrates how to use the lidar colorize function to update the color intensity fields of lidar points. The newly colorized points can be seen in Global Mapper if “Color by RGB/Elev” is selected in the Lidar Draw Mode dropdown list. This sample colorizes the lidar layer using 1 image layer, though many could be used. To do this, load more colored reference image layers into img_list.

img_list = [layers_list[1]]
colorize_setup = gm.GM_LidarColorizeSetup_t()
colorize_setup.mIntensityMethod = gm.GM_LidarColorizeInten_Default
colorize_setup.mFlags = gm.GM_LidarColorize_DefaultFlags
colorize_setup.mFlags += gm.GM_LidarColorize_HideProgress
err_LidarColorize = gm.LidarColorize(layers_list[0], img_list, colorize_setup, 0)

Thin and Compare Lidar

Note

Before running this sample, configure these variables:

• run_script_dir: the absolute path of the directory to be exported to

A sample that shows how the LidarThin and LidarCompare functions work. The first point cloud is thinned, this reduces the file size and can help improve processing speed for large datasets. The two point clouds are then compared, the differences between them are placed in a newly generated layer. LidarCompare can help identify geographic changes over time in a specific area. Afterwards, all of the loaded LAZ files are exported to a specified directory.

Learn how RunScript was used to run commands from the Global Mapper scripting language here.

# Thin the point clouds for both lidar layers
thin_setup = gm.GM_LidarSpatialThinSetup_t()
thin_setup.mCreateNewCloud = False
err_LidarThin = gm.LidarThin(layers_list[0], thin_setup, 0)

thin_setup = gm.GM_LidarSpatialThinSetup_t()
thin_setup.mCreateNewCloud = False
err_LidarThin = gm.LidarThin(layers_list[1], thin_setup, 0)

# This function can compare multiple lidar layers against each other,
# but for this sample only 2 layers are being compared
# Compare the Kilauea_2009.laz to the Kilauea_2018.laz to view changes
list1 = [layers_list[1]]  # Thinned
list2 = [layers_list[0]]  # Original
compare_setup = gm.GM_LidarCompareSetup_t()
compare_setup.mFlags = gm.GM_LidarCompare_CreateNewLayers
compare_setup.mNewLayerDesc = "Kilauea_changes.laz"
err_LidarCompare = gm.LidarCompare(list1, list2, compare_setup, 0)

# Get handle for the loaded layers again to include the layer created by LidarCompare()
arr_ptr, arr_size = gm.GetLoadedLayerList()
layers_array = gm.GM_LayerHandle_array_frompointer(arr_ptr)
layers_list = gm.carray_to_list(layers_array, arr_size)

# Exports the lidar data, reflecting new class changes
for i in range(len(layers_list)):
    # Get the name of the layer to be exported
    layer_name = gm.GetLayerInfo(layers_list[i]).mDescription
    if ".laz" in layer_name:
        # Use the RunScript func to export lidar data as a .laz file
        export_string = "EXPORT_VECTOR "
        export_string += "FILENAME=" + run_script_dir + "new_" + layer_name + " "
        export_string += "EXPORT_LAYER=" + layer_name + " "
        export_string += "TYPE=LIDAR_LAS "
        export_string += "ELEV_UNITS=METERS "
        export_string += "VERT_CS_CODE=NAVD88 "
        err_RunScript, handle, size = gm.RunScript(
            export_string, gm.GM_LoadFlags_UseDefaultLoadOpts, 0
        )

Lidar Query Samples

Classifying Query Points

Note

Before running this sample, configure these variables:

• run_script_dir: the absolute path of the directory to be exported to

This sample demonstrates a workflow using a lidar query. It loads multiple different lidar layers into a query. Query functions are then applied to classify all of the lidar layers in the query, demonstrating how multiple lidar layers can be classified with one function. The layers are then exported and some of their statistics are displayed.

Learn how RunScript was used to run commands from the Global Mapper scripting language here.

# Create the query of lidar layers (from a list of lidar layer handles)
err_CreateLidarQuery, query_pointer, query_size = gm.CreateLidarQuery(
    layers_list, None, None
)

# Gets general info from the query
print("\nInitial Lidar Query Data:")
err_GetLidarQueryInfo, query_info = gm.GetLidarQueryInfo(query_pointer)
print("Number of layers: " + str(query_info.mLayerCount))
print("Total points among layers: " + str(query_info.mNumPoints) + "\n")

# Classify all lidar layers in the query (Default Settings)
print("Classifying noise...")
noise_setup = gm.GM_ClassifyNoiseSetup_t()
# Find any noise and reclassify it
noise_setup.mFlags += gm.GM_ClassifyNoise_UnclassifiedOnly
err_LidarClassifyNoiseFromQuery = gm.LidarClassifyNoiseFromQuery(
    query_pointer, noise_setup, 0
)
print("Classifying ground...")
err_LidarClassifyGroundFromQuery = gm.LidarClassifyGroundFromQuery(
    query_pointer, None, 0
)
print("Classifying buildings and vegetation...")
err_LidarClassifyBuildingVegFromQuery = gm.LidarClassifyBuildingVegFromQuery(
    query_pointer, None, 0
)

# Free memory
err_FreeLidarQuery = gm.FreeLidarQuery(query_pointer)

# Export all of the newly classified lidar data
for i in range(len(layers_list)):
    # Get the name of the layer to be exported
    layer_name = gm.GetLayerInfo(layers_list[i]).mDescription
    # Use the RunScript func to export lidar data as a .laz file
    export_string = "EXPORT_VECTOR "
    export_string += "FILENAME=" + run_script_dir + "new_" + layer_name + " "
    export_string += "EXPORT_LAYER=" + layer_name + " "
    export_string += "TYPE=LIDAR_LAS "
    export_string += "ELEV_UNITS=METERS "
    export_string += "VERT_CS_CODE=NAVD88 "
    err_RunScript, handle, size = gm.RunScript(
        export_string, gm.GM_LoadFlags_UseDefaultLoadOpts, 0
    )

    # Displays stats about each layer
    err_GetLayerLidarStats, stats = gm.GetLayerLidarStats(layers_list[i], 0)
    print("\n" + layer_name + " Stats:")
    print("Average Elevation: " + str(stats.mAvgElev))
    print("Standard Deviation (Elevation): " + str(stats.mStdDevElev))
    print("Average Intensity: " + str(stats.mAvgIntensity))
    print("Standard Deviation (Intensity): " + str(stats.mStdDevIntensity))

Classifying a Graph and Extracting Query Features

Note

Before running this sample, configure these variables:

• run_script_dir: the absolute path of the directory to be exported to

This sample demonstrates how to classify a graph and extract vector features from a lidar query. Ground, building, and high vegetation points are classified from the query as a graph-classification. Then all building and tree features are extracted from the lidar layers. Finally, the generated vector layers and loaded LAZ layers are exported to a specified directory. The exported vector features will keep their previous attribute values.

Learn how RunScript was used to run commands from the Global Mapper scripting language here.

# Create the query of lidar layers (from a list of lidar layer handles)
err_CreateLidarQuery1, query_pointer1, query_size1 = gm.CreateLidarQuery(
    layers_list, None, None
)

# Classify ground points and then classify graph for building and vegetation points
err_LidarClassifyGroundFromQuery = gm.LidarClassifyGroundFromQuery(
    query_pointer1, None, 0
)
seg_opts = gm.GM_SegmentationSetup_t()
seg_opts.mMaxMahalanobisDistance = 1
seg_opts.mMinClusterSize = 6
build_veg_opts = gm.GM_ClassifyBuildingVegSetup_t()
build_veg_opts.mGridBinSizeUnits = gm.GM_LidarResMult_Meters
build_veg_opts.mMeanGridBinMult = 3.4239962
build_veg_opts.mVegMinHeightAboveGround = 0
build_veg_opts.mBuildingMinHeightAboveGround = 0.5
build_veg_opts.mPlanarAngleMaxDiff = 5
build_veg_opts.mPlanarThresholdM = 0.079999998
build_veg_opts.mRoofInvalidMaxAngle = 78
build_veg_opts.mRoofInvalidMinAngle = 65
build_veg_opts.mTreeThresholdM = 0.15
graph_opts = gm.GM_ClassifyGraphSetup_t()
graph_opts.mTypes = 6
graph_opts.mFlags = 0
graph_opts.mSegmentationOpts = seg_opts
graph_opts.mBuildingVegOpts = build_veg_opts
err_LidarClassifyGraphFromQuery = gm.LidarClassifyGraphFromQuery(
    query_pointer1, graph_opts, 0
)

# Extract all buildings and trees from the query
extract_opts = gm.GM_LidarExtractSetup_t()
extract_opts.mTypes = gm.GM_LidarExtract_Trees
extract_opts.mTypes += gm.GM_LidarExtract_Buildings
err_LidarExtractFeaturesFromQuery, generated_layers = gm.LidarExtractFeaturesFromQuery(
    query_pointer1, extract_opts, 0
)

# Free memory
err_FreeLidarQuery1 = gm.FreeLidarQuery(query_pointer1)

# Get handle for the loaded layers again to access generated layers
arr_ptr, arr_size = gm.GetLoadedLayerList()
layers_array = gm.GM_LayerHandle_array_frompointer(arr_ptr)
layers_list = gm.carray_to_list(layers_array, arr_size)

for i in range(len(layers_list)):
    # Get the name of the layer to be exported
    layer_name = gm.GetLayerInfo(layers_list[i]).mDescription
    # Export all of the newly classified lidar data
    if ".laz" in layer_name:
        # Use the RunScript func to export lidar data as a .laz file
        export_string = "EXPORT_VECTOR "
        export_string += "FILENAME=" + run_script_dir + "new_" + layer_name + " "
        export_string += "EXPORT_LAYER=" + layer_name + " "
        export_string += "TYPE=LIDAR_LAS "
        export_string += "ELEV_UNITS=METERS "
        export_string += "VERT_CS_CODE=NAVD88 "
        err_RunScript, handle, size = gm.RunScript(
            export_string, gm.GM_LoadFlags_UseDefaultLoadOpts, 0
        )
    else:
        # Exports generated feature layers as shapefiles
        # Ensures that the exported features keep their attributes
        err_ExportVector = gm.ExportVector(
            run_script_dir + layer_name,
            gm.GM_Export_Shapefile,
            layers_list[i],
            None,
            gm.GM_VectorExportFlags_HideProgress + gm.GM_VectorExportFlags_ExportAttrs,
            gm.NULL,
        )

Enhancing and Combining Lidar Queries

Note

Before running this sample, configure these variables:

• run_script_dir: the absolute path of the directory to be exported to

The code below imports pre-classified lidar data and separates any points classified as buildings and medium vegetation into two separate queries. Each query of points is then modified upon meeting specific criteria. Building-classified points that are over 70 meters tall were assigned a lime green color. Medium Vegetation-classified points that are over 60 meters tall were assigned an orange color. The two queries are then merged and exported as a new LAZ file. The changes made to point color can be seen in Global Mapper if “Color by RGB/Elev” is selected in the Lidar Draw Mode dropdown list. This can make points of specific criteria easier to view/interact with.

Learn how RunScript was used to run commands from the Global Mapper scripting language here.

# Query to hold building points
err_CreateEmptyLidarQuery1, query_handle1 = gm.CreateEmptyLidarQuery()
# Query to hold meduim vegetation points
err_CreateEmptyLidarQuery2, query_handle2 = gm.CreateEmptyLidarQuery()

# Iterates through all of the lidar points in a layer and adds the ones classified as
# buildings and medium vegetation to their respective queries.
num_pts = gm.GetLayerInfo(layers_list[0]).mNumLidarPoints
for i in range(num_pts):
    err_GetLidarPoint, point = gm.GetLidarPoint(
        layers_list[0], i, gm.GM_GetFeature_UseLayerProj, 0
    )
    if not err_GetLidarPoint:
        point_class = point.mClass
        # Get Building points
        if point_class == 6:
            # Criteria for a tall building (change point color)
            if point.mElevMeters >= 70:
                point.mRed = 130
                point.mGreen = 253
                point.Blue = 0
                # Using SetLidarPoint can be extremely slow if you are modifying very
                # large numbers of lidar points. (This sample uses a lidar file with ~92,000 points)
                err_SetLidarPoint1 = gm.SetLidarPoint(layers_list[0], i, point, 0)
            # Add point to query
            err_AddRemovePointToLidarQuery1 = gm.AddRemovePointToLidarQuery(
                query_handle1, layers_list[0], i, 1
            )
        # Get Medium Vegetation points
        elif point_class == 4:
            # Criteria for a tall tree (change point color)
            if point.mElevMeters >= 60:
                point.mRed = 252
                point.mGreen = 132
                point.Blue = 3
                err_SetLidarPoint2 = gm.SetLidarPoint(layers_list[0], i, point, 0)
            # Add point to query
            err_AddRemovePointToLidarQuery2 = gm.AddRemovePointToLidarQuery(
                query_handle2, layers_list[0], i, 1
            )
    else:
        print("Error: Couldn't retrieve lidar point")

# Merge the two queries together and create a layer based on the modified data
err_AddLidarQueryToQuery = gm.AddLidarQueryToQuery(query_handle1, query_handle2, 0)
err_CreateLayerFromLidarQuery, new_layer_handle = gm.CreateLayerFromLidarQuery(
    query_handle1, "class_layer.laz", None, 0
)

# Free memory
err_FreeLidarQuery1 = gm.FreeLidarQuery(query_handle1)
err_FreeLidarQuery2 = gm.FreeLidarQuery(query_handle2)

# Get handle for the loaded layers (This step is repeated to include the new layer)
arr_ptr, arr_size = gm.GetLoadedLayerList()
layers_array = gm.GM_LayerHandle_array_frompointer(arr_ptr)
layers_list = gm.carray_to_list(layers_array, arr_size)

# To view the modified colors for class_layer.laz in Global Mapper,
# make sure "Color by RGB/Elev" is selected in the lidar Draw Mode dropdown list
# Exports the lidar data, reflecting new class changes
for i in range(len(layers_list)):
    # Get the name of the layer to be exported
    layer_name = gm.GetLayerInfo(layers_list[i]).mDescription
    # Use the RunScript func to export lidar data as a .laz file
    export_string = "EXPORT_VECTOR "
    export_string += "FILENAME=" + run_script_dir + "new_" + layer_name + " "
    export_string += "EXPORT_LAYER=" + layer_name + " "
    export_string += "TYPE=LIDAR_LAS "
    export_string += "ELEV_UNITS=METERS "
    export_string += "VERT_CS_CODE=NAVD88 "
    err_RunScript, handle, size = gm.RunScript(
        export_string, gm.GM_LoadFlags_UseDefaultLoadOpts, 0
    )