Visualize Point Clouds On Maps Using Coordinate Reference System From LAS/LAZ Files

Since R2025a

This example shows how to use coordinate reference system (CRS) information in LAS/LAZ files to transform a point cloud from one coordinate reference system to another. Converting two point clouds collected from the same location but stored in different reference systems is important for data analysis, combining multiple point clouds and comparison.

A CRS defines point cloud coordinates in real-world locations. The types of CRS relevant for this example are:

1. Projected: Assigns cartesian x and y map coordinates to physical locations.

2. Geographic: Assigns latitude, longitude to physical locations.

A projected CRS consists of a geographic CRS and several parameters that are used to transform coordinates to and from the geographic coordinates. Point clouds can also be transformed from one projected CRS to another. To learn more about CRS, see Transform Coordinates to a Different Projected CRS and Project and Display Raster Data examples.

A LAS/LAZ file contains CRS information, specified by American Society for Photogrammetry and Remote Sensing (ASPRS), enabling to transform coordinates of point clouds stored in them.

Load Point Cloud

Read point cloud from a LAZ file for an area near Tuscaloosa, Alabama [1].

% Load point cloud from a LAZ file
fileName = fullfile(toolboxdir("lidar"),"lidardata","las","aerialLidarData.laz");
lasReader = lasFileReader(fileName);
ptCloud = readPointCloud(lasReader);

% Visualize
pcviewer(ptCloud);

Read Coordinate System Information

Read projected CRS information from the LAZ file using the readCRS function. A projected CRS consists of a geographic CRS and several parameters that are used to transform coordinates between Cartesian and geographic coordinates.

crs = readCRS(lasReader)

crs = 
  projcrs with properties:

                    Name: "NAD83 / UTM zone 16N"
           GeographicCRS: [1×1 geocrs]
        ProjectionMethod: "Transverse Mercator"
              LengthUnit: "meter"
    ProjectionParameters: [1×1 map.crs.ProjectionParameters]

The reference system is identified by the name "NAD83 / UTM zone 16N" and covers parts of North America, including state of Alabama. Regions are typically associated with multiple reference systems, which are published and managed by organizations such as the EPSG.

Visualize Area of Data Acquisition

CRS information can be used to visualize on a map, the area from which the point cloud was acquired. This is helpful for performing geospatial analysis and working with other data sources for the same area, such as, satellite imagery.

% Define an area of interest (AOI)
xLimits = ptCloud.XLimits;
yLimits = ptCloud.YLimits;
aoi = aoiquad(xLimits,yLimits,crs);

% Visualize AOI overlaid on a map
figure
geobasemap satellite
geoplot(aoi)

Visualize Point Cloud on 3-D Globe

CRS information can be used to plot a point cloud on a geographic globe. This map visualization aids in verifying data acquisition and is also useful for assessing map-building algorithms such as simultaneous localization and mapping (SLAM).

% Convert x-y coordinates to latitude-longitude.
[lat,lon] = projinv(crs,ptCloud.Location(:,1),ptCloud.Location(:,2));

% Use z-coordinate as height
height = ptCloud.Location(:,3);

% Visualize
fig = uifigure;
globe = geoglobe(fig);
gplot = geoplot3(globe,lat,lon,height,"o",MarkerSize=2);

% Adjust graphics camera position to bring it closer to the 
% point cloud for better visibility.
campos(globe, 33.25475, -87.75359,  89.24658)
camheight(globe, 89.2466)
camheading(globe, 356.3342)
campitch(globe, -17.8413)
camroll(globe, 360)

The tree indicated by the arrow in the image above is represented in point cloud on the map in the image below. To tilt and rotate the globe, hold the Ctrl key and drag.

Geographic globes displayed using the geoglobe function appear in a separate window. To display a screenshot of the geographic globe, use the following code.

im = getframe(fig);
figure
imshow(im.cdata,Border="tight")

Working with Geographic CRS

Point clouds can also be stored in geographic coordinates in a LAS/LAZ file. In this case, the readCRS function returns a geocrs (Mapping Toolbox) object. Read a new point cloud from the same location in Tuscaloosa, Alabama, stored in geographic coordinates (latitude and longitude).

