Create, Process, and Export Digital Surface Model from Lidar Data

Read Aerial Lidar Data

Read 3-D point cloud data for an area near Tuscaloosa, Alabama from a LAZ file [1]. The area includes roads, trees, and buildings.

lazFileName = fullfile(toolboxdir("lidar"),"lidardata","las","aerialLidarData.laz");
lasReader = lasFileReader(lazFileName);
ptCloud = readPointCloud(lasReader);

Display the data.

figure
pcshow(ptCloud.Location)

Create DSM

A DSM includes the elevations of ground points, the elevations of natural features such as trees, and the elevations of artificial features such as buildings. Create a DSM from the point cloud data by using the pc2dem function. Use the maximum point from each element of the point cloud, which corresponds to the first return pulse of the lidar data, by specifying the CornerFillMethod as "max". The function returns an array of elevation values and the x- and y-limits of the data.

gridRes = 1;
[Z,xlimits,ylimits] = pc2dem(ptCloud,gridRes,CornerFillMethod="max");

Spatially Reference DSM

Spatially reference the DSM by creating a map reference object.

R = maprefpostings(xlimits,ylimits,size(Z))

R = 
  MapPostingsReference with properties:

             XWorldLimits: [429745.02 430146.02]
             YWorldLimits: [3679830.75 3680114.75]
               RasterSize: [285 402]
     RasterInterpretation: 'postings'
         ColumnsStartFrom: 'south'
            RowsStartFrom: 'west'
    SampleSpacingInWorldX: 1
    SampleSpacingInWorldY: 1
     RasterExtentInWorldX: 401
     RasterExtentInWorldY: 284
         XIntrinsicLimits: [1 402]
         YIntrinsicLimits: [1 285]
       TransformationType: 'rectilinear'
     CoordinateSystemType: 'planar'
             ProjectedCRS: []

The reference object contains information such as the limits, the distance between the points, and the directions of the columns and rows. By default, the reference object assumes that columns start from the south and rows start from the west. These default values are consistent with the output of the pc2dem function, which creates the elevation array such that the first element represents the southwesternmost point.

The ProjectedCRS property of the reference object is empty, which means the DSM is not associated with a projected CRS. Read the CRS from the LAZ file and update the ProjectedCRS property.

p = readCRS(lasReader);
R.ProjectedCRS = p;
disp(p)

  projcrs with properties:

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

A projected CRS consists of a geographic CRS and several parameters that are used to transform coordinates to and from the geographic CRS. A geographic CRS consists of a datum (including a reference ellipsoid), a prime meridian, and an angular unit of measurement. View the geographic CRS and its reference ellipsoid.

g = p.GeographicCRS

g = 
  geocrs with properties:

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

g.Spheroid

ans = 
referenceEllipsoid with defining properties:

                 Code: 7019
                 Name: 'GRS 1980'
           LengthUnit: 'meter'
        SemimajorAxis: 6378137
        SemiminorAxis: 6356752.31414036
    InverseFlattening: 298.257222101
         Eccentricity: 0.0818191910428158

  and additional properties:

    Flattening
    ThirdFlattening
    MeanRadius
    SurfaceArea
    Volume

Display the spatially referenced DSM as an overhead surface by using the mapshow function.

figure
mapshow(Z,R,DisplayType="surface")
axis image
title("Digital Surface Model (DSM) from Aerial Lidar Data")

Crop DSM to Region of Interest

Represent the DSM region as a polygon by using a mappolyshape object. Update the ProjectedCRS property to match the CRS of the DSM.

bboxx = xlimits([1 1 2 2 1]);
bboxy = ylimits([1 2 2 1 1]);
bboxshape = mappolyshape(bboxx,bboxy);
bboxshape.ProjectedCRS = p;

View the region using satellite imagery. You can visually confirm that the satellite imagery aligns with the DSM visualization created using the mapshow function.

regioncolors = lines(2);
geoplot(bboxshape, ...
    EdgeColor=regioncolors(1,:), ...
    FaceAlpha=0.2, ...
    LineWidth=2, ...
    DisplayName="Aerial Lidar Data Region")
hold on
geobasemap satellite
legend

Select and display a region of interest. To use a predefined region that is bounded by roads on the east, north, and west, specify interactivelySelectPoints as false. Alternatively, you can interactively select four points that define a region by specifying interactivelySelectPoints as true.

interactivelySelectPoints = false;
if interactivelySelectPoints
    [cropbboxlat,cropbboxlon] = ginput(4); %#ok<UNRCH> 
