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,
...

alt2opt,actualradius)

vis = los2(Z,R,lat1,lon1,lat2,lon2,alt1,alt2,alt1opt,
...

alt2opt,actualradius,effectiveradius)

[vis,visprofile,dist,H,lattrk,lontrk]
= los2(...)

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.

alt2opt,actualradius)

```
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

alt2opt,actualradius,effectiveradius)`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.

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'

los2(Z,R,lat1,lon1,lat2,lon2,alt1,alt1);

Z = 500*peaks(100); refvec = [1000 0 0]; [lat1, lon1, lat2, lon2] = deal(-0.027, 0.05, -0.093, 0.042); los2(Z,refvec,lat1,lon1,lat2,lon2,100);

figure; axesm('globe','geoid',earthRadius('meters')) meshm(Z, refvec, size(Z), Z); axis tight camposm(-10,-10,1e6); camupm(0,0) demcmap('inc', Z, 1000); shading interp; camlight [vis,visprofile,dist,h,lattrk,lontrk] = ... los2(Z,refvec,lat1,lon1,lat2,lon2,100); plot3m(lattrk([1;end]),lontrk([1; end]),... h([1; end])+[100; 0],'r','linewidth',2) plotm(lattrk(~visprofile),lontrk(~visprofile),... h(~visprofile),'r.','markersize',10) plotm(lattrk(visprofile),lontrk(visprofile),... h(visprofile),'g.','markersize',10)