##// END OF EJS Templates
schain
RSS
merged branches
merged branches
jespinoza - - Load All Authors

File last commit:

r1361:1dcf28d5b762 merge
r1370:81f892b894eb merge
Show More
bltrproc_parameters.py
400 lines | 15.6 KiB | text/x-python | PythonLexer
/ schainpy / model / proc / bltrproc_parameters.py
History | Source | Raw |Copy content |Copy permalink
ebocanegra
first commit
 r965 '''
Created on Oct 24, 2016
@author: roj- LouVD
'''
import numpy
import datetime
import time
jespinoza
Fix h5py Dataset value attribute deprecation
 r1360 from schainpy.model.proc.jroproc_base import ProcessingUnit, Operation
Juan C. Espinoza
BLTR ok
 r1018 from schainpy.model.data.jrodata import Parameters
ebocanegra
first commit
 r965
Juan C. Espinoza
Rewrite controller, remove MPDecorator to units (keep for plots an writers) use of queues for interproc comm instead of zmq, self operations are no longer supported
 r1287
Juan C. Espinoza
BLTRParamreader ready
 r1010 class BLTRParametersProc(ProcessingUnit):
ebocanegra
first commit
 r965 '''
Juan C. Espinoza
BLTRParamreader ready
 r1010 Processing unit for BLTR parameters data (winds)
ebocanegra
first commit
 r965 Inputs:
self.dataOut.nmodes - Number of operation modes
self.dataOut.nchannels - Number of channels
self.dataOut.nranges - Number of ranges
Juan C. Espinoza
BLTRParamreader ready
 r1010
jespinoza
Add update method to plots to pass data (no more changes in jrodata)
 r1343 self.dataOut.data_snr - SNR array
ebocanegra
first commit
 r965 self.dataOut.data_output - Zonal, Vertical and Meridional velocity array
self.dataOut.height - Height array (km)
self.dataOut.time - Time array (seconds)
Juan C. Espinoza
BLTRParamreader ready
 r1010
ebocanegra
first commit
 r965 self.dataOut.fileIndex -Index of the file currently read
self.dataOut.lat - Latitude coordinate of BLTR location
Juan C. Espinoza
BLTRParamreader ready
 r1010
ebocanegra
first commit
 r965 self.dataOut.doy - Experiment doy (number of the day in the current year)
self.dataOut.month - Experiment month
self.dataOut.day - Experiment day
self.dataOut.year - Experiment year
'''
Juan C. Espinoza
BLTRParamreader ready
 r1010
Juan C. Espinoza
Multiprocessing for BLTR (sswma) data
 r1181 def __init__(self):
ebocanegra
first commit
 r965 '''
Juan C. Espinoza
BLTRParamreader ready
 r1010 Inputs: None
ebocanegra
first commit
 r965 '''
Juan C. Espinoza
Multiprocessing for BLTR (sswma) data
 r1181 ProcessingUnit.__init__(self)
ebocanegra
first commit
 r965 self.dataOut = Parameters()
Juan C. Espinoza
Operation MAD2Writer done Task #343
 r1021 def setup(self, mode):
'''
Juan C. Espinoza
BLTRParamreader ready
 r1010 '''
Juan C. Espinoza
Operation MAD2Writer done Task #343
 r1021 self.dataOut.mode = mode
Juan C. Espinoza
BLTR ok
 r1018
Juan C. Espinoza
Operation MAD2Writer done Task #343
 r1021 def run(self, mode, snr_threshold=None):
'''
Juan C. Espinoza
BLTR ok
 r1018 Inputs:
mode = High resolution (0) or Low resolution (1) data
snr_threshold = snr filter value
Juan C. Espinoza
BLTRParamreader ready
 r1010 '''
Juan C. Espinoza
Operation MAD2Writer done Task #343
 r1021
if not self.isConfig:
self.setup(mode)
self.isConfig = True
Juan C. Espinoza
BLTR ok
 r1018 if self.dataIn.type == 'Parameters':
ebocanegra
first commit
 r965 self.dataOut.copy(self.dataIn)
Juan C. Espinoza
Fix utc times and bugs in BLTR modules
 r1085
self.dataOut.data_param = self.dataOut.data[mode]
Juan C. Espinoza
BLTR ok
 r1018 self.dataOut.heightList = self.dataOut.height[0]
jespinoza
Add update method to plots to pass data (no more changes in jrodata)
 r1343 self.dataOut.data_snr = self.dataOut.data_snr[mode]
jespinoza
Fix h5py Dataset value attribute deprecation
 r1360 SNRavg = numpy.average(self.dataOut.data_snr, axis=0)
