''' Created on Set 10, 2017 @author: Juan C. Espinoza ''' import os import sys import time import glob import datetime import numpy from schainpy.model.proc.jroproc_base import ProcessingUnit from schainpy.model.data.jrodata import Parameters from schainpy.model.io.jroIO_base import JRODataReader, isNumber from schainpy.utils import log FILE_HEADER_STRUCTURE = numpy.dtype([ ('year', 'f'), ('doy', 'f'), ('nint', 'f'), ('navg', 'f'), ('fh', 'f'), ('dh', 'f'), ('nheights', 'f'), ('ipp', 'f') ]) REC_HEADER_STRUCTURE = numpy.dtype([ ('magic', 'f'), ('hours', 'f'), ('timedelta', 'f'), ('h0', 'f'), ('nheights', 'f'), ('snr1', 'f'), ('snr2', 'f'), ('snr', 'f'), ]) DATA_STRUCTURE = numpy.dtype([ ('range', ' dateFile) or (endDate < dateFile): continue self.fileList.append(thisFile) self.dateFileList.append(dateFile) return def setNextFile(self): file_id = self.fileIndex if file_id == len(self.fileList): log.success('No more files in the folder', self.name) self.flagNoMoreFiles = 1 return 0 log.success('Opening {}'.format(self.fileList[file_id]), self.name) filename = os.path.join(self.path, self.fileList[file_id]) dirname, name = os.path.split(filename) self.siteFile = name.split('.')[0] if self.filename is not None: self.fp.close() self.filename = filename self.fp = open(self.filename, 'rb') self.t0 = [7, 0, 0] self.tf = [18, 0, 0] self.header_file = numpy.fromfile(self.fp, FILE_HEADER_STRUCTURE, 1) yy = self.header_file['year'] - 1900 * (self.header_file['year'] > 3000) self.year = int(yy + 1900 * (yy < 1000)) self.doy = int(self.header_file['doy']) self.H0 = self.clockpulse*numpy.round(self.header_file['fh']/self.clockpulse) self.dH = self.clockpulse*numpy.round(self.header_file['dh']/self.clockpulse) self.ipp = self.clockpulse*numpy.round(self.header_file['ipp']/self.clockpulse) tau = self.ipp / 1.5E5 self.nheights = int(self.header_file['nheights']) self.heights = numpy.arange(self.nheights) * self.dH + self.H0 self.sizeOfFile = os.path.getsize(self.filename) self.counter_records = 0 self.flagIsNewFile = 0 self.fileIndex += 1 return 1 def readNextBlock(self): while True: if self.fp.tell() == self.sizeOfFile: self.flagIsNewFile = 1 if not self.setNextFile(): return 0 self.readBlock() if (self.datatime < datetime.datetime.combine(self.startDate, self.startTime)) or \ (self.datatime > datetime.datetime.combine(self.endDate, self.endTime)): log.warning( 'Reading Record No. {} -> {} [Skipping]'.format( self.counter_records, self.datatime.ctime()), self.name) continue break log.log('Reading Record No. {} -> {}'.format( self.counter_records, self.datatime.ctime()), self.name) return 1 def readBlock(self): header_rec = numpy.fromfile(self.fp, REC_HEADER_STRUCTURE, 1) self.timedelta = header_rec['timedelta'] if header_rec['magic'] == 888.: h0 = self.clockpulse * numpy.round(header_rec['h0'] / self.clockpulse) nheights = int(header_rec['nheights']) hours = float(header_rec['hours'][0]) self.datatime = datetime.datetime(self.year, 1, 1) + datetime.timedelta(days=self.doy-1, hours=hours) self.time = (self.datatime - datetime.datetime(1970,1,1)).total_seconds() buffer = numpy.fromfile(self.fp, 'f', 8*nheights).reshape(nheights, 8) idh0 = int(numpy.round((self.H0 - h0) / self.dH)) if idh0 == 0 and buffer[0,0] > 1E8: buffer[0,0] = buffer[1,0] pow0 = buffer[:, 0] pow1 = buffer[:, 1] acf0 = (buffer[:,2] + buffer[:,3]*1j) / pow0 acf1 = (buffer[:,4] + buffer[:,5]*1j) / pow1 dccf = (buffer[:,6] + buffer[:,7]*1j) / (pow0*pow1) ### SNR sno = (pow0 + pow1 - header_rec['snr']) / header_rec['snr'] sno10 = numpy.log10(sno) dsno = 1.0 / numpy.sqrt(self.header_file['nint'] * self.header_file['navg']) * (1 + (1 / sno)) ### Vertical Drift sp = numpy.sqrt(numpy.abs(acf0)*numpy.abs(acf1)) sp[numpy.where(numpy.abs(sp) >= 1.0)] = numpy.sqrt(0.9999) vzo = -numpy.arctan2(acf0.imag + acf1.imag,acf0.real + acf1.real)*1.5E5*1.5/(self.ipp*numpy.pi) dvzo = numpy.sqrt(1.0 - sp*sp)*0.338*1.5E5/(numpy.sqrt(self.header_file['nint']*self.header_file['navg'])*sp*self.ipp) err = numpy.where(dvzo <= 0.1) dvzo[err] = 0.1 #Zonal Drifts dt = self.header_file['nint']*self.ipp / 1.5E5 coh = numpy.sqrt(numpy.abs(dccf)) err = numpy.where(coh >= 1.0) coh[err] = numpy.sqrt(0.99999) err = numpy.where(coh <= 0.1) coh[err] = numpy.sqrt(0.1) vxo = numpy.arctan2(dccf.imag, dccf.real)*self.heights[idh0]*1.0E3/(self.kd*dt) 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) err = numpy.where(dvxo <= 0.1) dvxo[err] = 0.1 N = range(len(pow0)) self.buffer = numpy.zeros((6, self.nheights)) + numpy.nan self.buffer[0, idh0+numpy.array(N)] = sno10 self.buffer[1, idh0+numpy.array(N)] = vzo self.buffer[2, idh0+numpy.array(N)] = vxo self.buffer[3, idh0+numpy.array(N)] = dvzo self.buffer[4, idh0+numpy.array(N)] = dvxo self.buffer[5, idh0+numpy.array(N)] = dsno self.counter_records += 1 return def readHeader(self): ''' RecordHeader of BLTR rawdata file ''' header_structure = numpy.dtype( REC_HEADER_STRUCTURE.descr + [ ('antenna_coord', 'f4', (2, self.nchannels)), ('rx_gains', 'u4', (self.nchannels,)), ('rx_analysis', 'u4', (self.nchannels,)) ] ) self.header_rec = numpy.fromfile(self.fp, header_structure, 1) self.lat = self.header_rec['lat'][0] self.lon = self.header_rec['lon'][0] self.delta = self.header_rec['delta_r'][0] self.correction = self.header_rec['dmode_rngcorr'][0] self.imode = self.header_rec['dmode_index'][0] self.antenna = self.header_rec['antenna_coord'] self.rx_gains = self.header_rec['rx_gains'] self.time = self.header_rec['time'][0] dt = datetime.datetime.utcfromtimestamp(self.time) if dt.date()>self.datatime.date(): self.flagDiscontinuousBlock = 1 self.datatime = dt def readData(self): ''' Reading and filtering data block record of BLTR rawdata file, filtering is according to status_value. Input: status_value - Array data is set to NAN for values that are not equal to status_value ''' data_structure = numpy.dtype( DATA_STRUCTURE.descr + [ ('rx_saturation', 'u4', (self.nchannels,)), ('chan_offset', 'u4', (2 * self.nchannels,)), ('rx_amp', 'u4', (self.nchannels,)), ('rx_snr', 'f4', (self.nchannels,)), ('cross_snr', 'f4', (self.kchan,)), ('sea_power_relative', 'f4', (self.kchan,))] ) data = numpy.fromfile(self.fp, data_structure, self.nranges) height = data['range'] winds = numpy.array( (data['zonal'], data['meridional'], data['vertical'])) snr = data['rx_snr'].T winds[numpy.where(winds == -9999.)] = numpy.nan winds[:, numpy.where(data['status'] != self.status_value)] = numpy.nan snr[numpy.where(snr == -9999.)] = numpy.nan snr[:, numpy.where(data['status'] != self.status_value)] = numpy.nan snr = numpy.power(10, snr / 10) return height, winds, snr def set_output(self): ''' Storing data from databuffer to dataOut object ''' self.dataOut.data_SNR = self.buffer[0] self.dataOut.heights = self.heights self.dataOut.data_param = self.buffer[1:,] self.dataOut.utctimeInit = self.time self.dataOut.utctime = self.dataOut.utctimeInit self.dataOut.useLocalTime = False self.dataOut.paramInterval = self.timedelta self.dataOut.timezone = self.timezone self.dataOut.sizeOfFile = self.sizeOfFile self.dataOut.flagNoData = False self.dataOut.flagDiscontinuousBlock = self.flagDiscontinuousBlock def getData(self): ''' Storing data from databuffer to dataOut object ''' if self.flagNoMoreFiles: self.dataOut.flagNoData = True log.success('No file left to process', self.name) return 0 if not self.readNextBlock(): self.dataOut.flagNoData = True return 0 self.set_output() return 1