##// END OF EJS Templates
schain-jc
RSS
-Functional HDF5 Format Writing and Reading Unit...
Julio Valdez -
r804:4c41460abf8c
parent child
Show More
Show file before | Show file after | Hide comments
@@ -1150,7 +1150,7 class Parameters(JROData):
1150 1150
1151 1151 self.type = "Parameters"
1152 1152
1153 def getTimeRange1(self):
1153 def getTimeRange1(self, interval):
1154 1154
1155 1155 datatime = []
1156 1156
@@ -1162,7 +1162,7 class Parameters(JROData):
1162 1162 # datatime.append(self.utctimeInit)
1163 1163 # datatime.append(self.utctimeInit + self.outputInterval)
1164 1164 datatime.append(time1)
1165 datatime.append(time1 + self.outputInterval)
1165 datatime.append(time1 + interval)
1166 1166
1167 1167 datatime = numpy.array(datatime)
1168 1168
Show file before | Show file after | Hide comments
@@ -251,7 +251,7 class SkyMapPlot(Figure):
251 251 self.addAxes(1, 1, 0, 0, 1, 1, True)
252 252
253 253 def run(self, dataOut, id, wintitle="", channelList=None, showprofile=False,
254 tmin=None, tmax=None, timerange=None,
254 tmin=0, tmax=24, timerange=None,
255 255 save=False, figpath='./', figfile=None, show=True, ftp=False, wr_period=1,
256 256 server=None, folder=None, username=None, password=None,
257 257 ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0, realtime=False):
@@ -272,19 +272,19 class SkyMapPlot(Figure):
272 272 zmax : None
273 273 """
274 274
275 arrayParameters = dataOut.data_param[0,:]
275 arrayParameters = dataOut.data_param
276 276 error = arrayParameters[:,-1]
277 277 indValid = numpy.where(error == 0)[0]
278 278 finalMeteor = arrayParameters[indValid,:]
279 finalAzimuth = finalMeteor[:,4]
280 finalZenith = finalMeteor[:,5]
279 finalAzimuth = finalMeteor[:,3]
280 finalZenith = finalMeteor[:,4]
281 281
282 282 x = finalAzimuth*numpy.pi/180
283 283 y = finalZenith
284 x1 = dataOut.getTimeRange()
284 x1 = [dataOut.ltctime, dataOut.ltctime]
285 285
286 286 #thisDatetime = dataOut.datatime
287 thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[0])
287 thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.ltctime)
288 288 title = wintitle + " Parameters"
289 289 xlabel = "Zonal Zenith Angle (deg) "
290 290 ylabel = "Meridional Zenith Angle (deg)"
@@ -439,22 +439,13 class WindProfilerPlot(Figure):
439 439 zmax : None
440 440 """
441 441
442 if channelList == None:
443 channelIndexList = dataOut.channelIndexList
444 else:
445 channelIndexList = []
446 for channel in channelList:
447 if channel not in dataOut.channelList:
448 raise ValueError, "Channel %d is not in dataOut.channelList"
449 channelIndexList.append(dataOut.channelList.index(channel))
450
451 442 # if timerange is not None:
452 443 # self.timerange = timerange
453 444 #
454 445 # tmin = None
455 446 # tmax = None
456 447
457 x = dataOut.getTimeRange1()
448 x = dataOut.getTimeRange1(dataOut.outputInterval)
458 449 # y = dataOut.heightList
459 450 y = dataOut.heightList
460 451
@@ -480,7 +471,7 class WindProfilerPlot(Figure):
480 471
481 472 # showprofile = False
482 473 # thisDatetime = dataOut.datatime
483 thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[0])
474 thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.ltctime)
484 475 title = wintitle + "Wind"
485 476 xlabel = ""
486 477 ylabel = "Range (Km)"
@@ -707,7 +698,7 class ParametersPlot(Figure):
707 698 if parameterIndex == None:
708 699 parameterIndex = 1
709 700
710 x = dataOut.getTimeRange1()
701 x = dataOut.getTimeRange1(dataOut.paramInterval)
711 702 y = dataOut.heightList
712 703 z = data_param[channelIndexList,parameterIndex,:].copy()
713 704
@@ -1100,7 +1091,7 class EWDriftsPlot(Figure):
1100 1091 tmin = None
1101 1092 tmax = None
1102 1093
1103 x = dataOut.getTimeRange1()
1094 x = dataOut.getTimeRange1(dataOut.outputInterval)
1104 1095 # y = dataOut.heightList
1105 1096 y = dataOut.heightList
1106 1097
@@ -1280,7 +1271,7 class PhasePlot(Figure):
1280 1271
1281 1272 tmin = None
1282 1273 tmax = None
1283 x = dataOut.getTimeRange1()
1274 x = dataOut.getTimeRange1(dataOut.outputInterval)
1284 1275 y = dataOut.getHeiRange()
1285 1276
1286 1277
Show file before | Show file after | Hide comments
This diff has been collapsed as it changes many lines, (728 lines changed) Show them Hide them
@@ -3,13 +3,29 import time
3 3 import os
4 4 import h5py
5 5 import re
6 import datetime
6 7
7 8 from schainpy.model.data.jrodata import *
8 9 from schainpy.model.proc.jroproc_base import ProcessingUnit, Operation
10 # from jroIO_base import *
9 11 from schainpy.model.io.jroIO_base import *
12 import schainpy
10 13
11 14
12 15 class HDF5Reader(ProcessingUnit):
16 '''
17 Reads HDF5 format files
18
19 path
20
21 startDate
22
23 endDate
24
25 startTime
26
27 endTime
28 '''
13 29
14 30 ext = ".hdf5"
15 31
@@ -17,15 +33,19 class HDF5Reader(ProcessingUnit):
17 33
18 34 timezone = None
19 35
20 secStart = None
36 startTime = None
21 37
22 secEnd = None
38 endTime = None
23 39
24 40 fileIndex = None
25 41
26 blockIndex = None
42 utcList = None #To select data in the utctime list
27 43
28 blocksPerFile = None
44 blockList = None #List to blocks to be read from the file
45
46 blocksPerFile = None #Number of blocks to be read
47
48 blockIndex = None
29 49
30 50 path = None
31 51
@@ -37,10 +57,6 class HDF5Reader(ProcessingUnit):
37 57
38 58 #Hdf5 File
39 59
40 fpMetadata = None
41
42 pathMeta = None
43
44 60 listMetaname = None
45 61
46 62 listMeta = None
@@ -57,153 +73,80 class HDF5Reader(ProcessingUnit):
57 73
58 74 dataOut = None
59 75
60 nRecords = None
61
62 76
63 77 def __init__(self):
64 self.dataOut = self.__createObjByDefault()
78 self.dataOut = Parameters()
65 79 return
66
67 def __createObjByDefault(self):
68
69 dataObj = Parameters()
70
71 return dataObj
72
73 def setup(self,path=None,
74 startDate=None,
75 endDate=None,
76 startTime=datetime.time(0,0,0),
77 endTime=datetime.time(23,59,59),
78 walk=True,
79 timezone='ut',
80 all=0,
81 online=False,
82 ext=None):
83
84 if ext==None:
85 ext = self.ext
86 self.timezone = timezone
87 # self.all = all
88 # self.online = online
89 self.path = path
90
91 startDateTime = datetime.datetime.combine(startDate,startTime)
92 endDateTime = datetime.datetime.combine(endDate,endTime)
93 secStart = (startDateTime-datetime.datetime(1970,1,1)).total_seconds()
94 secEnd = (endDateTime-datetime.datetime(1970,1,1)).total_seconds()
95
96 self.secStart = secStart
97 self.secEnd = secEnd
98
99 if not(online):
100 #Busqueda de archivos offline
101 self.__searchFilesOffline(path, startDate, endDate, ext, startTime, endTime, secStart, secEnd, walk)
102 else:
103 self.__searchFilesOnline(path, walk)
104 80
105 if not(self.filenameList):
81 def setup(self, **kwargs):
82
83 path = kwargs['path']
84 startDate = kwargs['startDate']
85 endDate = kwargs['endDate']
86 startTime = kwargs['startTime']
87 endTime = kwargs['endTime']
88 walk = kwargs['walk']
89 if kwargs.has_key('ext'):
90 ext = kwargs['ext']
91 else:
92 ext = '.hdf5'
93
94 print "[Reading] Searching files in offline mode ..."
95 pathList, filenameList = self.__searchFilesOffLine(path, startDate=startDate, endDate=endDate,
96 startTime=startTime, endTime=endTime,
97 ext=ext, walk=walk)
98
99 if not(filenameList):
106 100 print "There is no files into the folder: %s"%(path)
107 101 sys.exit(-1)
108 102
109 # self.__getExpParameters()
110
111 103 self.fileIndex = -1
112
113 self.__setNextFileOffline()
104 self.startTime = startTime
105 self.endTime = endTime
114 106
115 107 self.__readMetadata()
116 108
117 self.blockIndex = 0
109 self.__setNextFileOffline()
118 110
119 111 return
120
121 def __searchFilesOffline(self,
112
113 def __searchFilesOffLine(self,
122 114 path,
123 startDate,
124 endDate,
125 ext,
115 startDate=None,
116 endDate=None,
126 117 startTime=datetime.time(0,0,0),
127 118 endTime=datetime.time(23,59,59),
128 secStart = 0,
129 secEnd = numpy.inf,
119 ext='.hdf5',
130 120 walk=True):
131 121
132 # self.__setParameters(path, startDate, endDate, startTime, endTime, walk)
133 #
134 # self.__checkPath()
135 #
136 # self.__findDataForDates()
137 #
138 # self.__selectDataForTimes()
139 #
140 # for i in range(len(self.filenameList)):
141 # print "%s" %(self.filenameList[i])
122 expLabel = ''
123 self.filenameList = []
124 self.datetimeList = []
142 125
143 126 pathList = []
144 127
145 if not walk:
146 #pathList.append(path)
147 multi_path = path.split(',')
148 for single_path in multi_path:
149 pathList.append(single_path)
150
151 else:
152 #dirList = []
153 multi_path = path.split(',')
154 for single_path in multi_path:
155 dirList = []
156 for thisPath in os.listdir(single_path):
157 if not os.path.isdir(os.path.join(single_path,thisPath)):
158 continue
159 if not isDoyFolder(thisPath):
160 continue
128 JRODataObj = JRODataReader()
129 dateList, pathList = JRODataObj.findDatafiles(path, startDate, endDate, expLabel, ext, walk, include_path=True)
130
131 if dateList == []:
132 print "[Reading] No *%s files in %s from %s to %s)"%(ext, path,
133 datetime.datetime.combine(startDate,startTime).ctime(),
134 datetime.datetime.combine(endDate,endTime).ctime())
161 135
162 dirList.append(thisPath)
163
164 if not(dirList):
165 return None, None
166
167 thisDate = startDate
168
169 while(thisDate <= endDate):
170 year = thisDate.timetuple().tm_year
171 doy = thisDate.timetuple().tm_yday
172
173 matchlist = fnmatch.filter(dirList, '?' + '%4.4d%3.3d' % (year,doy) + '*')
174 if len(matchlist) == 0:
175 thisDate += datetime.timedelta(1)
176 continue
177 for match in matchlist:
178 pathList.append(os.path.join(single_path,match))
179
180 thisDate += datetime.timedelta(1)
181
182 if pathList == []:
183 print "Any folder was found for the date range: %s-%s" %(startDate, endDate)
184 136 return None, None
185 137
186 print "%d folder(s) was(were) found for the date range: %s - %s" %(len(pathList), startDate, endDate)
187
138 if len(dateList) > 1:
139 print "[Reading] %d days were found in date range: %s - %s" %(len(dateList), startDate, endDate)
140 else:
141 print "[Reading] data was found for the date %s" %(dateList[0])
142
188 143 filenameList = []
189 144 datetimeList = []
190 pathDict = {}
191 filenameList_to_sort = []
192
193 for i in range(len(pathList)):
194
195 thisPath = pathList[i]
196
197 fileList = glob.glob1(thisPath, "*%s" %ext)
198 fileList.sort()
199 pathDict.setdefault(fileList[0])
200 pathDict[fileList[0]] = i
201 filenameList_to_sort.append(fileList[0])
202 145
203 filenameList_to_sort.sort()
146 #----------------------------------------------------------------------------------
204 147
205 for file in filenameList_to_sort:
206 thisPath = pathList[pathDict[file]]
148 for thisPath in pathList:
149 # thisPath = pathList[pathDict[file]]
207 150
208 151 fileList = glob.glob1(thisPath, "*%s" %ext)
209 152 fileList.sort()
@@ -211,7 +154,11 class HDF5Reader(ProcessingUnit):
211 154 for file in fileList:
212 155
213 156 filename = os.path.join(thisPath,file)
214 thisDatetime = self.__isFileinThisTime(filename, secStart, secEnd)
157
158 if not isFileInDateRange(filename, startDate, endDate):
159 continue
160
161 thisDatetime = self.__isFileInTimeRange(filename, startDate, endDate, startTime, endTime)
215 162
216 163 if not(thisDatetime):
217 164 continue
@@ -220,27 +167,32 class HDF5Reader(ProcessingUnit):
220 167 datetimeList.append(thisDatetime)
221 168
222 169 if not(filenameList):
223 print "Any file was found for the time range %s - %s" %(startTime, endTime)
170 print "[Reading] Any file was found int time range %s - %s" %(datetime.datetime.combine(startDate,startTime).ctime(), datetime.datetime.combine(endDate,endTime).ctime())
224 171 return None, None
225 172
226 print "%d file(s) was(were) found for the time range: %s - %s" %(len(filenameList), startTime, endTime)
173 print "[Reading] %d file(s) was(were) found in time range: %s - %s" %(len(filenameList), startTime, endTime)
227 174 print
228 175
229 176 for i in range(len(filenameList)):
230 print "%s -> [%s]" %(filenameList[i], datetimeList[i].ctime())
177 print "[Reading] %s -> [%s]" %(filenameList[i], datetimeList[i].ctime())
231 178
232 179 self.filenameList = filenameList
233 180 self.datetimeList = datetimeList
234
181
235 182 return pathList, filenameList
236
237 def __isFileinThisTime(self, filename, startSeconds, endSeconds):
183
184 def __isFileInTimeRange(self,filename, startDate, endDate, startTime, endTime):
185
238 186 """
239 187 Retorna 1 si el archivo de datos se encuentra dentro del rango de horas especificado.
240 188
241 189 Inputs:
242 190 filename : nombre completo del archivo de datos en formato Jicamarca (.r)
243 191
192 startDate : fecha inicial del rango seleccionado en formato datetime.date
193
194 endDate : fecha final del rango seleccionado en formato datetime.date
195
244 196 startTime : tiempo inicial del rango seleccionado en formato datetime.time
245 197
246 198 endTime : tiempo final del rango seleccionado en formato datetime.time
@@ -254,43 +206,61 class HDF5Reader(ProcessingUnit):
254 206 Si la cabecera no puede ser leida.
255 207
256 208 """
257
209
258 210 try:
259 fp = fp = h5py.File(filename,'r')
211 fp = h5py.File(filename,'r')
212 grp1 = fp['Data']
213
260 214 except IOError:
261 215 traceback.print_exc()
262 216 raise IOError, "The file %s can't be opened" %(filename)
217 #chino rata
218 #In case has utctime attribute
219 grp2 = grp1['utctime']
220 # thisUtcTime = grp2.value[0] - 5*3600 #To convert to local time
221 thisUtcTime = grp2.value[0]
222
223 fp.close()
263 224
264 grp = fp['Data']
265 timeAux = grp['time']
266 time0 = timeAux[:][0].astype(numpy.float) #Time Vector
225 thisDatetime = datetime.datetime.fromtimestamp(thisUtcTime[0])
226 thisDate = thisDatetime.date()
227 thisTime = thisDatetime.time()
267 228
268 fp.close()
229 startUtcTime = (datetime.datetime.combine(thisDate,startTime)- datetime.datetime(1970, 1, 1)).total_seconds()
230 endUtcTime = (datetime.datetime.combine(thisDate,endTime)- datetime.datetime(1970, 1, 1)).total_seconds()
269 231
270 if self.timezone == 'lt':
271 time0 -= 5*3600
272
273 boolTimer = numpy.logical_and(time0 >= startSeconds,time0 < endSeconds)
274
275 if not (numpy.any(boolTimer)):
232 #General case
233 # o>>>>>>>>>>>>>><<<<<<<<<<<<<<o
234 #-----------o----------------------------o-----------
235 # startTime endTime
236
237 if endTime >= startTime:
238 thisUtcLog = numpy.logical_and(thisUtcTime > startUtcTime, thisUtcTime < endUtcTime)
239 if numpy.any(thisUtcLog): #If there is one block between the hours mentioned
240 return thisDatetime
241 return None
242
243 #If endTime < startTime then endTime belongs to the next day
244 #<<<<<<<<<<<o o>>>>>>>>>>>
245 #-----------o----------------------------o-----------
246 # endTime startTime
247
248 if (thisDate == startDate) and numpy.all(thisUtcTime < startUtcTime):
249 return None
250
251 if (thisDate == endDate) and numpy.all(thisUtcTime > endUtcTime):
252 return None
253
254 if numpy.all(thisUtcTime < startUtcTime) and numpy.all(thisUtcTime > endUtcTime):
276 255 return None
277 256
278 thisDatetime = datetime.datetime.utcfromtimestamp(time0[0])
279 257 return thisDatetime
280
281 def __checkPath(self):
282 if os.path.exists(self.path):
283 self.status = 1
284 else:
285 self.status = 0
286 print 'Path:%s does not exists'%self.path
287
288 return
289
258
290 259 def __setNextFileOffline(self):
291 idFile = self.fileIndex
292 idFile += 1
293 260
261 self.fileIndex += 1
262 idFile = self.fileIndex
263
294 264 if not(idFile < len(self.filenameList)):
295 265 print "No more Files"
296 266 return 0
@@ -298,43 +268,51 class HDF5Reader(ProcessingUnit):
298 268 filename = self.filenameList[idFile]
299 269
300 270 filePointer = h5py.File(filename,'r')
301
302 self.flagIsNewFile = 1
303 self.fileIndex = idFile
271
304 272 self.filename = filename
305 273
306 274 self.fp = filePointer
307 275
308 276 print "Setting the file: %s"%self.filename
309 277
310 self.__readMetadata()
278 # self.__readMetadata()
311 279 self.__setBlockList()
280 self.__readData()
312 281 # self.nRecords = self.fp['Data'].attrs['blocksPerFile']
313 self.nRecords = self.fp['Data'].attrs['nRecords']
282 # self.nRecords = self.fp['Data'].attrs['nRecords']
314 283 self.blockIndex = 0
315 284 return 1
316 285
317 286 def __setBlockList(self):
318 287 '''
288 Selects the data within the times defined
289
319 290 self.fp
320 self.startDateTime
321 self.endDateTime
291 self.startTime
292 self.endTime
322 293
323 294 self.blockList
324 295 self.blocksPerFile
325 296
326 297 '''
327 filePointer = self.fp
328 secStart = self.secStart
329 secEnd = self.secEnd
298 fp = self.fp
299 startTime = self.startTime
300 endTime = self.endTime
330 301
331 grp = filePointer['Data']
332 timeVector = grp['time'].value.astype(numpy.float)[0]
302 grp = fp['Data']
303 thisUtcTime = grp['utctime'].value.astype(numpy.float)[0]
333 304
334 305 if self.timezone == 'lt':
335 timeVector -= 5*3600
306 thisUtcTime -= 5*3600
307
308 thisDatetime = datetime.datetime.fromtimestamp(thisUtcTime[0])
309 thisDate = thisDatetime.date()
310 thisTime = thisDatetime.time()
336 311
337 ind = numpy.where(numpy.logical_and(timeVector >= secStart , timeVector < secEnd))[0]
312 startUtcTime = (datetime.datetime.combine(thisDate,startTime) - datetime.datetime(1970, 1, 1)).total_seconds()
313 endUtcTime = (datetime.datetime.combine(thisDate,endTime) - datetime.datetime(1970, 1, 1)).total_seconds()
314
315 ind = numpy.where(numpy.logical_and(thisUtcTime >= startUtcTime, thisUtcTime < endUtcTime))[0]
338 316
339 317 self.blockList = ind
340 318 self.blocksPerFile = len(ind)
@@ -343,6 +321,8 class HDF5Reader(ProcessingUnit):
343 321
344 322 def __readMetadata(self):
345 323 '''
324 Reads Metadata
325
346 326 self.pathMeta
347 327
348 328 self.listShapes
@@ -351,41 +331,46 class HDF5Reader(ProcessingUnit):
351 331
352 332 '''
353 333
354 grp = self.fp['Data']
355 pathMeta = os.path.join(self.path, grp.attrs['metadata'])
356
357 if pathMeta == self.pathMeta:
358 return
359 else:
360 self.pathMeta = pathMeta
361
362 filePointer = h5py.File(self.pathMeta,'r')
363 groupPointer = filePointer['Metadata']
334 # grp = self.fp['Data']
335 # pathMeta = os.path.join(self.path, grp.attrs['metadata'])
336 #
337 # if pathMeta == self.pathMeta:
338 # return
339 # else:
340 # self.pathMeta = pathMeta
341 #
342 # filePointer = h5py.File(self.pathMeta,'r')
343 # groupPointer = filePointer['Metadata']
344
345 filename = self.filenameList[0]
346
347 fp = h5py.File(filename,'r')
348
349 gp = fp['Metadata']
364 350
365 351 listMetaname = []
366 352 listMetadata = []
367 for item in groupPointer.items():
353 for item in gp.items():
368 354 name = item[0]
369 355
370 356 if name=='array dimensions':
371 table = groupPointer[name][:]
357 table = gp[name][:]
372 358 listShapes = {}
373 359 for shapes in table:
374 listShapes[shapes[0]] = numpy.array([shapes[1],shapes[2],shapes[3],shapes[4]])
360 listShapes[shapes[0]] = numpy.array([shapes[1],shapes[2],shapes[3],shapes[4],shapes[5]])
375 361 else:
376 data = groupPointer[name].value
362 data = gp[name].value
377 363 listMetaname.append(name)
378 364 listMetadata.append(data)
379 365
380 if name=='type':
381 self.__initDataOut(data)
382
383 filePointer.close()
384
366 # if name=='type':
367 # self.__initDataOut(data)
368
385 369 self.listShapes = listShapes
386 370 self.listMetaname = listMetaname
387 371 self.listMeta = listMetadata
388 372
373 fp.close()
389 374 return
390 375
391 376 def __readData(self):
@@ -395,113 +380,115 class HDF5Reader(ProcessingUnit):
395 380
396 381 for item in grp.items():
397 382 name = item[0]
398
399 if name == 'time':
400 listdataname.append('utctime')
401 timeAux = grp[name].value.astype(numpy.float)[0]
402 listdata.append(timeAux)
403 continue
404
405 383 listdataname.append(name)
406 array = self.__setDataArray(self.nRecords, grp[name],self.listShapes[name])
384
385 array = self.__setDataArray(grp[name],self.listShapes[name])
407 386 listdata.append(array)
408 387
409 388 self.listDataname = listdataname
410 389 self.listData = listdata
411 390 return
412 391
413 def __setDataArray(self, nRecords, dataset, shapes):
392 def __setDataArray(self, dataset, shapes):
393
394 nDims = shapes[0]
414 395
415 nChannels = shapes[0] #Dimension 0
396 nDim2 = shapes[1] #Dimension 0
416 397
417 nPoints = shapes[1] #Dimension 1, number of Points or Parameters
398 nDim1 = shapes[2] #Dimension 1, number of Points or Parameters
418 399
419 nSamples = shapes[2] #Dimension 2, number of samples or ranges
400 nDim0 = shapes[3] #Dimension 2, number of samples or ranges
420 401
421 mode = shapes[3]
402 mode = shapes[4] #Mode of storing
422 403
423 # if nPoints>1:
424 # arrayData = numpy.zeros((nRecords,nChannels,nPoints,nSamples))
425 # else:
426 # arrayData = numpy.zeros((nRecords,nChannels,nSamples))
427 #
428 # chn = 'channel'
429 #
430 # for i in range(nChannels):
431 #
432 # data = dataset[chn + str(i)].value
433 #
434 # if nPoints>1:
435 # data = numpy.rollaxis(data,2)
436 #
437 # arrayData[:,i,:] = data
438
439 arrayData = numpy.zeros((nRecords,nChannels,nPoints,nSamples))
440 doSqueeze = False
441 if mode == 0:
442 strds = 'channel'
443 nDatas = nChannels
444 newShapes = (nRecords,nPoints,nSamples)
445 if nPoints == 1:
446 doSqueeze = True
447 axisSqueeze = 2
448 else:
404 blockList = self.blockList
405
406 blocksPerFile = self.blocksPerFile
407
408 #Depending on what mode the data was stored
409 # if mode == 0: #Divided in channels
410 # strds = 'channel'
411 # nDatas = nDim2
412 # newShapes = (blocksPerFile,nDim1,nDim0)
413 if mode == 1: #Divided in parameter
449 414 strds = 'param'
450 nDatas = nPoints
451 newShapes = (nRecords,nChannels,nSamples)
452 if nChannels == 1:
453 doSqueeze = True
454 axisSqueeze = 1
455
456 for i in range(nDatas):
415 nDatas = nDim1
416 newShapes = (blocksPerFile,nDim2,nDim0)
417 elif mode==2: #Concatenated in a table
418 strds = 'table0'
419 arrayData = dataset[strds].value
420 #Selecting part of the dataset
421 utctime = arrayData[:,0]
422 u, indices = numpy.unique(utctime, return_index=True)
457 423
458 data = dataset[strds + str(i)].value
459 data = data.reshape(newShapes)
460
461 if mode == 0:
462 arrayData[:,i,:,:] = data
463 else:
424 if blockList.size != indices.size:
425 indMin = indices[blockList[0]]
426 indMax = indices[blockList[-1] + 1]
427 arrayData = arrayData[indMin:indMax,:]
428 return arrayData
429
430 #------- One dimension ---------------
431 if nDims == 1:
432 arrayData = dataset.value.astype(numpy.float)[0][blockList]
433
434 #------- Two dimensions -----------
435 elif nDims == 2:
436 arrayData = numpy.zeros((blocksPerFile,nDim1,nDim0))
437 newShapes = (blocksPerFile,nDim0)
438 nDatas = nDim1
439
440 for i in range(nDatas):
441 data = dataset[strds + str(i)].value
442 arrayData[:,i,:] = data[blockList,:]
443
444 #------- Three dimensions ---------
445 else:
446 arrayData = numpy.zeros((blocksPerFile,nDim2,nDim1,nDim0))
447 for i in range(nDatas):
448
449 data = dataset[strds + str(i)].value
450 data = data[blockList,:,:]
451 data = data.reshape(newShapes)
452 # if mode == 0:
453 # arrayData[:,i,:,:] = data
454 # else:
464 455 arrayData[:,:,i,:] = data
465
466 if doSqueeze:
467 arrayData = numpy.squeeze(arrayData, axis=axisSqueeze)
468
456
469 457 return arrayData
470
471 def __initDataOut(self, type):
472 458
473 # if type =='Parameters':
474 # self.dataOut = Parameters()
475 # elif type =='Spectra':
476 # self.dataOut = Spectra()
477 # elif type =='Voltage':
478 # self.dataOut = Voltage()
479 # elif type =='Correlation':
480 # self.dataOut = Correlation()
481
482 return
483
484 459 def __setDataOut(self):
485 460 listMeta = self.listMeta
486 461 listMetaname = self.listMetaname
487 462 listDataname = self.listDataname
488 463 listData = self.listData
464 listShapes = self.listShapes
489 465
490 466 blockIndex = self.blockIndex
491 blockList = self.blockList
467 # blockList = self.blockList
492 468
493 469 for i in range(len(listMeta)):
494 470 setattr(self.dataOut,listMetaname[i],listMeta[i])
495 471
496 472 for j in range(len(listData)):
497 if listDataname[j]=='utctime':
498 # setattr(self.dataOut,listDataname[j],listData[j][blockList[blockIndex]])
499 setattr(self.dataOut,'utctimeInit',listData[j][blockList[blockIndex]])
500 continue
501
502 setattr(self.dataOut,listDataname[j],listData[j][blockList[blockIndex],:])
473 nShapes = listShapes[listDataname[j]][0]
474 mode = listShapes[listDataname[j]][4]
475 if nShapes == 1:
476 setattr(self.dataOut,listDataname[j],listData[j][blockIndex])
477 elif nShapes > 1:
478 setattr(self.dataOut,listDataname[j],listData[j][blockIndex,:])
479 #Mode Meteors
480 elif mode ==2:
481 selectedData = self.__selectDataMode2(listData[j], blockIndex)
482 setattr(self.dataOut, listDataname[j], selectedData)
483 return
484
485 def __selectDataMode2(self, data, blockIndex):
486 utctime = data[:,0]
487 aux, indices = numpy.unique(utctime, return_inverse=True)
488 selInd = numpy.where(indices == blockIndex)[0]
489 selData = data[selInd,:]
503 490
504 return self.dataOut.data_param
491 return selData
505 492
506 493 def getData(self):
507 494
@@ -514,13 +501,11 class HDF5Reader(ProcessingUnit):
514 501 if not( self.__setNextFileOffline() ):
515 502 self.dataOut.flagNoData = True
516 503 return 0
517
518 #
504
519 505 # if self.datablock == None: # setear esta condicion cuando no hayan datos por leers
520 506 # self.dataOut.flagNoData = True
521 507 # return 0
522
523 self.__readData()
508 # self.__readData()
524 509 self.__setDataOut()
525 510 self.dataOut.flagNoData = False
526 511
@@ -540,6 +525,21 class HDF5Reader(ProcessingUnit):
540 525 return
541 526
542 527 class HDF5Writer(Operation):
528 '''
529 HDF5 Writer, stores parameters data in HDF5 format files
530
531 path: path where the files will be stored
532
533 blocksPerFile: number of blocks that will be saved in per HDF5 format file
534
535 mode: selects the data stacking mode: '0' channels, '1' parameters, '3' table (for meteors)
536
537 metadataList: list of attributes that will be stored as metadata
538
539 dataList: list of attributes that will be stores as data
540
541 '''
542
543 543
544 544 ext = ".hdf5"
545 545
@@ -605,9 +605,6 class HDF5Writer(Operation):
605 605
606 606 self.path = kwargs['path']
607 607
608 if kwargs.has_key('ext'):
609 self.ext = kwargs['ext']
610
611 608 if kwargs.has_key('blocksPerFile'):
612 609 self.blocksPerFile = kwargs['blocksPerFile']
613 610 else:
@@ -625,7 +622,7 class HDF5Writer(Operation):
625 622 if type(mode) == int:
626 623 mode = numpy.zeros(len(self.dataList)) + mode
627 624 else:
628 mode = numpy.zeros(len(self.dataList))
625 mode = numpy.ones(len(self.dataList))
629 626
630 627 self.mode = mode
631 628
@@ -641,23 +638,40 class HDF5Writer(Operation):
641 638
642 639 dataAux = getattr(self.dataOut, self.dataList[i])
643 640
644 if type(dataAux)==float or type(dataAux)==int:
641 #--------------------- Conditionals ------------------------
642 #There is no data
643 if dataAux == None:
644 return 0
645
646 #Not array, just a number
647 if type(dataAux)==float or type(dataAux)==int:
645 648 arrayDim[i,0] = 1
649 mode[i] = 0
650
651 #Mode meteors
652 elif mode[i] == 2:
653 arrayDim[i,3] = dataAux.shape[-1]
654 arrayDim[i,4] = mode[i] #Mode the data was stored
655
656 #All the rest
646 657 else:
658 arrayDim0 = dataAux.shape #Data dimensions
659 arrayDim[i,0] = len(arrayDim0) #Number of array dimensions
660 arrayDim[i,4] = mode[i] #Mode the data was stored
647 661
648 if dataAux == None:
649 return 0
650
651 arrayDim0 = dataAux.shape
652 arrayDim[i,0] = len(arrayDim0)
653 arrayDim[i,4] = mode[i]
654
662 # Three-dimension arrays
655 663 if len(arrayDim0) == 3:
656 664 arrayDim[i,1:-1] = numpy.array(arrayDim0)
665
666 # Two-dimension arrays
657 667 elif len(arrayDim0) == 2:
658 arrayDim[i,2:-1] = numpy.array(arrayDim0) #nHeights
668 arrayDim[i,2:-1] = numpy.array(arrayDim0)
669
670 # One-dimension arrays
659 671 elif len(arrayDim0) == 1:
660 672 arrayDim[i,3] = arrayDim0
673
674 # No array, just a number
661 675 elif len(arrayDim0) == 0:
662 676 arrayDim[i,0] = 1
663 677 arrayDim[i,3] = 1
@@ -816,6 +830,7 class HDF5Writer(Operation):
816 830
817 831 for i in range(len(self.dataList)):
818 832
833 #One-dimension data
819 834 if nDims[i]==1:
820 835 # ds0 = grp.create_dataset(self.dataList[i], (1,1), maxshape=(1,self.blocksPerFile) , chunks = True, dtype='S20')
821 836 ds0 = grp.create_dataset(self.dataList[i], (1,1), maxshape=(1,self.blocksPerFile) , chunks = True, dtype=numpy.float64)
@@ -824,16 +839,26 class HDF5Writer(Operation):
824 839 nDimsForDs.append(nDims[i])
825 840 else:
826 841
827 if mode[i]==0:
828 strMode = "channel"
829 nDatas[i] = self.arrayDim[i,1]
830 else:
842 #Channel mode
843 # if mode[i] == 0:
844 # strMode = "channel"
845 #
846 # #nDatas is the number of arrays per variable
847 # if nDims[i] == 1:
848 # nDatas[i] = self.arrayDim[i,1]
849 # elif nDims[i] == 2:
850 # nDatas[i] = self.arrayDim[i,2]
851
852 #Parameters mode
853 if mode[i] == 1:
831 854 strMode = "param"
832 855 nDatas[i] = self.arrayDim[i,2]
833
834 if nDims[i]==2:
835 nDatas[i] = self.arrayDim[i,2]
836
856
857 #Meteors mode
858 elif mode[i] == 2:
859 strMode = "table"
860 nDatas[i] = 1
861
837 862 grp0 = grp.create_group(self.dataList[i])
838 863
839 864 for j in range(int(nDatas[i])):
@@ -841,43 +866,43 class HDF5Writer(Operation):
841 866
842 867 if nDims[i] == 3:
843 868 ds0 = grp0.create_dataset(tableName, (nDim1[i],nDim0[i],1) , data = numpy.zeros((nDim1[i],nDim0[i],1)) ,maxshape=(None,nDim0[i],None), chunks=True)
869
844 870 else:
845 871 ds0 = grp0.create_dataset(tableName, (1,nDim0[i]), data = numpy.zeros((1,nDim0[i])) , maxshape=(None,nDim0[i]), chunks=True)
846 872
847 873 ds.append(ds0)
848 874 data.append([])
849 875 nDimsForDs.append(nDims[i])
876
877 fp.flush()
878 fp.close()
879
850 880 self.nDatas = nDatas
851 881 self.nDims = nDims
852 882 self.nDimsForDs = nDimsForDs
853 883 #Saving variables
854 884 print 'Writing the file: %s'%filename
855 885 self.filename = filename
856 self.fp = fp
857 self.grp = grp
858 self.grp.attrs.modify('nRecords', 1)
886 # self.fp = fp
887 # self.grp = grp
888 # self.grp.attrs.modify('nRecords', 1)
859 889 self.ds = ds
860 890 self.data = data
861
862 self.setFile = setFile
891 #
892 # self.setFile = setFile
863 893 self.firsttime = True
864 894 self.blockIndex = 0
865 895 return
866 896
867 897 def putData(self):
868
869 if not self.firsttime:
870 self.readBlock()
871 898
872 899 if self.blockIndex == self.blocksPerFile or self.dateFlag():
873
874 900 self.setNextFile()
875 901
876 self.setBlock()
877 self.writeBlock()
878
879 self.fp.flush()
880 self.fp.close()
902 # if not self.firsttime:
903 self.readBlock()
904 self.setBlock() #Prepare data to be written
905 self.writeBlock() #Write data
881 906
882 907 return
883 908
@@ -903,11 +928,13 class HDF5Writer(Operation):
903 928 ds[ind] = ds0
904 929 ind += 1
905 930 else:
906 if self.mode[i]==0:
907 strMode = "channel"
908 else:
931 # if self.mode[i] == 0:
932 # strMode = "channel"
933 if self.mode[i] == 1:
909 934 strMode = "param"
910
935 elif self.mode[i] == 2:
936 strMode = "table"
937
911 938 grp0 = grp[self.dataList[i]]
912 939
913 940 for j in range(int(self.nDatas[i])):
@@ -915,8 +942,7 class HDF5Writer(Operation):
915 942 ds0 = grp0[tableName]
916 943 ds[ind] = ds0
917 944 ind += 1
918
919
945
920 946 self.fp = fp
921 947 self.grp = grp
922 948 self.ds = ds
@@ -940,32 +966,24 class HDF5Writer(Operation):
940 966 for i in range(len(self.dataList)):
941 967 dataAux = getattr(self.dataOut,self.dataList[i])
942 968
943 if nDims[i] == 1:
944 # data[ind] = numpy.array([str(dataAux)]).reshape((1,1))
969 if nDims[i] == 1 or mode[i] == 2:
945 970 data[ind] = dataAux
946 # if not self.firsttime:
947 # data[ind] = numpy.hstack((self.ds[ind][:], self.data[ind]))
948 ind += 1
949 else:
971 ind += 1
972
973 elif nDims[i] == 2:
950 974 for j in range(int(nDatas[i])):
951 if (mode[i] == 0) or (nDims[i] == 2): #In case division per channel or Dimensions is only 1
952 data[ind] = dataAux[j,:]
953 else:
954 data[ind] = dataAux[:,j,:]
975 data[ind] = dataAux[j,:]
976 ind += 1
955 977
956 # if nDims[i] == 3:
957 # data[ind] = data[ind].reshape((data[ind].shape[0],data[ind].shape[1],1))
958
959 # if not self.firsttime:
960 # data[ind] = numpy.dstack((self.ds[ind][:], data[ind]))
961
978 elif nDims[i] == 3:
979 for j in range(int(nDatas[i])):
980 # Extinct mode 0
981 # if (mode[i] == 0):
982 # data[ind] = dataAux[j,:,:]
962 983 # else:
963 # data[ind] = data[ind].reshape((1,data[ind].shape[0]))
964
965 # if not self.firsttime:
966 # data[ind] = numpy.vstack((self.ds[ind][:], data[ind]))
984 data[ind] = dataAux[:,j,:]
967 985 ind += 1
968
986
969 987 self.data = data
970 988 return
971 989
@@ -974,6 +992,8 class HDF5Writer(Operation):
974 992 Saves the block in the HDF5 file
975 993 '''
976 994 for i in range(len(self.ds)):
995
996 # First time
977 997 if self.firsttime:
978 998 # self.ds[i].resize(self.data[i].shape)
979 999 # self.ds[i][self.blockIndex,:] = self.data[i]
@@ -982,30 +1002,38 class HDF5Writer(Operation):
982 1002
983 1003 if nDims1 == 3:
984 1004 self.data[i] = self.data[i].reshape((self.data[i].shape[0],self.data[i].shape[1],1))
985
986 self.ds[i].resize(self.data[i].shape)
987 self.ds[i][:] = self.data[i]
1005
1006 self.ds[i].resize(self.data[i].shape)
1007
1008 self.ds[i][:] = self.data[i]
988 1009 else:
989 if self.nDimsForDs[i] == 1:
1010
1011 # From second time
1012 # Meteors!
1013 if self.mode[i] == 2:
1014 dataShape = self.data[i].shape
1015 dsShape = self.ds[i].shape
1016 self.ds[i].resize((self.ds[i].shape[0] + dataShape[0],self.ds[i].shape[1]))
1017 self.ds[i][dsShape[0]:,:] = self.data[i]
1018 # One dimension
1019 elif self.nDimsForDs[i] == 1:
990 1020 self.ds[i].resize((self.ds[i].shape[0], self.ds[i].shape[1] + 1))
991 1021 self.ds[i][0,-1] = self.data[i]
1022 # Two dimension
992 1023 elif self.nDimsForDs[i] == 2:
993 1024 self.ds[i].resize((self.ds[i].shape[0] + 1,self.ds[i].shape[1]))
994 1025 self.ds[i][self.blockIndex,:] = self.data[i]
1026 # Three dimensions
995 1027 elif self.nDimsForDs[i] == 3:
996
997 dataShape = self.data[i].shape
998 dsShape = self.ds[i].shape
999
1000 if dataShape[0]==dsShape[0]:
1001 self.ds[i].resize((self.ds[i].shape[0],self.ds[i].shape[1],self.ds[i].shape[2]+1))
1002 self.ds[i][:,:,-1] = self.data[i]
1003 else:
1004 self.ds[i].resize((self.ds[i].shape[0] + dataShape[0],self.ds[i].shape[1],self.ds[i].shape[2]))
1005 self.ds[i][dsShape[0]:,:,0] = self.data[i]
1028 self.ds[i].resize((self.ds[i].shape[0],self.ds[i].shape[1],self.ds[i].shape[2]+1))
1029 self.ds[i][:,:,-1] = self.data[i]
1006 1030
1031 self.firsttime = False
1007 1032 self.blockIndex += 1
1008 self.firsttime = False
1033
1034 #Close to save changes
1035 self.fp.flush()
1036 self.fp.close()
1009 1037 return
1010 1038
1011 1039 def run(self, dataOut, **kwargs):
Show file before | Show file after | Hide comments
@@ -365,7 +365,7 class ParametersProc(ProcessingUnit):
365 365 #------------------- Detect Meteors ------------------------------
366 366
367 367 def MeteorDetection(self, hei_ref = None, tauindex = 0,
368 predefinedPhaseShifts = None, centerReceiverIndex = 2, saveAll = False,
368 predefinedPhaseShifts = None, centerReceiverIndex = 2,
369 369 cohDetection = False, cohDet_timeStep = 1, cohDet_thresh = 25,
370 370 noise_timeStep = 4, noise_multiple = 4,
371 371 multDet_timeLimit = 1, multDet_rangeLimit = 3,
@@ -456,10 +456,10 class ParametersProc(ProcessingUnit):
456 456 if predefinedPhaseShifts != None:
457 457 hardwarePhaseShifts = numpy.array(predefinedPhaseShifts)*numpy.pi/180
458 458
459 elif beaconPhaseShifts:
460 #get hardware phase shifts using beacon signal
461 hardwarePhaseShifts = self.__getHardwarePhaseDiff(self.dataOut.data_pre, pairslist, newheis, 10)
462 hardwarePhaseShifts = numpy.insert(hardwarePhaseShifts,centerReceiverIndex,0)
459 # elif beaconPhaseShifts:
460 # #get hardware phase shifts using beacon signal
461 # hardwarePhaseShifts = self.__getHardwarePhaseDiff(self.dataOut.data_pre, pairslist, newheis, 10)
462 # hardwarePhaseShifts = numpy.insert(hardwarePhaseShifts,centerReceiverIndex,0)
463 463
464 464 else:
465 465 hardwarePhaseShifts = numpy.zeros(5)
@@ -529,6 +529,9 class ParametersProc(ProcessingUnit):
529 529 listMeteors4 = self.__getRadialVelocity(listMeteors3, listMeteorsVolts, radialStdThresh, pairslist, self.dataOut.timeInterval)
530 530
531 531 if len(listMeteors4) > 0:
532 #Setting New Array
533 date = self.dataOut.utctime
534 arrayParameters = self.__setNewArrays(listMeteors4, date, heiRang)
532 535
533 536 pairsList = []
534 537 pairx = (0,3)
@@ -536,10 +539,6 class ParametersProc(ProcessingUnit):
536 539 pairsList.append(pairx)
537 540 pairsList.append(pairy)
538 541
539 #Setting New Array
540 date = repr(self.dataOut.datatime)
541 arrayParameters = self.__setNewArrays(listMeteors4, date, heiRang)
542
543 542 meteorOps = MeteorOperations()
544 543 jph = numpy.array([0,0,0,0])
545 544 h = (hmin,hmax)
@@ -561,14 +560,30 class ParametersProc(ProcessingUnit):
561 560 #********************* END OF PARAMETERS CALCULATION **************************
562 561
563 562 #***************************+ PASS DATA TO NEXT STEP **********************
564 arrayFinal = arrayParameters.reshape((1,arrayParameters.shape[0],arrayParameters.shape[1]))
565 self.dataOut.data_param = arrayFinal
563 # arrayFinal = arrayParameters.reshape((1,arrayParameters.shape[0],arrayParameters.shape[1]))
564 self.dataOut.data_param = arrayParameters
566 565
567 if arrayFinal == None:
566 if arrayParameters == None:
568 567 self.dataOut.flagNoData = True
569 568
570 569 return
571 570
571 def correctMeteorPhases(self):
572
573 arrayParameters = self.dataOut.data_param
574 pairsList = []
575 pairx = (0,3)
576 pairy = (1,2)
577 pairsList.append(pairx)
578 pairsList.append(pairy)
579
580 meteorOps = MeteorOperations()
581 jph = numpy.array([0,0,0,0])
582 h = (hmin,hmax)
583 arrayParameters = meteorOps.getMeteorParams(arrayParameters, azimuth, h, pairsList, jph)
584 self.dataOut.data_param = arrayParameters
585 return
586
572 587 def __getHardwarePhaseDiff(self, voltage0, pairslist, newheis, n):
573 588
574 589 minIndex = min(newheis[0])
@@ -1055,133 +1070,134 class ParametersProc(ProcessingUnit):
1055 1070 arrayParameters = numpy.zeros((len(listMeteors), 14))
1056 1071
1057 1072 #Date inclusion
1058 date = re.findall(r'\((.*?)\)', date)
1059 date = date[0].split(',')
1060 date = map(int, date)
1061
1062 if len(date)<6:
1063 date.append(0)
1064
1065 date = [date[0]*10000 + date[1]*100 + date[2], date[3]*10000 + date[4]*100 + date[5]]
1066 arrayDate = numpy.tile(date, (len(listMeteors), 1))
1073 # date = re.findall(r'\((.*?)\)', date)
1074 # date = date[0].split(',')
1075 # date = map(int, date)
1076 #
1077 # if len(date)<6:
1078 # date.append(0)
1079 #
1080 # date = [date[0]*10000 + date[1]*100 + date[2], date[3]*10000 + date[4]*100 + date[5]]
1081 # arrayDate = numpy.tile(date, (len(listMeteors), 1))
1082 arrayDate = numpy.tile(date, (len(listMeteors)))
1067 1083
1068 1084 #Meteor array
1069 1085 # arrayMeteors[:,0] = heiRang[arrayMeteors[:,0].astype(int)]
1070 1086 # arrayMeteors = numpy.hstack((arrayDate, arrayMeteors))
1071 1087
1072 1088 #Parameters Array
1073 arrayParameters[:,:2] = arrayDate #Date
1074 arrayParameters[:,2] = heiRang[arrayMeteors[:,0].astype(int)] #Range
1075 arrayParameters[:,7:9] = arrayMeteors[:,-3:-1] #Radial velocity and its error
1076 arrayParameters[:,9:13] = arrayMeteors[:,7:11] #Phases
1089 arrayParameters[:,0] = arrayDate #Date
1090 arrayParameters[:,1] = heiRang[arrayMeteors[:,0].astype(int)] #Range
1091 arrayParameters[:,6:8] = arrayMeteors[:,-3:-1] #Radial velocity and its error
1092 arrayParameters[:,8:12] = arrayMeteors[:,7:11] #Phases
1077 1093 arrayParameters[:,-1] = arrayMeteors[:,-1] #Error
1078 1094
1079 1095
1080 1096 return arrayParameters
1081 1097
1082 def __getAOA(self, phases, pairsList, error, AOAthresh, azimuth):
1083
1084 arrayAOA = numpy.zeros((phases.shape[0],3))
1085 cosdir0, cosdir = self.__getDirectionCosines(phases, pairsList)
1086
1087 arrayAOA[:,:2] = self.__calculateAOA(cosdir, azimuth)
1088 cosDirError = numpy.sum(numpy.abs(cosdir0 - cosdir), axis = 1)
1089 arrayAOA[:,2] = cosDirError
1090
1091 azimuthAngle = arrayAOA[:,0]
1092 zenithAngle = arrayAOA[:,1]
1093
1094 #Setting Error
1095 #Number 3: AOA not fesible
1096 indInvalid = numpy.where(numpy.logical_and((numpy.logical_or(numpy.isnan(zenithAngle), numpy.isnan(azimuthAngle))),error == 0))[0]
1097 error[indInvalid] = 3
1098 #Number 4: Large difference in AOAs obtained from different antenna baselines
1099 indInvalid = numpy.where(numpy.logical_and(cosDirError > AOAthresh,error == 0))[0]
1100 error[indInvalid] = 4
1101 return arrayAOA, error
1102
1103 def __getDirectionCosines(self, arrayPhase, pairsList):
1104
1105 #Initializing some variables
1106 ang_aux = numpy.array([-8,-7,-6,-5,-4,-3,-2,-1,0,1,2,3,4,5,6,7,8])*2*numpy.pi
1107 ang_aux = ang_aux.reshape(1,ang_aux.size)
1108
1109 cosdir = numpy.zeros((arrayPhase.shape[0],2))
1110 cosdir0 = numpy.zeros((arrayPhase.shape[0],2))
1111
1112
1113 for i in range(2):
1114 #First Estimation
1115 phi0_aux = arrayPhase[:,pairsList[i][0]] + arrayPhase[:,pairsList[i][1]]
1116 #Dealias
1117 indcsi = numpy.where(phi0_aux > numpy.pi)
1118 phi0_aux[indcsi] -= 2*numpy.pi
1119 indcsi = numpy.where(phi0_aux < -numpy.pi)
1120 phi0_aux[indcsi] += 2*numpy.pi
1121 #Direction Cosine 0
1122 cosdir0[:,i] = -(phi0_aux)/(2*numpy.pi*0.5)
1123
1124 #Most-Accurate Second Estimation
1125 phi1_aux = arrayPhase[:,pairsList[i][0]] - arrayPhase[:,pairsList[i][1]]
1126 phi1_aux = phi1_aux.reshape(phi1_aux.size,1)
1127 #Direction Cosine 1
1128 cosdir1 = -(phi1_aux + ang_aux)/(2*numpy.pi*4.5)
1129
1130 #Searching the correct Direction Cosine
1131 cosdir0_aux = cosdir0[:,i]
1132 cosdir0_aux = cosdir0_aux.reshape(cosdir0_aux.size,1)
1133 #Minimum Distance
1134 cosDiff = (cosdir1 - cosdir0_aux)**2
1135 indcos = cosDiff.argmin(axis = 1)
1136 #Saving Value obtained
1137 cosdir[:,i] = cosdir1[numpy.arange(len(indcos)),indcos]
1138
1139 return cosdir0, cosdir
1140
1141 def __calculateAOA(self, cosdir, azimuth):
1142 cosdirX = cosdir[:,0]
1143 cosdirY = cosdir[:,1]
1144
1145 zenithAngle = numpy.arccos(numpy.sqrt(1 - cosdirX**2 - cosdirY**2))*180/numpy.pi
1146 azimuthAngle = numpy.arctan2(cosdirX,cosdirY)*180/numpy.pi + azimuth #0 deg north, 90 deg east
1147 angles = numpy.vstack((azimuthAngle, zenithAngle)).transpose()
1148
1149 return angles
1150
1151 def __getHeights(self, Ranges, zenith, error, minHeight, maxHeight):
1152
1153 Ramb = 375 #Ramb = c/(2*PRF)
1154 Re = 6371 #Earth Radius
1155 heights = numpy.zeros(Ranges.shape)
1156
1157 R_aux = numpy.array([0,1,2])*Ramb
1158 R_aux = R_aux.reshape(1,R_aux.size)
1159
1160 Ranges = Ranges.reshape(Ranges.size,1)
1161
1162 Ri = Ranges + R_aux
1163 hi = numpy.sqrt(Re**2 + Ri**2 + (2*Re*numpy.cos(zenith*numpy.pi/180)*Ri.transpose()).transpose()) - Re
1164
1165 #Check if there is a height between 70 and 110 km
1166 h_bool = numpy.sum(numpy.logical_and(hi > minHeight, hi < maxHeight), axis = 1)
1167 ind_h = numpy.where(h_bool == 1)[0]
1168
1169 hCorr = hi[ind_h, :]
1170 ind_hCorr = numpy.where(numpy.logical_and(hi > minHeight, hi < maxHeight))
1171
1172 hCorr = hi[ind_hCorr]
1173 heights[ind_h] = hCorr
1174
1175 #Setting Error
1176 #Number 13: Height unresolvable echo: not valid height within 70 to 110 km
1177 #Number 14: Height ambiguous echo: more than one possible height within 70 to 110 km
1178
1179 indInvalid2 = numpy.where(numpy.logical_and(h_bool > 1, error == 0))[0]
1180 error[indInvalid2] = 14
1181 indInvalid1 = numpy.where(numpy.logical_and(h_bool == 0, error == 0))[0]
1182 error[indInvalid1] = 13
1183
1184 return heights, error
1098 # def __getAOA(self, phases, pairsList, error, AOAthresh, azimuth):
1099 #
1100 # arrayAOA = numpy.zeros((phases.shape[0],3))
1101 # cosdir0, cosdir = self.__getDirectionCosines(phases, pairsList)
1102 #
1103 # arrayAOA[:,:2] = self.__calculateAOA(cosdir, azimuth)
1104 # cosDirError = numpy.sum(numpy.abs(cosdir0 - cosdir), axis = 1)
1105 # arrayAOA[:,2] = cosDirError
1106 #
1107 # azimuthAngle = arrayAOA[:,0]
1108 # zenithAngle = arrayAOA[:,1]
1109 #
1110 # #Setting Error
1111 # #Number 3: AOA not fesible
1112 # indInvalid = numpy.where(numpy.logical_and((numpy.logical_or(numpy.isnan(zenithAngle), numpy.isnan(azimuthAngle))),error == 0))[0]
1113 # error[indInvalid] = 3
1114 # #Number 4: Large difference in AOAs obtained from different antenna baselines
1115 # indInvalid = numpy.where(numpy.logical_and(cosDirError > AOAthresh,error == 0))[0]
1116 # error[indInvalid] = 4
1117 # return arrayAOA, error
1118 #
1119 # def __getDirectionCosines(self, arrayPhase, pairsList):
1120 #
1121 # #Initializing some variables
1122 # ang_aux = numpy.array([-8,-7,-6,-5,-4,-3,-2,-1,0,1,2,3,4,5,6,7,8])*2*numpy.pi
1123 # ang_aux = ang_aux.reshape(1,ang_aux.size)
1124 #
1125 # cosdir = numpy.zeros((arrayPhase.shape[0],2))
1126 # cosdir0 = numpy.zeros((arrayPhase.shape[0],2))
1127 #
1128 #
1129 # for i in range(2):
1130 # #First Estimation
1131 # phi0_aux = arrayPhase[:,pairsList[i][0]] + arrayPhase[:,pairsList[i][1]]
1132 # #Dealias
1133 # indcsi = numpy.where(phi0_aux > numpy.pi)
1134 # phi0_aux[indcsi] -= 2*numpy.pi
1135 # indcsi = numpy.where(phi0_aux < -numpy.pi)
1136 # phi0_aux[indcsi] += 2*numpy.pi
1137 # #Direction Cosine 0
1138 # cosdir0[:,i] = -(phi0_aux)/(2*numpy.pi*0.5)
1139 #
1140 # #Most-Accurate Second Estimation
1141 # phi1_aux = arrayPhase[:,pairsList[i][0]] - arrayPhase[:,pairsList[i][1]]
1142 # phi1_aux = phi1_aux.reshape(phi1_aux.size,1)
1143 # #Direction Cosine 1
1144 # cosdir1 = -(phi1_aux + ang_aux)/(2*numpy.pi*4.5)
1145 #
1146 # #Searching the correct Direction Cosine
1147 # cosdir0_aux = cosdir0[:,i]
1148 # cosdir0_aux = cosdir0_aux.reshape(cosdir0_aux.size,1)
1149 # #Minimum Distance
1150 # cosDiff = (cosdir1 - cosdir0_aux)**2
1151 # indcos = cosDiff.argmin(axis = 1)
1152 # #Saving Value obtained
1153 # cosdir[:,i] = cosdir1[numpy.arange(len(indcos)),indcos]
1154 #
1155 # return cosdir0, cosdir
1156 #
1157 # def __calculateAOA(self, cosdir, azimuth):
1158 # cosdirX = cosdir[:,0]
1159 # cosdirY = cosdir[:,1]
1160 #
1161 # zenithAngle = numpy.arccos(numpy.sqrt(1 - cosdirX**2 - cosdirY**2))*180/numpy.pi
1162 # azimuthAngle = numpy.arctan2(cosdirX,cosdirY)*180/numpy.pi + azimuth #0 deg north, 90 deg east
1163 # angles = numpy.vstack((azimuthAngle, zenithAngle)).transpose()
1164 #
1165 # return angles
1166 #
1167 # def __getHeights(self, Ranges, zenith, error, minHeight, maxHeight):
1168 #
1169 # Ramb = 375 #Ramb = c/(2*PRF)
1170 # Re = 6371 #Earth Radius
1171 # heights = numpy.zeros(Ranges.shape)
1172 #
1173 # R_aux = numpy.array([0,1,2])*Ramb
1174 # R_aux = R_aux.reshape(1,R_aux.size)
1175 #
1176 # Ranges = Ranges.reshape(Ranges.size,1)
1177 #
1178 # Ri = Ranges + R_aux
1179 # hi = numpy.sqrt(Re**2 + Ri**2 + (2*Re*numpy.cos(zenith*numpy.pi/180)*Ri.transpose()).transpose()) - Re
1180 #
1181 # #Check if there is a height between 70 and 110 km
1182 # h_bool = numpy.sum(numpy.logical_and(hi > minHeight, hi < maxHeight), axis = 1)
1183 # ind_h = numpy.where(h_bool == 1)[0]
1184 #
1185 # hCorr = hi[ind_h, :]
1186 # ind_hCorr = numpy.where(numpy.logical_and(hi > minHeight, hi < maxHeight))
1187 #
1188 # hCorr = hi[ind_hCorr]
1189 # heights[ind_h] = hCorr
1190 #
1191 # #Setting Error
1192 # #Number 13: Height unresolvable echo: not valid height within 70 to 110 km
1193 # #Number 14: Height ambiguous echo: more than one possible height within 70 to 110 km
1194 #
1195 # indInvalid2 = numpy.where(numpy.logical_and(h_bool > 1, error == 0))[0]
1196 # error[indInvalid2] = 14
1197 # indInvalid1 = numpy.where(numpy.logical_and(h_bool == 0, error == 0))[0]
1198 # error[indInvalid1] = 13
1199 #
1200 # return heights, error
1185 1201
1186 1202 def SpectralFitting(self, getSNR = True, path=None, file=None, groupList=None):
1187 1203
@@ -1321,7 +1337,6 class ParametersProc(ProcessingUnit):
1321 1337 chisq=numpy.dot((dp-fmp).T,(dp-fmp))
1322 1338 return chisq
1323 1339
1324
1325 1340
1326 1341 class WindProfiler(Operation):
1327 1342
@@ -1600,11 +1615,11 class WindProfiler(Operation):
1600 1615 winds = numpy.zeros((2,nInt))*numpy.nan
1601 1616
1602 1617 #Filter errors
1603 error = numpy.where(arrayMeteor[0,:,-1] == 0)[0]
1604 finalMeteor = arrayMeteor[0,error,:]
1618 error = numpy.where(arrayMeteor[:,-1] == 0)[0]
1619 finalMeteor = arrayMeteor[error,:]
1605 1620
1606 1621 #Meteor Histogram
1607 finalHeights = finalMeteor[:,3]
1622 finalHeights = finalMeteor[:,2]
1608 1623 hist = numpy.histogram(finalHeights, bins = nInt, range = (heightMin,heightMax))
1609 1624 nMeteorsPerI = hist[0]
1610 1625 heightPerI = hist[1]
@@ -1625,9 +1640,9 class WindProfiler(Operation):
1625 1640 meteorAux = finalMeteor2[ind1:ind2,:]
1626 1641
1627 1642 if meteorAux.shape[0] >= meteorThresh:
1628 vel = meteorAux[:, 7]
1629 zen = meteorAux[:, 5]*numpy.pi/180
1630 azim = meteorAux[:, 4]*numpy.pi/180
1643 vel = meteorAux[:, 6]
1644 zen = meteorAux[:, 4]*numpy.pi/180
1645 azim = meteorAux[:, 3]*numpy.pi/180
1631 1646
1632 1647 n = numpy.cos(zen)
1633 1648 # m = (1 - n**2)/(1 - numpy.tan(azim)**2)
@@ -1745,7 +1760,7 class WindProfiler(Operation):
1745 1760 self.__firstdata = copy.copy(dataOut)
1746 1761
1747 1762 else:
1748 self.__buffer = numpy.hstack((self.__buffer, dataOut.data_param))
1763 self.__buffer = numpy.vstack((self.__buffer, dataOut.data_param))
1749 1764
1750 1765 self.__checkTime(dataOut.utctime, dataOut.paramInterval, dataOut.outputInterval) #Check if the buffer is ready
1751 1766
@@ -1996,13 +2011,13 class PhaseCalibration(Operation):
1996 2011 h = (hmin, hmax)
1997 2012 pairsList = ((0,3),(1,2))
1998 2013
1999 meteorsArray = self.__buffer[0,:,:]
2014 meteorsArray = self.__buffer
2000 2015 error = meteorsArray[:,-1]
2001 2016 boolError = (error==0)|(error==3)|(error==4)|(error==13)|(error==14)
2002 2017 ind1 = numpy.where(boolError)[0]
2003 2018 meteorsArray = meteorsArray[ind1,:]
2004 2019 meteorsArray[:,-1] = 0
2005 phases = meteorsArray[:,9:13]
2020 phases = meteorsArray[:,8:12]
2006 2021
2007 2022 #Calculate Gammas
2008 2023 gammas = self.__getGammas(pairs, k, distances, phases)
@@ -2033,14 +2048,14 class MeteorOperations():
2033 2048 #JONES ET AL. 1998
2034 2049 AOAthresh = numpy.pi/8
2035 2050 error = arrayParameters[:,-1]
2036 phases = -arrayParameters[:,9:13] + jph
2037 arrayParameters[:,4:7], arrayParameters[:,-1] = self.__getAOA(phases, pairsList, error, AOAthresh, azimuth)
2051 phases = -arrayParameters[:,8:12] + jph
2052 arrayParameters[:,3:6], arrayParameters[:,-1] = self.__getAOA(phases, pairsList, error, AOAthresh, azimuth)
2038 2053
2039 2054 #Calculate Heights (Error N 13 and 14)
2040 2055 error = arrayParameters[:,-1]
2041 Ranges = arrayParameters[:,2]
2042 zenith = arrayParameters[:,5]
2043 arrayParameters[:,3], arrayParameters[:,-1] = self.__getHeights(Ranges, zenith, error, hmin, hmax)
2056 Ranges = arrayParameters[:,1]
2057 zenith = arrayParameters[:,4]
2058 arrayParameters[:,2], arrayParameters[:,-1] = self.__getHeights(Ranges, zenith, error, hmin, hmax)
2044 2059
2045 2060 #----------------------- Get Final data ------------------------------------
2046 2061 # error = arrayParameters[:,-1]
General Comments 0
You need to be logged in to leave comments. Login now