在二维和三维地图上可视化无人机飞行路径

本示例直观展示了从夏威夷的莫纳罗亚基线观测站飞往莫纳罗亚火山顶部的仿真无人机 (UAV) 飞行过程。首先，将轨迹显示在地理坐标区和地理球体上。然后，通过使用相机导航函数，同步视图并可视化飞行路径。最后，欣赏毛纳洛亚火山顶部的全景图。

在二维空间中可视化兴趣面

利用无人机监测火山周边不断变化的地形特征、气体及火山灰云的轨迹，正逐渐成为科学家们的一个重要研究领域 [1]。无人机可以在对火山学家而言存在危险的兴趣面 (AOIs) 内飞行。在派遣无人机执行任务之前，仿真其飞行路径有助于了解地形和地貌。要全面了解感兴趣区域 (AOI) 并从二维视角观察，请在地理坐标区中查看莫纳罗亚基线观测站和莫纳罗亚火山的位置。

获取莫纳罗亚基线观测站的坐标

请指定莫纳罗亚基线观测站的坐标 [2]。观测站的高度以米为单位，以平均海平面 (MSL) 为基准。

obslat = 19.5362;
obslon = -155.5763;
obsH = 3397.00;

获取莫纳罗亚火山的坐标

请指定莫纳罗亚火山顶部的坐标 [3]。该火山的高度为正高，单位为米。

mllat = 19.475;
mllon = -155.608;
mlH = 4169;

以二维视图查看莫纳罗亚基线观测站和莫纳罗亚火山

若要从二维视角观察感兴趣区域 (AOI)，请使用 geoaxes 和 geoplot 绘制观测站和火山顶部的位置。

figure
geoaxes(Basemap="satellite",ZoomLevel=12)
hold on
geoplot(obslat,obslon,"ow",MarkerSize=10,MarkerFaceColor="magenta",DisplayName="Mauna Loa Observatory")
geoplot(mllat,mllon,"ow",MarkerSize=10,MarkerFaceColor="blue",DisplayName="Mauna Loa Volcano")
legend

同步显示莫纳罗亚基线观测站的二维和三维视图

使用地理坐标区以二维视角查看天文台，使用地理球体以三维视角查看天文台。

在同一图中创建地理坐标区和地理球体

通过在同一 UI 图形中创建地理坐标区和地理球体，设置 2D 和三维地图显示。若要查看更多二维地图内容，请将地理坐标区的 InnerPosition 设置为其 OuterPosition。若要使两个地图显示使用同一底图，请将地理坐标区的底图设置为 "satellite"。

figpos = [1000 500 800 400];
uif = uifigure(Position=figpos);
ug = uigridlayout(uif,[1,2]);
p1 = uipanel(ug);
p2 = uipanel(ug);
gx = geoaxes(p1,Basemap="satellite"); 
gg = geoglobe(p2); 
gx.InnerPosition = gx.OuterPosition;
gg.Position = [0 0 1 1];

2-D map and 3-D globe in the same figure

以二维视图查看观测站

从海拔 200 米的地形高处眺望观测站。通过调整地图中心和缩放级别，控制地理坐标区的显示范围。您可以通过将球体的摄像机高度转换为地理坐标区的缩放级别，从而使地理坐标区的视图与球体的视图保持同步。使用 heightToZoomLevel 本地函数，根据地形高度计算近似缩放级别。

heightAboveTerrain = 200;
gx.MapCenter = [obslat obslon];
zoomLevel = heightToZoomLevel(heightAboveTerrain,obslat);
gx.ZoomLevel = zoomLevel;

以三维形式查看观测站

通过调整摄像机的位置来控制地理球体的视角。campos 函数要求您指定椭球高（相对于 WGS84 椭球），而不是正高（相对于平均海平面）。将观测站的高度转换为椭球高。所有高度均以米为单位。

N = egm96geoid(obslat,obslon);
obsh = obsH + N;
ellipsoidalHeight = obsh + heightAboveTerrain;
campos(gg,obslat,obslon,ellipsoidalHeight)
drawnow

2-D map and 3-D globe showing satellite images of the observatory

导入飞行轨迹数据并计算航向和三维距离

将从莫纳罗亚基线观测站到莫纳罗亚火山顶部的仿真飞行轨迹导入。该文件包含无人机路径的纬度、经度和高度数据，这些数据均以平均海平面为基准。

T = readgeotable("sample_uavtrack.gpx",Layer="track_points");
tlat = T.Shape.Latitude';
tlon = T.Shape.Longitude';
talt = T.Elevation';

