Skip to content

Commit 6e2e246

Browse files
brendancolbrendancol
committed
Add visibility analysis user guide notebook (#1145)
1 parent 4a110da commit 6e2e246

File tree

1 file changed

+222
-0
lines changed

1 file changed

+222
-0
lines changed

‎examples/user_guide/48_Visibility_Analysis.ipynb‎

Lines changed: 222 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,222 @@
1+
{
2+
"cells": [
3+
{
4+
"cell_type": "markdown",
5+
"metadata": {},
6+
"source": [
7+
"# Visibility Analysis\n",
8+
"\n",
9+
"This notebook demonstrates multi-observer viewshed analysis and point-to-point\n",
10+
"line-of-sight profiling.\n",
11+
"\n",
12+
"**Functions covered:**\n",
13+
"- `cumulative_viewshed` -- count of observers with line-of-sight to each cell\n",
14+
"- `visibility_frequency` -- normalized version (fraction of observers)\n",
15+
"- `line_of_sight` -- elevation profile, visibility, and optional Fresnel zone clearance"
16+
]
17+
},
18+
{
19+
"cell_type": "code",
20+
"execution_count": null,
21+
"metadata": {},
22+
"outputs": [],
23+
"source": [
24+
"import numpy as np\n",
25+
"import xarray as xr\n",
26+
"import matplotlib.pyplot as plt\n",
27+
"\n",
28+
"from xrspatial.terrain import generate_terrain\n",
29+
"from xrspatial.hillshade import hillshade\n",
30+
"from xrspatial.visibility import (\n",
31+
" cumulative_viewshed,\n",
32+
" visibility_frequency,\n",
33+
" line_of_sight,\n",
34+
")"
35+
]
36+
},
37+
{
38+
"cell_type": "markdown",
39+
"metadata": {},
40+
"source": [
41+
"## Generate synthetic terrain"
42+
]
43+
},
44+
{
45+
"cell_type": "code",
46+
"execution_count": null,
47+
"metadata": {},
48+
"outputs": [],
49+
"source": [
50+
"terrain = generate_terrain(width=200, height=200, seed=42)\n",
51+
"hs = hillshade(terrain)\n",
52+
"\n",
53+
"fig, ax = plt.subplots(figsize=(8, 8))\n",
54+
"hs.plot.imshow(ax=ax, cmap='gray', add_colorbar=False)\n",
55+
"terrain.plot.imshow(ax=ax, alpha=0.4, cmap='terrain', add_colorbar=True)\n",
56+
"ax.set_title('Synthetic terrain')\n",
57+
"plt.tight_layout()"
58+
]
59+
},
60+
{
61+
"cell_type": "markdown",
62+
"metadata": {},
63+
"source": [
64+
"## Cumulative viewshed\n",
65+
"\n",
66+
"Place three observers on the terrain and count how many can see each cell."
67+
]
68+
},
69+
{
70+
"cell_type": "code",
71+
"execution_count": null,
72+
"metadata": {},
73+
"outputs": [],
74+
"source": [
75+
"xs = terrain.coords['x'].values\n",
76+
"ys = terrain.coords['y'].values\n",
77+
"\n",
78+
"observers = [\n",
79+
" {'x': float(xs[40]), 'y': float(ys[40]), 'observer_elev': 20},\n",
80+
" {'x': float(xs[160]), 'y': float(ys[80]), 'observer_elev': 20},\n",
81+
" {'x': float(xs[100]), 'y': float(ys[160]), 'observer_elev': 20},\n",
82+
"]\n",
83+
"\n",
84+
"cum = cumulative_viewshed(terrain, observers)\n",
85+
"\n",
86+
"fig, ax = plt.subplots(figsize=(8, 8))\n",
87+
"hs.plot.imshow(ax=ax, cmap='gray', add_colorbar=False)\n",
88+
"cum.plot.imshow(ax=ax, alpha=0.6, cmap='YlOrRd',\n",
89+
" add_colorbar=True, vmin=0, vmax=3)\n",
90+
"for obs in observers:\n",
91+
" ax.plot(obs['x'], obs['y'], 'k^', markersize=10)\n",
92+
"ax.set_title('Cumulative viewshed (observer count)')\n",
93+
"plt.tight_layout()"
94+
]
95+
},
96+
{
97+
"cell_type": "markdown",
98+
"metadata": {},
99+
"source": [
100+
"## Visibility frequency\n",
101+
"\n",
102+
"Same data, but normalized to [0, 1]."
103+
]
104+
},
105+
{
106+
"cell_type": "code",
107+
"execution_count": null,
108+
"metadata": {},
109+
"outputs": [],
110+
"source": [
111+
"freq = visibility_frequency(terrain, observers)\n",
112+
"\n",
113+
"fig, ax = plt.subplots(figsize=(8, 8))\n",
114+
"hs.plot.imshow(ax=ax, cmap='gray', add_colorbar=False)\n",
115+
"freq.plot.imshow(ax=ax, alpha=0.6, cmap='RdYlGn',\n",
116+
" add_colorbar=True, vmin=0, vmax=1)\n",
117+
"for obs in observers:\n",
118+
" ax.plot(obs['x'], obs['y'], 'k^', markersize=10)\n",
119+
"ax.set_title('Visibility frequency')\n",
120+
"plt.tight_layout()"
121+
]
122+
},
123+
{
124+
"cell_type": "markdown",
125+
"metadata": {},
126+
"source": [
127+
"## Line of sight\n",
128+
"\n",
129+
"Draw a transect between two points and check visibility along it."
130+
]
131+
},
132+
{
133+
"cell_type": "code",
134+
"execution_count": null,
135+
"metadata": {},
136+
"outputs": [],
137+
"source": [
138+
"ox, oy = float(xs[20]), float(ys[100])\n",
139+
"tx, ty = float(xs[180]), float(ys[100])\n",
140+
"\n",
141+
"los = line_of_sight(terrain, x0=ox, y0=oy, x1=tx, y1=ty,\n",
142+
" observer_elev=15, target_elev=5)\n",
143+
"\n",
144+
"fig, ax = plt.subplots(figsize=(12, 4))\n",
145+
"ax.fill_between(los['distance'].values, 0, los['elevation'].values,\n",
146+
" color='sienna', alpha=0.4, label='Terrain')\n",
147+
"ax.plot(los['distance'].values, los['los_height'].values,\n",
148+
" 'r--', label='LOS ray')\n",
149+
"\n",
150+
"vis = los['visible'].values\n",
151+
"d = los['distance'].values\n",
152+
"ax.scatter(d[vis], los['elevation'].values[vis],\n",
153+
" c='green', s=8, label='Visible', zorder=3)\n",
154+
"ax.scatter(d[~vis], los['elevation'].values[~vis],\n",
155+
" c='red', s=8, label='Blocked', zorder=3)\n",
156+
"\n",
157+
"ax.set_xlabel('Distance')\n",
158+
"ax.set_ylabel('Elevation')\n",
159+
"ax.set_title('Line-of-sight profile')\n",
160+
"ax.legend()\n",
161+
"plt.tight_layout()"
162+
]
163+
},
164+
{
165+
"cell_type": "markdown",
166+
"metadata": {},
167+
"source": [
168+
"## Fresnel zone clearance\n",
169+
"\n",
170+
"For radio link planning, check whether the first Fresnel zone is clear at 900 MHz."
171+
]
172+
},
173+
{
174+
"cell_type": "code",
175+
"execution_count": null,
176+
"metadata": {},
177+
"outputs": [],
178+
"source": [
179+
"los_f = line_of_sight(terrain, x0=ox, y0=oy, x1=tx, y1=ty,\n",
180+
" observer_elev=50, target_elev=50,\n",
181+
" frequency_mhz=900)\n",
182+
"\n",
183+
"fig, ax = plt.subplots(figsize=(12, 4))\n",
184+
"d = los_f['distance'].values\n",
185+
"ax.fill_between(d, 0, los_f['elevation'].values,\n",
186+
" color='sienna', alpha=0.4, label='Terrain')\n",
187+
"ax.plot(d, los_f['los_height'].values, 'r--', label='LOS ray')\n",
188+
"\n",
189+
"# Fresnel zone envelope\n",
190+
"lh = los_f['los_height'].values\n",
191+
"fr = los_f['fresnel_radius'].values\n",
192+
"ax.fill_between(d, lh - fr, lh + fr, alpha=0.15, color='blue',\n",
193+
" label='1st Fresnel zone')\n",
194+
"\n",
195+
"fc = los_f['fresnel_clear'].values\n",
196+
"ax.scatter(d[fc], los_f['elevation'].values[fc],\n",
197+
" c='green', s=8, label='Fresnel clear', zorder=3)\n",
198+
"ax.scatter(d[~fc], los_f['elevation'].values[~fc],\n",
199+
" c='red', s=8, label='Fresnel blocked', zorder=3)\n",
200+
"\n",
201+
"ax.set_xlabel('Distance')\n",
202+
"ax.set_ylabel('Elevation')\n",
203+
"ax.set_title('Fresnel zone clearance (900 MHz)')\n",
204+
"ax.legend()\n",
205+
"plt.tight_layout()"
206+
]
207+
}
208+
],
209+
"metadata": {
210+
"kernelspec": {
211+
"display_name": "Python 3",
212+
"language": "python",
213+
"name": "python3"
214+
},
215+
"language_info": {
216+
"name": "python",
217+
"version": "3.10.0"
218+
}
219+
},
220+
"nbformat": 4,
221+
"nbformat_minor": 4
222+
}

0 commit comments

Comments
 (0)