Georeference Aerial Point Cloud for Scene Generation

This example uses:

This example shows how to georeference an aerial point cloud by using the coordinate reference system (CRS) information. This example also uses the georeferenced point cloud to add elevation information to maps such as an OpenStreetMap® or RoadRunner HD Map.

You can use a combination of sensors to create a virtual scene that includes both a road network and static objects. For example, you can use a camera and GPS data to reconstruct a road network, while, using lidar data to recreate static scene objects such as buildings and trees. Alternatively, you can use the OpenStreetMap® live map service to reconstruct road networks. However you create your virtual road network, you must ensure that the data from various sensors and maps align within a consistent spatial coordinate system. A georeferenced point cloud ensures that the point locations within the point cloud correspond accurately with maps, such as an OpenStreetMap or RoadRunner HD Map.

You can mount lidar sensor on an aerial vehicle to record aerial lidar data. In this example, you georeference a point cloud obtained from aerial lidar data. The aerial lidar point cloud contains CRS information, which you can use to georeference the point cloud.

In this example, you:

Load and visualize an aerial point cloud
Georeference the aerial point cloud
Validate the georeferenced point cloud

This example requires the Scenario Builder for Automated Driving Toolbox™ support package. Check if the support package is installed and, if not, install it using the Add-On Manager. For more information on installing support packages, see Get and Manage Add-Ons.

checkIfScenarioBuilderIsInstalled

Load and Visualize Data

This example uses aerial point cloud data downloaded from the U.S. Geological Survey (USGS) 3DEP LidarExplorer [1]. This is an online platform that provides access to recorded lidar data in LAZ format for various regions of the United States. To download lidar data, in the 3DEP LidarExplorer, hold Ctrl and drag to draw a rectangular area of interest on the map. The USGS Lidar Explorer shows a list of available lidar data in the selected region. Select a data set from the list, and download lidar point cloud as a LAZ file.

For this example, download a ZIP file containing lidar data from the 3DEP LidarExplorer, and then unzip the file. Load the point cloud data into the workspace using the readPointCloud (Lidar Toolbox) function of the lasFileReader (Lidar Toolbox) object, and visualize the point cloud.

dataFolder = tempdir;
lidarDataFileName = "USGSLidarExplorerData.zip"; 
url = "https://ssd.mathworks.com/supportfiles/driving/data/" + lidarDataFileName; 
filePath = fullfile(dataFolder,lidarDataFileName); 
if ~isfile(filePath)
    websave(filePath,url);
end
unzip(filePath,dataFolder)

Warning: Cannot overwrite file "C:\Users\rampatro\AppData\Local\Temp\USGS_LPC_CA_NoCAL_3DEP_Supp_Funding_2018_D18_w2276n1958.laz". The file is already open.

% Specify the path of the downloaded LAZ file.
lasfile = fullfile(dataFolder,"USGS_LPC_CA_NoCAL_3DEP_Supp_Funding_2018_D18_w2276n1958.laz");

% Create a lasFileReader object.
lasReader = lasFileReader(lasfile);

% Read the point cloud.
pc = readPointCloud(lasReader);

Visualize the central region of the point cloud using the helperViewPointCloudCenterRegion helper function. To visualize the complete point cloud, you can use the pcshow function.

Observe that the x- and y-locations of the point cloud are in the projected coordinate system.

helperViewPointCloudCenterRegion(pc)
title("USGS Aerial Lidar Data: Central region")

Display the range of locations along the x-axis.

xLimits = pc.XLimits

xLimits = 1×2
10⁶ ×

   -2.2760   -2.2750

Display the range of locations along the y-axis.

yLimits = pc.YLimits

yLimits = 1×2
10⁶ ×

    1.9580    1.9590

Georeference Aerial Point Cloud

Georeference the aerial point cloud data by using the helperGeoreferencePointCloud helper function, attached to this example as a supporting file. The function takes a lasFileReader (Lidar Toolbox) object as input and returns a georeferenced point cloud along with its georeference point in the form [latitude longitude altitude]. By default, the function sets the georeference as the center of the point cloud. If the input lasFileReader (Lidar Toolbox) object does not contain CRS information, you can specify an EPSG code, which defines the coordinate reference system of the point cloud. The aerial point cloud data used in this example has a projected CRS EPSG code of 6350. For information on valid EPSG codes, see the EPSG home page.

if hasCRSData(lasReader)
    % Obtain a georeferenced point cloud along with its georeference point.
    [georeferencedPointCloud,georeference] = helperGeoreferencePointCloud(lasReader);

