Surface Analysis Samples

Surface analysis involves examining terrain data to extract geographic information, such as slope, aspect, and elevation. Using digital elevation models (DEMs), watersheds, viewshed maps, contour lines, and elevation grids can be generated. These tools can help better understand landscape characteristics and assess the environmental impact of enironmental changes like natural disasters.

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

Elevation Grids from Multiple Lidar Files

Note

Before running this sample, configure these variables:

• import_dir: the directory containing sample data

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

A sample that iterates through multiple lidar files, generates elevation grids for each one, determines the appropriate dimensions for the exported image, then exports the elevation grid as a JPEG. This demonstrates how multiple files in a directory can be easily accessed and modified.

# Retrieve all lidar files in the directory
filenames = []
for f in dir_files:
    file = f.lower()
    if file.endswith(".las") or file.endswith(".laz"):
        filenames.append(file)

# Modify the display options
vert_display_opts = gm.GM_VerticalDisplayOptions_t()
vert_display_opts.mSize = 50
vert_display_opts.mShaderName = "Atlas Shader"
err_SetVerticalDisplayOptions = gm.SetVerticalDisplayOptions(vert_display_opts)

# Loop through all the lidar layers
for file in filenames:
    err_LoadLayerList, arr_ptr, arr_size = gm.LoadLayerList(
        import_dir + file, gm.GM_LoadFlags_UseDefaultLoadOpts
    )
    layers_arr = gm.GM_LayerHandle_array_frompointer(arr_ptr)

    # Generate an elevation grid with all the loaded layers
    grid_setup = gm.GM_GridGenSetup_t()
    grid_setup.mGridAlg = gm.GM_GridAlg_BinAverage
    grid_setup.mElevUnits = gm.GM_ElevUnit_Meters
    err_GenerateElevationGrid2, grid_layer = gm.GenerateElevationGrid2(
        [layers_arr[0]], grid_setup
    )

    # Get the right dimensions for the image
    bounds = gm.GetLayerInfo(grid_layer).mGlobalRect
    ratio = abs(bounds.mMaxX - bounds.mMinX) / abs(bounds.mMaxY - bounds.mMinY)
    width = 1920
    height = width * (1 / ratio)
    # Make sure the exported layers are at least HD
    if height < 1080: 
        height = 1080
        width = height * ratio
    width = int(width)
    height = int(height)

    # Export the raster image, maintaining file names
    export_name = file.split('.')
    err_ExportRaster = gm.ExportRaster(
        export_dir + export_name[0] + "_elevation_grid.jpg",
        gm.GM_Export_JPG,
        grid_layer,
        None,
        width,
        height,
        0x0,
    )
    err_CloseLayer1 = gm.CloseLayer(layers_arr[0])
    err_CloseLayer2 = gm.CloseLayer(grid_layer)

Contour maps from Multiple USGS DEMs

Note

Before running this sample, configure these variables:

• import_dir: the directory containing sample data

• DXF_dir: the directory for DXF files to be exported to

• SHP_dir: the directory for SHP files to be exported to

This script generates a vector contour map based on elevation data and exports the contours into multiple vector formats. All of the data is retrieved from a provided directory, filtering only the desired extensions. The DEMs are then loaded and contour lines are generated at intervals of 50ft. The resulting contours are then exported as DXFs and shapefiles.

# Retrieve all files in the directory with the desired extension(s)
filenames = []
for f in dir_files:
    file = f.lower()
    if (
        file.endswith(".dem")
        or file.endswith(".stds")
        or file.endswith(".tar")
        or file.endswith(".gz")
    ):
        filenames.append(file)

# Loop over and load in all the found files
for filename in filenames:
    err_LoadLayerList, arr_ptr, arr_size = gm.LoadLayerList(
        import_dir + filename, gm.GM_LoadFlags_UseDefaultLoadOpts
    )
    layer_arr = gm.GM_LayerHandle_array_frompointer(arr_ptr)
    layer = layer_arr[0]
    # Generate the contour lines at intervals of 50 feet
    contour_params = gm.GM_ContourParams_t()
    contour_params.mIntervalInFeet = 50
    contour_params.mShowProgress = False
    err_GenerateContours, contour_layer = gm.GenerateContours(layer, contour_params)

    # Export the generate contours file to this path
    filename_wo_ext = filename[: filename.index(".")]
    err_ExportVector1 = gm.ExportVector(
        DXF_dir + filename_wo_ext + "_contours.dxf",
        gm.GM_Export_DXF,
        contour_layer,
        None,
        gm.GM_VectorExportFlags_HideProgress + gm.GM_VectorExportFlags_GenPRJFile,
        0x0,
    )
    err_ExportVector2 = gm.ExportVector(
        SHP_dir + filename_wo_ext + "_contours.shp",
        gm.GM_Export_Shapefile,
        contour_layer,
        None,
        gm.GM_VectorExportFlags_HideProgress
        + gm.GM_VectorExportFlags_GenPRJFile
        + gm.GM_VectorExportFlags_ExportLines
        + gm.GM_VectorExportFlags_Export3D,
        0x0,
    )

    # Close the layers after we're done with them
    err_CloseLayer1 = gm.CloseLayer(layer)
    err_CloseLayer2 = gm.CloseLayer(contour_layer)

