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 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");
pc2dem function returns values of type
single, but many functions require inputs to be of type
double. Convert the values to
Z = double(Z); xlimits = double(xlimits); ylimits = double(ylimits);
Spatially reference the DSM by creating a map reference object.
R = maprefpostings(xlimits,ylimits,size(Z))
R = MapPostingsReference with properties: XWorldLimits: [429745.03125 430146.03125] 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. The metadata for the LAZ file  indicates the projected CRS is UTM Zone 16N, specified by EPSG authority code
26916. Create a projected CRS object.
epsgCode = 26916; p = projcrs(epsgCode)
p = projcrs with properties: Name: "NAD83 / UTM zone 16N" GeographicCRS: [1×1 geocrs] ProjectionMethod: "Transverse Mercator" LengthUnit: "meter" ProjectionParameters: [1×1 map.crs.ProjectionParameters]
ProjectedCRS property of the reference object.
R.ProjectedCRS = p;
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")
Transform the projected x- and y-limit coordinates into geographic latitude and longitude limit coordinates.
[latlim,lonlim] = projinv(p,xlimits,ylimits);
View the region using satellite imagery. You can visually confirm that the satellite imagery aligns with the DSM visualization created using the
bboxlat = latlim([1 1 2 2 1]); bboxlon = lonlim([1 2 2 1 1]); geoplot(bboxlat,bboxlon,"r", ... LineWidth=2, ... DisplayName="Aerial Lidar Data Region") hold on geobasemap satellite legend
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 between the satellite imagery and DSM visualization.
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 geoplot(cropbboxlat([1:4 1]),cropbboxlon([1:4 1]),"y", ... LineWidth=2, ... DisplayName="Selected Region of Interest")
Transform the latitude and longitude limit coordinates for the region to x- and y-limit coordinates.
[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);
Write the cropped DSM to a GeoTIFF file called
lidardsm.tif. Specify the projected CRS by using the
datafile = "lidardsm.tif"; geotiffwrite(datafile,Zcrop,Rcrop,CoordRefSysCode=epsgCode)
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.