##// END OF EJS Templates
schain-jc
RSS
Review BLTRSpectraReader Unit
Juan C. Espinoza -
r1199:13d983fa985c
parent child
Show More
Side by Side Unified
Context file:
r1199:13d983fa985c Loading diff...:
Show file before | Show file after | Hide comments
@@ -0,0 +1,466
1 import os
2 import sys
3 import glob
4 import fnmatch
5 import datetime
6 import time
7 import re
8 import h5py
9 import numpy
10
11 import pylab as plb
12 from scipy.optimize import curve_fit
13 from scipy import asarray as ar, exp
14
15 SPEED_OF_LIGHT = 299792458
16 SPEED_OF_LIGHT = 3e8
17
18 try:
19 from gevent import sleep
20 except:
21 from time import sleep
22
23 from .utils import folder_in_range
24
25 from schainpy.model.data.jrodata import Spectra
26 from schainpy.model.proc.jroproc_base import ProcessingUnit, Operation, MPDecorator
27 from schainpy.utils import log
28 from schainpy.model.io.jroIO_base import JRODataReader
29
30 def pol2cart(rho, phi):
31 x = rho * numpy.cos(phi)
32 y = rho * numpy.sin(phi)
33 return(x, y)
34
35 FILE_STRUCTURE = numpy.dtype([ # HEADER 48bytes
36 ('FileMgcNumber', '<u4'), # 0x23020100
37 ('nFDTdataRecors', '<u4'),
38 ('OffsetStartHeader', '<u4'),
39 ('RadarUnitId', '<u4'),
40 ('SiteName', 'S32'), # Null terminated
41 ])
42
43
44 class FileHeaderBLTR():
45
46 def __init__(self, fo):
47
48 self.fo = fo
49 self.size = 48
50 self.read()
51
52 def read(self):
53
54 header = numpy.fromfile(self.fo, FILE_STRUCTURE, 1)
55 self.FileMgcNumber = hex(header['FileMgcNumber'][0])
56 self.nFDTdataRecors = int(header['nFDTdataRecors'][0])
57 self.RadarUnitId = int(header['RadarUnitId'][0])
58 self.OffsetStartHeader = int(header['OffsetStartHeader'][0])
59 self.SiteName = header['SiteName'][0]
60
61 def write(self, fp):
62
63 headerTuple = (self.FileMgcNumber,
64 self.nFDTdataRecors,
65 self.RadarUnitId,
66 self.SiteName,
67 self.size)
68
69 header = numpy.array(headerTuple, FILE_STRUCTURE)
70 header.tofile(fp)
71 ''' ndarray.tofile(fid, sep, format) Write array to a file as text or binary (default).
72
73 fid : file or str
74 An open file object, or a string containing a filename.
75
76 sep : str
77 Separator between array items for text output. If "" (empty), a binary file is written,
78 equivalent to file.write(a.tobytes()).
79
80 format : str
81 Format string for text file output. Each entry in the array is formatted to text by
82 first converting it to the closest Python type, and then using "format" % item.
83
84 '''
85
86 return 1
87
88
89 RECORD_STRUCTURE = numpy.dtype([ # RECORD HEADER 180+20N bytes
90 ('RecMgcNumber', '<u4'), # 0x23030001
91 ('RecCounter', '<u4'), # Record counter(0,1, ...)
92 # Offset to start of next record form start of this record
93 ('Off2StartNxtRec', '<u4'),
94 # Offset to start of data from start of this record
95 ('Off2StartData', '<u4'),
96 # Epoch time stamp of start of acquisition (seconds)
97 ('nUtime', '<i4'),
98 # Millisecond component of time stamp (0,...,999)
99 ('nMilisec', '<u4'),
100 # Experiment tag name (null terminated)
101 ('ExpTagName', 'S32'),
102 # Experiment comment (null terminated)
103 ('ExpComment', 'S32'),
104 # Site latitude (from GPS) in degrees (positive implies North)
105 ('SiteLatDegrees', '<f4'),
106 # Site longitude (from GPS) in degrees (positive implies East)
107 ('SiteLongDegrees', '<f4'),
108 # RTC GPS engine status (0=SEEK, 1=LOCK, 2=NOT FITTED, 3=UNAVAILABLE)
109 ('RTCgpsStatus', '<u4'),
110 ('TransmitFrec', '<u4'), # Transmit frequency (Hz)
111 ('ReceiveFrec', '<u4'), # Receive frequency
112 # First local oscillator frequency (Hz)
113 ('FirstOsciFrec', '<u4'),
114 # (0="O", 1="E", 2="linear 1", 3="linear2")
115 ('Polarisation', '<u4'),
116 # Receiver filter settings (0,1,2,3)
117 ('ReceiverFiltSett', '<u4'),
118 # Number of modes in use (1 or 2)
119 ('nModesInUse', '<u4'),
120 # Dual Mode index number for these data (0 or 1)
121 ('DualModeIndex', '<u4'),
122 # Dual Mode range correction for these data (m)
123 ('DualModeRange', '<u4'),
124 # Number of digital channels acquired (2*N)
125 ('nDigChannels', '<u4'),
126 # Sampling resolution (meters)
127 ('SampResolution', '<u4'),
128 # Number of range gates sampled
129 ('nHeights', '<u4'),
130 # Start range of sampling (meters)
131 ('StartRangeSamp', '<u4'),
132 ('PRFhz', '<u4'), # PRF (Hz)
133 ('nCohInt', '<u4'), # Integrations
134 # Number of data points transformed
135 ('nProfiles', '<u4'),
136 # Number of receive beams stored in file (1 or N)
137 ('nChannels', '<u4'),
138 ('nIncohInt', '<u4'), # Number of spectral averages
139 # FFT windowing index (0 = no window)
140 ('FFTwindowingInd', '<u4'),
141 # Beam steer angle (azimuth) in degrees (clockwise from true North)
142 ('BeamAngleAzim', '<f4'),
143 # Beam steer angle (zenith) in degrees (0=> vertical)
144 ('BeamAngleZen', '<f4'),
145 # Antenna coordinates (Range(meters), Bearing(degrees)) - N pairs
146 ('AntennaCoord0', '<f4'),
147 # Antenna coordinates (Range(meters), Bearing(degrees)) - N pairs
148 ('AntennaAngl0', '<f4'),
149 # Antenna coordinates (Range(meters), Bearing(degrees)) - N pairs
150 ('AntennaCoord1', '<f4'),
151 # Antenna coordinates (Range(meters), Bearing(degrees)) - N pairs
152 ('AntennaAngl1', '<f4'),
153 # Antenna coordinates (Range(meters), Bearing(degrees)) - N pairs
154 ('AntennaCoord2', '<f4'),
155 # Antenna coordinates (Range(meters), Bearing(degrees)) - N pairs
156 ('AntennaAngl2', '<f4'),
157 # Receiver phase calibration (degrees) - N values
158 ('RecPhaseCalibr0', '<f4'),
159 # Receiver phase calibration (degrees) - N values
160 ('RecPhaseCalibr1', '<f4'),
161 # Receiver phase calibration (degrees) - N values
162 ('RecPhaseCalibr2', '<f4'),
163 # Receiver amplitude calibration (ratio relative to receiver one) - N values
164 ('RecAmpCalibr0', '<f4'),
165 # Receiver amplitude calibration (ratio relative to receiver one) - N values
166 ('RecAmpCalibr1', '<f4'),
167 # Receiver amplitude calibration (ratio relative to receiver one) - N values
168 ('RecAmpCalibr2', '<f4'),
169 # Receiver gains in dB - N values
170 ('ReceiverGaindB0', '<i4'),
171 # Receiver gains in dB - N values
172 ('ReceiverGaindB1', '<i4'),
173 # Receiver gains in dB - N values
174 ('ReceiverGaindB2', '<i4'),
175 ])
176
177
178 class RecordHeaderBLTR():
179
180 def __init__(self, fo):
181
182 self.fo = fo
183 self.OffsetStartHeader = 48
184 self.Off2StartNxtRec = 811248
185
186 def read(self, block):
187 OffRHeader = self.OffsetStartHeader + block * self.Off2StartNxtRec
188 self.fo.seek(OffRHeader, os.SEEK_SET)
189 header = numpy.fromfile(self.fo, RECORD_STRUCTURE, 1)
190 self.RecMgcNumber = hex(header['RecMgcNumber'][0]) # 0x23030001
191 self.RecCounter = int(header['RecCounter'][0])
192 self.Off2StartNxtRec = int(header['Off2StartNxtRec'][0])
193 self.Off2StartData = int(header['Off2StartData'][0])
194 self.nUtime = header['nUtime'][0]
195 self.nMilisec = header['nMilisec'][0]
196 self.ExpTagName = '' # str(header['ExpTagName'][0])
197 self.ExpComment = '' # str(header['ExpComment'][0])
198 self.SiteLatDegrees = header['SiteLatDegrees'][0]
199 self.SiteLongDegrees = header['SiteLongDegrees'][0]
200 self.RTCgpsStatus = header['RTCgpsStatus'][0]
201 self.TransmitFrec = header['TransmitFrec'][0]
202 self.ReceiveFrec = header['ReceiveFrec'][0]
203 self.FirstOsciFrec = header['FirstOsciFrec'][0]
204 self.Polarisation = header['Polarisation'][0]
205 self.ReceiverFiltSett = header['ReceiverFiltSett'][0]
206 self.nModesInUse = header['nModesInUse'][0]
207 self.DualModeIndex = header['DualModeIndex'][0]
208 self.DualModeRange = header['DualModeRange'][0]
209 self.nDigChannels = header['nDigChannels'][0]
210 self.SampResolution = header['SampResolution'][0]
211 self.nHeights = header['nHeights'][0]
212 self.StartRangeSamp = header['StartRangeSamp'][0]
213 self.PRFhz = header['PRFhz'][0]
214 self.nCohInt = header['nCohInt'][0]
215 self.nProfiles = header['nProfiles'][0]
216 self.nChannels = header['nChannels'][0]
217 self.nIncohInt = header['nIncohInt'][0]
218 self.FFTwindowingInd = header['FFTwindowingInd'][0]
219 self.BeamAngleAzim = header['BeamAngleAzim'][0]
220 self.BeamAngleZen = header['BeamAngleZen'][0]
221 self.AntennaCoord0 = header['AntennaCoord0'][0]
222 self.AntennaAngl0 = header['AntennaAngl0'][0]
223 self.AntennaCoord1 = header['AntennaCoord1'][0]
224 self.AntennaAngl1 = header['AntennaAngl1'][0]
225 self.AntennaCoord2 = header['AntennaCoord2'][0]
226 self.AntennaAngl2 = header['AntennaAngl2'][0]
227 self.RecPhaseCalibr0 = header['RecPhaseCalibr0'][0]
228 self.RecPhaseCalibr1 = header['RecPhaseCalibr1'][0]
229 self.RecPhaseCalibr2 = header['RecPhaseCalibr2'][0]
230 self.RecAmpCalibr0 = header['RecAmpCalibr0'][0]
231 self.RecAmpCalibr1 = header['RecAmpCalibr1'][0]
232 self.RecAmpCalibr2 = header['RecAmpCalibr2'][0]
233 self.ReceiverGaindB0 = header['ReceiverGaindB0'][0]
234 self.ReceiverGaindB1 = header['ReceiverGaindB1'][0]
235 self.ReceiverGaindB2 = header['ReceiverGaindB2'][0]
236 self.ipp = 0.5 * (SPEED_OF_LIGHT / self.PRFhz)
237 self.RHsize = 180 + 20 * self.nChannels
238 self.Datasize = self.nProfiles * self.nChannels * self.nHeights * 2 * 4
239 endFp = self.OffsetStartHeader + self.RecCounter * self.Off2StartNxtRec
240
241
242 if OffRHeader > endFp:
243 sys.stderr.write(
244 "Warning %s: Size value read from System Header is lower than it has to be\n" % fp)
245 return 0
246
247 if OffRHeader < endFp:
248 sys.stderr.write(
249 "Warning %s: Size value read from System Header size is greater than it has to be\n" % fp)
250 return 0
251
252 return 1
253
254 @MPDecorator
255 class BLTRSpectraReader (ProcessingUnit):
256
257 def __init__(self):
258
259 ProcessingUnit.__init__(self)
260
261 self.ext = ".fdt"
262 self.optchar = "P"
263 self.fpFile = None
264 self.fp = None
265 self.BlockCounter = 0
266 self.fileSizeByHeader = None
267 self.filenameList = []
268 self.fileSelector = 0
269 self.Off2StartNxtRec = 0
270 self.RecCounter = 0
271 self.flagNoMoreFiles = 0
272 self.data_spc = None
273 self.data_cspc = None
274 self.path = None
275 self.OffsetStartHeader = 0
276 self.Off2StartData = 0
277 self.ipp = 0
278 self.nFDTdataRecors = 0
279 self.blocksize = 0
280 self.dataOut = Spectra()
281 self.dataOut.flagNoData = False
282
283 def search_files(self):
284 '''
285 Function that indicates the number of .fdt files that exist in the folder to be read.
286 It also creates an organized list with the names of the files to read.
287 '''
288
289 files = glob.glob(os.path.join(self.path, '*{}'.format(self.ext)))
290 files = sorted(files)
291 for f in files:
292 filename = f.split('/')[-1]
293 if folder_in_range(filename.split('.')[1], self.startDate, self.endDate, '%Y%m%d'):
294 self.filenameList.append(f)
295
296 def run(self, **kwargs):
297 '''
298 This method will be the one that will initiate the data entry, will be called constantly.
299 You should first verify that your Setup () is set up and then continue to acquire
300 the data to be processed with getData ().
301 '''
302 if not self.isConfig:
303 self.setup(**kwargs)
304 self.isConfig = True
305
306 self.getData()
307
308 def setup(self,
309 path=None,
310 startDate=None,
311 endDate=None,
312 startTime=None,
313 endTime=None,
314 walk=True,
315 code=None,
316 online=False,
317 mode=None,
318 **kwargs):
319
320 self.isConfig = True
321
322 self.path = path
323 self.startDate = startDate
324 self.endDate = endDate
325 self.startTime = startTime
326 self.endTime = endTime
327 self.walk = walk
328 self.mode = int(mode)
329 self.search_files()
330 if self.filenameList:
331 self.readFile()
332
333 def getData(self):
334 '''
335 Before starting this function, you should check that there is still an unread file,
336 If there are still blocks to read or if the data block is empty.
337
338 You should call the file "read".
339
340 '''
341
342 if self.flagNoMoreFiles:
343 self.dataOut.flagNoData = True
344 self.dataOut.error = 'No more files'
345
346 self.readBlock()
347
348 def readFile(self):
349 '''
350 You must indicate if you are reading in Online or Offline mode and load the
351 The parameters for this file reading mode.
352
353 Then you must do 2 actions:
354
355 1. Get the BLTR FileHeader.
356 2. Start reading the first block.
357 '''
358
359 if self.fileSelector < len(self.filenameList):
360 log.success('Opening file: {}'.format(self.filenameList[self.fileSelector]), self.name)
361 self.fp = open(self.filenameList[self.fileSelector])
362 self.fheader = FileHeaderBLTR(self.fp)
363 self.rheader = RecordHeaderBLTR(self.fp)
364 self.nFDTdataRecors = self.fheader.nFDTdataRecors
365 self.fileSelector += 1
366 self.BlockCounter = 0
367 return 1
368 else:
369 self.flagNoMoreFiles = True
370 self.dataOut.flagNoData = True
371 return 0
372
373 def readBlock(self):
374 '''
375 It should be checked if the block has data, if it is not passed to the next file.
376
377 Then the following is done:
378
379 1. Read the RecordHeader
380 2. Fill the buffer with the current block number.
381
382 '''
383
384 if self.BlockCounter == self.nFDTdataRecors:
385 if not self.readFile():
386 return
387
388 if self.mode == 1:
389 self.rheader.read(self.BlockCounter+1)
390 elif self.mode == 0:
391 self.rheader.read(self.BlockCounter)
392
393 self.RecCounter = self.rheader.RecCounter
394 self.OffsetStartHeader = self.rheader.OffsetStartHeader
395 self.Off2StartNxtRec = self.rheader.Off2StartNxtRec
396 self.Off2StartData = self.rheader.Off2StartData
397 self.nProfiles = self.rheader.nProfiles
398 self.nChannels = self.rheader.nChannels
399 self.nHeights = self.rheader.nHeights
400 self.frequency = self.rheader.TransmitFrec
401 self.DualModeIndex = self.rheader.DualModeIndex
402 self.pairsList = [(0, 1), (0, 2), (1, 2)]
403 self.dataOut.pairsList = self.pairsList
404 self.nRdPairs = len(self.dataOut.pairsList)
405 self.dataOut.nRdPairs = self.nRdPairs
406 self.dataOut.heightList = (self.rheader.StartRangeSamp + numpy.arange(self.nHeights) * self.rheader.SampResolution) / 1000.
407 self.dataOut.channelList = range(self.nChannels)
408 self.dataOut.nProfiles=self.rheader.nProfiles
409 self.dataOut.nIncohInt=self.rheader.nIncohInt
410 self.dataOut.nCohInt=self.rheader.nCohInt
411 self.dataOut.ippSeconds= 1/float(self.rheader.PRFhz)
412 self.dataOut.PRF=self.rheader.PRFhz
413 self.dataOut.nFFTPoints=self.rheader.nProfiles
414 self.dataOut.utctime=self.rheader.nUtime
415 self.dataOut.timeZone = 0
416 self.dataOut.useLocalTime = False
417 self.dataOut.nmodes = 2
418 log.log('Reading block {} - {}'.format(self.BlockCounter, self.dataOut.datatime), self.name)
419 OffDATA = self.OffsetStartHeader + self.RecCounter * \
420 self.Off2StartNxtRec + self.Off2StartData
421 self.fp.seek(OffDATA, os.SEEK_SET)
422
423 self.data_fft = numpy.fromfile(self.fp, [('complex','<c8')], self.nProfiles*self.nChannels*self.nHeights )
424 self.data_fft = self.data_fft.astype(numpy.dtype('complex'))
425 self.data_block = numpy.reshape(self.data_fft,(self.nHeights, self.nChannels, self.nProfiles))
426 self.data_block = numpy.transpose(self.data_block, (1,2,0))
427 copy = self.data_block.copy()
428 spc = copy * numpy.conjugate(copy)
429 self.data_spc = numpy.absolute(spc) # valor absoluto o magnitud
430 self.dataOut.data_spc = self.data_spc
431
432 cspc = self.data_block.copy()
433 self.data_cspc = self.data_block.copy()
434
435 for i in range(self.nRdPairs):
436
437 chan_index0 = self.dataOut.pairsList[i][0]
438 chan_index1 = self.dataOut.pairsList[i][1]
439
440 self.data_cspc[i, :, :] = cspc[chan_index0, :, :] * numpy.conjugate(cspc[chan_index1, :, :])
441
442 '''Getting Eij and Nij'''
443 (AntennaX0, AntennaY0) = pol2cart(
444 self.rheader.AntennaCoord0, self.rheader.AntennaAngl0 * numpy.pi / 180)
445 (AntennaX1, AntennaY1) = pol2cart(
446 self.rheader.AntennaCoord1, self.rheader.AntennaAngl1 * numpy.pi / 180)
447 (AntennaX2, AntennaY2) = pol2cart(
448 self.rheader.AntennaCoord2, self.rheader.AntennaAngl2 * numpy.pi / 180)
449
450 E01 = AntennaX0 - AntennaX1
451 N01 = AntennaY0 - AntennaY1
452
453 E02 = AntennaX0 - AntennaX2
454 N02 = AntennaY0 - AntennaY2
455
456 E12 = AntennaX1 - AntennaX2
457 N12 = AntennaY1 - AntennaY2
458
459 self.ChanDist = numpy.array(
460 [[E01, N01], [E02, N02], [E12, N12]])
461
462 self.dataOut.ChanDist = self.ChanDist
463
464 self.BlockCounter += 2
465 self.dataOut.data_spc = self.data_spc
466 self.dataOut.data_cspc =self.data_cspc
Show file before | Show file after | Hide comments
@@ -1,23 +1,23
1 '''
1 '''
2
2
3 $Author: murco $
3 $Author: murco $
4 $Id: JRODataIO.py 169 2012-11-19 21:57:03Z murco $
4 $Id: JRODataIO.py 169 2012-11-19 21:57:03Z murco $
5 '''
5 '''
6
6
7 from .jroIO_voltage import *
7 from .jroIO_voltage import *
8 from .jroIO_spectra import *
8 from .jroIO_spectra import *
9 from .jroIO_heispectra import *
9 from .jroIO_heispectra import *
10 from .jroIO_usrp import *
10 from .jroIO_usrp import *
11 from .jroIO_digitalRF import *
11 from .jroIO_digitalRF import *
12 from .jroIO_kamisr import *
12 from .jroIO_kamisr import *
13 from .jroIO_param import *
13 from .jroIO_param import *
14 from .jroIO_hf import *
14 from .jroIO_hf import *
15
15
16 from .jroIO_madrigal import *
16 from .jroIO_madrigal import *
17
17
18 from .bltrIO_param import *
18 from .bltrIO_param import *
19 from .jroIO_bltr import *
19 from .bltrIO_spectra import *
20 from .jroIO_mira35c import *
20 from .jroIO_mira35c import *
21 from .julIO_param import *
21 from .julIO_param import *
22
22
23 from .pxIO_param import * No newline at end of file
23 from .pxIO_param import *
Show file before | Show file after | Hide comments
1 NO CONTENT: file was removed
NO CONTENT: file was removed
This diff has been collapsed as it changes many lines, (697 lines changed) Show them Hide them
General Comments 0
You need to be logged in to leave comments. Login now