% Load point cloud
geolasreader = lasFileReader("aerialLidarDataGeo.laz");

Each point in this point cloud represents the same geographic location as the points in aerialLidarData.laz. The Location property of the point cloud object stores the latitudes, longitudes and heights in three dimensions, respectively. The point cloud is in geographic coordinates, so the readPointCloud function issues a warning.

% Display point cloud. 
geoPtCloud = readPointCloud(geolasreader)

geoPtCloud = 
  pointCloud with properties:

     Location: [1018047×3 double]
        Count: 1018047
      XLimits: [-87.7543 -87.7499]
      YLimits: [33.2552 33.2578]
      ZLimits: [72.7900 125.8200]
        Color: []
       Normal: []
    Intensity: [1018047×1 uint16]

% Read geographic CRS.
crsGeo = readCRS(geolasreader)

crsGeo = 
  geocrs with properties:

             Name: "NAD83"
            Datum: "North American Datum 1983"
         Spheroid: [1×1 referenceEllipsoid]
    PrimeMeridian: 0
        AngleUnit: "degree"

The geographic CRS here is the North American Datum of 1983,as identified by the name "NAD83" in the object. Another example of a geographic CRS commonly found in LAS/LAZ file is the World Geodetic System of 1984 (WGS84), a global reference system used by GPS and widely adopted in mapping and navigation.

Incorrect Visualization with Geographic Coordinates

Point cloud functions and applications require Cartesian coordinates to visualize and process point clouds. As demonstrated by this figure, point clouds in geographic coordinates are visualized incorrectly.

figure
pcshow(geoPtCloud)
title("Incorrect Visualization")

Convert to Cartesian Coordinates

To visualize or process this point cloud, convert the geographic coordinates to Cartesian coordinates. Projected CRS information is necessary for this conversion. You can create a projected CRS object in one of the following ways:

• Using a known code specified by an authority like EPSG

• Using a well-known text (WKT) string.

• From a corresponding data source, for example, satellite imagery or vector data stored in a file such as shapefile or KML.

• Find a code by searching, browsing, and exploring the Spatial Reference home page. Use geographic region information to narrow down the search. If there are multiple choices, codes with a "UTM Zone" identifier are typically reasonable choices.

For this example, EPSG:26916 was selected based on the name and the geographic region of data collection, using the Spatial Reference home page. Create a projcrs (Mapping Toolbox) object.

proj = projcrs(26916)

proj = 
  projcrs with properties:

                    Name: "NAD83 / UTM zone 16N"
           GeographicCRS: [1×1 geocrs]
        ProjectionMethod: "Transverse Mercator"
              LengthUnit: "meter"
    ProjectionParameters: [1×1 map.crs.ProjectionParameters]

Verify that the geographic CRS in the object matches the CRS obtained from LAZ file. If the CRSs do not match, then the projected Cartesian coordinates might be misaligned from the geographic locations they represent.

isequal(proj.GeographicCRS,crsGeo)

ans = logical
   1

Project latitude-longitude coordinates to x-y coordinates using projfwd (Mapping Toolbox) function, and use the height from the point cloud object as the z-coordinate.

[x,y] = projfwd(proj, geoPtCloud.Location(:,2), geoPtCloud.Location(:,1));
z = geoPtCloud.Location(:,3);
tformedPtCloud = pointCloud([x y z]);

Visualize the point cloud converted to Cartesian coordinates.

pcviewer(tformedPtCloud);

To verify the conversion from geographic coordinates to Cartesian, you can compare the point cloud from aerialLidarData.laz which is stored in Cartesian coordinates and tformedPtCloud, which is a converted point cloud from geographic to Cartesian coordinates. Visualizing the two point clouds using pcshowpair, shows that there is no difference between them, verifying the accuracy of the conversion. If there was a difference, the point clouds would be misaligned.

figure
pcshowpair(ptCloud, tformedPtCloud)

References

[1] OpenTopography. Tuscaloosa, AL: Seasonal Inundation Dynamics And Invertebrate Communities. OpenTopography, 2011. DOI.org (Datacite), https://doi.org/10.5069/G9SF2T3K