计算飞行航向

使用 azimuth 函数计算每个航迹点处的无人机航向。azimuth 函数采用顺时针为正的约定来计算值。

wgs84 = wgs84Ellipsoid;
theading = azimuth(tlat(1:end-1),tlon(1:end-1),tlat(2:end),tlon(2:end),wgs84);
theading = [theading(1);theading(:)];

计算三维距离

计算无人机飞行轨迹的累计距离。distance 函数不会考虑高程或海拔的变化。为了计算无人机在三维空间中从一点移动到另一点的距离，需要使用地心笛卡尔坐标系（X、Y、Z）。使用 ecefOffset 函数计算点对点偏移分量（单位：米）。该无人机飞行的高度数据是以平均海平面为基准的。要使用 ecefOffset 函数，高度值必须以椭球体为基准。将飞行轨迹的正高转换为椭球高（以 WGS84 椭球为基准）。所有高度均以米为单位。

N = egm96geoid(tlat,tlon);
h = talt + N;

计算距离偏移量。

lat1 = tlat(1:end-1);
lat2 = tlat(2:end);
lon1 = tlon(1:end-1);
lon2 = tlon(2:end);
h1 = h(1:end-1);
h2 = h(2:end);
[dx,dy,dz] = ecefOffset(wgs84,lat1,lon1,h1,lat2,lon2,h2);

使用 hypot 函数计算每对相邻点之间的欧几里得距离。距离单位是米。

distanceIncrementIn3D = hypot(hypot(dx, dy), dz);

计算三维空间中的累计距离以及总距离（单位：米）。

cumulativeDistanceIn3D = cumsum(distanceIncrementIn3D);
totalDistanceIn3D = sum(distanceIncrementIn3D);
fprintf("Total UAV track distance is %f meters.\n",totalDistanceIn3D)

Total UAV track distance is 8931.072120 meters.

为用于绘制动画的累计距离分配一个变量。

tdist = [0 cumulativeDistanceIn3D];

绘制从莫纳罗亚基线观测站到莫纳罗亚火山顶部的飞行路线

绘制从莫纳罗亚基线观测站到莫纳罗亚火山顶部的仿真飞行路线。

绘制飞行轨迹。默认情况下，地理球体会将该线置于显示区域的中心。由于您之前已设置了 MapCenter 和 ZoomLevel，因此地理坐标区的显示不会发生变化。

geoplot3(gg,tlat,tlon,talt,"c",LineWidth=2,HeightReference="geoid")
ptrack = geoplot(gx,tlat,tlon,"c",LineWidth=2);

通过将球体的摄像机高度转换为坐标轴的缩放级别，使地图中心和缩放级别与三维视图保持一致。

[clat,clon,cheight] = campos(gg);
gx.MapCenter = [clat clon];
gx.ZoomLevel = heightToZoomLevel(cheight,clat);
drawnow

2-D map and 3-D globe showing flight line plotted over satellite imagery

将初始视点从莫纳罗亚基线观测站设置为莫纳罗亚火山顶部

将摄像机位置设置为轨迹的第一个坐标，即可从起始位置查看飞行轨迹。为了获得更好的视角，请将摄像机高度设置为 75 米，这大约是轨迹的高度。将摄像机俯仰角设置为-90，即可俯视观测站。由于轨迹的前两个点位于同一位置，且这些位置的计算航向均为 0，因此请将航向设置为计算航向数组的第三个元素来查看该轨迹。

campos(gg,tlat(1),tlon(1))
camheight(gg,talt(1) + 75)
campitch(gg,-90)
camheading(gg,theading(3))

在二维地图上，使用标记显示飞行轨迹的起始和终点位置，并使用图标显示无人机的位置。为无人机轨迹、标记和图标创建图例。

hold(gx,"on")
mstart = geoplot(gx,tlat(1),tlon(1),"ow",MarkerSize=10,MarkerFaceColor="magenta");
mend = geoplot(gx,tlat(end),tlon(end),"ow",MarkerSize=10,MarkerFaceColor="blue");
icon = geoiconchart(gx,tlat(1),tlon(1),"uav.png",SizeData=30);

icon.DisplayName = "Current Location";
mstart.DisplayName = "Start Location";
mend.DisplayName = "End Location";
ptrack.DisplayName = "UAV Track";
legend(gx)

从美国地质调查局 (USGS) 国家地图中添加 USGS 地形图底图。然后，通过更改底图和缩放级别来查看 AOI 的拓扑结构。

