Считайте данные о лидаре антенны

Считайте 3-D данные об облаке точек для области под Таскалусой, Алабама из файла [1] LAZ. Область включает дороги, деревья и создания.

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

Отобразите данные.

figure
pcshow(ptCloud.Location)

Создайте DSM

DSM включает вертикальные изменения наземных точек, вертикальные изменения природных объектов, такие как деревья и вертикальные изменения искусственных функций, такие как создания. Создайте DSM из данных об облаке точек при помощи pc2dem функция. Используйте максимальную точку от каждого элемента облака точек, которое соответствует первому импульсу возврата данных о лидаре путем определения CornerFillMethod как "max". Функция возвращает массив значений вертикального изменения и x-и y-пределов данных.

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

pc2dem функция возвращает значения типа single, но много функций требуют, чтобы входные параметры имели тип double. Преобразуйте значения в double.

Z = double(Z);
xlimits = double(xlimits);
ylimits = double(ylimits);

Пространственно ссылочный DSM

Пространственно сошлитесь на DSM путем создания объекта ссылки карты.

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: []

Ссылочный объект содержит информацию, такую как пределы, расстояние между точками и направления столбцов и строк. По умолчанию ссылочный объект принимает, что столбцы запускаются с юга, и строки запускаются с запада. Эти значения по умолчанию сопоставимы с выходом pc2dem функция, которая создает вертикальное изменение, выстраивает таким образом, что первый элемент представляет точку southwesternmost.

ProjectedCRS свойство ссылочного объекта пусто, что означает, что DSM не сопоставлен со спроектированным CRS. Метаданные для файла [1] LAZ указывают, что спроектированный CRS является Зоной UTM 16 Н, заданных кодом полномочий EPSG 26916. Создайте спроектированный объект CRS.

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 свойство ссылочного объекта.

R.ProjectedCRS = p;

Спроектированный CRS состоит из географического CRS и нескольких параметров, которые используются, чтобы преобразовать координаты к и от географического CRS. Географический CRS состоит из данной величины (включая опорный эллипсоид), нулевой меридиан и угловая единица измерения. Просмотрите географический CRS и его опорный эллипсоид.

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

Отобразите DSM, на который пространственно ссылаются, как служебную поверхность при помощи mapshow функция.

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

Обрежьте DSM к необходимой области

Преобразуйте спроектированный x-и y-предельные координаты в географические предельные координаты широты и долготы.

[latlim,lonlim] = projinv(p,xlimits,ylimits);

Просмотрите область с помощью спутниковых снимков. Можно визуально подтвердить, что спутниковые снимки выравниваются с визуализацией DSM, созданной с помощью mapshow функция.

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

Географический CRS, лежащий в основе satellite основная карта является WGS84, в то время как географический CRS, лежащий в основе данных DSM, является NAD83. NAD83 и WGS84 подобны, но не идентичны. В результате могут быть несоответствия между спутниковыми снимками и визуализацией DSM.

Выберите и отобразите необходимую область. Чтобы использовать предопределенную область, которая ограничена дорогами на востоке, север и запад, задают interactivelySelectPoints как false. В качестве альтернативы можно в интерактивном режиме выбрать четыре точки, которые задают область путем определения interactivelySelectPoints как 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
geoplot(cropbboxlat([1:4 1]),cropbboxlon([1:4 1]),"y", ...
    LineWidth=2, ...
    DisplayName="Selected Region of Interest")

Преобразуйте предельные координаты широты и долготы для области к y-предельные координаты и x-.

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

Создайте пределы обрезки путем нахождения границ x-и y-координат.

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

Создайте новое, пространственно сослался на DSM, который содержит данные в необходимой области.

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

Экспортируйте DSM в файл GeoTIFF

Запишите обрезанный DSM в файл GeoTIFF под названием lidardsm.tif. Задайте спроектированный CRS при помощи CoordRefSysCode аргумент.

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

Один способ подтвердить файл GeoTIFF состоит в том, чтобы возвратить информацию о файле как RasterInfo объект. Например, проверьте, что спроектированный CRS находится в файле путем запроса CoordinateReferenceSystem свойство RasterInfo объект.

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]

Другой способ подтвердить файл GeoTIFF путем отображения его. Считайте новый DSM как массив и ссылочный объект при помощи readgeoraster функция. Затем отобразите DSM.

[Z2,R2] = readgeoraster(datafile);

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

Можно использовать файл GeoTIFF в других приложениях, которые импортируют данные о GIS. Например, RoadRunner позволяет вам добавить данные о вертикальном изменении от файлов GeoTIFF до сцен.

Ссылки

[1] OpenTopography. “Таскалуса, AL: Сезонная Динамика Наплыва И Бесхарактерные Сообщества”, 2011. https://doi.org/10.5069/G9SF2T3K.