jroproc_bltr.py
575 lines
| 22.2 KiB
| text/x-python
|
PythonLexer
ebocanegra
|
r965 | ''' | |
| Created on Oct 24, 2016 | |||
| @author: roj- LouVD | |||
| ''' | |||
| import numpy | |||
| import copy | |||
| import datetime | |||
| import time | |||
| from time import gmtime | |||
| from jroproc_base import ProcessingUnit | |||
| from schainpy.model.data.jrodata import Parameters | |||
| from numpy import transpose | |||
| from matplotlib import cm | |||
| import matplotlib.pyplot as plt | |||
| from matplotlib.mlab import griddata | |||
| class BLTRProcess(ProcessingUnit): | |||
| isConfig = False | |||
| ''' | |||
| Processing unit for BLTR rawdata | |||
| Inputs: | |||
| self.dataOut.nmodes - Number of operation modes | |||
| self.dataOut.nchannels - Number of channels | |||
| self.dataOut.nranges - Number of ranges | |||
| self.dataOut.data_SNR - SNR array | |||
| self.dataOut.data_output - Zonal, Vertical and Meridional velocity array | |||
| self.dataOut.height - Height array (km) | |||
| self.dataOut.time - Time array (seconds) | |||
| self.dataOut.fileIndex -Index of the file currently read | |||
| self.dataOut.lat - Latitude coordinate of BLTR location | |||
| 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 | |||
| ''' | |||
| def __init__(self): | |||
| ''' | |||
| Inputs: None | |||
| ''' | |||
| ProcessingUnit.__init__(self) | |||
| self.dataOut = Parameters() | |||
| # Filters | |||
| snr_val = None | |||
| value = None | |||
| svalue2 = None | |||
| method = None | |||
| factor = None | |||
| filter = None | |||
| npoints = None | |||
| status_value = None | |||
| width = None | |||
| self.flagfirstmode = 0 | |||
| def run (self): | |||
| if self.dataIn.type == "Parameters": | |||
| self.dataOut.copy(self.dataIn) | |||
| def TimeSelect(self): | |||
| ''' | |||
| Selecting the time array according to the day of the experiment with a duration of 24 hours | |||
| ''' | |||
| k1 = datetime.datetime(self.dataOut.year, self.dataOut.month, self.dataOut.day) - datetime.timedelta(hours=5) | |||
| k2 = datetime.datetime(self.dataOut.year, self.dataOut.month, self.dataOut.day) + datetime.timedelta(hours=25) - datetime.timedelta(hours=5) | |||
| limit_sec1 = time.mktime(k1.timetuple()) | |||
| limit_sec2 = time.mktime(k2.timetuple()) | |||
| valid_data = 0 | |||
| doy = self.dataOut.doy | |||
| t1 = numpy.where(self.dataOut.time[0, :] >= limit_sec1) | |||
| t2 = numpy.where(self.dataOut.time[0, :] < 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') | |||
| valid_data = valid_data + len(time_select) | |||
| if len(time_select) > 0: | |||
| self.f_timesec = self.dataOut.time[:, time_select] | |||
| snr = self.dataOut.data_SNR[time_select, :, :, :] | |||
| zon = self.dataOut.data_output[0][time_select, :, :] | |||
| mer = self.dataOut.data_output[1][time_select, :, :] | |||
| ver = self.dataOut.data_output[2][time_select, :, :] | |||
| if valid_data > 0: | |||
| self.timesec1 = self.f_timesec[0, :] | |||
| self.f_height = self.dataOut.height | |||
| self.f_zon = zon | |||
| self.f_mer = mer | |||
| self.f_ver = ver | |||
| self.f_snr = snr | |||
| self.f_timedate = [] | |||
| self.f_time = [] | |||
| for valuet in self.timesec1: | |||
| time_t = time.gmtime(valuet) | |||
| year = time_t.tm_year | |||
| month = time_t.tm_mon | |||
| day = time_t.tm_mday | |||
| hour = time_t.tm_hour | |||
| minute = time_t.tm_min | |||
| second = time_t.tm_sec | |||
| f_timedate_0 = datetime.datetime(year, month, day, hour, minute, second) | |||
| self.f_timedate.append(f_timedate_0) | |||
| return self.f_timedate, self.f_timesec, self.f_height, self.f_zon, self.f_mer, self.f_ver, self.f_snr | |||
| else: | |||
| self.f_timesec = None | |||
| self.f_timedate = None | |||
| self.f_height = None | |||
| self.f_zon = None | |||
| self.f_mer = None | |||
| self.f_ver = None | |||
| self.f_snr = None | |||
| print 'Invalid time' | |||
| return self.f_timedate, self.f_height, self.f_zon, self.f_mer, self.f_ver, self.f_snr | |||
| def SnrFilter(self, snr_val,modetofilter): | |||
| ''' | |||
| Inputs: snr_val - Threshold value | |||
| ''' | |||
| if modetofilter!=2 and modetofilter!=1 : | |||
| raise ValueError,'Mode to filter should be "1" or "2". {} is not valid, check "Modetofilter" value.'.format(modetofilter) | |||
| m = modetofilter-1 | |||
| print ' SNR filter [mode {}]: SNR <= {}: data_output = NA'.format(modetofilter,snr_val) | |||
| for k in range(self.dataOut.nchannels): | |||
| for r in range(self.dataOut.nranges): | |||
| if self.dataOut.data_SNR[r,k,m] <= snr_val: | |||
| self.dataOut.data_output[2][r,m] = numpy.nan | |||
| self.dataOut.data_output[1][r,m] = numpy.nan | |||
| self.dataOut.data_output[0][r,m] = numpy.nan | |||
| def OutliersFilter(self,modetofilter,svalue,svalue2,method,factor,filter,npoints): | |||
| ''' | |||
| Inputs: | |||
| svalue - string to select array velocity | |||
| svalue2 - string to choose axis filtering | |||
| method - 0 for SMOOTH or 1 for MEDIAN | |||
| factor - number used to set threshold | |||
| filter - 1 for data filtering using the standard deviation criteria else 0 | |||
| npoints - number of points for mask filter | |||
| ''' | |||
| if modetofilter!=2 and modetofilter!=1 : | |||
| raise ValueError,'Mode to filter should be "1" or "2". {} is not valid, check "Modetofilter" value.'.format(modetofilter) | |||
| m = modetofilter-1 | |||
| print ' Outliers Filter [mode {}]: {} {} / threshold = {}'.format(modetofilter,svalue,svalue,factor) | |||
| npoints = 9 | |||
| novalid = 0.1 | |||
| if svalue == 'zonal': | |||
| value = self.dataOut.data_output[0] | |||
| elif svalue == 'meridional': | |||
| value = self.dataOut.data_output[1] | |||
| elif svalue == 'vertical': | |||
| value = self.dataOut.data_output[2] | |||
| else: | |||
| print 'value is not defined' | |||
| return | |||
| if svalue2 == 'inTime': | |||
| yaxis = self.dataOut.height | |||
| xaxis = numpy.array([[self.dataOut.time1],[self.dataOut.time1]]) | |||
| elif svalue2 == 'inHeight': | |||
| yaxis = numpy.array([[self.dataOut.time1],[self.dataOut.time1]]) | |||
| xaxis = self.dataOut.height | |||
| else: | |||
| print 'svalue2 is required, either inHeight or inTime' | |||
| return | |||
| output_array = value | |||
| value_temp = value[:,m] | |||
| error = numpy.zeros(len(self.dataOut.time[m,:])) | |||
| if svalue2 == 'inHeight': | |||
| value_temp = numpy.transpose(value_temp) | |||
| error = numpy.zeros(len(self.dataOut.height)) | |||
| htemp = yaxis[m,:] | |||
| std = value_temp | |||
| for h in range(len(htemp)): | |||
| if filter: #standard deviation filtering | |||
| std[h] = numpy.std(value_temp[h],ddof = npoints) | |||
| value_temp[numpy.where(std[h] > 5),h] = numpy.nan | |||
| error[numpy.where(std[h] > 5)] = error[numpy.where(std[h] > 5)] + 1 | |||
| nvalues_valid = len(numpy.where(numpy.isfinite(value_temp[h]))[0]) | |||
| minvalid = novalid*len(xaxis[m,:]) | |||
| if minvalid <= npoints: | |||
| minvalid = npoints | |||
| #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 | |||
| def prePlot(self,modeselect=None): | |||
| ''' | |||
| Inputs: | |||
| self.dataOut.data_output - Zonal, Meridional and Vertical velocity array | |||
| self.dataOut.height - height array | |||
| self.dataOut.time - Time array (seconds) | |||
| self.dataOut.data_SNR - SNR array | |||
| ''' | |||
| m = modeselect -1 | |||
| print ' [Plotting mode {}]'.format(modeselect) | |||
| if not (m ==1 or m==0): | |||
| raise IndexError("'Mode' must be egual to : 1 or 2") | |||
| # | |||
| if self.flagfirstmode==0: | |||
| #copy of the data | |||
| self.data_output_copy = self.dataOut.data_output.copy() | |||
| self.data_height_copy = self.dataOut.height.copy() | |||
| self.data_time_copy = self.dataOut.time.copy() | |||
| self.data_SNR_copy = self.dataOut.data_SNR.copy() | |||
| self.flagfirstmode = 1 | |||
| else: | |||
| self.dataOut.data_output = self.data_output_copy | |||
| self.dataOut.height = self.data_height_copy | |||
| self.dataOut.time = self.data_time_copy | |||
| self.dataOut.data_SNR = self.data_SNR_copy | |||
| self.flagfirstmode = 0 | |||
| #select data for mode m | |||
| #self.dataOut.data_output = self.dataOut.data_output[:,:,m] | |||
| self.dataOut.heightList = self.dataOut.height[0,:] | |||
| data_SNR = self.dataOut.data_SNR[:,:,m] | |||
| self.dataOut.data_SNR= transpose(data_SNR) | |||
| if m==1 and self.dataOut.counter_records%2==0: | |||
| print '*********' | |||
| print 'MODO 2' | |||
| #print 'Zonal', self.dataOut.data_output[0] | |||
| #print 'Meridional', self.dataOut.data_output[1] | |||
| #print 'Vertical', self.dataOut.data_output[2] | |||
| print '*********' | |||
| Vx=self.dataOut.data_output[0,:,m] | |||
| Vy=self.dataOut.data_output[1,:,m] | |||
| Vmag=numpy.sqrt(Vx**2+Vy**2) | |||
| Vang=numpy.arctan2(Vy,Vx) | |||
| #print 'Vmag', Vmag | |||
| #print 'Vang', Vang | |||
| self.dataOut.data_output[0,:,m]=Vmag | |||
| self.dataOut.data_output[1,:,m]=Vang | |||
| prin= self.dataOut.data_output[0,:,m][~numpy.isnan(self.dataOut.data_output[0,:,m])] | |||
| print ' ' | |||
| print 'VmagAverage',numpy.mean(prin) | |||
| print ' ' | |||
| self.dataOut.data_output = self.dataOut.data_output[:,:,m] | |||