los2
Line-of-sight visibility between two points in terrain
Syntax
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(...)
Description
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.
Examples
Find Elevation Angle of Point 90 Degrees from Observer
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'
AngleUnit: 'degree'
los2(Z,R,lat1,lon1,lat2,lon2,alt1,alt1);
Determine Line-of-Sight
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)
