-
Notifications
You must be signed in to change notification settings - Fork 1
/
test_horizon.py
71 lines (52 loc) · 1.73 KB
/
test_horizon.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
# -*- coding: utf-8 -*-
"""
Created on Tue Jun 12 16:56:36 2018
@author: slauniai
"""
import numpy as np
import matplotlib.pyplot as plt
import dem_horizon as h
def run_demo():
import numpy as np
import matplotlib.pyplot as plt
import dem_horizon as h
# import dem example
demfile = r'sve_1_dem_16m_aggr.asc'
dem, lat, lon, _ = h.read_AsciiGrid(demfile)
#plt.figure(1)
#plt.imshow(dem)
#plt.xlabel('Lon ix') # column index
#plt.ylabel('Lat ix') # row index
#cb = plt.colorbar()
#cb.set_label('Elev [m]')
# set parameters
dx = 16.0 # dem resolution (m)
view_dist = 500.0 # viewing distance
sectors = 36 # sectors = 360 (deg) / azimuth interval (deg)
# create Horizon object
H = h.Horizon(dem, dx=16)
print('rows, cols', H.shape)
# let's compute horizon for point Lat_ix = 40, Lon_ix 50
lat_ix = 140
lon_ix = 60
P0 = h.Point(lat_ix, lon_ix, dem[lat_ix, lon_ix])
# calculate horizon height angle (deg)
Az = np.linspace(0, 360, sectors)
Az, elev_angle, terrain_profile = H.calc_horizon(P0, R=view_dist, Az_deg=Az, figs=True)
# plot elevation profiles for few directions
Az_ix = np.arange(0, len(Az), 5)
fig = plt.figure()
for k in range(len(Az_ix)):
y = (terrain_profile['elev'][Az_ix[k]])
x = np.arange(0, len(y)) * dx
s = '%.1f deg' % terrain_profile['Az'][Az_ix[k]]
plt.plot(x, y, '-', label=s)
plt.legend()
plt.title('terrain profiles')
plt.ylabel('elev [m]')
plt.xlabel('distance from target [m]')
fig.savefig('horizon.png')
fig.show()
return terrain_profile
if __name__ == "__main__":
run_demo()