Main Content


Line-of-sight visibility between two points in terrain


vis = los2(Z,R,lat1,lon1,lat2,lon2)
vis = los2(Z,R,lat1,lon1,lat2,lon2,alt1)
vis = los2(Z,R,lat1,lon1,lat2,lon2,alt1,alt2)
vis = los2(Z,R,lat1,lon1,lat2,lon2,alt1,alt2,alt1opt)
vis = los2(Z,R,lat1,lon1,lat2,lon2,alt1,alt2,alt1opt,alt2opt)
vis = los2(Z,R,lat1,lon1,lat2,lon2,alt1,alt2,alt1opt, ...
vis = los2(Z,R,lat1,lon1,lat2,lon2,alt1,alt2,alt1opt, ...
[vis,visprofile,dist,H,lattrk,lontrk] = los2(...)


los2 computes the mutual visibility between two points on a displayed digital elevation map. los2 uses the current object if it is a regular data grid, or the first regular data grid found on the current axes. The grid's zdata is used for the profile. The color data is used in the absence of data in z. The two points are selected by clicking on the map. The result is displayed in a new figure. Markers indicate visible and obscured points along the profile. The profile is shown in a Cartesian coordinate system with the origin at the observer's location. The displayed z coordinate accounts for the elevation of the terrain and the curvature of the body.

vis = los2(Z,R,lat1,lon1,lat2,lon2) computes the mutual visibility between pairs of points on a digital elevation map. The elevations are provided as a regular data grid Z containing elevations in units of meters. The two points are provided as vectors of latitudes and longitudes in units of degrees. The resulting logical variable vis is equal to one when the two points are visible to each other, and zero when the line of sight is obscured by terrain. If any of the input arguments are empty, los2 attempts to gather the data from the current axes. With one or more output arguments, no figures are created and only the data is returned.

R can be a geographic raster reference object, a referencing vector, or a referencing matrix. If R is a geographic raster reference object, its RasterSize property must be consistent with size(Z).

If R is a referencing vector, it must be a 1-by-3 with elements:

[cells/degree northern_latitude_limit western_longitude_limit]

If R is a referencing matrix, it must be 3-by-2 and transform raster row and column indices to or from geographic coordinates according to:

[lon lat] = [row col 1] * R

If R is a referencing matrix, it must define a (non-rotational, non-skewed) relationship in which each column of the data grid falls along a meridian and each row falls along a parallel. Nearest-neighbor interpolation is used by default. NaN is returned for points outside the grid limits or for which lat or lon contain NaN. All angles are in units of degrees.

vis = los2(Z,R,lat1,lon1,lat2,lon2,alt1) places the first point at the specified altitude in meters above the surface (on a tower, for instance). This is equivalent to putting the point on a tower. If omitted, point 1 is assumed to be on the surface. alt1 may be either a scalar or a vector with the same length as lat1, lon1, lat2, and lon2.

vis = los2(Z,R,lat1,lon1,lat2,lon2,alt1,alt2) places both points at a specified altitudes in meters above the surface. alt2 may be either a scalar or a vector with the same length as lat1, lon1, lat2, and lon2. If alt2 is omitted, point 2 is assumed to be on the surface.

vis = los2(Z,R,lat1,lon1,lat2,lon2,alt1,alt2,alt1opt) controls the interpretation of alt1 as either a relative altitude (alt1opt equals 'AGL', the default) or an absolute altitude (alt1opt equals 'MSL'). If the altitude option is 'AGL', alt1 is interpreted as the altitude of point 1 in meters above the terrain (“above ground level”). If alt1opt is 'MSL', alt1 is interpreted as altitude above zero (“mean sea level”).

vis = los2(Z,R,lat1,lon1,lat2,lon2,alt1,alt2,alt1opt,alt2opt) controls the interpretation ALT2.

vis = los2(Z,R,lat1,lon1,lat2,lon2,alt1,alt2,alt1opt, ...
does the visibility calculation on a sphere with the specified radius. If omitted, the radius of the earth in meters is assumed. The altitudes, elevations and the radius should be in the same units. This calling form is most useful for computations on bodies other than the earth.

vis = los2(Z,R,lat1,lon1,lat2,lon2,alt1,alt2,alt1opt, ...
assumes a larger radius for propagation of the line of sight. This can account for the curvature of the signal path due to refraction in the atmosphere. For example, radio propagation in the atmosphere is commonly treated as straight line propagation on a sphere with 4/3 the radius of the earth. In that case the last two arguments would be R_e and 4/3*R_e, where R_e is the radius of the earth. Use Inf as the effective radius for flat earth visibility calculations. The altitudes, elevations and radii should be in the same units.

[vis,visprofile,dist,H,lattrk,lontrk] = los2(...), for scalar inputs (lat1, lon1, etc.), returns vectors of points along the path between the two points. visprofile is a logical vector containing true (logical(1)) where the intermediate points are visible and false (logical(0)) otherwise. dist is the distance along the path (in meters or the units of the radius). H contains the terrain profile relative to the vertical datum along the path. lattrk and lontrk are the latitudes and longitudes of the points along the path. For vector inputs los2 returns visprofile, dist, H, lattrk, and lontrk as cell arrays, with one cell per element of lat1,lon1, etc.

los2(...), with no output arguments, displays the visibility profile between the two points in a new figure.


collapse all

Determine the line-of-sight for two points. Display the visibility profile in 2-D and 3-D.

Import terrain data for an area around South Boulder Peak in Colorado. Prepare the data for analysis by specifying the output type as 'double'.

[Z,R] = readgeoraster("n39_w106_3arc_v2.dt1","OutputType","double");

Specify the coordinates of South Boulder Peak and a nearby point. Determine whether the nearby point is visible from a location 100 meters above the peak by assigning the result of los2 to a variable.

latPeak = 39.95384;
lonPeak = -105.29916;
latNearby = 39.96955;
lonNearby = -105.29867;
vis = los2(Z,R,latPeak,lonPeak,latNearby,lonNearby,100)
vis = logical

The value of vis is 0, which means the point is not visible from the peak.

To display the visibility profile in 2-D, call los2 without an output argument.


To display the visibility profile in 3-D, return visibility information about points between the coordinates you specified. The outputs lattrk and lontrk contain the coordinates of the points. The output visprofile contains logical values indicating whether the points are visible.

[~,visprofile,~,~,lattrk,lontrk] = los2(Z,R,latPeak,lonPeak,latNearby,lonNearby,100);

Create a geographic globe that includes custom terrain for the area around South Boulder Peak.

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

Plot the line-of-sight between the coordinates you specified. Use a magenta line to indicate the point is not visible from above the peak.

geoplot3(g,[latPeak latNearby],[lonPeak lonNearby],[100 0],"m","LineWidth",2, ...

Plot the intermediate points. Use magenta markers for points that are not visible and green markers for points that are visible. Plot each value 10 meters above the terrain so the markers are not obscured by the terrain.

geoplot3(g,lattrk(~visprofile),lontrk(~visprofile),10, ...
geoplot3(g,lattrk(visprofile),lontrk(visprofile),10,"go", ...

Adjust the view by setting the camera position, pitch, roll, and heading.


Close the figure containing the geographic globe and remove the custom terrain.


Find the elevation angle of a point 90 degrees from an observer assuming that the observer and the target are both 1000 km above the Earth.

lat1 = 0; 
lon1 = 0; 
alt1 = 1000*1000;
lat2 = 0; 
lon2 = 90; 
alt2 = 1000*1000;
[az, elev, r] = geodetic2aer(lat2,lon2,alt2,lat1,lon1,alt1,referenceEllipsoid('grs 80'))
az = 90
elev = -45
r = 1.0434e+07

Visually check the result using the los2 line of sight function. Construct a data grid of zeros to represent the Earth's surface. The los2 function with no output arguments creates a figure displaying the geometry.

Z = zeros(181,361);
R = georefpostings([-90 90],[-180 180], size(Z))
R = 
  GeographicPostingsReference with properties:

              LatitudeLimits: [-90 90]
             LongitudeLimits: [-180 180]
                  RasterSize: [181 361]
        RasterInterpretation: 'postings'
            ColumnsStartFrom: 'south'
               RowsStartFrom: 'west'
     SampleSpacingInLatitude: 1
    SampleSpacingInLongitude: 1
      RasterExtentInLatitude: 180
     RasterExtentInLongitude: 360
            XIntrinsicLimits: [1 361]
            YIntrinsicLimits: [1 181]
        CoordinateSystemType: 'geographic'
               GeographicCRS: []
                   AngleUnit: 'degree'


Figure contains an axes object. The axes object contains 5 objects of type line. These objects represent Terrain, Line of Sight, Visible, Obscured, Observer.

Version History

Introduced before R2006a