Create, Process, and Export Digital Surface Model from Lidar Data
This example shows how to process aerial lidar data received from an airborne lidar system into a GeoTIFF file. Import a LAZ file containing aerial lidar data, create a spatially referenced digital surface model (DSM) from the data, crop the DSM to an area of interest, and export the cropped DSM to a GeoTIFF file.
When you export a DSM to a GeoTIFF file, you also export the projected coordinate reference system (CRS) for the data. Projected CRSs associate x- and y-coordinates to locations on Earth. Specifying the projected CRS is important when creating a model because the same coordinates in different projected CRSs can refer to different locations.
Read Aerial Lidar Data
Read 3-D point cloud data for an area near Tuscaloosa, Alabama from a LAZ file . The area includes roads, trees, and buildings.
lazFileName = fullfile(toolboxdir("lidar"),"lidardata","las","aerialLidarData.laz"); lasReader = lasFileReader(lazFileName); ptCloud = readPointCloud(lasReader);
Display the data.
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
"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.
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
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"
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
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
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
false. Alternatively, you can interactively select four points that define a region by specifying
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  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.
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, LENGTHUNIT["metre",1]], AXIS["northing",north, ORDER, 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
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.
 OpenTopography. “Tuscaloosa, AL: Seasonal Inundation Dynamics And Invertebrate Communities,” 2011. https://doi.org/10.5069/G9SF2T3K.
pc2dem(Lidar Toolbox) |