viewshed

Create sample terrain data by using the peaks and georefcells functions.

Z = 500*peaks(100);
R = georefcells([-0.1 0],[0 0.1],size(Z));

Calculate the viewshed for a sample point that is 200 meters above the terrain.

latObs = -0.027;
lonObs = 0.05;
[vis,visR] = viewshed(Z,R,latObs,lonObs);

Display the viewshed as a surface on an axesm-based map. Plot the sample point so that it appears above the terrain.

figure
axesm("globe","geoid",earthRadius)
meshm(vis,visR,size(Z),Z)
axis tight
camposm(-10,-10,1e6)
camupm(0,0)

[row,col] = geographicToDiscrete(R,latObs,lonObs);
h = Z(row,col);
plot3m(latObs,lonObs,h+200,"ko","MarkerSize",7,"MarkerFaceColor","r")

Adjust the colormap and lighting. Add a colorbar.

colormap(flipud(summer(2))) 
brighten(0.75)
shading interp
camlight

cb = lcolorbar(["Obscured" "Visible"]);
cb.Position = [0.85 0.5 0.02 0.1];

Display the area that is visible by radar for an aircraft flying above the terrain. You can use the viewshed function to model radio wave propagation in the atmosphere as straight lines on a sphere.

Load terrain elevation data for an area around the Korean peninsula. To make the water areas appear consistent on the plot, set all elevations below sea level (korea5c<0) to the same value.

load korea5c
korea5c(korea5c<0) = -1;

Specify the location of the aircraft. Then, calculate visible areas for the aircraft, at 3000 meters above mean sea level, of the surface. Model radio propagation in the atmosphere by specifying an effective radius for the reference sphere.

latObs = 34.0931;
lonObs = 125.6578; 
hObs = 3000;
hTarget = 0;
rad = earthRadius;
[vis1,vis1R] = viewshed(korea5c,korea5cR,latObs,lonObs,hObs,hTarget,"MSL","AGL",rad,4/3*rad);

Set up a relief map for the area by adjusting the data aspect ratio, the plot box aspect ratio, and the lighting.

figure
worldmap(korea5c,korea5cR)
setm(gca,"geoid",[1 0])

da = daspect;
pba = pbaspect;
da(3) = 7.5*pba(3)/da(3);
daspect(da);

camlight(90,5)
camlight(0,5)
material([0.25 0.8 0])

Display the terrain data, the viewshed, and the aircraft on the map. Adjust the colormap so that the visible areas appear in blue and the obscured areas appear in red. Note that some ocean areas are obscured by an island.

meshm(vis1,vis1R,size(korea5c),korea5c)
plotm(latObs,lonObs,"ko","MarkerFaceColor","k")
colormap([1 0 0; 0 0 1])

Calculate the viewshed again, this time calculating visible areas for the aircraft at 3000 meters of other aircraft flying above it at 5000 meters. The visible area is larger, with some of the area obscured by the island and mountains.

hTarget = 5000;
[vis2,vis2R] = viewshed(korea5c,korea5cR,latObs,lonObs,hObs,hTarget,"MSL","MSL",rad,4/3*rad);
clmo surface
meshm(vis2,vis2R,size(korea5c),korea5c)

Since R2024a

Calculate a viewshed using terrain from multiple adjacent terrain tiles by merging the tiles into one raster.

Read elevation data from two DTED files into the workspace. Merge the data by using the mergetiles function.

[Z1,R1] = readgeoraster("n39_w106_3arc_v2.dt1",OutputType="double");
[Z2,R2] = readgeoraster("n40_w106_3arc_v2.dt1",OutputType="double");
[Z,R] = mergetiles(Z1,R1,Z2,R2);

Specify the location of an observation point that is 100 meters above the terrain.

lat = 39.9597;
lon = -105.6883;
h = 100;

Calculate the viewshed for the observation point. The output indicates visible areas using 1 values and obscured areas using 0 values. To avoid plotting the obscured areas, replace the 0 values with NaN values.

v = viewshed(Z,R,lat,lon,h);

idx = (v == 0);
v(idx) = NaN;

Create a map of the region.

figure
usamap(Z,R)

Display the elevation data as a surface. To make the elevation data appear beneath the viewshed, specify the surface heights using 0s and specify the surface colors using the elevations. Change the colormap.

zrs = zeros(R.RasterSize);
geoshow(zrs,R,DisplayType="surface",CData=Z)
colormap gray

Display the viewshed as a blue surface.

geoshow(v,R,DisplayType="surface",FaceColor="#0072BD")

Display the observation point using a yellow circle marker.

geoshow(lat,lon,ZData=h,DisplayType="point", ...
    Marker="o",MarkerFaceColor="#EDB120",MarkerEdgeColor="k")

Add a labeled color bar.

c = colorbar;
c.Label.String = "Elevation (m)";

The elevation data used in this example is from the US Geological Survey.