else
    % If the LAZ file does not contain CRS data, specify an EPSG code which defines the CRS of the point cloud.
    code = '6350';
    
    % Obtain a georeferenced point cloud along with its georeference point.
    [georeferencedPointCloud, georeference] = helperGeoreferencePointCloud(lasReader,code);
end

Visualize the central region of the point cloud by using the helperViewPointCloudCenterRegion helper function.

Observe that the xy-coordinates of the point cloud are in the local coordinate system, whose origin is at the center of the point cloud, specified by the returned georeference location georeference.

helperViewPointCloudCenterRegion(georeferencedPointCloud)
title("Georeferenced Point Cloud")

Display the range of locations along the x-axis.

xLimits = georeferencedPointCloud.XLimits

xLimits = 1×2

 -621.3181  626.3163

Display the range of locations along the y-axis.

yLimits = georeferencedPointCloud.YLimits

yLimits = 1×2

 -611.2522  612.3971

Validate Georeferenced Point Cloud

To validate the georeferenced point cloud, you must verify whether the locations of points in the georeferenced point cloud align with the locations of road boundaries in the associated OpenStreetMap (OSM) [2].

Import the OSM road network by using the getMapROI function. Use the georeference value of the point cloud to specify the latitude and longitude values for the OSM road network, and the mean of the limits of the georeferenced point cloud to specify the extent of the road network. To get the road network in the RoadRunner HD Map format, import the roads from the OSM web service into a driving scenario and use the getRoadRunnerHDMap function.

osmExtent = mean(abs([xLimits yLimits]));
mapParameters = getMapROI(georeference(1),georeference(2),Extent=osmExtent);
osmFileName = websave("RoadScene.osm",mapParameters.osmUrl,weboptions(ContentType="xml"));
scenario = drivingScenario;
roadNetwork(scenario,OpenStreetMap=osmFileName)

% Get the roadrunnerHDMap object from the drivingScenario.
rrMap = getRoadRunnerHDMap(scenario);

View the georeferenced point cloud and the overlaid RoadRunner HD Map. Observe that the xy-locations of the georeferenced point cloud align with the xy-locations of the road network in the RoadRunner HD Map. However, the OSM might not contain accurate elevation information, because of which the z-locations of the georeferenced point cloud might not align with the z-locations of the road network.

ax = plot(rrMap);
hold on
helperViewPointCloudCenterRegion(georeferencedPointCloud,ax)
title("Georeferenced Point Cloud Overlaid on OSM")
hold off

To refine the elevation in the map using the elevation information from the georeferenced point cloud, use the addElevation function.

elevatedSceneMap = addElevation(rrMap,georeferencedPointCloud);

View the georeferenced point cloud and the overlaid, elevation-adjusted RoadRunner HD Map. Observe that the z-locations of the road network in the elevated RoadRunner HD Map now align with the z-locations of the georeferenced point cloud.

This validates that the point cloud is georeferenced, and that the locations of the points in the georeferenced point cloud align with the locations of the boundaries in the road network.

ax2 = plot(elevatedSceneMap);
hold on
helperViewPointCloudCenterRegion(georeferencedPointCloud,ax2)
hold off

Summary

In this example, you learned how to georeference a point cloud obtained from an aerial lidar sensor by using CRS information from the EPSG. You also used the georeferenced point cloud to add elevation information to maps such as an OpenStreetMap or RoadRunner HD Map.

You can also extract static objects from a georeferenced point cloud to create a scene that spatially aligns these objects with a road network that has been reconstructed using the OSM web service.

References

[1] 3DEP LidarExplorer courtesy of the U.S. Geological Survey https://apps.nationalmap.gov/lidar-explorer/.

[2] You can download OpenStreetMap files from https://www.openstreetmap.org, which provides access to crowd-sourced map data all over the world. The data is licensed under the Open Data Commons Open Database License (ODbL), https://opendatacommons.org/licenses/odbl/.

Helper Function

helperViewPointCloudCenterRegion — Visualizes the central region of a registered point cloud.

function helperViewPointCloudCenterRegion(varargin)
    ptCld = varargin{1};
    if nargin == 1
        ax = figure;
    else
        ax = varargin{2};
    end

    % Define the ROI.
    cent = [sum(ptCld.XLimits)/2 sum(ptCld.YLimits)/2];
    roi = [cent(1)-250 cent(1)+200 cent(2)-250 cent(2)+200 ptCld.ZLimits];
    % Visualize the point cloud.
    pcshow(ptCld)
    xlim(roi(1:2))
    ylim(roi(3:4))
    zlim(roi(5:6))
    xlabel("X")
    ylabel("Y")
    zlabel("Z")
end