##// END OF EJS Templates
schain-jc
RSS
Julia reader
Juan C. Espinoza -
r1101:f8230ce3f7eb
parent child
Show More
Side by Side Unified
Context file:
r1101:f8230ce3f7eb Loading diff...:
Show file before | Show file after | Hide comments
@@ -0,0 +1,385
1 '''
2 Created on Set 10, 2017
3
4 @author: Juan C. Espinoza
5 '''
6
7
8 import os
9 import sys
10 import time
11 import glob
12 import datetime
13
14 import numpy
15
16 from schainpy.model.proc.jroproc_base import ProcessingUnit
17 from schainpy.model.data.jrodata import Parameters
18 from schainpy.model.io.jroIO_base import JRODataReader, isNumber
19 from schainpy.utils import log
20
21 FILE_HEADER_STRUCTURE = numpy.dtype([
22 ('year', 'f'),
23 ('doy', 'f'),
24 ('nint', 'f'),
25 ('navg', 'f'),
26 ('fh', 'f'),
27 ('dh', 'f'),
28 ('nheights', 'f'),
29 ('ipp', 'f')
30 ])
31
32 REC_HEADER_STRUCTURE = numpy.dtype([
33 ('magic', 'f'),
34 ('hours', 'f'),
35 ('timedelta', 'f'),
36 ('h0', 'f'),
37 ('nheights', 'f'),
38 ('snr1', 'f'),
39 ('snr2', 'f'),
40 ('snr', 'f'),
41 ])
42
43 DATA_STRUCTURE = numpy.dtype([
44 ('range', '<u4'),
45 ('status', '<u4'),
46 ('zonal', '<f4'),
47 ('meridional', '<f4'),
48 ('vertical', '<f4'),
49 ('zonal_a', '<f4'),
50 ('meridional_a', '<f4'),
51 ('corrected_fading', '<f4'), # seconds
52 ('uncorrected_fading', '<f4'), # seconds
53 ('time_diff', '<f4'),
54 ('major_axis', '<f4'),
55 ('axial_ratio', '<f4'),
56 ('orientation', '<f4'),
57 ('sea_power', '<u4'),
58 ('sea_algorithm', '<u4')
59 ])
60
61
62 class JULIAParamReader(JRODataReader, ProcessingUnit):
63 '''
64 Julia data (eej, spf, 150km) *.dat files
65 '''
66
67 ext = '.dat'
68
69 def __init__(self, **kwargs):
70
71 ProcessingUnit.__init__(self, **kwargs)
72
73 self.dataOut = Parameters()
74 self.counter_records = 0
75 self.flagNoMoreFiles = 0
76 self.isConfig = False
77 self.filename = None
78 self.clockpulse = 0.15
79 self.kd = 213.6
80
81 def setup(self,
82 path=None,
83 startDate=None,
84 endDate=None,
85 ext=None,
86 startTime=datetime.time(0, 0, 0),
87 endTime=datetime.time(23, 59, 59),
88 timezone=0,
89 format=None,
90 **kwargs):
91
92 self.path = path
93 self.startDate = startDate
94 self.endDate = endDate
95 self.startTime = startTime
96 self.endTime = endTime
97 self.datatime = datetime.datetime(1900, 1, 1)
98 self.format = format
99
100 if self.path is None:
101 raise ValueError, "The path is not valid"
102
103 if ext is None:
104 ext = self.ext
105
106 self.search_files(self.path, startDate, endDate, ext)
107 self.timezone = timezone
108 self.fileIndex = 0
109
110 if not self.fileList:
111 log.warning('There is no files matching these date in the folder: {}'.format(
112 path), self.name)
113
114 self.setNextFile()
115
116 def search_files(self, path, startDate, endDate, ext):
117 '''
118 Searching for BLTR rawdata file in path
119 Creating a list of file to proces included in [startDate,endDate]
120
121 Input:
122 path - Path to find BLTR rawdata files
123 startDate - Select file from this date
124 enDate - Select file until this date
125 ext - Extension of the file to read
126 '''
127
128 log.success('Searching files in {} '.format(path), self.name)
129 fileList0 = glob.glob1(path, '{}*{}'.format(self.format.upper(), ext))
130 fileList0.sort()
131
132 self.fileList = []
133 self.dateFileList = []
134
135 for thisFile in fileList0:
136 year = thisFile[2:4]
137 if not isNumber(year):
138 continue
139
140 month = thisFile[4:6]
141 if not isNumber(month):
142 continue
143
144 day = thisFile[6:8]
145 if not isNumber(day):
146 continue
147
148 year, month, day = int(year), int(month), int(day)
149 dateFile = datetime.date(year+2000, month, day)
150
151 if (startDate > dateFile) or (endDate < dateFile):
152 continue
153
154 self.fileList.append(thisFile)
155 self.dateFileList.append(dateFile)
156
157 return
158
159 def setNextFile(self):
160
161 file_id = self.fileIndex
162
163 if file_id == len(self.fileList):
164 log.success('No more files in the folder', self.name)
165 self.flagNoMoreFiles = 1
166 return 0
167
168 log.success('Opening {}'.format(self.fileList[file_id]), self.name)
169 filename = os.path.join(self.path, self.fileList[file_id])
170
171 dirname, name = os.path.split(filename)
172 self.siteFile = name.split('.')[0]
173 if self.filename is not None:
174 self.fp.close()
175 self.filename = filename
176 self.fp = open(self.filename, 'rb')
177 self.t0 = [7, 0, 0]
178 self.tf = [18, 0, 0]
179
180 self.header_file = numpy.fromfile(self.fp, FILE_HEADER_STRUCTURE, 1)
181 yy = self.header_file['year'] - 1900 * (self.header_file['year'] > 3000)
182 self.year = int(yy + 1900 * (yy < 1000))
183 self.doy = int(self.header_file['doy'])
184 self.H0 = self.clockpulse*numpy.round(self.header_file['fh']/self.clockpulse)
185 self.dH = self.clockpulse*numpy.round(self.header_file['dh']/self.clockpulse)
186 self.ipp = self.clockpulse*numpy.round(self.header_file['ipp']/self.clockpulse)
187
188 tau = self.ipp / 1.5E5
189
190 self.nheights = int(self.header_file['nheights'])
191 self.heights = numpy.arange(self.nheights) * self.dH + self.H0
192
193 self.sizeOfFile = os.path.getsize(self.filename)
194 self.counter_records = 0
195 self.flagIsNewFile = 0
196 self.fileIndex += 1
197
198 return 1
199
200 def readNextBlock(self):
201
202 while True:
203 if self.fp.tell() == self.sizeOfFile:
204 self.flagIsNewFile = 1
205 if not self.setNextFile():
206 return 0
207
208 self.readBlock()
209
210 if (self.datatime < datetime.datetime.combine(self.startDate, self.startTime)) or \
211 (self.datatime > datetime.datetime.combine(self.endDate, self.endTime)):
212 log.warning(
213 'Reading Record No. {} -> {} [Skipping]'.format(
214 self.counter_records,
215 self.datatime.ctime()),
216 self.name)
217 continue
218 break
219
220 log.log('Reading Record No. {} -> {}'.format(
221 self.counter_records,
222 self.datatime.ctime()), self.name)
223
224 return 1
225
226 def readBlock(self):
227
228 header_rec = numpy.fromfile(self.fp, REC_HEADER_STRUCTURE, 1)
229 self.timedelta = header_rec['timedelta']
230 if header_rec['magic'] == 888.:
231 h0 = self.clockpulse * numpy.round(header_rec['h0'] / self.clockpulse)
232 nheights = int(header_rec['nheights'])
233 hours = float(header_rec['hours'][0])
234 self.datatime = datetime.datetime(self.year, 1, 1) + datetime.timedelta(days=self.doy-1, hours=hours)
235 self.time = (self.datatime - datetime.datetime(1970,1,1)).total_seconds()
236
237 buffer = numpy.fromfile(self.fp, 'f', 8*nheights).reshape(nheights, 8)
238 idh0 = int(numpy.round((self.H0 - h0) / self.dH))
239 if idh0 == 0 and buffer[0,0] > 1E8:
240 buffer[0,0] = buffer[1,0]
241
242 pow0 = buffer[:, 0]
243 pow1 = buffer[:, 1]
244 acf0 = (buffer[:,2] + buffer[:,3]*1j) / pow0
245 acf1 = (buffer[:,4] + buffer[:,5]*1j) / pow1
246 dccf = (buffer[:,6] + buffer[:,7]*1j) / (pow0*pow1)
247
248 ### SNR
249 sno = (pow0 + pow1 - header_rec['snr']) / header_rec['snr']
250 sno10 = numpy.log10(sno)
251 dsno = 1.0 / numpy.sqrt(self.header_file['nint'] * self.header_file['navg']) * (1 + (1 / sno))
252
253 ### Vertical Drift
254 sp = numpy.sqrt(numpy.abs(acf0)*numpy.abs(acf1))
255 sp[numpy.where(numpy.abs(sp) >= 1.0)] = numpy.sqrt(0.9999)
256
257 vzo = -numpy.arctan2(acf0.imag + acf1.imag,acf0.real + acf1.real)*1.5E5*1.5/(self.ipp*numpy.pi)
258 dvzo = numpy.sqrt(1.0 - sp*sp)*0.338*1.5E5/(numpy.sqrt(self.header_file['nint']*self.header_file['navg'])*sp*self.ipp)
259 err = numpy.where(dvzo <= 0.1)
260 dvzo[err] = 0.1
261
262 #Zonal Drifts
263 dt = self.header_file['nint']*self.ipp / 1.5E5
264 coh = numpy.sqrt(numpy.abs(dccf))
265 err = numpy.where(coh >= 1.0)
266 coh[err] = numpy.sqrt(0.99999)
267
268 err = numpy.where(coh <= 0.1)
269 coh[err] = numpy.sqrt(0.1)
270
271 vxo = numpy.arctan2(dccf.imag, dccf.real)*self.heights[idh0]*1.0E3/(self.kd*dt)
272 dvxo = numpy.sqrt(1.0 - coh*coh)*self.heights[idh0]*1.0E3/(numpy.sqrt(self.header_file['nint']*self.header_file['navg'])*coh*self.kd*dt)
273
274 err = numpy.where(dvxo <= 0.1)
275 dvxo[err] = 0.1
276
277 N = range(len(pow0))
278
279 self.buffer = numpy.zeros((6, self.nheights)) + numpy.nan
280
281 self.buffer[0, idh0+numpy.array(N)] = sno10
282 self.buffer[1, idh0+numpy.array(N)] = vzo
283 self.buffer[2, idh0+numpy.array(N)] = vxo
284 self.buffer[3, idh0+numpy.array(N)] = dvzo
285 self.buffer[4, idh0+numpy.array(N)] = dvxo
286 self.buffer[5, idh0+numpy.array(N)] = dsno
287
288 self.counter_records += 1
289
290 return
291
292 def readHeader(self):
293 '''
294 RecordHeader of BLTR rawdata file
295 '''
296
297 header_structure = numpy.dtype(
298 REC_HEADER_STRUCTURE.descr + [
299 ('antenna_coord', 'f4', (2, self.nchannels)),
300 ('rx_gains', 'u4', (self.nchannels,)),
301 ('rx_analysis', 'u4', (self.nchannels,))
302 ]
303 )
304
305 self.header_rec = numpy.fromfile(self.fp, header_structure, 1)
306 self.lat = self.header_rec['lat'][0]
307 self.lon = self.header_rec['lon'][0]
308 self.delta = self.header_rec['delta_r'][0]
309 self.correction = self.header_rec['dmode_rngcorr'][0]
310 self.imode = self.header_rec['dmode_index'][0]
311 self.antenna = self.header_rec['antenna_coord']
312 self.rx_gains = self.header_rec['rx_gains']
313 self.time = self.header_rec['time'][0]
314 dt = datetime.datetime.utcfromtimestamp(self.time)
315 if dt.date()>self.datatime.date():
316 self.flagDiscontinuousBlock = 1
317 self.datatime = dt
318
319 def readData(self):
320 '''
321 Reading and filtering data block record of BLTR rawdata file, filtering is according to status_value.
322
323 Input:
324 status_value - Array data is set to NAN for values that are not equal to status_value
325
326 '''
327
328 data_structure = numpy.dtype(
329 DATA_STRUCTURE.descr + [
330 ('rx_saturation', 'u4', (self.nchannels,)),
331 ('chan_offset', 'u4', (2 * self.nchannels,)),
332 ('rx_amp', 'u4', (self.nchannels,)),
333 ('rx_snr', 'f4', (self.nchannels,)),
334 ('cross_snr', 'f4', (self.kchan,)),
335 ('sea_power_relative', 'f4', (self.kchan,))]
336 )
337
338 data = numpy.fromfile(self.fp, data_structure, self.nranges)
339
340 height = data['range']
341 winds = numpy.array(
342 (data['zonal'], data['meridional'], data['vertical']))
343 snr = data['rx_snr'].T
344
345 winds[numpy.where(winds == -9999.)] = numpy.nan
346 winds[:, numpy.where(data['status'] != self.status_value)] = numpy.nan
347 snr[numpy.where(snr == -9999.)] = numpy.nan
348 snr[:, numpy.where(data['status'] != self.status_value)] = numpy.nan
349 snr = numpy.power(10, snr / 10)
350
351 return height, winds, snr
352
353 def set_output(self):
354 '''
355 Storing data from databuffer to dataOut object
356 '''
357
358 self.dataOut.data_SNR = self.buffer[0]
359 self.dataOut.heights = self.heights
360 self.dataOut.data_param = self.buffer[1:,]
361 self.dataOut.utctimeInit = self.time
362 self.dataOut.utctime = self.dataOut.utctimeInit
363 self.dataOut.useLocalTime = False
364 self.dataOut.paramInterval = self.timedelta
365 self.dataOut.timezone = self.timezone
366 self.dataOut.sizeOfFile = self.sizeOfFile
367 self.dataOut.flagNoData = False
368 self.dataOut.flagDiscontinuousBlock = self.flagDiscontinuousBlock
369
370 def getData(self):
371 '''
372 Storing data from databuffer to dataOut object
373 '''
374 if self.flagNoMoreFiles:
375 self.dataOut.flagNoData = True
376 log.success('No file left to process', self.name)
377 return 0
378
379 if not self.readNextBlock():
380 self.dataOut.flagNoData = True
381 return 0
382
383 self.set_output()
384
385 return 1
General Comments 0
You need to be logged in to leave comments. Login now