Visualize Aircraft Line-of-Sight over Terrain

This example shows how to compute and visualize the line-of-sight visibility of an aircraft from a ground location over terrain. First, import terrain data for a region and apply it to a 3-D geographic globe. Then, perform point-to-point visibility analysis from a ground location to a simulated flight path, and display the results on a 3-D geographic globe. Finally, perform point-to-area visibility analysis from the ground location corresponding to aircraft flying at constant altitude, and display the results on a 2-D geographic axes.

Use line-of-sight analysis for ground-to-air scenarios where unobstructed visibility is important, such as for radar surveillance, communications, and UAV path planning. This example applies the analysis to radar surveillance for an airport.

Import Terrain Data

Specify a DTED-format terrain file to use for data analysis and 3-D visualization. The terrain file was downloaded from the "SRTM Void Filled" data set available from the United States Geological Survey (USGS).

dtedfile = "n39_w106_3arc_v2.dt1";
attribution = "SRTM 3 arc-second resolution. Data available from the U.S. Geological Survey.";

Import DTED file data into the workspace as an array and a geographic raster reference object, specifying the return type as double so that the data works with all analysis functions.

[Zterrain,Rterrain] = readgeoraster(dtedfile,"OutputType","double");

View the geographic limits and sample resolution of the terrain data by accessing properties of the geographic raster reference object. The limits for the file correspond to the region around Boulder, Colorado, US, and the resolution corresponds to the DTED level-1 format, which has sample resolution of 3 arc seconds, or about 90 meters.

latlim = Rterrain.LatitudeLimits;
lonlim = Rterrain.LongitudeLimits;
latspc = Rterrain.SampleSpacingInLatitude;
lonspc = Rterrain.SampleSpacingInLongitude;
disp("Latitude limits of terrain: " + mat2str(latlim) + newline + ...
    "Longitude limits of terrain: " + mat2str(lonlim) + newline + ...
    "Terrain resolution in latitude: " + latspc*3600 + " arc seconds" + newline + ...
    "Terrain resolution in longitude: " + lonspc*3600 + " arc seconds")

Latitude limits of terrain: [39 40]
Longitude limits of terrain: [-106 -105]
Terrain resolution in latitude: 3 arc seconds
Terrain resolution in longitude: 3 arc seconds

Visualize Aircraft Trajectory Line-of-Sight on a 3-D Map

Create Geographic Globe with Custom Terrain

Add custom data with the DTED file for use with 3-D visualization.

addCustomTerrain("southboulder",dtedfile,"Attribution",attribution)

Specify the custom terrain with a new geographic globe.

fig = uifigure;
g = geoglobe(fig,"Terrain","southboulder");

View Radar Ground Location

Define a radar ground location at Rocky Mountain Metropolitan Airport. The radar is mounted on a tower 10 meters above the ground. The radar altitude is the sum of the ground elevation and the radar tower height, referenced to mean sea level.

rdrlat = 39.913756;
rdrlon = -105.118062;
rdrtowerht = 10;
rdralt = 1717 + rdrtowerht;

Plot the radar location on the geographic globe.

geoplot3(g,rdrlat,rdrlon,rdralt,"co", ...
    "LineWidth",6, ...
    "MarkerSize",1)

Simulate Aircraft Trajectory

Simulate the trajectory of an aircraft circling over the mountains.

Define the center location of a circling aircraft.

tlat0 = 39.80384;
tlon0 = -105.49916;
tht0 = 3000;

Define trajectory waypoints for the aircraft using east-north-up (ENU) Cartesian coordinates. Specify a curve with a radius of 5 km (5000 m) and a vertical offset of 1 km (1000 m) over 1.5 revolutions. Then, convert the ENU coordinates to geodetic coordinates that are referenced to the WGS84 ellipsoid.

azs = 1:540;
r = 5000;
[X,Y] = pol2cart(deg2rad(azs),r);
Z = linspace(0,1000,numel(azs));
wgs84 = wgs84Ellipsoid;
[tlat,tlon,tht] = enu2geodetic(X,Y,Z,tlat0,tlon0,tht0,wgs84);

View Aircraft Trajectory over Terrain

Plot the aircraft trajectory on the geographic globe. The default view, or camera position, is overhead and oriented down.

hold(g,"on")
traj = geoplot3(g,tlat,tlon,tht,"y", ...
    "HeightReference","ellipsoid", ...
    "LineWidth",3);

View the 3-D terrain and radar location from a distance by changing the camera position and rotation angles.

campos(g,39.77114,-105.62662,6670)
camheading(g,70)
campitch(g,-12)

Compute Line-of-Sight Visibility with Aircraft Trajectory

Compute line-of-sight visibility with the los2 function and the DTED data.

The los2 function supports either orthometric height (height above mean sea level) or height above ground level. Convert the aircraft trajectory heights from ellipsoidal height to orthometric height. Then, compute the line-of-sight from the airport radar location to each aircraft trajectory waypoint and convert the results to a logical array.

numwaypts = numel(tlat);
isvis = zeros(1,numwaypts);
talt = tht - egm96geoid(tlat,tlon);
for k = 1:numwaypts
    isvis(k) = los2(Zterrain,Rterrain,rdrlat,rdrlon,tlat(k),tlon(k),rdralt,talt(k),"MSL","MSL");
