''' Created on Nov 9, 2016 @author: roj- LouVD ''' import numpy import os.path import sys import time import datetime from sys import path from os.path import dirname from mimify import HeaderFile from numpy import size, asarray from datetime import datetime from schainpy.model.proc.jroproc_base import ProcessingUnit, Operation from schainpy.model.data.jrodata import Parameters from schainpy.model.data.jroheaderIO import RadarControllerHeader, SystemHeader from schainpy.model.graphics.jroplot_parameters import WindProfilerPlot from schainpy.model.io.jroIO_base import * import schainpy #import madrigal #import madrigal.cedar #from madrigal.cedar import MadrigalCatalogRecord import warnings from time import gmtime from math import floor warnings.simplefilter("error") from numpy.lib.nanfunctions import nansum warnings.simplefilter('ignore', FutureWarning) class testBLTRReader(ProcessingUnit): def __init__(self): path = None startDate = None endDate = None startTime = None endTime = None startTime = None endTime = None isConfig = False dataOut = None walk = None ext = 'swwma' fileList = [] fileIndex = -1 timezone = None filename = None timearray = None height = None snr_ref = None zon_ref = None ver_ref = None mer_ref = None nmodes = None nchannels = None nranges = None year = None month = None day = None lat = None lon = None siteFile = None ProcessingUnit.__init__(self) self.dataOut = self.createObjByDefault() self.imode = 0 self.counter_records = 0 self.isConfig = False self.flagNoMoreFiles = 0 self.buffer = None def createObjByDefault(self): dataObj = Parameters() return dataObj def info(self): ''' Experience information ''' self.hoy = datetime.datetime.now() place = 'Jicamarca Radio Observatory' signalchainweb='http://jro-dev.igp.gob.pe:3000/projects/signal-chain/wiki/Manual_de_Desarrollador' print '{} at {}'.format(self.hoy,place) print 'Boundary Layer and Tropospheric Radar (BLTR) script, Wind velocities and SNR from *.sswma files' print '{} \n'.format(signalchainweb) def run(self, path, startDate, endDate, ext, startTime, endTime): if not(self.isConfig): self.setup(path, startDate, endDate, ext) self.isConfig = True self.getData() def setup(self, path=None, startDate=None, endDate=None, ext=None, startTime=datetime.time(0, 0, 0), endTime=datetime.time(23, 59, 59), timezone=0): self.info() self.path = path if self.path == None: raise ValueError, "The path is not valid" if ext == None: ext = self.ext self.searchFiles(self.path, startDate, endDate, ext) self.timezone = timezone self.ext = ext self.fileIndex = -1 if not(self.fileList): raise Warning, "There is no files matching these date in the folder: %s. \n Check 'startDate' and 'endDate' "%(path) if not(self.setNextFile()): print 'not next file' if (startDate!=None) and (endDate!=None): print "No files in range: %s - %s" %(datetime.datetime.combine(startDate,startTime).ctime(), datetime.datetime.combine(endDate,endTime).ctime()) elif startDate != None: print "No files in range: %s" %(datetime.datetime.combine(startDate,startTime).ctime()) else: print "No files" sys.exit(-1) def searchFiles(self, path, startDate, endDate, ext=None): ''' Searching for BLTR rawdata file in path Creating a list of file to proces included in [startDate,endDate] Input: path - Path to find BLTR rawdata files startDate - Select file from this date enDate - Select file until this date ext - Extension of the file to read ''' fullpath = path foldercounter = 0 print 'Searching file in %s ' % (fullpath) fileList0 = glob.glob1(fullpath, "*%s" % ext) fileList0.sort() self.fileList = [] self.dateFileList = [] for thisFile in fileList0: year = thisFile[-14:-10] if not isNumber(year): continue month = thisFile[-10:-8] if not isNumber(month): continue day = thisFile[-8:-6] if not isNumber(day): continue year, month, day = int(year), int(month), int(day) dateFile = datetime.date(year, month, day) if not ((startDate <= dateFile) and (endDate > dateFile)): continue self.fileList.append(thisFile) self.dateFileList.append(dateFile) return 1 def setNextFile(self): idFile = self.fileIndex while (True): idFile += 1 if idFile >= len(self.fileList): print '\nNo more files in the folder' print 'Total number of file(s) read : {}'.format(self.fileIndex + 1) print 'Time of processing : {}'.format(datetime.datetime.now()- self.hoy) self.flagNoMoreFiles = 1 return 0 if self.isConfig: print '------------------------[Next File]---------------------------' filename = os.path.join(self.path, self.fileList[idFile]) self.Open(filename) print '\n[Setting file] (%s) ...' % self.fileList[idFile] break self.flagIsNewFile =0 self.fileIndex = idFile self.filename = filename print 'File:',self.filename return 1 def readDataBlock(self): self.readHeader() self.dataRecords(0) print '[New Record] record: {} /{} // file {}/{}'.format(self.counter_records,self.nrecords,self.fileIndex+1,len(self.fileList)) self.setDataBuffer() self.flagIsNewBlock = 1 if self.counter_records > self.nrecords: self.flagIsNewFile = 1 return 0 return 1 def setDataBuffer(self): ''' Storing data from one block ''' self.t = datetime.datetime(self.year, self.month, self.day) self.doy = time.localtime(time.mktime(self.t.timetuple())).tm_yday self.buffer = numpy.squeeze(numpy.array([[self.one_snr],[self.one_zonal],[self.one_vertical],[self.one_meridional], [self.time],[self.height],[self.fileIndex], [self.year],[self.month],[self.day],[self.t],[self.doy]])) self.dataOut.time1 = self.time1 def Open(self, filename): ''' Opening BLTR rawdata file defined by filename Inputs: filename - Full path name of BLTR rawdata file ''' [dir, name] = os.path.split(filename) strFile = name.split('.') self.siteFile = strFile[0] # 'peru2' ---> Piura - 'peru1' ---> Huancayo or Porcuya self.filename = filename if os.path.isfile(self.filename) == False: print 'File do not exist. Check "filename"' sys.exit(0) self.h_file = numpy.dtype([ ('FMN', '= 0., self.dataOut.height[0, :] <= self.hcm, numpy.isfinite(self.dataOut.height[0, :]))) else: h_select = numpy.where(numpy.bitwise_and(self.dataOut.height[0, :] >= 0., self.dataOut.height[0, :] < 20, numpy.isfinite(self.dataOut.height[0, :]))) ht = h_select[0] self.o_height = self.dataOut.height[self.im, ht] self.o_zon = self.z_zon[ht, self.im] self.o_mer = self.z_mer[ht, self.im] self.o_ver = self.z_ver[ht, self.im] o_snr = self.dataOut.data_SNR[ :, :, self.im] o_snr = o_snr[ht, :] ndiv = numpy.nansum((numpy.isfinite(o_snr)), 1) ndiv = ndiv.astype(float) sel_div = numpy.where(ndiv == 0.) ndiv[sel_div] = numpy.nan if self.nchannels > 1: msnr = numpy.nansum(o_snr, axis=1) else: msnr = o_snr try: self.msnr = 10 * numpy.log10(msnr / ndiv) except ZeroDivisionError: self.msnr = 10 * numpy.log10(msnr /1) print 'Number of division (ndiv) egal to 1 by default. Check SNR' time_t = time.gmtime(self.dataOut.time1) 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 timedate_0 = datetime.datetime(year, month, day, hour, minute, second) # 1d parameters GDLATR = self.lat GDLONR = self.lon GDLAT2 = self.lat GLON2 = self.lon # 2d parameters GDALT = self.o_height SNL = self.msnr VN1P2 = self.o_zon VN2P2 = self.o_mer EL2 = self.o_ver NROW = len(self.o_height) startTime = timedate_0 endTime = startTime self.dataRec = madrigal.cedar.MadrigalDataRecord(self.kinst, self.kindat, startTime.year, startTime.month, startTime.day, startTime.hour, startTime.minute, startTime.second, 0, endTime.year, endTime.month, endTime.day, endTime.hour, endTime.minute, endTime.second, 0, ('gdlatr', 'gdlonr', 'gdlat2', 'glon2'), ('gdalt', 'snl', 'vn1p2', 'vn2p2', 'el2'), NROW, ind2DList=['gdalt']) # Setting 1d values self.dataRec.set1D('gdlatr', GDLATR) self.dataRec.set1D('gdlonr', GDLONR) self.dataRec.set1D('gdlat2', GDLAT2) self.dataRec.set1D('glon2', GLON2) # Setting 2d values for n in range(self.o_height.shape[0]): self.dataRec.set2D('gdalt', n, GDALT[n]) self.dataRec.set2D('snl', n, SNL[n]) self.dataRec.set2D('vn1p2', n, VN1P2[n]) self.dataRec.set2D('vn2p2', n, VN2P2[n]) self.dataRec.set2D('el2', n, EL2[n]) # Appending new data record ''' [MADRIGAL3]There are two ways to write to a MadrigalCedarFile. Either this method (write) is called after all the records have been appended to the MadrigalCedarFile, or dump is called after a certain number of records are appended, and then at the end dump is called a final time if there were any records not yet dumped, followed by addArray. ''' self.cedarObj.append(self.dataRec) print ' [Writing] records {} (mode {}).'.format(self.dataOut.counter_records,self.im+1) self.cedarObj.dump() def setHeader(self): ''' - Creating self.catHeadObj - Adding information catalog - Writing file header ''' self.catHeadObj = madrigal.cedar.CatalogHeaderCreator(self.fullname) kindatDesc, comments, analyst, history, principleInvestigator = self._info_BLTR() self.catHeadObj.createCatalog(principleInvestigator="Jarjar", expPurpose='characterize the atmospheric dynamics in this region where frequently it happens the El Nino', sciRemarks="http://madrigal3.haystack.mit.edu/static/CEDARMadrigalHdf5Format.pdf") self.catHeadObj.createHeader(kindatDesc, analyst, comments, history) self.catHeadObj.write() print '[File created] path: %s' % (self.fullname) def putData(self): if self.dataOut.flagNoData: return 0 if self.dataOut.counter_records == 1: self.setFile() print '[Writing] Setting new hdf5 file for the mode {}'.format(self.im+1) if self.dataOut.counter_records <= self.dataOut.nrecords: self.writeBlock() if self.dataOut.counter_records == self.dataOut.nrecords: self.cedarObj.addArray() self.setHeader() self.flagIsNewFile = 1 def _info_BLTR(self): kindatDesc = '''--This header is for KINDAT = %d''' % self.kindat history = None analyst = '''Jarjar''' principleInvestigator = ''' Jarjar Radio Observatorio de Jicamarca Instituto Geofisico del Peru ''' if self.type == 1: comments = ''' --These data are provided by two Boundary Layer and Tropospheric Radar (BLTR) deployed at two different locations at Peru(GMT-5), one of them at Piura(5.17 S, 80.64W) and another located at Huancayo (12.04 S, 75.32 W). --The purpose of conducting these observations is to measure wind in the differents levels of height, this radar makes measurements the Zonal(U), Meridional(V) and Vertical(W) wind velocities component in northcoast from Peru. And the main purpose of these mensurations is to characterize the atmospheric dynamics in this region where frequently it happens the 'El Nino Phenomenon' --In Kindat = 1600, contains information of wind velocities component since 0 Km to 3 Km. --In Kindat = 1601, contains information of wind velocities component since 0 Km to 10 Km. --The Huancayo-BLTR is a VHF Profiler Radar System is a 3 channel coherent receiver pulsed radar utilising state-of-the-art software and computing techniques to acquire, decode, and translate signals obtained from partial reflection echoes in the troposphere, lower stratosphere and mesosphere. It uses an array of three horizontal spaced and vertically directed receiving antennas. The data is recorded thirty seconds, averaged to one minute mean values of Height, Zonal, Meridional and Vertical wind. --The Huancayo-BLTR was installed in January 2010. This instrument was designed and constructed by Genesis Soft Pty. Ltd. Is constituted by three groups of spaced antennas (distributed) forming an isosceles triangle. Station _______ Geographic Coord ______ Geomagnetic Coord _______________ Latitude _ Longitude __ Latitude _ Longitude Huancayo (HUA) __12.04 S ___ 75.32 W _____ -12.05 ____ 352.85 Piura (PIU) _____ 5.17 S ___ 80.64 W ______ 5.18 ____ 350.93 WIND OBSERVATIONS --To obtain wind the BLTR uses Spaced Antenna technique (e.g., Briggs 1984). The scatter and reflection it still provided by variations in the refractive index as in the Doppler method(Gage and Basley,1978; Balsley and Gage 1982; Larsen and Rottger 1982), but instead of using the Doppler shift to derive the velocity components, the cross-correlation between signals in an array of three horizontally spaced and vertically directed receiving antennas is used. ...................................................................... For more information, consult the following references: - Balsley, B. B., and K. S. Gage., On the use of radars for operational wind profiling, Bull. Amer. Meteor.Soc.,63, 1009-1018, 1982. - Briggs, B. H., The analysis of spaced sensor data by correations techniques, Handbook for MAP, Vol. 13, SCOTEP Secretariat, University of Illinois, Urbana, 166-186, 1984. - Gage, K. S., and B.B. Balsley., Doppler radar probing of the clear atmosphere, Bull. Amer. Meteor.Soc., 59, 1074-1093, 1978. - Larsen, M. F., The Spaced Antenna Technique for Radar Wind Profiling, Journal of Atm. and Ocean. Technology. , Vol.6, 920-937, 1989. - Larsen, M. F., A method for single radar voracity measurements?, Handbook for MAP,SCOSTEP Secretariat, University of the Illinois, Urban, in press, 1989. ...................................................................... ACKNOWLEDGEMENTS: --The Piura and Huancayo BLTR are part of the network of instruments operated by the Jicamarca Radio Observatory. --The Jicamarca Radio Observatory is a facility of the Instituto Geofisico del Peru operated with support from the NSF Cooperative Agreement ATM-0432565 through Cornell University ...................................................................... Further questions and comments should be addressed to: Radio Observatorio de Jicamarca Instituto Geofisico del Peru Lima, Peru Web URL: http://jro.igp.gob.pe ...................................................................... ''' return kindatDesc, comments, analyst, history, principleInvestigator