##// END OF EJS Templates
schain
RSS
Class to filter data by snr threshold. To calculate vertical and zonal velocities, vertical and zonal errors and SNR.
imanay -
r1662:30ff981335bb
parent child
Show More
Side by Side Unified
Context file:
r1662:30ff981335bb Loading diff...:
Show file before | Show file after | Hide comments
@@ -1899,7 +1899,158 class SpectralMoments(Operation):
1899 if talk:
1899 if talk:
1900 print('noise =', noise)
1900 print('noise =', noise)
1901 return noise
1901 return noise
1902
1903 class JULIADriftsEstimation(Operation):
1904
1905 def __init__(self):
1906 Operation.__init__(self)
1907
1908
1909 def newtotal(self, data):
1910 return numpy.nansum(data)
1911
1912 #def data_filter(self, parm, snrth=-19.5, swth=20, wErrth=500):
1913 def data_filter(self, parm, snrth=-20, swth=20, wErrth=500):
1914
1915 Sz0 = parm.shape # Sz0: h,p
1916 drift = parm[:,0]
1917 sw = 2*parm[:,1]
1918 snr = 10*numpy.log10(parm[:,2])
1919 Sz = drift.shape # Sz: h
1920 mask = numpy.ones((Sz[0]))
1921 th=0
1922 valid=numpy.where(numpy.isfinite(snr))
1923 cvalid = len(valid[0])
1924 if cvalid >= 1:
1925 # CΓ‘lculo del ruido promedio de snr para el i-Γ©simo grupo de alturas
1926 nbins = int(numpy.max(snr)-numpy.min(snr))+1 # bin size = 1, similar to IDL
1927 h = numpy.histogram(snr,bins=nbins)
1928 hist = h[0]
1929 values = numpy.round_(h[1])
1930 moda = values[numpy.where(hist == numpy.max(hist))]
1931 indNoise = numpy.where(numpy.abs(snr - numpy.min(moda)) < 3)[0]
1932
1933 noise = snr[indNoise]
1934 noise_mean = numpy.sum(noise)/len(noise)
1935 # CΓ‘lculo de media de snr
1936 med = numpy.median(snr)
1937 # Establece el umbral de snr
1938 if noise_mean > med + 3:
1939 th = med
1940 else:
1941 th = noise_mean + 3
1942 # Establece mΓ‘scara
1943 novalid = numpy.where(snr <= th)[0]
1944 mask[novalid] = numpy.nan
1945 # Elimina datos que no sobrepasen el umbral: PARAMETRO
1946 novalid = numpy.where(snr <= snrth)
1947 cnovalid = len(novalid[0])
1948 if cnovalid > 0:
1949 mask[novalid] = numpy.nan
1950 novalid = numpy.where(numpy.isnan(snr))
1951 cnovalid = len(novalid[0])
1952 if cnovalid > 0:
1953 mask[novalid] = numpy.nan
1954 new_parm = numpy.zeros((Sz0[0],Sz0[1]))
1955 for h in range(Sz0[0]):
1956 for p in range(Sz0[1]):
1957 if numpy.isnan(mask[h]):
1958 new_parm[h,p]=numpy.nan
1959 else:
1960 new_parm[h,p]=parm[h,p]
1961
1962 return new_parm, th
1963
1964 def run(self, dataOut, zenith, zenithCorrection,heights=None, statistics=0, otype=0):
1965
1966 nCh=dataOut.spcpar.shape[0]
1967
1968 nHei=dataOut.spcpar.shape[1]
1969 nParam=dataOut.spcpar.shape[2]
1970 # Solo las alturas de interes
1971 hei=dataOut.heightList
1972 hvalid=numpy.where([hei >= heights[0]][0] & [hei <= heights[1]][0])[0]
1973 nhvalid=len(hvalid)
1974 parm = numpy.zeros((nCh,nhvalid,nParam))
1975 parm = dataOut.spcpar[:,hvalid,:]
1976
1977 # Primer filtrado: Umbral de SNR
1978 #snrth=-19
1979 for i in range(nCh):
1980 #print('snr:',parm[i,:,2])
1981 #dataOut.spcpar[i,hvalid,:] = self.data_filter(parm[i,:,:],snrth)[0]
1982 dataOut.spcpar[i,hvalid,:] = self.data_filter(parm[i,:,:])[0]
1983 #print('dataOut.spcpar[0,:,2]',dataOut.spcpar[0,:,2])
1984 #print('dataOut.spcpar[1,:,2]',dataOut.spcpar[1,:,2])
1985 zenith = numpy.array(zenith)
1986 zenith -= zenithCorrection
1987 zenith *= numpy.pi/180
1988 alpha = zenith[0]
1989 beta = zenith[1]
1990
1991 dopplerCH0 = dataOut.spcpar[0,:,0]
1992 dopplerCH1 = dataOut.spcpar[1,:,0]
1993 swCH0 = dataOut.spcpar[0,:,1]
1994 swCH1 = dataOut.spcpar[1,:,1]
1995 snrCH0 = 10*numpy.log10(dataOut.spcpar[0,:,2])
1996 snrCH1 = 10*numpy.log10(dataOut.spcpar[1,:,2])
1997 noiseCH0 = dataOut.spcpar[0,:,3]
1998 noiseCH1 = dataOut.spcpar[1,:,3]
1999 wErrCH0 = dataOut.spcpar[0,:,5]
2000 wErrCH1 = dataOut.spcpar[1,:,5]
2001
2002 # Vertical and zonal calculation according to geometry
2003 sinB_A = numpy.sin(beta)*numpy.cos(alpha) - numpy.sin(alpha)* numpy.cos(beta)
2004 drift = -(dopplerCH0 * numpy.sin(beta) - dopplerCH1 * numpy.sin(alpha))/ sinB_A
2005 '''
2006 print('drift.shape:',drift.shape)
2007 print('drift min:', numpy.nanmin(drift))
2008 print('drift max:', numpy.nanmax(drift))
2009 '''
1902
2010
2011 '''
2012 print('shape:', dopplerCH0[hvalid].shape)
2013 print('dopplerCH0:', dopplerCH0[hvalid])
2014 print('dopplerCH1:', dopplerCH1[hvalid])
2015 print('drift:', drift[hvalid])
2016 '''
2017 zonal = (dopplerCH0 * numpy.cos(beta) - dopplerCH1 * numpy.cos(alpha))/ sinB_A
2018 '''
2019 print('zonal min:', numpy.nanmin(zonal))
2020 print('zonal max:', numpy.nanmax(zonal))
2021 '''
2022 #print('zonal:', zonal[hvalid])
2023 snr = (snrCH0 + snrCH1)/2
2024 '''
2025 print('snr min:', 10*numpy.log10(numpy.nanmin(snr)))
2026 print('snr max:', 10*numpy.log10(numpy.nanmax(snr)))
2027 '''
2028 noise = (noiseCH0 + noiseCH1)/2
2029 sw = (swCH0 + swCH1)/2
2030 w_w_err= numpy.sqrt(numpy.power(wErrCH0 * numpy.sin(beta)/numpy.abs(sinB_A),2) + numpy.power(wErrCH1 * numpy.sin(alpha)/numpy.abs(sinB_A),2))
2031 w_e_err= numpy.sqrt(numpy.power(wErrCH0 * numpy.cos(beta)/numpy.abs(-1*sinB_A),2) + numpy.power(wErrCH1 * numpy.cos(alpha)/numpy.abs(-1*sinB_A),2))
2032
2033 # for statistics150km
2034 if statistics:
2035 print('Implemented offline.')
2036
2037 if otype == 0:
2038 winds = numpy.vstack((snr, drift, zonal, noise, sw, w_w_err, w_e_err)) # to process statistics drifts
2039 elif otype == 3:
2040 winds = numpy.vstack((snr, drift, zonal)) # to generic plot: 3 RTI's
2041 elif otype == 4:
2042 winds = numpy.vstack((snrCH0, drift, snrCH1, zonal)) # to generic plot: 4 RTI's
2043
2044 snr1 = numpy.vstack((snrCH0, snrCH1))
2045
2046 dataOut.data_output = winds
2047 dataOut.data_snr = snr1
2048
2049 dataOut.utctimeInit = dataOut.utctime
2050 dataOut.outputInterval = dataOut.timeInterval
2051
2052 return dataOut
2053
1903 class SALags(Operation):
2054 class SALags(Operation):
1904 '''
2055 '''
1905 Function GetMoments()
2056 Function GetMoments()
General Comments 0
You need to be logged in to leave comments. Login now