SNRavgdB = 10*numpy.log10(SNRavg)
self.dataOut.data_snr_avg_db = SNRavgdB.reshape(1, *SNRavgdB.shape)
ebocanegra
first commit
 r965
Danny Scipión
Actualizacion de SNR como parte de data_param
 r1357 # Censoring Data
Juan C. Espinoza
BLTR ok
 r1018 if snr_threshold is not None:
for i in range(3):
Juan C. Espinoza
Fix utc times and bugs in BLTR modules
 r1085 self.dataOut.data_param[i][SNRavgdB <= snr_threshold] = numpy.nan
Juan C. Espinoza
BLTR ok
 r1018
# TODO
Juan C. Espinoza
Rewrite controller, remove MPDecorator to units (keep for plots an writers) use of queues for interproc comm instead of zmq, self operations are no longer supported
 r1287
Juan C. Espinoza
BLTR ok
 r1018 class OutliersFilter(Operation):
George Yong
Multiprocessing for BLTR (all operations) working
 r1185 def __init__(self):
ebocanegra
first commit
 r965 '''
'''
George Yong
Multiprocessing for BLTR (all operations) working
 r1185 Operation.__init__(self)
ebocanegra
first commit
 r965
Juan C. Espinoza
BLTR ok
 r1018 def run(self, svalue2, method, factor, filter, npoints=9):
ebocanegra
first commit
 r965 '''
Inputs:
svalue - string to select array velocity
svalue2 - string to choose axis filtering
method - 0 for SMOOTH or 1 for MEDIAN
Juan C. Espinoza
BLTR ok
 r1018 factor - number used to set threshold
ebocanegra
first commit
 r965 filter - 1 for data filtering using the standard deviation criteria else 0
npoints - number of points for mask filter
Juan C. Espinoza
BLTR ok
 r1018 '''
George Yong
Python 2to3, Spectra (all operations) working
 r1167 print(' Outliers Filter {} {} / threshold = {}'.format(svalue, svalue, factor))
Juan C. Espinoza
BLTR ok
 r1018
ebocanegra
first commit
 r965
Juan C. Espinoza
BLTR ok
 r1018 yaxis = self.dataOut.heightList
xaxis = numpy.array([[self.dataOut.utctime]])
ebocanegra
first commit
 r965
Juan C. Espinoza
BLTR ok
 r1018 # Zonal
value_temp = self.dataOut.data_output[0]
ebocanegra
first commit
 r965
Juan C. Espinoza
BLTR ok
 r1018 # Zonal
value_temp = self.dataOut.data_output[1]
ebocanegra
first commit
 r965
Juan C. Espinoza
BLTR ok
 r1018 # Vertical
value_temp = numpy.transpose(self.dataOut.data_output[2])
htemp = yaxis
ebocanegra
first commit
 r965 std = value_temp
for h in range(len(htemp)):
nvalues_valid = len(numpy.where(numpy.isfinite(value_temp[h]))[0])
Juan C. Espinoza
BLTR ok
 r1018 minvalid = npoints
ebocanegra
first commit
 r965