else
    cropbboxlat = [33.2571550; 33.2551982; 33.2551982; 33.2571125];
    cropbboxlon = [-87.7530648; -87.7530139; -87.7509086; -87.7509086];
end
cropbboxlat(end+1) = cropbboxlat(1);
cropbboxlon(end+1) = cropbboxlon(1);
cropbboxshape = geopolyshape(cropbboxlat,cropbboxlon);

geoplot(cropbboxshape, ...
    EdgeColor=regioncolors(2,:), ...
    FaceAlpha=0.2, ...
    LineWidth=2, ...
    DisplayName="Selected Region of Interest")

Transform the latitude and longitude limit coordinates for the region to x- and y-limit coordinates. The geographic CRS underlying the satellite basemap is WGS84, while the geographic CRS underlying the DSM data is NAD83. NAD83 and WGS84 are similar, but not identical. As a result, there can be discrepancies in coordinates between the satellite imagery and DSM.

[cropbboxx,cropbboxy] = projfwd(p,cropbboxlat(:),cropbboxlon(:));

Create the crop limits by finding the bounds of the x- and y-coordinates.

[cropxlimmin,cropxlimmax] = bounds(cropbboxx);
cropxlimits = [cropxlimmin cropxlimmax];
[cropylimmin,cropylimmax] = bounds(cropbboxy);
cropylimits = [cropylimmin cropylimmax];

Create a new spatially referenced DSM that contains data within the region of interest.

[Zcrop,Rcrop] = mapcrop(Z,R,cropxlimits,cropylimits);

Export DSM to GeoTIFF File

Write the cropped DSM to a GeoTIFF file called lidardsm.tif. Specify the projected CRS by using the CoordRefSysCode argument. The metadata for the LAZ file [1] indicates the projected CRS is UTM Zone 16N, specified by EPSG authority code 26916.

datafile = "lidardsm.tif";
epsgCode = 26916;
geotiffwrite(datafile,Zcrop,Rcrop,CoordRefSysCode=epsgCode)

You can also find the authority code by displaying the well-known text (WKT) string for the projected CRS. For this WKT, the authority code is in the last line.

wktstring(p,"Format","formatted")

ans = 
    "PROJCRS["NAD83 / UTM zone 16N",
         BASEGEOGCRS["NAD83",
             DATUM["North American Datum 1983",
                 ELLIPSOID["GRS 1980",6378137,298.257222101,
                     LENGTHUNIT["metre",1]]],
             PRIMEM["Greenwich",0,
                 ANGLEUNIT["degree",0.0174532925199433]],
             ID["EPSG",4269]],
         CONVERSION["UTM zone 16N",
             METHOD["Transverse Mercator",
                 ID["EPSG",9807]],
             PARAMETER["Latitude of natural origin",0,
                 ANGLEUNIT["degree",0.0174532925199433],
                 ID["EPSG",8801]],
             PARAMETER["Longitude of natural origin",-87,
                 ANGLEUNIT["degree",0.0174532925199433],
                 ID["EPSG",8802]],
             PARAMETER["Scale factor at natural origin",0.9996,
                 SCALEUNIT["unity",1],
                 ID["EPSG",8805]],
             PARAMETER["False easting",500000,
                 LENGTHUNIT["metre",1],
                 ID["EPSG",8806]],
             PARAMETER["False northing",0,
                 LENGTHUNIT["metre",1],
                 ID["EPSG",8807]]],
         CS[Cartesian,2],
             AXIS["easting",east,
                 ORDER[1],
                 LENGTHUNIT["metre",1]],
             AXIS["northing",north,
                 ORDER[2],
                 LENGTHUNIT["metre",1]],
         ID["EPSG",26916]]"

One way to validate the GeoTIFF file is to return information about the file as a RasterInfo object. For example, verify that the projected CRS is in the file by querying the CoordinateReferenceSystem property of the RasterInfo object.

info = georasterinfo(datafile);
info.CoordinateReferenceSystem

ans = 
  projcrs with properties:

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

Another way to validate the GeoTIFF file is by displaying it. Read the new DSM as an array and a reference object by using the readgeoraster function. Then, display the DSM.

[Z2,R2] = readgeoraster(datafile);

figure
mapshow(Z2,R2,DisplayType="surface")
axis image
title("Cropped DSM from Aerial Lidar Data")

You can use the GeoTIFF file in other applications that import GIS data. For example, RoadRunner enables you to add elevation data from GeoTIFF files to scenes.

References

[1] OpenTopography. “Tuscaloosa, AL: Seasonal Inundation Dynamics And Invertebrate Communities,” 2011. https://doi.org/10.5069/G9SF2T3K.