Generating a Viewshed and Watershed

Note

Before running this sample, configure these variables:

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

A script that generates a viewshed centered on a specified (X,Y) coordinate and creates a watershed based on the provided GMG file. It calculates key watershed attributes, such as drainage area size, inflow, and outflow values. The resulting viewshed, stream network, and drainage areas are then exported to the specified export directory.

# Generate a viewshed at a specific coordinate-location
center_point = gm.GM_Point_t(435986.9, 4893585.6)
viewshed_opts = gm.GM_ViewShedParams_t()
viewshed_opts.mAtmosphericCorr = 1.33333
viewshed_opts.mCenterPoint = center_point
viewshed_opts.mColor = gm.toRGB(255, 255, 255)
viewshed_opts.mDesc = "Viewshed Analysis"
viewshed_opts.mGenCoveragePolys = True
viewshed_opts.mRadius = 3000
viewshed_opts.mTransmitterHeight = 0
viewshed_opts.mReceiverHeight = 2
err_CalcViewShed, viewshed_handle = gm.CalcViewShed(layers_list[0], viewshed_opts)

# Generate a watershed
watershed_param = gm.GM_WatershedParams_t()
watershed_param.mFlags = gm.GM_Watershed_CreateAreas
watershed_param.mFlags += gm.GM_Watershed_KeepZeroAtZero
watershed_param.mFlags += gm.GM_Watershed_SmoothStreams
watershed_param.mDesc = "DRAINAGE NETWORK"
watershed_param.mMaxDepressionDepth = 12
watershed_param.mStreamThreshold = 10000
watershed_param.mBounds = gm.GetLayerInfo(layers_list[0]).mGlobalRect
watershed_layer = [layers_list[0]]
err_GenerateWatershed, stream_layer, drainage_area_layer = gm.GenerateWatershed(
    watershed_layer, watershed_param
)

# Gather watershed data
watershed_data = []
i = 0
stream_feature = gm.GetLineFeature(stream_layer, i)
while stream_feature:
    # Get the STREAM_ID, IN_FLOW, OUT_FLOW, and DRAIN_AREA
    attr_array = gm.GM_AttrValue_array_frompointer(
        stream_feature.mFeatureInfo.mAttrList
    )
    attribute_values = gm.carray_to_list(
        attr_array, stream_feature.mFeatureInfo.mNumAttrs
    )

    # Create a dictionary to hold the stream data
    stream_data = {
        attribute_values[0].mName: attribute_values[0].mVal,
        attribute_values[1].mName: attribute_values[1].mVal,
        attribute_values[2].mName: attribute_values[2].mVal,
        attribute_values[3].mName: attribute_values[3].mVal,
    }

    watershed_data.append(stream_data)

    # Frees deep copies from memory
    # Without this step the script causes the program to freeze
    gm.FreeLineFeature(stream_feature)
    i += 1
    stream_feature = gm.GetLineFeature(stream_layer, i)

# Display watershed data
for data in watershed_data:
    print("STREAM_ID", data["STREAM_ID"] + ":")
    print("  IN_FLOW:", data["IN_FLOW"])
    print("  OUT_FLOW:", data["OUT_FLOW"])
    print("  DRAIN_AREA:", data["DRAIN_AREA"])
    print("-" * 30)

# Export all created layers
err_ExportVector1 = gm.ExportVector(
    export_dir + "view_shed.kmz",
    gm.GM_Export_KML,
    viewshed_handle,
    None,
    gm.GM_ExportFlags_HideProgress,
    gm.NULL,
)
err_ExportVector2 = gm.ExportVector(
    export_dir + "streams.shp",
    gm.GM_Export_Shapefile,
    stream_layer,
    None,
    gm.GM_ExportFlags_HideProgress
    + gm.GM_ExportFlags_GenPRJFile
    + gm.GM_ExportFlags_GenWorldFile,
    gm.NULL,
)
err_ExportVector3 = gm.ExportVector(
    export_dir + "drainage_areas.shp",
    gm.GM_Export_Shapefile,
    drainage_area_layer,
    None,
    gm.GM_ExportFlags_HideProgress
    + gm.GM_ExportFlags_GenPRJFile
    + gm.GM_ExportFlags_GenWorldFile,
    gm.NULL,
)

err_CloseLayer1 = gm.CloseLayer(viewshed_handle)
err_CloseLayer2 = gm.CloseLayer(stream_layer)
err_CloseLayer3 = gm.CloseLayer(drainage_area_layer)