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.
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
Add custom data with the DTED file for use with 3-D visualization.
Specify the custom terrain with a new geographic globe. Preserve the custom terrain on the globe when data is added by setting hold on.
fig = uifigure; g = geoglobe(fig,"Terrain","southboulder"); hold(g,"on")
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 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);
Plot the aircraft trajectory on the geographic globe. The default view, or camera position, is overhead and oriented down.
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 the
los2 function and the DTED data.
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);
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.
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)
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 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") hold on 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); 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.
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.
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 by closing the geographic globe and removing the imported terrain data.
if isvalid(fig) close(fig) end removeCustomTerrain("southboulder")