##// END OF EJS Templates
schain-jc
RSS
Review BLTRSpectraReader Unit
Review BLTRSpectraReader Unit
Juan C. Espinoza - - Load All Authors

File last commit:

r1199:13d983fa985c
r1199:13d983fa985c
Show More
bltrIO_spectra.py
466 lines | 18.7 KiB | text/x-python | PythonLexer
/ schainpy / model / io / bltrIO_spectra.py
History | Source | Raw |Copy content |Copy permalink
Juan C. Espinoza
Review BLTRSpectraReader Unit
 r1199 import os
import sys
import glob
import fnmatch
import datetime
import time
import re
import h5py
import numpy
import pylab as plb
from scipy.optimize import curve_fit
from scipy import asarray as ar, exp
SPEED_OF_LIGHT = 299792458
SPEED_OF_LIGHT = 3e8
try:
from gevent import sleep
except:
from time import sleep
from .utils import folder_in_range
from schainpy.model.data.jrodata import Spectra
from schainpy.model.proc.jroproc_base import ProcessingUnit, Operation, MPDecorator
from schainpy.utils import log
from schainpy.model.io.jroIO_base import JRODataReader
def pol2cart(rho, phi):
x = rho * numpy.cos(phi)
y = rho * numpy.sin(phi)
return(x, y)
FILE_STRUCTURE = numpy.dtype([ # HEADER 48bytes
('FileMgcNumber', '<u4'), # 0x23020100
('nFDTdataRecors', '<u4'),
('OffsetStartHeader', '<u4'),
('RadarUnitId', '<u4'),
('SiteName', 'S32'), # Null terminated
])
class FileHeaderBLTR():
def __init__(self, fo):
self.fo = fo
self.size = 48
self.read()
def read(self):
header = numpy.fromfile(self.fo, FILE_STRUCTURE, 1)
self.FileMgcNumber = hex(header['FileMgcNumber'][0])
self.nFDTdataRecors = int(header['nFDTdataRecors'][0])
self.RadarUnitId = int(header['RadarUnitId'][0])
self.OffsetStartHeader = int(header['OffsetStartHeader'][0])
self.SiteName = header['SiteName'][0]
def write(self, fp):
headerTuple = (self.FileMgcNumber,
self.nFDTdataRecors,
self.RadarUnitId,
self.SiteName,
self.size)
header = numpy.array(headerTuple, FILE_STRUCTURE)
header.tofile(fp)
''' ndarray.tofile(fid, sep, format) Write array to a file as text or binary (default).
fid : file or str
An open file object, or a string containing a filename.
sep : str
Separator between array items for text output. If "" (empty), a binary file is written,
equivalent to file.write(a.tobytes()).
format : str
Format string for text file output. Each entry in the array is formatted to text by
first converting it to the closest Python type, and then using "format" % item.
'''
return 1
RECORD_STRUCTURE = numpy.dtype([ # RECORD HEADER 180+20N bytes
('RecMgcNumber', '<u4'), # 0x23030001
('RecCounter', '<u4'), # Record counter(0,1, ...)
# Offset to start of next record form start of this record
('Off2StartNxtRec', '<u4'),
# Offset to start of data from start of this record
('Off2StartData', '<u4'),
# Epoch time stamp of start of acquisition (seconds)
('nUtime', '<i4'),
# Millisecond component of time stamp (0,...,999)
('nMilisec', '<u4'),
# Experiment tag name (null terminated)
('ExpTagName', 'S32'),
# Experiment comment (null terminated)
('ExpComment', 'S32'),
# Site latitude (from GPS) in degrees (positive implies North)
('SiteLatDegrees', '<f4'),
# Site longitude (from GPS) in degrees (positive implies East)
('SiteLongDegrees', '<f4'),
# RTC GPS engine status (0=SEEK, 1=LOCK, 2=NOT FITTED, 3=UNAVAILABLE)
('RTCgpsStatus', '<u4'),
('TransmitFrec', '<u4'), # Transmit frequency (Hz)
('ReceiveFrec', '<u4'), # Receive frequency
# First local oscillator frequency (Hz)
('FirstOsciFrec', '<u4'),
# (0="O", 1="E", 2="linear 1", 3="linear2")
('Polarisation', '<u4'),
# Receiver filter settings (0,1,2,3)
('ReceiverFiltSett', '<u4'),
# Number of modes in use (1 or 2)
('nModesInUse', '<u4'),
# Dual Mode index number for these data (0 or 1)
('DualModeIndex', '<u4'),
# Dual Mode range correction for these data (m)
('DualModeRange', '<u4'),
# Number of digital channels acquired (2*N)
('nDigChannels', '<u4'),
# Sampling resolution (meters)
('SampResolution', '<u4'),
# Number of range gates sampled
('nHeights', '<u4'),
# Start range of sampling (meters)
('StartRangeSamp', '<u4'),
('PRFhz', '<u4'), # PRF (Hz)
('nCohInt', '<u4'), # Integrations
# Number of data points transformed
('nProfiles', '<u4'),
# Number of receive beams stored in file (1 or N)
('nChannels', '<u4'),
('nIncohInt', '<u4'), # Number of spectral averages
# FFT windowing index (0 = no window)
('FFTwindowingInd', '<u4'),
# Beam steer angle (azimuth) in degrees (clockwise from true North)
('BeamAngleAzim', '<f4'),
# Beam steer angle (zenith) in degrees (0=> vertical)
('BeamAngleZen', '<f4'),
# Antenna coordinates (Range(meters), Bearing(degrees)) - N pairs
('AntennaCoord0', '<f4'),
# Antenna coordinates (Range(meters), Bearing(degrees)) - N pairs
('AntennaAngl0', '<f4'),
# Antenna coordinates (Range(meters), Bearing(degrees)) - N pairs
('AntennaCoord1', '<f4'),
# Antenna coordinates (Range(meters), Bearing(degrees)) - N pairs
('AntennaAngl1', '<f4'),
# Antenna coordinates (Range(meters), Bearing(degrees)) - N pairs
('AntennaCoord2', '<f4'),
# Antenna coordinates (Range(meters), Bearing(degrees)) - N pairs
('AntennaAngl2', '<f4'),
# Receiver phase calibration (degrees) - N values
('RecPhaseCalibr0', '<f4'),
# Receiver phase calibration (degrees) - N values
('RecPhaseCalibr1', '<f4'),
# Receiver phase calibration (degrees) - N values
('RecPhaseCalibr2', '<f4'),
# Receiver amplitude calibration (ratio relative to receiver one) - N values
('RecAmpCalibr0', '<f4'),
# Receiver amplitude calibration (ratio relative to receiver one) - N values
('RecAmpCalibr1', '<f4'),
# Receiver amplitude calibration (ratio relative to receiver one) - N values
('RecAmpCalibr2', '<f4'),
# Receiver gains in dB - N values
('ReceiverGaindB0', '<i4'),
# Receiver gains in dB - N values
('ReceiverGaindB1', '<i4'),
# Receiver gains in dB - N values
('ReceiverGaindB2', '<i4'),
])
class RecordHeaderBLTR():
def __init__(self, fo):
self.fo = fo
self.OffsetStartHeader = 48
self.Off2StartNxtRec = 811248
def read(self, block):
OffRHeader = self.OffsetStartHeader + block * self.Off2StartNxtRec
self.fo.seek(OffRHeader, os.SEEK_SET)
header = numpy.fromfile(self.fo, RECORD_STRUCTURE, 1)
self.RecMgcNumber = hex(header['RecMgcNumber'][0]) # 0x23030001
self.RecCounter = int(header['RecCounter'][0])
self.Off2StartNxtRec = int(header['Off2StartNxtRec'][0])
self.Off2StartData = int(header['Off2StartData'][0])
self.nUtime = header['nUtime'][0]
self.nMilisec = header['nMilisec'][0]
self.ExpTagName = '' # str(header['ExpTagName'][0])
self.ExpComment = '' # str(header['ExpComment'][0])
self.SiteLatDegrees = header['SiteLatDegrees'][0]
self.SiteLongDegrees = header['SiteLongDegrees'][0]
self.RTCgpsStatus = header['RTCgpsStatus'][0]
self.TransmitFrec = header['TransmitFrec'][0]
self.ReceiveFrec = header['ReceiveFrec'][0]
self.FirstOsciFrec = header['FirstOsciFrec'][0]
self.Polarisation = header['Polarisation'][0]
self.ReceiverFiltSett = header['ReceiverFiltSett'][0]
self.nModesInUse = header['nModesInUse'][0]
self.DualModeIndex = header['DualModeIndex'][0]
self.DualModeRange = header['DualModeRange'][0]
self.nDigChannels = header['nDigChannels'][0]
self.SampResolution = header['SampResolution'][0]
self.nHeights = header['nHeights'][0]
self.StartRangeSamp = header['StartRangeSamp'][0]
self.PRFhz = header['PRFhz'][0]
self.nCohInt = header['nCohInt'][0]
self.nProfiles = header['nProfiles'][0]
self.nChannels = header['nChannels'][0]
self.nIncohInt = header['nIncohInt'][0]
self.FFTwindowingInd = header['FFTwindowingInd'][0]
self.BeamAngleAzim = header['BeamAngleAzim'][0]
self.BeamAngleZen = header['BeamAngleZen'][0]
self.AntennaCoord0 = header['AntennaCoord0'][0]
self.AntennaAngl0 = header['AntennaAngl0'][0]
self.AntennaCoord1 = header['AntennaCoord1'][0]
self.AntennaAngl1 = header['AntennaAngl1'][0]
self.AntennaCoord2 = header['AntennaCoord2'][0]
self.AntennaAngl2 = header['AntennaAngl2'][0]
self.RecPhaseCalibr0 = header['RecPhaseCalibr0'][0]
self.RecPhaseCalibr1 = header['RecPhaseCalibr1'][0]
self.RecPhaseCalibr2 = header['RecPhaseCalibr2'][0]
self.RecAmpCalibr0 = header['RecAmpCalibr0'][0]
self.RecAmpCalibr1 = header['RecAmpCalibr1'][0]
self.RecAmpCalibr2 = header['RecAmpCalibr2'][0]
self.ReceiverGaindB0 = header['ReceiverGaindB0'][0]
self.ReceiverGaindB1 = header['ReceiverGaindB1'][0]
self.ReceiverGaindB2 = header['ReceiverGaindB2'][0]
self.ipp = 0.5 * (SPEED_OF_LIGHT / self.PRFhz)
self.RHsize = 180 + 20 * self.nChannels
self.Datasize = self.nProfiles * self.nChannels * self.nHeights * 2 * 4
endFp = self.OffsetStartHeader + self.RecCounter * self.Off2StartNxtRec
if OffRHeader > endFp:
sys.stderr.write(
"Warning %s: Size value read from System Header is lower than it has to be\n" % fp)
return 0
if OffRHeader < endFp:
sys.stderr.write(
"Warning %s: Size value read from System Header size is greater than it has to be\n" % fp)
return 0
return 1
@MPDecorator
class BLTRSpectraReader (ProcessingUnit):
def __init__(self):
ProcessingUnit.__init__(self)
self.ext = ".fdt"
self.optchar = "P"
self.fpFile = None
self.fp = None
self.BlockCounter = 0
self.fileSizeByHeader = None
self.filenameList = []
self.fileSelector = 0
self.Off2StartNxtRec = 0
self.RecCounter = 0
self.flagNoMoreFiles = 0
self.data_spc = None
self.data_cspc = None
self.path = None
self.OffsetStartHeader = 0
self.Off2StartData = 0
self.ipp = 0
self.nFDTdataRecors = 0
self.blocksize = 0
self.dataOut = Spectra()
self.dataOut.flagNoData = False
def search_files(self):
'''
Function that indicates the number of .fdt files that exist in the folder to be read.
It also creates an organized list with the names of the files to read.
'''
files = glob.glob(os.path.join(self.path, '*{}'.format(self.ext)))
files = sorted(files)
for f in files:
filename = f.split('/')[-1]
if folder_in_range(filename.split('.')[1], self.startDate, self.endDate, '%Y%m%d'):
self.filenameList.append(f)
def run(self, **kwargs):
'''
This method will be the one that will initiate the data entry, will be called constantly.
You should first verify that your Setup () is set up and then continue to acquire
the data to be processed with getData ().
'''
if not self.isConfig:
self.setup(**kwargs)
self.isConfig = True
self.getData()
def setup(self,
path=None,
startDate=None,
endDate=None,
startTime=None,
endTime=None,
walk=True,
code=None,
online=False,
mode=None,
**kwargs):
self.isConfig = True
self.path = path
self.startDate = startDate
self.endDate = endDate
self.startTime = startTime
self.endTime = endTime
self.walk = walk
self.mode = int(mode)
self.search_files()
if self.filenameList:
self.readFile()
def getData(self):
'''
Before starting this function, you should check that there is still an unread file,
If there are still blocks to read or if the data block is empty.
You should call the file "read".
'''
if self.flagNoMoreFiles:
self.dataOut.flagNoData = True
self.dataOut.error = 'No more files'
self.readBlock()
def readFile(self):
'''
You must indicate if you are reading in Online or Offline mode and load the
The parameters for this file reading mode.
Then you must do 2 actions:
1. Get the BLTR FileHeader.
2. Start reading the first block.
'''
if self.fileSelector < len(self.filenameList):
log.success('Opening file: {}'.format(self.filenameList[self.fileSelector]), self.name)
self.fp = open(self.filenameList[self.fileSelector])
self.fheader = FileHeaderBLTR(self.fp)
self.rheader = RecordHeaderBLTR(self.fp)
self.nFDTdataRecors = self.fheader.nFDTdataRecors
self.fileSelector += 1
self.BlockCounter = 0
return 1
else:
self.flagNoMoreFiles = True
self.dataOut.flagNoData = True
return 0
def readBlock(self):
'''
It should be checked if the block has data, if it is not passed to the next file.
Then the following is done:
1. Read the RecordHeader
2. Fill the buffer with the current block number.
'''
if self.BlockCounter == self.nFDTdataRecors:
if not self.readFile():
return
if self.mode == 1:
self.rheader.read(self.BlockCounter+1)
elif self.mode == 0:
self.rheader.read(self.BlockCounter)
self.RecCounter = self.rheader.RecCounter
self.OffsetStartHeader = self.rheader.OffsetStartHeader
self.Off2StartNxtRec = self.rheader.Off2StartNxtRec
self.Off2StartData = self.rheader.Off2StartData
self.nProfiles = self.rheader.nProfiles
self.nChannels = self.rheader.nChannels
self.nHeights = self.rheader.nHeights
self.frequency = self.rheader.TransmitFrec
self.DualModeIndex = self.rheader.DualModeIndex
self.pairsList = [(0, 1), (0, 2), (1, 2)]
self.dataOut.pairsList = self.pairsList
self.nRdPairs = len(self.dataOut.pairsList)
self.dataOut.nRdPairs = self.nRdPairs
self.dataOut.heightList = (self.rheader.StartRangeSamp + numpy.arange(self.nHeights) * self.rheader.SampResolution) / 1000.
self.dataOut.channelList = range(self.nChannels)
self.dataOut.nProfiles=self.rheader.nProfiles
self.dataOut.nIncohInt=self.rheader.nIncohInt
self.dataOut.nCohInt=self.rheader.nCohInt
self.dataOut.ippSeconds= 1/float(self.rheader.PRFhz)
self.dataOut.PRF=self.rheader.PRFhz
self.dataOut.nFFTPoints=self.rheader.nProfiles
self.dataOut.utctime=self.rheader.nUtime
self.dataOut.timeZone = 0
self.dataOut.useLocalTime = False
self.dataOut.nmodes = 2
log.log('Reading block {} - {}'.format(self.BlockCounter, self.dataOut.datatime), self.name)
OffDATA = self.OffsetStartHeader + self.RecCounter * \
self.Off2StartNxtRec + self.Off2StartData
self.fp.seek(OffDATA, os.SEEK_SET)
self.data_fft = numpy.fromfile(self.fp, [('complex','<c8')], self.nProfiles*self.nChannels*self.nHeights )
self.data_fft = self.data_fft.astype(numpy.dtype('complex'))
self.data_block = numpy.reshape(self.data_fft,(self.nHeights, self.nChannels, self.nProfiles))
self.data_block = numpy.transpose(self.data_block, (1,2,0))
copy = self.data_block.copy()
spc = copy * numpy.conjugate(copy)
self.data_spc = numpy.absolute(spc) # valor absoluto o magnitud
self.dataOut.data_spc = self.data_spc
cspc = self.data_block.copy()
self.data_cspc = self.data_block.copy()
for i in range(self.nRdPairs):
chan_index0 = self.dataOut.pairsList[i][0]
chan_index1 = self.dataOut.pairsList[i][1]
self.data_cspc[i, :, :] = cspc[chan_index0, :, :] * numpy.conjugate(cspc[chan_index1, :, :])
'''Getting Eij and Nij'''
(AntennaX0, AntennaY0) = pol2cart(
self.rheader.AntennaCoord0, self.rheader.AntennaAngl0 * numpy.pi / 180)
(AntennaX1, AntennaY1) = pol2cart(
self.rheader.AntennaCoord1, self.rheader.AntennaAngl1 * numpy.pi / 180)
(AntennaX2, AntennaY2) = pol2cart(
self.rheader.AntennaCoord2, self.rheader.AntennaAngl2 * numpy.pi / 180)
E01 = AntennaX0 - AntennaX1
N01 = AntennaY0 - AntennaY1
E02 = AntennaX0 - AntennaX2
N02 = AntennaY0 - AntennaY2
E12 = AntennaX1 - AntennaX2
N12 = AntennaY1 - AntennaY2
self.ChanDist = numpy.array(
[[E01, N01], [E02, N02], [E12, N12]])
self.dataOut.ChanDist = self.ChanDist
self.BlockCounter += 2
self.dataOut.data_spc = self.data_spc
self.dataOut.data_cspc =self.data_cspc