jroIO_bltr.py
1179 lines
| 48.4 KiB
| text/x-python
|
PythonLexer
|
r1075 | import os | |
import sys | |||
|
r965 | import glob | |
import fnmatch | |||
import datetime | |||
import time | |||
import re | |||
import h5py | |||
import numpy | |||
import pylab as plb | |||
from scipy.optimize import curve_fit | |||
|
r1010 | from scipy import asarray as ar, exp | |
|
r965 | from scipy import stats | |
from numpy.ma.core import getdata | |||
SPEED_OF_LIGHT = 299792458 | |||
SPEED_OF_LIGHT = 3e8 | |||
try: | |||
from gevent import sleep | |||
except: | |||
from time import sleep | |||
from schainpy.model.data.jrodata import Spectra | |||
#from schainpy.model.data.BLTRheaderIO import FileHeader, RecordHeader | |||
from schainpy.model.proc.jroproc_base import ProcessingUnit, Operation | |||
#from schainpy.model.io.jroIO_bltr import BLTRReader | |||
from numpy import imag, shape, NaN | |||
|
r1167 | from .jroIO_base import JRODataReader | |
|
r965 | ||
class Header(object): | |||
|
r1075 | ||
|
r965 | def __init__(self): | |
raise NotImplementedError | |||
|
r1075 | ||
|
r965 | def read(self): | |
|
r1075 | ||
|
r965 | raise NotImplementedError | |
|
r1075 | ||
|
r965 | def write(self): | |
|
r1075 | ||
|
r965 | raise NotImplementedError | |
|
r1075 | ||
|
r965 | def printInfo(self): | |
|
r1075 | ||
message = "#" * 50 + "\n" | |||
|
r965 | message += self.__class__.__name__.upper() + "\n" | |
|
r1075 | message += "#" * 50 + "\n" | |
|
r1167 | keyList = list(self.__dict__.keys()) | |
|
r965 | keyList.sort() | |
|
r1075 | ||
|
r965 | for key in keyList: | |
|
r1075 | message += "%s = %s" % (key, self.__dict__[key]) + "\n" | |
|
r965 | if "size" not in keyList: | |
attr = getattr(self, "size") | |||
|
r1075 | if attr: | |
message += "%s = %s" % ("size", attr) + "\n" | |||
|
r965 | ||
|
r1075 | # print message | |
|
r965 | ||
|
r1075 | FILE_STRUCTURE = numpy.dtype([ # HEADER 48bytes | |
('FileMgcNumber', '<u4'), # 0x23020100 | |||
# No Of FDT data records in this file (0 or more) | |||
('nFDTdataRecors', '<u4'), | |||
('OffsetStartHeader', '<u4'), | |||
('RadarUnitId', '<u4'), | |||
('SiteName', numpy.str_, 32), # Null terminated | |||
]) | |||
|
r965 | ||
class FileHeaderBLTR(Header): | |||
|
r1075 | ||
|
r965 | def __init__(self): | |
|
r1075 | ||
self.FileMgcNumber = 0 # 0x23020100 | |||
# No Of FDT data records in this file (0 or more) | |||
self.nFDTdataRecors = 0 | |||
self.RadarUnitId = 0 | |||
self.OffsetStartHeader = 0 | |||
self.SiteName = "" | |||
|
r965 | self.size = 48 | |
|
r1075 | ||
|
r965 | def FHread(self, fp): | |
|
r1075 | # try: | |
startFp = open(fp, "rb") | |||
header = numpy.fromfile(startFp, FILE_STRUCTURE, 1) | |||
|
r1167 | print(' ') | |
print('puntero file header', startFp.tell()) | |||
print(' ') | |||
|
r1075 | ||
|
r965 | ''' numpy.fromfile(file, dtype, count, sep='') | |
file : file or str | |||
Open file object or filename. | |||
dtype : data-type | |||
Data type of the returned array. For binary files, it is used to determine | |||
the size and byte-order of the items in the file. | |||
count : int | |||
Number of items to read. -1 means all items (i.e., the complete file). | |||
sep : str | |||
Separator between items if file is a text file. Empty ("") separator means | |||
the file should be treated as binary. Spaces (" ") in the separator match zero | |||
or more whitespace characters. A separator consisting only of spaces must match | |||
at least one whitespace. | |||
''' | |||
|
r1075 | self.FileMgcNumber = hex(header['FileMgcNumber'][0]) | |
# No Of FDT data records in this file (0 or more) | |||
self.nFDTdataRecors = int(header['nFDTdataRecors'][0]) | |||
self.RadarUnitId = int(header['RadarUnitId'][0]) | |||
self.OffsetStartHeader = int(header['OffsetStartHeader'][0]) | |||
self.SiteName = str(header['SiteName'][0]) | |||
# print 'Numero de bloques', self.nFDTdataRecors | |||
|
r965 | ||
|
r1075 | if self.size < 48: | |
|
r965 | return 0 | |
|
r1075 | ||
|
r965 | return 1 | |
|
r1075 | ||
|
r965 | def write(self, fp): | |
|
r1075 | ||
|
r965 | headerTuple = (self.FileMgcNumber, | |
self.nFDTdataRecors, | |||
self.RadarUnitId, | |||
self.SiteName, | |||
self.size) | |||
|
r1075 | ||
|
r965 | header = numpy.array(headerTuple, FILE_STRUCTURE) | |
# numpy.array(object, dtype=None, copy=True, order=None, subok=False, ndmin=0) | |||
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. | |||
''' | |||
|
r1075 | return 1 | |
|
r965 | ||
|
r1075 | 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', numpy.str_, 32), | |||
# Experiment comment (null terminated) | |||
('ExpComment', numpy.str_, 32), | |||
# 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'), | |||
]) | |||
|
r965 | ||
|
r1075 | class RecordHeaderBLTR(Header): | |
|
r965 | ||
|
r1075 | def __init__(self, RecMgcNumber=None, RecCounter=0, Off2StartNxtRec=811248, | |
nUtime=0, nMilisec=0, ExpTagName=None, | |||
ExpComment=None, SiteLatDegrees=0, SiteLongDegrees=0, | |||
RTCgpsStatus=0, TransmitFrec=0, ReceiveFrec=0, | |||
FirstOsciFrec=0, Polarisation=0, ReceiverFiltSett=0, | |||
nModesInUse=0, DualModeIndex=0, DualModeRange=0, | |||
nDigChannels=0, SampResolution=0, nHeights=0, | |||
StartRangeSamp=0, PRFhz=0, nCohInt=0, | |||
nProfiles=0, nChannels=0, nIncohInt=0, | |||
FFTwindowingInd=0, BeamAngleAzim=0, BeamAngleZen=0, | |||
AntennaCoord0=0, AntennaCoord1=0, AntennaCoord2=0, | |||
RecPhaseCalibr0=0, RecPhaseCalibr1=0, RecPhaseCalibr2=0, | |||
RecAmpCalibr0=0, RecAmpCalibr1=0, RecAmpCalibr2=0, | |||
AntennaAngl0=0, AntennaAngl1=0, AntennaAngl2=0, | |||
ReceiverGaindB0=0, ReceiverGaindB1=0, ReceiverGaindB2=0, Off2StartData=0, OffsetStartHeader=0): | |||
|
r965 | ||
|
r1075 | self.RecMgcNumber = RecMgcNumber # 0x23030001 | |
self.RecCounter = RecCounter | |||
self.Off2StartNxtRec = Off2StartNxtRec | |||
|
r965 | self.Off2StartData = Off2StartData | |
|
r1075 | self.nUtime = nUtime | |
self.nMilisec = nMilisec | |||
self.ExpTagName = ExpTagName | |||
self.ExpComment = ExpComment | |||
self.SiteLatDegrees = SiteLatDegrees | |||
self.SiteLongDegrees = SiteLongDegrees | |||
self.RTCgpsStatus = RTCgpsStatus | |||
self.TransmitFrec = TransmitFrec | |||
self.ReceiveFrec = ReceiveFrec | |||
self.FirstOsciFrec = FirstOsciFrec | |||
self.Polarisation = Polarisation | |||
self.ReceiverFiltSett = ReceiverFiltSett | |||
self.nModesInUse = nModesInUse | |||
self.DualModeIndex = DualModeIndex | |||
self.DualModeRange = DualModeRange | |||
|
r965 | self.nDigChannels = nDigChannels | |
|
r1075 | self.SampResolution = SampResolution | |
self.nHeights = nHeights | |||
self.StartRangeSamp = StartRangeSamp | |||
self.PRFhz = PRFhz | |||
self.nCohInt = nCohInt | |||
self.nProfiles = nProfiles | |||
self.nChannels = nChannels | |||
self.nIncohInt = nIncohInt | |||
self.FFTwindowingInd = FFTwindowingInd | |||
self.BeamAngleAzim = BeamAngleAzim | |||
self.BeamAngleZen = BeamAngleZen | |||
|
r965 | self.AntennaCoord0 = AntennaCoord0 | |
self.AntennaAngl0 = AntennaAngl0 | |||
self.AntennaAngl1 = AntennaAngl1 | |||
self.AntennaAngl2 = AntennaAngl2 | |||
self.AntennaCoord1 = AntennaCoord1 | |||
|
r1075 | self.AntennaCoord2 = AntennaCoord2 | |
|
r965 | self.RecPhaseCalibr0 = RecPhaseCalibr0 | |
self.RecPhaseCalibr1 = RecPhaseCalibr1 | |||
|
r1075 | self.RecPhaseCalibr2 = RecPhaseCalibr2 | |
|
r965 | self.RecAmpCalibr0 = RecAmpCalibr0 | |
self.RecAmpCalibr1 = RecAmpCalibr1 | |||
|
r1075 | self.RecAmpCalibr2 = RecAmpCalibr2 | |
|
r965 | self.ReceiverGaindB0 = ReceiverGaindB0 | |
self.ReceiverGaindB1 = ReceiverGaindB1 | |||
|
r1075 | self.ReceiverGaindB2 = ReceiverGaindB2 | |
self.OffsetStartHeader = 48 | |||
|
r965 | def RHread(self, fp): | |
|
r1075 | # print fp | |
# startFp = open('/home/erick/Documents/Data/huancayo.20161019.22.fdt',"rb") #The method tell() returns the current position of the file read/write pointer within the file. | |||
# The method tell() returns the current position of the file read/write pointer within the file. | |||
startFp = open(fp, "rb") | |||
# RecCounter=0 | |||
# Off2StartNxtRec=811248 | |||
OffRHeader = self.OffsetStartHeader + self.RecCounter * self.Off2StartNxtRec | |||
|
r1167 | print(' ') | |
print('puntero Record Header', startFp.tell()) | |||
print(' ') | |||
|
r1075 | ||
|
r965 | startFp.seek(OffRHeader, os.SEEK_SET) | |
|
r1075 | ||
|
r1167 | print(' ') | |
print('puntero Record Header con seek', startFp.tell()) | |||
print(' ') | |||
|
r1075 | ||
# print 'Posicion del bloque: ',OffRHeader | |||
header = numpy.fromfile(startFp, RECORD_STRUCTURE, 1) | |||
|
r1167 | print(' ') | |
print('puntero Record Header con seek', startFp.tell()) | |||
print(' ') | |||
|
r1075 | ||
|
r1167 | print(' ') | |
|
r965 | # | |
|
r1075 | # print 'puntero Record Header despues de seek', header.tell() | |
|
r1167 | print(' ') | |
|
r1075 | ||
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] | |||
|
r965 | self.nDigChannels = header['nDigChannels'][0] | |
|
r1075 | 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] | |||
|
r965 | self.AntennaCoord0 = header['AntennaCoord0'][0] | |
self.AntennaAngl0 = header['AntennaAngl0'][0] | |||
self.AntennaCoord1 = header['AntennaCoord1'][0] | |||
|
r1075 | self.AntennaAngl1 = header['AntennaAngl1'][0] | |
|
r965 | self.AntennaCoord2 = header['AntennaCoord2'][0] | |
|
r1075 | self.AntennaAngl2 = header['AntennaAngl2'][0] | |
|
r965 | self.RecPhaseCalibr0 = header['RecPhaseCalibr0'][0] | |
self.RecPhaseCalibr1 = header['RecPhaseCalibr1'][0] | |||
|
r1075 | self.RecPhaseCalibr2 = header['RecPhaseCalibr2'][0] | |
|
r965 | self.RecAmpCalibr0 = header['RecAmpCalibr0'][0] | |
self.RecAmpCalibr1 = header['RecAmpCalibr1'][0] | |||
|
r1075 | self.RecAmpCalibr2 = header['RecAmpCalibr2'][0] | |
|
r965 | self.ReceiverGaindB0 = header['ReceiverGaindB0'][0] | |
self.ReceiverGaindB1 = header['ReceiverGaindB1'][0] | |||
self.ReceiverGaindB2 = header['ReceiverGaindB2'][0] | |||
|
r1075 | 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 | |||
# print 'Datasize',self.Datasize | |||
endFp = self.OffsetStartHeader + self.RecCounter * self.Off2StartNxtRec | |||
|
r1167 | print('==============================================') | |
print('RecMgcNumber ', self.RecMgcNumber) | |||
print('RecCounter ', self.RecCounter) | |||
print('Off2StartNxtRec ', self.Off2StartNxtRec) | |||
print('Off2StartData ', self.Off2StartData) | |||
print('Range Resolution ', self.SampResolution) | |||
print('First Height ', self.StartRangeSamp) | |||
print('PRF (Hz) ', self.PRFhz) | |||
print('Heights (K) ', self.nHeights) | |||
print('Channels (N) ', self.nChannels) | |||
print('Profiles (J) ', self.nProfiles) | |||
print('iCoh ', self.nCohInt) | |||
print('iInCoh ', self.nIncohInt) | |||
print('BeamAngleAzim ', self.BeamAngleAzim) | |||
print('BeamAngleZen ', self.BeamAngleZen) | |||
|
r1075 | ||
# print 'ModoEnUso ',self.DualModeIndex | |||
# print 'UtcTime ',self.nUtime | |||
# print 'MiliSec ',self.nMilisec | |||
# print 'Exp TagName ',self.ExpTagName | |||
# print 'Exp Comment ',self.ExpComment | |||
# print 'FFT Window Index ',self.FFTwindowingInd | |||
# print 'N Dig. Channels ',self.nDigChannels | |||
|
r1167 | print('Size de bloque ', self.RHsize) | |
print('DataSize ', self.Datasize) | |||
print('BeamAngleAzim ', self.BeamAngleAzim) | |||
|
r1075 | # print 'AntennaCoord0 ',self.AntennaCoord0 | |
# print 'AntennaAngl0 ',self.AntennaAngl0 | |||
# print 'AntennaCoord1 ',self.AntennaCoord1 | |||
# print 'AntennaAngl1 ',self.AntennaAngl1 | |||
# print 'AntennaCoord2 ',self.AntennaCoord2 | |||
# print 'AntennaAngl2 ',self.AntennaAngl2 | |||
|
r1167 | print('RecPhaseCalibr0 ', self.RecPhaseCalibr0) | |
print('RecPhaseCalibr1 ', self.RecPhaseCalibr1) | |||
print('RecPhaseCalibr2 ', self.RecPhaseCalibr2) | |||
print('RecAmpCalibr0 ', self.RecAmpCalibr0) | |||
print('RecAmpCalibr1 ', self.RecAmpCalibr1) | |||
print('RecAmpCalibr2 ', self.RecAmpCalibr2) | |||
print('ReceiverGaindB0 ', self.ReceiverGaindB0) | |||
print('ReceiverGaindB1 ', self.ReceiverGaindB1) | |||
print('ReceiverGaindB2 ', self.ReceiverGaindB2) | |||
print('==============================================') | |||
|
r1075 | ||
|
r965 | if OffRHeader > endFp: | |
|
r1075 | sys.stderr.write( | |
"Warning %s: Size value read from System Header is lower than it has to be\n" % fp) | |||
|
r965 | return 0 | |
|
r1075 | ||
|
r965 | if OffRHeader < endFp: | |
|
r1075 | sys.stderr.write( | |
"Warning %s: Size value read from System Header size is greater than it has to be\n" % fp) | |||
|
r965 | return 0 | |
|
r1075 | ||
|
r965 | return 1 | |
|
r1075 | ||
|
r1010 | class BLTRSpectraReader (ProcessingUnit, FileHeaderBLTR, RecordHeaderBLTR, JRODataReader): | |
|
r1075 | ||
|
r965 | path = None | |
startDate = None | |||
endDate = None | |||
startTime = None | |||
endTime = None | |||
walk = None | |||
isConfig = False | |||
|
r1075 | ||
fileList = None | |||
# metadata | |||
TimeZone = None | |||
Interval = None | |||
heightList = None | |||
# data | |||
data = None | |||
utctime = None | |||
|
r1001 | def __init__(self, **kwargs): | |
|
r1075 | ||
# Eliminar de la base la herencia | |||
|
r1001 | ProcessingUnit.__init__(self, **kwargs) | |
|
r1075 | ||
|
r1010 | #self.isConfig = False | |
|
r1075 | ||
|
r965 | #self.pts2read_SelfSpectra = 0 | |
#self.pts2read_CrossSpectra = 0 | |||
#self.pts2read_DCchannels = 0 | |||
#self.datablock = None | |||
self.utc = None | |||
self.ext = ".fdt" | |||
self.optchar = "P" | |||
|
r1075 | self.fpFile = None | |
|
r965 | self.fp = None | |
|
r1075 | self.BlockCounter = 0 | |
|
r965 | self.dtype = None | |
self.fileSizeByHeader = None | |||
self.filenameList = [] | |||
self.fileSelector = 0 | |||
|
r1075 | self.Off2StartNxtRec = 0 | |
self.RecCounter = 0 | |||
|
r965 | self.flagNoMoreFiles = 0 | |
|
r1075 | self.data_spc = None | |
self.data_cspc = None | |||
self.data_output = None | |||
|
r965 | self.path = None | |
|
r1075 | self.OffsetStartHeader = 0 | |
self.Off2StartData = 0 | |||
|
r965 | self.ipp = 0 | |
|
r1075 | self.nFDTdataRecors = 0 | |
|
r965 | self.blocksize = 0 | |
self.dataOut = Spectra() | |||
|
r1075 | self.profileIndex = 1 # Always | |
self.dataOut.flagNoData = False | |||
|
r1097 | self.dataOut.nRdPairs = 0 | |
self.dataOut.data_spc = None | |||
|
r1075 | self.dataOut.velocityX = [] | |
self.dataOut.velocityY = [] | |||
self.dataOut.velocityV = [] | |||
|
r965 | ||
def Files2Read(self, fp): | |||
''' | |||
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. | |||
''' | |||
|
r1075 | # self.__checkPath() | |
# Gets the list of files within the fp address | |||
ListaData = os.listdir(fp) | |||
# Sort the list of files from least to largest by names | |||
ListaData = sorted(ListaData) | |||
nFiles = 0 # File Counter | |||
FileList = [] # A list is created that will contain the .fdt files | |||
for IndexFile in ListaData: | |||
if '.fdt' in IndexFile: | |||
|
r965 | FileList.append(IndexFile) | |
|
r1075 | nFiles += 1 | |
# print 'Files2Read' | |||
# print 'Existen '+str(nFiles)+' archivos .fdt' | |||
self.filenameList = FileList # List of files from least to largest by names | |||
|
r965 | 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 | |||
|
r1075 | ||
|
r965 | self.getData() | |
|
r1075 | # print 'running' | |
|
r965 | def setup(self, path=None, | |
|
r1075 | startDate=None, | |
endDate=None, | |||
startTime=None, | |||
endTime=None, | |||
walk=True, | |||
timezone='utc', | |||
code=None, | |||
online=False, | |||
ReadMode=None, | |||
**kwargs): | |||
|
r965 | self.isConfig = True | |
|
r1075 | ||
self.path = path | |||
self.startDate = startDate | |||
self.endDate = endDate | |||
self.startTime = startTime | |||
self.endTime = endTime | |||
self.walk = walk | |||
self.ReadMode = int(ReadMode) | |||
|
r965 | pass | |
|
r1075 | ||
|
r965 | 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. | |||
|
r1075 | ||
|
r965 | You should call the file "read". | |
|
r1075 | ||
|
r965 | ''' | |
|
r1075 | ||
|
r965 | if self.flagNoMoreFiles: | |
self.dataOut.flagNoData = True | |||
|
r1167 | print('NoData se vuelve true') | |
|
r965 | return 0 | |
|
r1075 | self.fp = self.path | |
|
r965 | self.Files2Read(self.fp) | |
self.readFile(self.fp) | |||
self.dataOut.data_spc = self.data_spc | |||
|
r1075 | self.dataOut.data_cspc = self.data_cspc | |
self.dataOut.data_output = self.data_output | |||
|
r1167 | print('self.dataOut.data_output', shape(self.dataOut.data_output)) | |
|
r1075 | ||
# self.removeDC() | |||
return self.dataOut.data_spc | |||
def readFile(self, fp): | |||
|
r965 | ''' | |
You must indicate if you are reading in Online or Offline mode and load the | |||
The parameters for this file reading mode. | |||
|
r1075 | ||
|
r965 | Then you must do 2 actions: | |
|
r1075 | ||
|
r965 | 1. Get the BLTR FileHeader. | |
2. Start reading the first block. | |||
''' | |||
|
r1075 | ||
# The address of the folder is generated the name of the .fdt file that will be read | |||
|
r1167 | print("File: ", self.fileSelector + 1) | |
|
r1075 | ||
|
r965 | if self.fileSelector < len(self.filenameList): | |
|
r1075 | ||
self.fpFile = str(fp) + '/' + \ | |||
str(self.filenameList[self.fileSelector]) | |||
# print self.fpFile | |||
|
r965 | fheader = FileHeaderBLTR() | |
|
r1075 | fheader.FHread(self.fpFile) # Bltr FileHeader Reading | |
self.nFDTdataRecors = fheader.nFDTdataRecors | |||
self.readBlock() # Block reading | |||
|
r965 | else: | |
|
r1167 | print('readFile FlagNoData becomes true') | |
|
r1075 | self.flagNoMoreFiles = True | |
|
r965 | self.dataOut.flagNoData = True | |
|
r1075 | return 0 | |
|
r965 | def getVelRange(self, extrapoints=0): | |
|
r1075 | Lambda = SPEED_OF_LIGHT / 50000000 | |
# 1./(self.dataOut.ippSeconds * self.dataOut.nCohInt) | |||
PRF = self.dataOut.PRF | |||
Vmax = -Lambda / (4. * (1. / PRF) * self.dataOut.nCohInt * 2.) | |||
deltafreq = PRF / (self.nProfiles) | |||
deltavel = (Vmax * 2) / (self.nProfiles) | |||
freqrange = deltafreq * \ | |||
(numpy.arange(self.nProfiles) - self.nProfiles / 2.) - deltafreq / 2 | |||
velrange = deltavel * \ | |||
(numpy.arange(self.nProfiles) - self.nProfiles / 2.) | |||
return velrange | |||
def readBlock(self): | |||
|
r965 | ''' | |
It should be checked if the block has data, if it is not passed to the next file. | |||
|
r1075 | ||
|
r965 | Then the following is done: | |
|
r1075 | ||
|
r965 | 1. Read the RecordHeader | |
2. Fill the buffer with the current block number. | |||
|
r1075 | ||
|
r965 | ''' | |
|
r1075 | ||
if self.BlockCounter < self.nFDTdataRecors - 2: | |||
|
r1167 | print(self.nFDTdataRecors, 'CONDICION!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') | |
|
r1075 | if self.ReadMode == 1: | |
rheader = RecordHeaderBLTR(RecCounter=self.BlockCounter + 1) | |||
elif self.ReadMode == 0: | |||
|
r965 | rheader = RecordHeaderBLTR(RecCounter=self.BlockCounter) | |
|
r1075 | ||
rheader.RHread(self.fpFile) # Bltr FileHeader Reading | |||
self.OffsetStartHeader = rheader.OffsetStartHeader | |||
self.RecCounter = rheader.RecCounter | |||
self.Off2StartNxtRec = rheader.Off2StartNxtRec | |||
self.Off2StartData = rheader.Off2StartData | |||
self.nProfiles = rheader.nProfiles | |||
self.nChannels = rheader.nChannels | |||
self.nHeights = rheader.nHeights | |||
self.frequency = rheader.TransmitFrec | |||
self.DualModeIndex = rheader.DualModeIndex | |||
self.pairsList = [(0, 1), (0, 2), (1, 2)] | |||
|
r965 | self.dataOut.pairsList = self.pairsList | |
|
r1075 | ||
self.nRdPairs = len(self.dataOut.pairsList) | |||
|
r965 | self.dataOut.nRdPairs = self.nRdPairs | |
|
r1075 | ||
self.__firstHeigth = rheader.StartRangeSamp | |||
self.__deltaHeigth = rheader.SampResolution | |||
self.dataOut.heightList = self.__firstHeigth + \ | |||
|
r1167 | numpy.array(list(range(self.nHeights))) * self.__deltaHeigth | |
self.dataOut.channelList = list(range(self.nChannels)) | |||
|
r1075 | self.dataOut.nProfiles = rheader.nProfiles | |
self.dataOut.nIncohInt = rheader.nIncohInt | |||
self.dataOut.nCohInt = rheader.nCohInt | |||
self.dataOut.ippSeconds = 1 / float(rheader.PRFhz) | |||
self.dataOut.PRF = rheader.PRFhz | |||
self.dataOut.nFFTPoints = rheader.nProfiles | |||
self.dataOut.utctime = rheader.nUtime | |||
self.dataOut.timeZone = 0 | |||
self.dataOut.normFactor = self.dataOut.nProfiles * \ | |||
self.dataOut.nIncohInt * self.dataOut.nCohInt | |||
self.dataOut.outputInterval = self.dataOut.ippSeconds * \ | |||
self.dataOut.nCohInt * self.dataOut.nIncohInt * self.nProfiles | |||
self.data_output = numpy.ones([3, rheader.nHeights]) * numpy.NaN | |||
|
r1167 | print('self.data_output', shape(self.data_output)) | |
|
r1075 | self.dataOut.velocityX = [] | |
self.dataOut.velocityY = [] | |||
self.dataOut.velocityV = [] | |||
|
r965 | '''Block Reading, the Block Data is received and Reshape is used to give it | |
shape. | |||
|
r1075 | ''' | |
# Procedure to take the pointer to where the date block starts | |||
startDATA = open(self.fpFile, "rb") | |||
OffDATA = self.OffsetStartHeader + self.RecCounter * \ | |||
self.Off2StartNxtRec + self.Off2StartData | |||
|
r965 | startDATA.seek(OffDATA, os.SEEK_SET) | |
|
r1075 | ||
|
r965 | def moving_average(x, N=2): | |
|
r1075 | return numpy.convolve(x, numpy.ones((N,)) / N)[(N - 1):] | |
def gaus(xSamples, a, x0, sigma): | |||
return a * exp(-(xSamples - x0)**2 / (2 * sigma**2)) | |||
def Find(x, value): | |||
|
r965 | for index in range(len(x)): | |
|
r1075 | if x[index] == value: | |
return index | |||
|
r965 | def pol2cart(rho, phi): | |
x = rho * numpy.cos(phi) | |||
y = rho * numpy.sin(phi) | |||
return(x, y) | |||
|
r1075 | ||
if self.DualModeIndex == self.ReadMode: | |||
self.data_fft = numpy.fromfile( | |||
startDATA, [('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)) | |||
|
r1001 | copy = self.data_block.copy() | |
|
r1075 | spc = copy * numpy.conjugate(copy) | |
self.data_spc = numpy.absolute( | |||
spc) # valor absoluto o magnitud | |||
|
r965 | factor = self.dataOut.normFactor | |
|
r1075 | ||
z = self.data_spc.copy() # /factor | |||
|
r965 | z = numpy.where(numpy.isfinite(z), z, numpy.NAN) | |
#zdB = 10*numpy.log10(z) | |||
|
r1167 | print(' ') | |
print('Z: ') | |||
print(shape(z)) | |||
print(' ') | |||
print(' ') | |||
|
r1075 | ||
self.dataOut.data_spc = self.data_spc | |||
self.noise = self.dataOut.getNoise( | |||
ymin_index=80, ymax_index=132) # /factor | |||
|
r965 | #noisedB = 10*numpy.log10(self.noise) | |
|
r1075 | ||
ySamples = numpy.ones([3, self.nProfiles]) | |||
phase = numpy.ones([3, self.nProfiles]) | |||
CSPCSamples = numpy.ones( | |||
[3, self.nProfiles], dtype=numpy.complex_) | |||
coherence = numpy.ones([3, self.nProfiles]) | |||
PhaseSlope = numpy.ones(3) | |||
PhaseInter = numpy.ones(3) | |||
|
r965 | '''****** Getting CrossSpectra ******''' | |
|
r1075 | cspc = self.data_block.copy() | |
self.data_cspc = self.data_block.copy() | |||
xFrec = self.getVelRange(1) | |||
VelRange = self.getVelRange(1) | |||
self.dataOut.VelRange = VelRange | |||
# print ' ' | |||
# print ' ' | |||
# print 'xFrec',xFrec | |||
# print ' ' | |||
# print ' ' | |||
# Height=35 | |||
for i in range(self.nRdPairs): | |||
|
r965 | chan_index0 = self.dataOut.pairsList[i][0] | |
chan_index1 = self.dataOut.pairsList[i][1] | |||
|
r1075 | ||
self.data_cspc[i, :, :] = cspc[chan_index0, :, | |||
:] * numpy.conjugate(cspc[chan_index1, :, :]) | |||
|
r965 | '''Getting Eij and Nij''' | |
|
r1075 | (AntennaX0, AntennaY0) = pol2cart( | |
rheader.AntennaCoord0, rheader.AntennaAngl0 * numpy.pi / 180) | |||
(AntennaX1, AntennaY1) = pol2cart( | |||
rheader.AntennaCoord1, rheader.AntennaAngl1 * numpy.pi / 180) | |||
(AntennaX2, AntennaY2) = pol2cart( | |||
rheader.AntennaCoord2, 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]]) | |||
|
r965 | self.dataOut.ChanDist = self.ChanDist | |
|
r1075 | ||
|
r965 | # for Height in range(self.nHeights): | |
|
r1075 | # | |
|
r965 | # for i in range(self.nRdPairs): | |
|
r1075 | # | |
|
r965 | # '''****** Line of Data SPC ******''' | |
# zline=z[i,:,Height] | |||
|
r1075 | # | |
|
r965 | # '''****** DC is removed ******''' | |
# DC=Find(zline,numpy.amax(zline)) | |||
# zline[DC]=(zline[DC-1]+zline[DC+1])/2 | |||
|
r1075 | # | |
# | |||
|
r965 | # '''****** SPC is normalized ******''' | |
# FactNorm= zline.copy() / numpy.sum(zline.copy()) | |||
# FactNorm= FactNorm/numpy.sum(FactNorm) | |||
|
r1075 | # | |
|
r965 | # SmoothSPC=moving_average(FactNorm,N=3) | |
|
r1075 | # | |
|
r965 | # xSamples = ar(range(len(SmoothSPC))) | |
# ySamples[i] = SmoothSPC-self.noise[i] | |||
|
r1075 | # | |
|
r965 | # for i in range(self.nRdPairs): | |
|
r1075 | # | |
|
r965 | # '''****** Line of Data CSPC ******''' | |
# cspcLine=self.data_cspc[i,:,Height].copy() | |||
|
r1075 | # | |
# | |||
# | |||
|
r965 | # '''****** CSPC is normalized ******''' | |
# chan_index0 = self.dataOut.pairsList[i][0] | |||
|
r1075 | # chan_index1 = self.dataOut.pairsList[i][1] | |
|
r965 | # CSPCFactor= numpy.sum(ySamples[chan_index0]) * numpy.sum(ySamples[chan_index1]) | |
|
r1075 | # | |
# | |||
|
r965 | # CSPCNorm= cspcLine.copy() / numpy.sqrt(CSPCFactor) | |
|
r1075 | # | |
# | |||
|
r965 | # CSPCSamples[i] = CSPCNorm-self.noise[i] | |
# coherence[i] = numpy.abs(CSPCSamples[i]) / numpy.sqrt(CSPCFactor) | |||
|
r1075 | # | |
|
r965 | # '''****** DC is removed ******''' | |
# DC=Find(coherence[i],numpy.amax(coherence[i])) | |||
# coherence[i][DC]=(coherence[i][DC-1]+coherence[i][DC+1])/2 | |||
# coherence[i]= moving_average(coherence[i],N=2) | |||
|
r1075 | # | |
|
r965 | # phase[i] = moving_average( numpy.arctan2(CSPCSamples[i].imag, CSPCSamples[i].real),N=1)#*180/numpy.pi | |
|
r1075 | # | |
# | |||
|
r965 | # '''****** Getting fij width ******''' | |
|
r1075 | # | |
|
r965 | # yMean=[] | |
|
r1075 | # yMean2=[] | |
# | |||
|
r965 | # for j in range(len(ySamples[1])): | |
# yMean=numpy.append(yMean,numpy.average([ySamples[0,j],ySamples[1,j],ySamples[2,j]])) | |||
|
r1075 | # | |
|
r965 | # '''******* Getting fitting Gaussian ******''' | |
# meanGauss=sum(xSamples*yMean) / len(xSamples) | |||
# sigma=sum(yMean*(xSamples-meanGauss)**2) / len(xSamples) | |||
# #print 'Height',Height,'SNR', meanGauss/sigma**2 | |||
|
r1075 | # | |
|
r965 | # if (abs(meanGauss/sigma**2) > 0.0001) : | |
|
r1075 | # | |
# try: | |||
|
r965 | # popt,pcov = curve_fit(gaus,xSamples,yMean,p0=[1,meanGauss,sigma]) | |
|
r1075 | # | |
|
r965 | # if numpy.amax(popt)>numpy.amax(yMean)*0.3: | |
# FitGauss=gaus(xSamples,*popt) | |||
|
r1075 | # | |
# else: | |||
|
r965 | # FitGauss=numpy.ones(len(xSamples))*numpy.mean(yMean) | |
# print 'Verificador: Dentro', Height | |||
# except RuntimeError: | |||
|
r1075 | # | |
|
r965 | # try: | |
# for j in range(len(ySamples[1])): | |||
# yMean2=numpy.append(yMean2,numpy.average([ySamples[1,j],ySamples[2,j]])) | |||
# popt,pcov = curve_fit(gaus,xSamples,yMean2,p0=[1,meanGauss,sigma]) | |||
# FitGauss=gaus(xSamples,*popt) | |||
# print 'Verificador: Exepcion1', Height | |||
# except RuntimeError: | |||
|
r1075 | # | |
|
r965 | # try: | |
# popt,pcov = curve_fit(gaus,xSamples,ySamples[1],p0=[1,meanGauss,sigma]) | |||
# FitGauss=gaus(xSamples,*popt) | |||
# print 'Verificador: Exepcion2', Height | |||
# except RuntimeError: | |||
# FitGauss=numpy.ones(len(xSamples))*numpy.mean(yMean) | |||
# print 'Verificador: Exepcion3', Height | |||
# else: | |||
# FitGauss=numpy.ones(len(xSamples))*numpy.mean(yMean) | |||
# #print 'Verificador: Fuera', Height | |||
|
r1075 | # | |
# | |||
# | |||
|
r965 | # Maximun=numpy.amax(yMean) | |
# eMinus1=Maximun*numpy.exp(-1) | |||
|
r1075 | # | |
|
r965 | # HWpos=Find(FitGauss,min(FitGauss, key=lambda value:abs(value-eMinus1))) | |
# HalfWidth= xFrec[HWpos] | |||
# GCpos=Find(FitGauss, numpy.amax(FitGauss)) | |||
# Vpos=Find(FactNorm, numpy.amax(FactNorm)) | |||
# #Vpos=numpy.sum(FactNorm)/len(FactNorm) | |||
# #Vpos=Find(FactNorm, min(FactNorm, key=lambda value:abs(value- numpy.mean(FactNorm) ))) | |||
# #print 'GCpos',GCpos, numpy.amax(FitGauss), 'HWpos',HWpos | |||
# '''****** Getting Fij ******''' | |||
|
r1075 | # | |
|
r965 | # GaussCenter=xFrec[GCpos] | |
# if (GaussCenter<0 and HalfWidth>0) or (GaussCenter>0 and HalfWidth<0): | |||
# Fij=abs(GaussCenter)+abs(HalfWidth)+0.0000001 | |||
# else: | |||
# Fij=abs(GaussCenter-HalfWidth)+0.0000001 | |||
|
r1075 | # | |
|
r965 | # '''****** Getting Frecuency range of significant data ******''' | |
|
r1075 | # | |
|
r965 | # Rangpos=Find(FitGauss,min(FitGauss, key=lambda value:abs(value-Maximun*0.10))) | |
|
r1075 | # | |
|
r965 | # if Rangpos<GCpos: | |
# Range=numpy.array([Rangpos,2*GCpos-Rangpos]) | |||
# else: | |||
# Range=numpy.array([2*GCpos-Rangpos,Rangpos]) | |||
|
r1075 | # | |
|
r965 | # FrecRange=xFrec[Range[0]:Range[1]] | |
|
r1075 | # | |
|
r965 | # #print 'FrecRange', FrecRange | |
# '''****** Getting SCPC Slope ******''' | |||
|
r1075 | # | |
|
r965 | # for i in range(self.nRdPairs): | |
|
r1075 | # | |
|
r965 | # if len(FrecRange)>5 and len(FrecRange)<self.nProfiles*0.5: | |
|
r1075 | # PhaseRange=moving_average(phase[i,Range[0]:Range[1]],N=3) | |
# | |||
|
r965 | # slope, intercept, r_value, p_value, std_err = stats.linregress(FrecRange,PhaseRange) | |
# PhaseSlope[i]=slope | |||
# PhaseInter[i]=intercept | |||
# else: | |||
# PhaseSlope[i]=0 | |||
# PhaseInter[i]=0 | |||
|
r1075 | # | |
|
r965 | # # plt.figure(i+15) | |
# # plt.title('FASE ( CH%s*CH%s )' %(self.dataOut.pairsList[i][0],self.dataOut.pairsList[i][1])) | |||
# # plt.xlabel('Frecuencia (KHz)') | |||
# # plt.ylabel('Magnitud') | |||
# # #plt.subplot(311+i) | |||
# # plt.plot(FrecRange,PhaseRange,'b') | |||
# # plt.plot(FrecRange,FrecRange*PhaseSlope[i]+PhaseInter[i],'r') | |||
|
r1075 | # | |
|
r965 | # #plt.axis([-0.6, 0.2, -3.2, 3.2]) | |
|
r1075 | # | |
# | |||
|
r965 | # '''Getting constant C''' | |
# cC=(Fij*numpy.pi)**2 | |||
|
r1075 | # | |
|
r965 | # # '''Getting Eij and Nij''' | |
# # (AntennaX0,AntennaY0)=pol2cart(rheader.AntennaCoord0, rheader.AntennaAngl0*numpy.pi/180) | |||
# # (AntennaX1,AntennaY1)=pol2cart(rheader.AntennaCoord1, rheader.AntennaAngl1*numpy.pi/180) | |||
# # (AntennaX2,AntennaY2)=pol2cart(rheader.AntennaCoord2, rheader.AntennaAngl2*numpy.pi/180) | |||
|
r1075 | # # | |
|
r965 | # # E01=AntennaX0-AntennaX1 | |
# # N01=AntennaY0-AntennaY1 | |||
|
r1075 | # # | |
|
r965 | # # E02=AntennaX0-AntennaX2 | |
# # N02=AntennaY0-AntennaY2 | |||
|
r1075 | # # | |
|
r965 | # # E12=AntennaX1-AntennaX2 | |
# # N12=AntennaY1-AntennaY2 | |||
|
r1075 | # | |
|
r965 | # '''****** Getting constants F and G ******''' | |
# MijEijNij=numpy.array([[E02,N02], [E12,N12]]) | |||
# MijResult0=(-PhaseSlope[1]*cC) / (2*numpy.pi) | |||
|
r1075 | # MijResult1=(-PhaseSlope[2]*cC) / (2*numpy.pi) | |
|
r965 | # MijResults=numpy.array([MijResult0,MijResult1]) | |
# (cF,cG) = numpy.linalg.solve(MijEijNij, MijResults) | |||
|
r1075 | # | |
|
r965 | # '''****** Getting constants A, B and H ******''' | |
# W01=numpy.amax(coherence[0]) | |||
# W02=numpy.amax(coherence[1]) | |||
# W12=numpy.amax(coherence[2]) | |||
|
r1075 | # | |
|
r965 | # WijResult0=((cF*E01+cG*N01)**2)/cC - numpy.log(W01 / numpy.sqrt(numpy.pi/cC)) | |
# WijResult1=((cF*E02+cG*N02)**2)/cC - numpy.log(W02 / numpy.sqrt(numpy.pi/cC)) | |||
# WijResult2=((cF*E12+cG*N12)**2)/cC - numpy.log(W12 / numpy.sqrt(numpy.pi/cC)) | |||
|
r1075 | # | |
|
r965 | # WijResults=numpy.array([WijResult0, WijResult1, WijResult2]) | |
|
r1075 | # | |
# WijEijNij=numpy.array([ [E01**2, N01**2, 2*E01*N01] , [E02**2, N02**2, 2*E02*N02] , [E12**2, N12**2, 2*E12*N12] ]) | |||
|
r965 | # (cA,cB,cH) = numpy.linalg.solve(WijEijNij, WijResults) | |
|
r1075 | # | |
|
r965 | # VxVy=numpy.array([[cA,cH],[cH,cB]]) | |
|
r1075 | # | |
|
r965 | # VxVyResults=numpy.array([-cF,-cG]) | |
# (Vx,Vy) = numpy.linalg.solve(VxVy, VxVyResults) | |||
# Vzon = Vy | |||
# Vmer = Vx | |||
# Vmag=numpy.sqrt(Vzon**2+Vmer**2) | |||
# Vang=numpy.arctan2(Vmer,Vzon) | |||
|
r1075 | # | |
|
r965 | # if abs(Vy)<100 and abs(Vy)> 0.: | |
# self.dataOut.velocityX=numpy.append(self.dataOut.velocityX, Vzon) #Vmag | |||
# #print 'Vmag',Vmag | |||
# else: | |||
# self.dataOut.velocityX=numpy.append(self.dataOut.velocityX, NaN) | |||
|
r1075 | # | |
|
r965 | # if abs(Vx)<100 and abs(Vx) > 0.: | |
# self.dataOut.velocityY=numpy.append(self.dataOut.velocityY, Vmer) #Vang | |||
|
r1075 | # #print 'Vang',Vang | |
|
r965 | # else: | |
# self.dataOut.velocityY=numpy.append(self.dataOut.velocityY, NaN) | |||
|
r1075 | # | |
|
r965 | # if abs(GaussCenter)<2: | |
# self.dataOut.velocityV=numpy.append(self.dataOut.velocityV, xFrec[Vpos]) | |||
|
r1075 | # | |
|
r965 | # else: | |
# self.dataOut.velocityV=numpy.append(self.dataOut.velocityV, NaN) | |||
|
r1075 | # | |
# | |||
|
r965 | # # print '********************************************' | |
# # print 'HalfWidth ', HalfWidth | |||
# # print 'Maximun ', Maximun | |||
# # print 'eMinus1 ', eMinus1 | |||
# # print 'Rangpos ', Rangpos | |||
# # print 'GaussCenter ',GaussCenter | |||
# # print 'E01 ',E01 | |||
# # print 'N01 ',N01 | |||
# # print 'E02 ',E02 | |||
# # print 'N02 ',N02 | |||
# # print 'E12 ',E12 | |||
# # print 'N12 ',N12 | |||
# #print 'self.dataOut.velocityX ', self.dataOut.velocityX | |||
# # print 'Fij ', Fij | |||
# # print 'cC ', cC | |||
# # print 'cF ', cF | |||
# # print 'cG ', cG | |||
# # print 'cA ', cA | |||
# # print 'cB ', cB | |||
# # print 'cH ', cH | |||
# # print 'Vx ', Vx | |||
# # print 'Vy ', Vy | |||
# # print 'Vmag ', Vmag | |||
# # print 'Vang ', Vang*180/numpy.pi | |||
# # print 'PhaseSlope ',PhaseSlope[0] | |||
# # print 'PhaseSlope ',PhaseSlope[1] | |||
# # print 'PhaseSlope ',PhaseSlope[2] | |||
# # print '********************************************' | |||
# #print 'data_output',shape(self.dataOut.velocityX), shape(self.dataOut.velocityY) | |||
|
r1075 | # | |
|
r965 | # #print 'self.dataOut.velocityX', len(self.dataOut.velocityX) | |
# #print 'self.dataOut.velocityY', len(self.dataOut.velocityY) | |||
# #print 'self.dataOut.velocityV', self.dataOut.velocityV | |||
|
r1075 | # | |
|
r965 | # self.data_output[0]=numpy.array(self.dataOut.velocityX) | |
# self.data_output[1]=numpy.array(self.dataOut.velocityY) | |||
# self.data_output[2]=numpy.array(self.dataOut.velocityV) | |||
|
r1075 | # | |
|
r965 | # prin= self.data_output[0][~numpy.isnan(self.data_output[0])] | |
# print ' ' | |||
|
r1075 | # print 'VmagAverage',numpy.mean(prin) | |
|
r965 | # print ' ' | |
# # plt.figure(5) | |||
# # plt.subplot(211) | |||
# # plt.plot(self.dataOut.velocityX,'yo:') | |||
# # plt.subplot(212) | |||
# # plt.plot(self.dataOut.velocityY,'yo:') | |||
|
r1075 | # | |
|
r965 | # # plt.figure(1) | |
# # # plt.subplot(121) | |||
# # # plt.plot(xFrec,ySamples[0],'k',label='Ch0') | |||
# # # plt.plot(xFrec,ySamples[1],'g',label='Ch1') | |||
# # # plt.plot(xFrec,ySamples[2],'r',label='Ch2') | |||
# # # plt.plot(xFrec,FitGauss,'yo:',label='fit') | |||
# # # plt.legend() | |||
# # plt.title('DATOS A ALTURA DE 2850 METROS') | |||
|
r1075 | # # | |
|
r965 | # # plt.xlabel('Frecuencia (KHz)') | |
# # plt.ylabel('Magnitud') | |||
# # # plt.subplot(122) | |||
# # # plt.title('Fit for Time Constant') | |||
# # #plt.plot(xFrec,zline) | |||
# # #plt.plot(xFrec,SmoothSPC,'g') | |||
# # plt.plot(xFrec,FactNorm) | |||
# # plt.axis([-4, 4, 0, 0.15]) | |||
# # # plt.xlabel('SelfSpectra KHz') | |||
|
r1075 | # # | |
|
r965 | # # plt.figure(10) | |
# # # plt.subplot(121) | |||
# # plt.plot(xFrec,ySamples[0],'b',label='Ch0') | |||
# # plt.plot(xFrec,ySamples[1],'y',label='Ch1') | |||
# # plt.plot(xFrec,ySamples[2],'r',label='Ch2') | |||
# # # plt.plot(xFrec,FitGauss,'yo:',label='fit') | |||
# # plt.legend() | |||
# # plt.title('SELFSPECTRA EN CANALES') | |||
|
r1075 | # # | |
|
r965 | # # plt.xlabel('Frecuencia (KHz)') | |
# # plt.ylabel('Magnitud') | |||
# # # plt.subplot(122) | |||
# # # plt.title('Fit for Time Constant') | |||
# # #plt.plot(xFrec,zline) | |||
# # #plt.plot(xFrec,SmoothSPC,'g') | |||
# # # plt.plot(xFrec,FactNorm) | |||
# # # plt.axis([-4, 4, 0, 0.15]) | |||
# # # plt.xlabel('SelfSpectra KHz') | |||
|
r1075 | # # | |
|
r965 | # # plt.figure(9) | |
|
r1075 | # # | |
# # | |||
|
r965 | # # plt.title('DATOS SUAVIZADOS') | |
# # plt.xlabel('Frecuencia (KHz)') | |||
# # plt.ylabel('Magnitud') | |||
# # plt.plot(xFrec,SmoothSPC,'g') | |||
|
r1075 | # # | |
|
r965 | # # #plt.plot(xFrec,FactNorm) | |
# # plt.axis([-4, 4, 0, 0.15]) | |||
# # # plt.xlabel('SelfSpectra KHz') | |||
|
r1075 | # # # | |
|
r965 | # # plt.figure(2) | |
# # # #plt.subplot(121) | |||
# # plt.plot(xFrec,yMean,'r',label='Mean SelfSpectra') | |||
# # plt.plot(xFrec,FitGauss,'yo:',label='Ajuste Gaussiano') | |||
# # # plt.plot(xFrec[Rangpos],FitGauss[Find(FitGauss,min(FitGauss, key=lambda value:abs(value-Maximun*0.1)))],'bo') | |||
# # # #plt.plot(xFrec,phase) | |||
# # # plt.xlabel('Suavizado, promediado KHz') | |||
# # plt.title('SELFSPECTRA PROMEDIADO') | |||
# # # #plt.subplot(122) | |||
# # # #plt.plot(xSamples,zline) | |||
# # plt.xlabel('Frecuencia (KHz)') | |||
# # plt.ylabel('Magnitud') | |||
# # plt.legend() | |||
|
r1075 | # # # | |
|
r965 | # # # plt.figure(3) | |
# # # plt.subplot(311) | |||
# # # #plt.plot(xFrec,phase[0]) | |||
# # # plt.plot(xFrec,phase[0],'g') | |||
# # # plt.subplot(312) | |||
# # # plt.plot(xFrec,phase[1],'g') | |||
# # # plt.subplot(313) | |||
# # # plt.plot(xFrec,phase[2],'g') | |||
# # # #plt.plot(xFrec,phase[2]) | |||
|
r1075 | # # # | |
|
r965 | # # # plt.figure(4) | |
|
r1075 | # # # | |
|
r965 | # # # plt.plot(xSamples,coherence[0],'b') | |
# # # plt.plot(xSamples,coherence[1],'r') | |||
# # # plt.plot(xSamples,coherence[2],'g') | |||
# # plt.show() | |||
|
r1075 | # # # | |
|
r965 | # # # plt.clf() | |
# # # plt.cla() | |||
|
r1075 | # # # plt.close() | |
# | |||
# print ' ' | |||
self.BlockCounter += 2 | |||
|
r965 | else: | |
|
r1075 | self.fileSelector += 1 | |
self.BlockCounter = 0 | |||
|
r1167 | print("Next File") |