-
Notifications
You must be signed in to change notification settings - Fork 3
/
srtm_path.py
131 lines (103 loc) · 5.29 KB
/
srtm_path.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
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
import numpy as np
from pathlib import Path
from scipy.interpolate import RectBivariateSpline
from tkinter.messagebox import showerror, showwarning
# folder name where hgt files are located:
hgtfolder = 'hgt'
hgtpath = Path.joinpath(Path.cwd(), hgtfolder)
x = np.linspace(0, 1, 3601, dtype=np.dtype('float32'))
# rewrited with showerror instead of print
class MissingHgtError(Exception):
pass
def get_path_profile(phi_t, psi_t, phi_r, psi_r, coord_samples=500, fill_missing=False):
[dist, atrdeg, di, delta, phi_values, psi_values] = get_path_geometry(phi_t, psi_t, phi_r, psi_r, coord_samples)
subpaths, ufnames, flagsmissing = get_path_sections(phi_values, psi_values, fill_missing)
hi = get_elevations(subpaths, ufnames, flagsmissing)
if not fill_missing:
if flagsmissing.any():
missing = np.array2string(ufnames[np.nonzero(flagsmissing)])
showwarning('Missing HGT files', f"Missing hgt files and filled with 0 values:\n{missing}")
# print(f'Distance: {dist:.5f}km\nBearing angle: {atrdeg:.1f}°\nMax elevation:{np.max(hi):.0f}m')
return dist, atrdeg, di, hi, delta
def get_elevations(subpaths, ufnames, flagsmissing):
interpvalues = np.empty(0, dtype=int)
for i in range(len(subpaths)):
# open hgt filename according to subpath
path = Path.joinpath(hgtpath, ufnames[i])
if not flagsmissing[i]:
with open(path, 'rb') as data:
elevationdata = np.fromfile(data, np.dtype('>i2'), 3601 ** 2).reshape(3601, 3601)
# interpolate elevation values with scipy.interpolate.griddata
f1 = RectBivariateSpline(x, x, elevationdata, kx=1, ky=1)
latcorr = int(ufnames[i][1:3]) + 1
loncorr = int(ufnames[i][4:7])
interpvalues = np.hstack((interpvalues, f1.ev(latcorr - subpaths[i][0], subpaths[i][1] - loncorr)))
else:
interpvalues = np.hstack((interpvalues, np.zeros_like(subpaths[i][0])))
return interpvalues
def get_path_sections(phi_values, psi_values, fill_missing):
fnames = get_hgt_names(phi_values, psi_values)
# check if all hgt files exist in hgt directory
_, indices = np.unique(fnames, return_index=True)
indices = np.sort(indices)
ufnames = fnames[indices]
flagsmissing = np.full_like(indices, 0)
for i in range(len(indices)):
file = fnames[indices[i]]
path = Path.joinpath(hgtpath, file)
flagsmissing[i] = not path.exists()
if not fill_missing:
if flagsmissing.any():
missing = np.array2string(ufnames[np.nonzero(flagsmissing)])
showerror('Missing HGT files', f"Missing hgt files from hgt directory: {missing}")
raise MissingHgtError(f"Missing hgt files from hgt directory:\n{missing}")
# create sub-paths for each difference hgt file:
subpaths = [None] * len(indices)
for i in range(len(indices) - 1):
startidx = indices[i]
endidx = indices[i+1]
subpaths[i] = (phi_values[startidx:endidx], psi_values[startidx:endidx])
# subpaths[i] = list(zip(phi_values[startidx:endidx], psi_values[startidx:endidx]))
subpaths[-1] = (phi_values[indices[-1]:], psi_values[indices[-1]:])
# subpaths[-1] = list(zip(phi_values[indices[-1]:], psi_values[indices[-1]:]))
return subpaths, ufnames, flagsmissing
def get_hgt_names(phi_values, psi_values):
n_s = np.where(phi_values >= 0, 'N', 'S')
e_w = np.where(psi_values >= 0, 'E', 'W')
lat = abs(phi_values).astype(int).astype(str)
lat = np.char.zfill(lat, 2)
lon = abs(psi_values).astype(int).astype(str)
lon = np.char.zfill(lon, 3)
fnames = np.char.add(np.char.add(n_s, lat), np.char.add(e_w, lon))
fnames = np.char.add(fnames, '.hgt')
return fnames
def get_path_geometry(phi_t, psi_t, phi_r, psi_r, coord_samples):
phi_t = np.deg2rad(phi_t)
psi_t = np.deg2rad(psi_t)
phi_r = np.deg2rad(phi_r)
psi_r = np.deg2rad(psi_r)
# The angle subtended by the path at the centre of the Earth, δ, from the stations’ geographic
# coordinates using eq (65):
delta = np.arccos(np.sin(phi_t) * np.sin(phi_r) +
np.cos(phi_t) * np.cos(phi_r) *
np.cos(psi_t - psi_r))
# great circle distance:
R = 6371.009
dist = R * delta # km
# bearing (azimuthal direction clockwise from true North) from station t to station r
# using eq 67:
atr = np.arccos((np.sin(phi_r) - np.sin(phi_t) * np.cos(delta)) /
(np.sin(delta) * np.cos(phi_t)))
if psi_t > psi_r:
atr = 2 * np.pi - atr
atrdeg = np.rad2deg(atr)
if coord_samples == 0:
return dist, atrdeg, None, None, None, None
di = np.linspace(0, dist, coord_samples)
phi_values = np.arcsin(np.sin(phi_t) * np.cos(di / R) + np.cos(phi_t) * np.sin(di / R) * np.cos(atr))
psi_values = psi_t + np.arctan2(np.sin(atr) * np.sin(di / R) * np.cos(phi_t),
np.cos(di / R) - np.sin(phi_t) * np.sin(phi_values))
psi_values = np.remainder(psi_values + 3 * np.pi, (2 * np.pi) - np.pi)
phi_values = np.rad2deg(phi_values)
psi_values = np.rad2deg(psi_values)
return [dist, atrdeg, di, delta, phi_values, psi_values]