#only if valid values greater than the minimum required (10%)
if nvalues_valid > minvalid:
if method == 0:
#SMOOTH
w = value_temp[h] - self.Smooth(input=value_temp[h], width=npoints, edge_truncate=1)
if method == 1:
#MEDIAN
w = value_temp[h] - self.Median(input=value_temp[h], width = npoints)
dw = numpy.std(w[numpy.where(numpy.isfinite(w))],ddof = 1)
threshold = dw*factor
value_temp[numpy.where(w > threshold),h] = numpy.nan
value_temp[numpy.where(w < -1*threshold),h] = numpy.nan
#At the end
if svalue2 == 'inHeight':
value_temp = numpy.transpose(value_temp)
output_array[:,m] = value_temp
if svalue == 'zonal':
self.dataOut.data_output[0] = output_array
elif svalue == 'meridional':
self.dataOut.data_output[1] = output_array
elif svalue == 'vertical':
self.dataOut.data_output[2] = output_array
return self.dataOut.data_output
def Median(self,input,width):
'''
Inputs:
input - Velocity array
width - Number of points for mask filter
'''
if numpy.mod(width,2) == 1:
pc = int((width - 1) / 2)
cont = 0
output = []
for i in range(len(input)):
if i >= pc and i < len(input) - pc:
new2 = input[i-pc:i+pc+1]
temp = numpy.where(numpy.isfinite(new2))
new = new2[temp]
value = numpy.median(new)
output.append(value)
output = numpy.array(output)
output = numpy.hstack((input[0:pc],output))
output = numpy.hstack((output,input[-pc:len(input)]))
return output
def Smooth(self,input,width,edge_truncate = None):
'''
Inputs:
input - Velocity array
width - Number of points for mask filter
edge_truncate - 1 for truncate the convolution product else
'''
if numpy.mod(width,2) == 0:
real_width = width + 1
nzeros = width / 2
else:
real_width = width
nzeros = (width - 1) / 2
half_width = int(real_width)/2
length = len(input)
gate = numpy.ones(real_width,dtype='float')
norm_of_gate = numpy.sum(gate)
nan_process = 0
nan_id = numpy.where(numpy.isnan(input))
if len(nan_id[0]) > 0:
nan_process = 1
pb = numpy.zeros(len(input))
pb[nan_id] = 1.
input[nan_id] = 0.
if edge_truncate == True:
output = numpy.convolve(input/norm_of_gate,gate,mode='same')
elif edge_truncate == False or edge_truncate == None:
output = numpy.convolve(input/norm_of_gate,gate,mode='valid')
output = numpy.hstack((input[0:half_width],output))
output = numpy.hstack((output,input[len(input)-half_width:len(input)]))
if nan_process:
pb = numpy.convolve(pb/norm_of_gate,gate,mode='valid')
pb = numpy.hstack((numpy.zeros(half_width),pb))
pb = numpy.hstack((pb,numpy.zeros(half_width)))
output[numpy.where(pb > 0.9999)] = numpy.nan
input[nan_id] = numpy.nan
return output
def Average(self,aver=0,nhaver=1):
'''
Inputs:
aver - Indicates the time period over which is averaged or consensus data
nhaver - Indicates the decimation factor in heights
'''
nhpoints = 48
lat_piura = -5.17
lat_huancayo = -12.04
lat_porcuya = -5.8
if '%2.2f'%self.dataOut.lat == '%2.2f'%lat_piura:
hcm = 3.
if self.dataOut.year == 2003 :
if self.dataOut.doy >= 25 and self.dataOut.doy < 64:
nhpoints = 12
elif '%2.2f'%self.dataOut.lat == '%2.2f'%lat_huancayo:
hcm = 3.
if self.dataOut.year == 2003 :
if self.dataOut.doy >= 25 and self.dataOut.doy < 64:
nhpoints = 12
elif '%2.2f'%self.dataOut.lat == '%2.2f'%lat_porcuya:
hcm = 5.#2
pdata = 0.2
taver = [1,2,3,4,6,8,12,24]
t0 = 0
tf = 24
ntime =(tf-t0)/taver[aver]
ti = numpy.arange(ntime)
tf = numpy.arange(ntime) + taver[aver]
old_height = self.dataOut.heightList
if nhaver > 1:
num_hei = len(self.dataOut.heightList)/nhaver/self.dataOut.nmodes
deltha = 0.05*nhaver
minhvalid = pdata*nhaver
for im in range(self.dataOut.nmodes):
new_height = numpy.arange(num_hei)*deltha + self.dataOut.height[im,0] + deltha/2.
data_fHeigths_List = []
data_fZonal_List = []
data_fMeridional_List = []
data_fVertical_List = []
startDTList = []
for i in range(ntime):
height = old_height
start = datetime.datetime(self.dataOut.year,self.dataOut.month,self.dataOut.day) + datetime.timedelta(hours = int(ti[i])) - datetime.timedelta(hours = 5)
stop = datetime.datetime(self.dataOut.year,self.dataOut.month,self.dataOut.day) + datetime.timedelta(hours = int(tf[i])) - datetime.timedelta(hours = 5)
limit_sec1 = time.mktime(start.timetuple())
limit_sec2 = time.mktime(stop.timetuple())
t1 = numpy.where(self.f_timesec >= limit_sec1)
t2 = numpy.where(self.f_timesec < limit_sec2)
time_select = []
for val_sec in t1[0]:
if val_sec in t2[0]:
time_select.append(val_sec)
time_select = numpy.array(time_select,dtype = 'int')
minvalid = numpy.ceil(pdata*nhpoints)
zon_aver = numpy.zeros([self.dataOut.nranges,self.dataOut.nmodes],dtype='f4') + numpy.nan
mer_aver = numpy.zeros([self.dataOut.nranges,self.dataOut.nmodes],dtype='f4') + numpy.nan
ver_aver = numpy.zeros([self.dataOut.nranges,self.dataOut.nmodes],dtype='f4') + numpy.nan
if nhaver > 1:
new_zon_aver = numpy.zeros([num_hei,self.dataOut.nmodes],dtype='f4') + numpy.nan
new_mer_aver = numpy.zeros([num_hei,self.dataOut.nmodes],dtype='f4') + numpy.nan
new_ver_aver = numpy.zeros([num_hei,self.dataOut.nmodes],dtype='f4') + numpy.nan
if len(time_select) > minvalid:
time_average = self.f_timesec[time_select]
for im in range(self.dataOut.nmodes):
for ih in range(self.dataOut.nranges):
if numpy.sum(numpy.isfinite(self.f_zon[time_select,ih,im])) >= minvalid:
zon_aver[ih,im] = numpy.nansum(self.f_zon[time_select,ih,im]) / numpy.sum(numpy.isfinite(self.f_zon[time_select,ih,im]))
if numpy.sum(numpy.isfinite(self.f_mer[time_select,ih,im])) >= minvalid:
mer_aver[ih,im] = numpy.nansum(self.f_mer[time_select,ih,im]) / numpy.sum(numpy.isfinite(self.f_mer[time_select,ih,im]))
if numpy.sum(numpy.isfinite(self.f_ver[time_select,ih,im])) >= minvalid:
ver_aver[ih,im] = numpy.nansum(self.f_ver[time_select,ih,im]) / numpy.sum(numpy.isfinite(self.f_ver[time_select,ih,im]))
if nhaver > 1:
for ih in range(num_hei):
hvalid = numpy.arange(nhaver) + nhaver*ih
if numpy.sum(numpy.isfinite(zon_aver[hvalid,im])) >= minvalid:
new_zon_aver[ih,im] = numpy.nansum(zon_aver[hvalid,im]) / numpy.sum(numpy.isfinite(zon_aver[hvalid,im]))
if numpy.sum(numpy.isfinite(mer_aver[hvalid,im])) >= minvalid:
new_mer_aver[ih,im] = numpy.nansum(mer_aver[hvalid,im]) / numpy.sum(numpy.isfinite(mer_aver[hvalid,im]))
if numpy.sum(numpy.isfinite(ver_aver[hvalid,im])) >= minvalid:
new_ver_aver[ih,im] = numpy.nansum(ver_aver[hvalid,im]) / numpy.sum(numpy.isfinite(ver_aver[hvalid,im]))
if nhaver > 1:
zon_aver = new_zon_aver
mer_aver = new_mer_aver
ver_aver = new_ver_aver
height = new_height
tstart = time_average[0]
tend = time_average[-1]
startTime = time.gmtime(tstart)
year = startTime.tm_year
month = startTime.tm_mon
day = startTime.tm_mday
hour = startTime.tm_hour
minute = startTime.tm_min
second = startTime.tm_sec
startDTList.append(datetime.datetime(year,month,day,hour,minute,second))
o_height = numpy.array([])
o_zon_aver = numpy.array([])
o_mer_aver = numpy.array([])
o_ver_aver = numpy.array([])
if self.dataOut.nmodes > 1:
for im in range(self.dataOut.nmodes):
if im == 0:
h_select = numpy.where(numpy.bitwise_and(height[0,:] >=0,height[0,:] <= hcm,numpy.isfinite(height[0,:])))
else:
h_select = numpy.where(numpy.bitwise_and(height[1,:] > hcm,height[1,:] < 20,numpy.isfinite(height[1,:])))
ht = h_select[0]
o_height = numpy.hstack((o_height,height[im,ht]))
o_zon_aver = numpy.hstack((o_zon_aver,zon_aver[ht,im]))
o_mer_aver = numpy.hstack((o_mer_aver,mer_aver[ht,im]))
o_ver_aver = numpy.hstack((o_ver_aver,ver_aver[ht,im]))
data_fHeigths_List.append(o_height)
data_fZonal_List.append(o_zon_aver)
data_fMeridional_List.append(o_mer_aver)
data_fVertical_List.append(o_ver_aver)
else:
h_select = numpy.where(numpy.bitwise_and(height[0,:] <= hcm,numpy.isfinite(height[0,:])))
ht = h_select[0]
o_height = numpy.hstack((o_height,height[im,ht]))
o_zon_aver = numpy.hstack((o_zon_aver,zon_aver[ht,im]))
o_mer_aver = numpy.hstack((o_mer_aver,mer_aver[ht,im]))
o_ver_aver = numpy.hstack((o_ver_aver,ver_aver[ht,im]))
data_fHeigths_List.append(o_height)
data_fZonal_List.append(o_zon_aver)
data_fMeridional_List.append(o_mer_aver)
data_fVertical_List.append(o_ver_aver)
return startDTList, data_fHeigths_List, data_fZonal_List, data_fMeridional_List, data_fVertical_List