url = "https://basemap.nationalmap.gov/ArcGIS/rest/services/USGSTopo/MapServer/tile/${z}/${y}/${x}/png";
addCustomBasemap("usgstopo",url,Attribution="USGS The National Map")

gx.Basemap = "usgstopo";
gx.ZoomLevel = 11;

创建一个与无人机当前纬度和经度相对应的数据提示。

dt = datatip(ptrack,DataIndex=1,Location="northwest");

通过设置纬度和经度行格式的样式，并添加显示无人机海拔高度、航向以及距观测站距离的行，自定义数据提示。

ptrack.DataTipTemplate.DataTipRows(1).Format = "%.3f";
ptrack.DataTipTemplate.DataTipRows(2).Format = "%.3f";

dtrow = dataTipTextRow("Distance",tdist,"%.2f");
dtrow(end+1) = dataTipTextRow("Altitude",talt,"%.2f");
dtrow(end+1) = dataTipTextRow("Heading",theading,"%.2f");
ptrack.DataTipTemplate.DataTipRows(end+1:end+3) = dtrow;

2-D map and 3-D globe. The 2-D map shows the flight path plotted over a topographic basemap. The 3-D globe shows the beginning of the flight path plotted over satellite imagery.

从莫纳罗亚基线观测站飞往莫纳罗亚火山顶峰

制作一段从莫纳罗亚基线观测站飞往莫纳罗亚火山顶部的动画。

为了更清晰地观察无人机飞往火山顶部的飞行轨迹，请更新摄像机的俯仰角值。

pitch = -2.7689;
campitch(gg,pitch)

在二维地图和三维球体上显示飞行轨迹。

通过更新图标的纬度、经度和旋转角度，在二维地图上查看无人机的位置。IconRotation 属性以逆时针方向（正方向）旋转角度，因此请使用 azimuth 函数返回值的负数。
通过更新数据提示以使用当前索引，查看无人机的当前位置、高度和航向。
通过设置摄像机位置，为三维飞行场景添加动画效果。

for k = 2:(length(tlat)-1)    
    % update icon on 2-D map
    set(icon,"LatitudeData",tlat(k),"LongitudeData",tlon(k),"IconRotation",-theading(k))
    
    % update data tip on 2-D map
    dt.DataIndex = k;
    
    % update camera position for 3-D globe
    campos(gg,tlat(k),tlon(k))
    camheight(gg,talt(k)+100)
    camheading(gg,theading(k))
   
    drawnow
    pause(.25)
end

campos(gg,tlat(end),tlon(end),talt(end)+100)
dt.DataIndex = length(tlat);

Animation of the flight path in 2-D and 3-D

欣赏从莫纳罗亚火山顶端拍摄的 360 度全景图

通过将相机方向和图标旋转 360 度，欣赏从莫纳罗亚火山顶端拍摄的 360 度全景图。以 5 度的步长顺时针旋转，并从下一个 5 度的步长开始。更新航向数据提示。

initialHeading = camheading(gg);
increment = 5;
initialHeading = initialHeading + (increment - mod(initialHeading,increment));

for degree = initialHeading:increment:initialHeading+360
    heading = mod(degree,360);
    camheading(gg,heading);
    icon.IconRotation = -heading;
    ptrack.DataTipTemplate.DataTipRows(end).Value(dt.DataIndex) = heading;
    drawnow
end

The same 2-D map and 3-D globe. The camera for the globe faces the horizon.

局部函数

将高度（以米为单位，相对于 WGS84 椭球体）转换为缩放级别

function zoomLevel = heightToZoomLevel(height, lat)
    earthCircumference = 2 * pi * 6378137;
    zoomLevel = log2((earthCircumference *cosd(lat)) / height) + 1;
    zoomLevel = max(0, zoomLevel);
    zoomLevel = min(19, zoomLevel);
end

参考

[1] Williams, Sarah C. P. “Studying Volcanic Eruptions with Aerial Drones.” Proceedings of the National Academy of Sciences of the United States of America 110, no. 27 (July 2, 2013): 10881. https://doi.org/10.1073/pnas.1309922110.

[2] NOAA. “Mauna Loa Baseline Observatory.” Global Monitoring Laboratory. Accessed June 16, 2020. https://gml.noaa.gov/obop/mlo/.

[3] USGS. “Mauna Loa.” Hawaiian Volcano Observatory. Accessed June 16, 2020. https://www.usgs.gov/volcanoes/mauna-loa.

另请参阅