-
Notifications
You must be signed in to change notification settings - Fork 35
/
convert5gmv1ToChannels.py
183 lines (166 loc) · 8.27 KB
/
convert5gmv1ToChannels.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
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
'''
Will parse all database and create numpy arrays that represent all channels in the database.
Specificities: some episodes do not have all scenes. And some scenes do not have all receivers.
Assuming Ne episodes, with Ns scenes each, and Nr receivers (given there was only one transmitter),
there are Ne x Ns x Nr channel matrices and each must represent L=25 rays.
With Ne=119, Ns=50, Nr=10, we have 59500 matrices with 25 rays. It is better to save
each episode in one file, with the matrix given by
scene 1:Ns x Tx_index x Rx_index x numberRays and 7 numbers, the following for each ray
path_gain
timeOfArrival
departure_elevation
departure_azimuth
arrival_elevation
arrival_azimuth
isLOS
to simplify we assume that all episodes have the same number of scenes (e.g. 50) and receivers (e.g. 10).
Episodes, scenes, etc, start counting from 0 (not 1).
'''
import datetime
import numpy as np
from shapely import geometry
#from matplotlib import pyplot as plt
import h5py
import matplotlib.pyplot as plt
from rwisimulation.positionmatrix import position_matrix_per_object_shape, calc_position_matrix
#from rwisimulation.calcrxpower import calc_rx_power
from rwisimulation.datamodel import save5gmdata as fgdb
class c:
#analysis_area = (648, 348, 850, 685)
analysis_area = (744, 429, 767, 679) #rosslyn
#analysis_area = (700, 600, 30, 18) #china
analysis_area_resolution = 0.5
analysis_polygon = geometry.Polygon([(c.analysis_area[0], c.analysis_area[1]),
(c.analysis_area[2], c.analysis_area[1]),
(c.analysis_area[2], c.analysis_area[3]),
(c.analysis_area[0], c.analysis_area[3])])
session = fgdb.Session()
totalNumEpisodes = session.query(fgdb.Episode).count()
#pm_per_object_shape = position_matrix_per_object_shape(c.analysis_area, c.analysis_area_resolution)
#print(pm_per_object_shape)
# just to report time
start = datetime.datetime.today()
perc_done = None
#Active if will restrict the analysis area
Use_analysis_polygon = False
#Generate npz file
use_npz = False
#if needed, manually create the output folder
fileNamePrefix = './insite_data_s007_carrie2.8GHz/beijing_mobile_2.8GHz_ts1s_VP' #prefix of output files
pythonExtension = '.npz'
matlabExtension = '.hdf5'
# assume 40 scenes per episode, 10 receivers per scene
numScenesPerEpisode = 40
numTxRxPairsPerScene = 10
numRaysPerTxRxPair = 100
numVariablePerRay = 7+1 #has the ray angle now
#plt.ion()
numEpisode = 0
numLOS = 0
numNLOS = 0
for ep in session.query(fgdb.Episode): #go over all episodes
print('Processing ', ep.number_of_scenes, ' scenes in episode ', ep.insite_pah,)
print('Start time = ', ep.simulation_time_begin, ' and sampling period = ', ep.sampling_time, ' seconds')
print('Episode: ' + str(numEpisode) + ' out of ' + str(totalNumEpisodes))
#initialization
#Ns x [Tx_index x Rx_index x numberRays] and 7 numbers, the following for each ray
allEpisodeData = np.zeros((numScenesPerEpisode, numTxRxPairsPerScene, numRaysPerTxRxPair,
numVariablePerRay), np.float32)
allEpisodeData.fill(np.nan)
#from the first scene, get all receiver names
rec_name_to_array_idx_map = [obj.name for obj in ep.scenes[0].objects if len(obj.receivers) > 0]
print(rec_name_to_array_idx_map)
#process each scene in this episode
#count # of ep.scenes
for sc_i, sc in enumerate(ep.scenes):
#print('Processing scene # ', sc_i)
polygon_list = []
polygon_z = []
polygons_of_interest_idx_list = []
rec_present = []
for obj in sc.objects:
if len(obj.receivers) == 0:
continue #do not process objects that are not receivers
obj_polygon = geometry.asMultiPoint(obj.vertice_array[:,(0,1)]).convex_hull
# check if object is inside the analysis_area
if not Use_analysis_polygon:
# if the object is a receiver and is within the analysis area
if len(obj.receivers) > 0:
rec_array_idx = rec_name_to_array_idx_map.index(obj.name)
for rec in obj.receivers: #for all receivers
ray_i = 0
isLOSChannel = 0
for ray in rec.rays: #for all rays
#gather all info
thisRayInfo = np.zeros(numVariablePerRay)
thisRayInfo[0] = ray.path_gain
thisRayInfo[1] = ray.time_of_arrival
thisRayInfo[2] = ray.departure_elevation
thisRayInfo[3] = ray.departure_azimuth
thisRayInfo[4] = ray.arrival_elevation
thisRayInfo[5] = ray.arrival_azimuth
thisRayInfo[6] = ray.is_los
thisRayInfo[7] = ray.phaseInDegrees
#allEpisodeData = np.zeros((numScenesPerEpisode, numTxRxPairsPerScene,
# numRaysPerTxRxPair, numVariablePerRay), np.float32)
allEpisodeData[sc_i][rec_array_idx][ray_i]=thisRayInfo
ray_i += 1
if ray.is_los == 1:
isLOSChannel = True #if one ray is LOS, the channel is
#print('AK:',sc_i, rec_array_idx)
if isLOSChannel == True:
numLOS += 1
else:
numNLOS += 1
# just for reporting spent time
elif obj_polygon.within(analysis_polygon):
if len(obj.receivers) > 0:
rec_array_idx = rec_name_to_array_idx_map.index(obj.name)
for rec in obj.receivers: #for all receivers
ray_i = 0
isLOSChannel = 0
for ray in rec.rays: #for all rays
#gather all info
thisRayInfo = np.zeros(numVariablePerRay)
thisRayInfo[0] = ray.path_gain
thisRayInfo[1] = ray.time_of_arrival
thisRayInfo[2] = ray.departure_elevation
thisRayInfo[3] = ray.departure_azimuth
thisRayInfo[4] = ray.arrival_elevation
thisRayInfo[5] = ray.arrival_azimuth
thisRayInfo[6] = ray.is_los
thisRayInfo[7] = ray.phaseInDegrees
#allEpisodeData = np.zeros((numScenesPerEpisode, numTxRxPairsPerScene,
# numRaysPerTxRxPair, numVariablePerRay), np.float32)
allEpisodeData[sc_i][rec_array_idx][ray_i]=thisRayInfo
ray_i += 1
if ray.is_los == 1:
isLOSChannel = True #if one ray is LOS, the channel is
#print('AK:',sc_i, rec_array_idx)
if isLOSChannel == True:
numLOS += 1
else:
numNLOS += 1
# just for reporting spent time
perc_done = ((sc_i + 1) / ep.number_of_scenes) * 100
elapsed_time = datetime.datetime.today() - start
time_p_perc = elapsed_time / perc_done
print('\r Done: {:.2f}% Scene: {} time per scene: {} time to finish: {}'.format(
perc_done,
sc_i + 1,
elapsed_time / (sc_i + 1),
time_p_perc * (100 - perc_done)), end='')
if use_npz:
print()
outputFileName = fileNamePrefix + '_e' + str(numEpisode) + pythonExtension
np.savez(outputFileName, allEpisodeData=allEpisodeData)
print('==> Wrote file ' + outputFileName)
outputFileName = fileNamePrefix + '_e' + str(numEpisode) + matlabExtension
print('==> Wrote file ' + outputFileName)
f = h5py.File(outputFileName, 'w')
f['allEpisodeData'] = allEpisodeData
f.close()
numEpisode += 1 #increment episode counter
print('numLOS = ', numLOS)
print('numNLOS = ', numNLOS)
print('Sum = ', numLOS + numNLOS)