end
isvis = logical(isvis);

Note that los2 calculates line-of-sight visibility assuming the data is referenced to a spherical Earth, whereas the data is actually referenced to the WGS84 ellipsoid, and as a result there may be minor discrepancies. The line-of-sight calculation also corresponds to optical line-of-sight and does not account for refraction through the atmosphere.

Visualize Line-of-Sight Visibility over Terrain

Plot the line-of-sight visibility. Use green markers where the aircraft is visible from the airport and magenta markers where it is not visible.

delete(traj)
geoplot3(g,tlat(isvis),tlon(isvis),tht(isvis),"og", ...
    "HeightReference","ellipsoid", ...
    "LineWidth",2, ...
    "MarkerSize",1)
geoplot3(g,tlat(~isvis),tlon(~isvis),tht(~isvis),"om", ...
    "HeightReference","ellipsoid", ...
    "LineWidth",2, ...
    "MarkerSize",1)

View the line-of-sight plot from the perspective of the airport. Get the geodetic coordinates of the position that is 900 meters east, 200 meters north, and 100 meters up from the radar location. Then, set the camera position and rotation angles. The green markers appear in view, but the magenta markers are either completely or partially obstructed by terrain.

rdrht = rdralt + egm96geoid(rdrlat,rdrlon);
[camlat,camlon,camht] = enu2geodetic(900,200,100,rdrlat,rdrlon,rdrht,wgs84);
campos(g,camlat,camlon,camht)
camheading(g,-110)
campitch(g,0)

Visualize Aircraft Line-of-Sight Visibility Contours on a 2-D Map

The previous sections performed point-to-point line-of-sight analysis and visualization from a radar location to an aircraft trajectory. Now perform point-to-area line-of-sight analysis and visualization from the same radar location over the terrain region. The visualization displays the edge of visibility for an aircraft flying at constant altitude.

Plot Radar Location and Terrain Limits on a 2-D Map

Plot the radar location in a new figure with a topographic 2-D map.

figure
geoplot(rdrlat,rdrlon,"co", ...
    "LineWidth",6, ...
    "MarkerSize",3, ...
    "DisplayName","Radar location")
geobasemap topographic
gx = gca;
gx.InnerPosition = gx.OuterPosition;

Display the limits of the custom terrain as a rectangle on the map.

latmin = latlim(1);
latmax = latlim(2);
lonmin = lonlim(1);
lonmax = lonlim(2);
hold on
geoplot([latmin latmin latmax latmax latmin],[lonmin lonmax lonmax lonmin lonmin], ...
    "LineWidth",1, ...
    "Color","k", ...
    "DisplayName","Terrain limits")

Display a legend in the northwest corner.

legend("Location","northwest")

Plot Visibility Contours for Aircraft Flying at Constant Altitude

Specify three altitudes above mean sea level for the aircraft. For each altitude:

• Compute the viewshed using the radar location as the observer. The viewshed defines the area that has line-of-sight visibility.

• Find the edge of aircraft visibility by computing contours from the viewshed data.

• Remove small contour segments.

• Plot the contours on the geographic axes.

tgtalts = [3000 4000 5000];

minVertices = 10;
cfig = figure("Visible","off"); % Suppress contour plot using invisible figure
cax = axes("Parent",cfig);
for tgtalt = tgtalts
    vis = viewshed(Zterrain,Rterrain,rdrlat,rdrlon,rdralt,tgtalt,"MSL","MSL");
    
    C = contourm(vis,Rterrain,"LevelList",1,"Parent",cax);
    clat = C(2,:);
    clon = C(1,:);
    
    clats = [];
    clons = []; 
    k = 1;   
    while k < size(C,2)
        numVertices = clat(k);
        if numVertices > minVertices % Do not plot small segments 
            clats = [clats clat(k+1:k+numVertices) NaN]; %#ok<AGROW> 
            clons = [clons clon(k+1:k+numVertices) NaN]; %#ok<AGROW> 
        end
        k = k + numVertices + 1;
    end
    
    geoplot(gx,clats,clons,"LineWidth",2, ...
        "DisplayName", "Aircraft: " + string(tgtalt) + " m");
end

The contours appear primarily to the west of the radar location over the mountains. The contours do not appear in other directions because visibility is not constrained by the terrain in those directions within the terrain data limits.

If the radar is constrained by line-of-sight visibility, then the contours correspond to radar coverage regions for varying altitude, where the nearest contour to the radar corresponds to radar coverage for an aircraft flying at 3000 meters and the furthest contour corresponds to radar coverage for an aircraft flying at 5000 meters.

As with los2, the viewshed function calculates line-of-sight visibility assuming the data is referenced to a spherical Earth, whereas the data is actually referenced to the WGS84 ellipsoid, and as a result there may be minor discrepancies. The line-of-sight calculation also corresponds to optical line-of-sight and does not account for refraction through the atmosphere.

Clean Up

Clean up by closing the geographic globe and removing the imported terrain data.

if isvalid(fig)
    close(fig)
end
removeCustomTerrain("southboulder")