Skip to content

Commit 1c64643

Browse files
brendancolbrendancol
committed
Add line_of_sight with Fresnel zone support (#1145)
1 parent f0257ed commit 1c64643

File tree

2 files changed

+207
-1
lines changed

2 files changed

+207
-1
lines changed

‎xrspatial/tests/test_visibility.py‎

Lines changed: 89 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -66,3 +66,92 @@ def test_dask_matches_numpy(self):
6666
elev_np, _, _ = _extract_transect(raster_np, cells)
6767
elev_da, _, _ = _extract_transect(raster_dask, cells)
6868
np.testing.assert_array_equal(elev_np, elev_da)
69+
70+
71+
from xrspatial.visibility import line_of_sight
72+
73+
74+
class TestLineOfSight:
75+
def test_flat_terrain_all_visible(self):
76+
data = np.zeros((5, 10), dtype=float)
77+
raster = _make_raster(data)
78+
result = line_of_sight(raster, x0=0, y0=2, x1=9, y1=2,
79+
observer_elev=10, target_elev=10)
80+
assert isinstance(result, xr.Dataset)
81+
assert 'visible' in result
82+
assert 'elevation' in result
83+
assert 'los_height' in result
84+
assert 'distance' in result
85+
assert result['visible'].all()
86+
87+
def test_obstruction_blocks_view(self):
88+
data = np.zeros((1, 10), dtype=float)
89+
data[0, 5] = 100 # tall wall in the middle
90+
raster = _make_raster(data)
91+
result = line_of_sight(raster, x0=0, y0=0, x1=9, y1=0,
92+
observer_elev=1, target_elev=0)
93+
vis = result['visible'].values
94+
# observer cell is visible
95+
assert vis[0]
96+
# cells before the wall are visible
97+
assert all(vis[:6])
98+
# at least some cells after the wall are blocked
99+
assert not all(vis[6:])
100+
101+
def test_observer_equals_target(self):
102+
data = np.ones((5, 5), dtype=float)
103+
raster = _make_raster(data)
104+
result = line_of_sight(raster, x0=2, y0=2, x1=2, y1=2)
105+
assert len(result['sample']) == 1
106+
assert result['visible'].values[0]
107+
108+
def test_elevation_offsets(self):
109+
data = np.zeros((1, 5), dtype=float)
110+
raster = _make_raster(data)
111+
result = line_of_sight(raster, x0=0, y0=0, x1=4, y1=0,
112+
observer_elev=10, target_elev=20)
113+
los = result['los_height'].values
114+
# LOS starts at 10, ends at 20
115+
assert abs(los[0] - 10.0) < 1e-10
116+
assert abs(los[-1] - 20.0) < 1e-10
117+
118+
def test_distance_monotonic(self):
119+
data = np.zeros((5, 10), dtype=float)
120+
raster = _make_raster(data)
121+
result = line_of_sight(raster, x0=0, y0=0, x1=9, y1=4)
122+
d = result['distance'].values
123+
assert all(d[i] <= d[i + 1] for i in range(len(d) - 1))
124+
125+
def test_fresnel_zone(self):
126+
data = np.zeros((1, 11), dtype=float)
127+
raster = _make_raster(data)
128+
result = line_of_sight(raster, x0=0, y0=0, x1=10, y1=0,
129+
observer_elev=50, target_elev=50,
130+
frequency_mhz=900)
131+
assert 'fresnel_radius' in result
132+
assert 'fresnel_clear' in result
133+
# midpoint has largest Fresnel radius
134+
fr = result['fresnel_radius'].values
135+
mid = len(fr) // 2
136+
assert fr[mid] >= fr[1]
137+
assert fr[mid] >= fr[-2]
138+
# with 50m clearance and flat terrain, Fresnel should be clear
139+
assert result['fresnel_clear'].all()
140+
141+
def test_no_fresnel_by_default(self):
142+
data = np.zeros((5, 5), dtype=float)
143+
raster = _make_raster(data)
144+
result = line_of_sight(raster, x0=0, y0=0, x1=4, y1=4)
145+
assert 'fresnel_radius' not in result
146+
assert 'fresnel_clear' not in result
147+
148+
def test_xy_coords_in_output(self):
149+
data = np.zeros((5, 10), dtype=float)
150+
raster = _make_raster(data)
151+
result = line_of_sight(raster, x0=0, y0=2, x1=9, y1=2)
152+
# first point should match observer
153+
assert abs(result['x'].values[0] - 0.0) < 1e-10
154+
assert abs(result['y'].values[0] - 2.0) < 1e-10
155+
# last point should match target
156+
assert abs(result['x'].values[-1] - 9.0) < 1e-10
157+
assert abs(result['y'].values[-1] - 2.0) < 1e-10

‎xrspatial/visibility.py‎

Lines changed: 118 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -12,8 +12,11 @@
1212
"""
1313

1414
import numpy as np
15+
import xarray
1516

16-
from .utils import has_cuda_and_cupy, has_dask_array, is_cupy_array
17+
from .utils import _validate_raster, has_cuda_and_cupy, has_dask_array, is_cupy_array
18+
19+
SPEED_OF_LIGHT = 299_792_458.0 # m/s
1720

1821

1922
def _bresenham_line(r0, c0, r1, c1):
@@ -64,3 +67,117 @@ def _extract_transect(raster, cells):
6467

6568
elevations = data[rows, cols].astype(np.float64)
6669
return elevations, x_coords, y_coords
70+
71+
72+
def _fresnel_radius_1(d1, d2, freq_hz):
73+
"""First Fresnel zone radius at a point d1 from transmitter, d2 from receiver."""
74+
D = d1 + d2
75+
if D == 0 or freq_hz == 0:
76+
return 0.0
77+
wavelength = SPEED_OF_LIGHT / freq_hz
78+
return np.sqrt(wavelength * d1 * d2 / D)
79+
80+
81+
def line_of_sight(
82+
raster: xarray.DataArray,
83+
x0: float, y0: float,
84+
x1: float, y1: float,
85+
observer_elev: float = 0,
86+
target_elev: float = 0,
87+
frequency_mhz: float = None,
88+
) -> xarray.Dataset:
89+
"""Compute elevation profile and visibility along a straight line.
90+
91+
Parameters
92+
----------
93+
raster : xarray.DataArray
94+
Elevation raster.
95+
x0, y0 : float
96+
Observer location in data-space coordinates.
97+
x1, y1 : float
98+
Target location in data-space coordinates.
99+
observer_elev : float
100+
Height above terrain at the observer.
101+
target_elev : float
102+
Height above terrain at the target.
103+
frequency_mhz : float, optional
104+
Radio frequency in MHz. When set, first Fresnel zone clearance
105+
is computed at each sample point.
106+
107+
Returns
108+
-------
109+
xarray.Dataset
110+
Dataset with dimension ``sample`` containing variables
111+
``distance``, ``elevation``, ``los_height``, ``visible``,
112+
``x``, ``y``, and optionally ``fresnel_radius`` and
113+
``fresnel_clear``.
114+
"""
115+
_validate_raster(raster, func_name='line_of_sight', name='raster')
116+
117+
x_coords = raster.coords['x'].values
118+
y_coords = raster.coords['y'].values
119+
120+
# snap to nearest grid cell
121+
c0 = int(np.argmin(np.abs(x_coords - x0)))
122+
r0 = int(np.argmin(np.abs(y_coords - y0)))
123+
c1 = int(np.argmin(np.abs(x_coords - x1)))
124+
r1 = int(np.argmin(np.abs(y_coords - y1)))
125+
126+
cells = _bresenham_line(r0, c0, r1, c1)
127+
elevations, xs, ys = _extract_transect(raster, cells)
128+
129+
n = len(cells)
130+
131+
# cumulative distance along the transect
132+
distance = np.zeros(n, dtype=np.float64)
133+
for i in range(1, n):
134+
dx = xs[i] - xs[i - 1]
135+
dy = ys[i] - ys[i - 1]
136+
distance[i] = distance[i - 1] + np.sqrt(dx * dx + dy * dy)
137+
138+
total_dist = distance[-1] if n > 1 else 0.0
139+
140+
# LOS height: linear interpolation from observer to target
141+
obs_h = elevations[0] + observer_elev
142+
tgt_h = elevations[-1] + target_elev if n > 1 else obs_h
143+
if total_dist > 0:
144+
los_height = obs_h + (tgt_h - obs_h) * (distance / total_dist)
145+
else:
146+
los_height = np.array([obs_h])
147+
148+
# visibility: track max elevation angle from observer
149+
visible = np.ones(n, dtype=bool)
150+
max_angle = -np.inf
151+
for i in range(1, n):
152+
if distance[i] == 0:
153+
continue
154+
angle = (elevations[i] - obs_h) / distance[i]
155+
if angle >= max_angle:
156+
max_angle = angle
157+
else:
158+
visible[i] = False
159+
160+
data_vars = {
161+
'distance': ('sample', distance),
162+
'elevation': ('sample', elevations),
163+
'los_height': ('sample', los_height),
164+
'visible': ('sample', visible),
165+
'x': ('sample', xs),
166+
'y': ('sample', ys),
167+
}
168+
169+
if frequency_mhz is not None:
170+
freq_hz = frequency_mhz * 1e6
171+
fresnel = np.zeros(n, dtype=np.float64)
172+
fresnel_clear = np.ones(n, dtype=bool)
173+
for i in range(n):
174+
d1 = distance[i]
175+
d2 = total_dist - d1
176+
fresnel[i] = _fresnel_radius_1(d1, d2, freq_hz)
177+
clearance = los_height[i] - elevations[i]
178+
if clearance < fresnel[i]:
179+
fresnel_clear[i] = False
180+
data_vars['fresnel_radius'] = ('sample', fresnel)
181+
data_vars['fresnel_clear'] = ('sample', fresnel_clear)
182+
183+
return xarray.Dataset(data_vars)

0 commit comments

Comments
 (0)