##// END OF EJS Templates
schain
RSS
classes for ESF and 150kM Processing
imanay -
r1765:827e92f42420
parent child
Show More
Side by Side Unified
Context file:
r1765:827e92f42420 Loading diff...:
Show file before | Show file after | Hide comments
This diff has been collapsed as it changes many lines, (800 lines changed) Show them Hide them
@@ -13,6 +13,7 from multiprocessing import Pool, TimeoutError
13 from multiprocessing.pool import ThreadPool
13 from multiprocessing.pool import ThreadPool
14 import time
14 import time
15
15
16 import matplotlib.pyplot as plt
16 from scipy.optimize import fmin_l_bfgs_b #optimize with bounds on state papameters
17 from scipy.optimize import fmin_l_bfgs_b #optimize with bounds on state papameters
17 from .jroproc_base import ProcessingUnit, Operation, MPDecorator
18 from .jroproc_base import ProcessingUnit, Operation, MPDecorator
18 from schainpy.model.data.jrodata import Parameters, hildebrand_sekhon
19 from schainpy.model.data.jrodata import Parameters, hildebrand_sekhon
@@ -128,6 +129,10 class ParametersProc(ProcessingUnit):
128 self.dataOut.data_pre = [self.dataIn.data_spc, self.dataIn.data_cspc]
129 self.dataOut.data_pre = [self.dataIn.data_spc, self.dataIn.data_cspc]
129 self.dataOut.data_spc = self.dataIn.data_spc
130 self.dataOut.data_spc = self.dataIn.data_spc
130 self.dataOut.data_cspc = self.dataIn.data_cspc
131 self.dataOut.data_cspc = self.dataIn.data_cspc
132 # for JULIA processing
133 self.dataOut.data_diffcspc = self.dataIn.data_diffcspc
134 self.dataOut.nDiffIncohInt = self.dataIn.nDiffIncohInt
135 # for JULIA processing
131 self.dataOut.nProfiles = self.dataIn.nProfiles
136 self.dataOut.nProfiles = self.dataIn.nProfiles
132 self.dataOut.nIncohInt = self.dataIn.nIncohInt
137 self.dataOut.nIncohInt = self.dataIn.nIncohInt
133 self.dataOut.nFFTPoints = self.dataIn.nFFTPoints
138 self.dataOut.nFFTPoints = self.dataIn.nFFTPoints
@@ -1421,67 +1426,10 class SpectralMoments(Operation):
1421
1426
1422 '''
1427 '''
1423
1428
1424 def run(self, dataOut, proc_type=0, mode_fit=0, exp='150EEJ'):
1425
1426 absc = dataOut.abscissaList[:-1]
1427 #noise = dataOut.noise
1428 nChannel = dataOut.data_pre[0].shape[0]
1429 nHei = dataOut.data_pre[0].shape[2]
1430 data_param = numpy.zeros((nChannel, 4 + proc_type*3, nHei))
1431
1432 if proc_type == 1:
1433 type1 = mode_fit
1434 fwindow = numpy.zeros(absc.size) + 1
1435 if exp == '150EEJ':
1436 b=64
1437 fwindow[0:absc.size//2 - b] = 0
1438 fwindow[absc.size//2 + b:] = 0
1439 vers = 1 # new
1440 nProfiles = dataOut.nProfiles
1441 nCohInt = dataOut.nCohInt
1442 nIncohInt = dataOut.nIncohInt
1443 M = numpy.power(numpy.array(1/(nProfiles * nCohInt) ,dtype='float32'),2)
1444 N = numpy.array(M / nIncohInt,dtype='float32')
1445 data = dataOut.data_pre[0] * N
1446 #noise = dataOut.noise * N
1447 noise = numpy.zeros(nChannel)
1448 for ind in range(nChannel):
1449 noise[ind] = self.__NoiseByChannel(nProfiles, nIncohInt, data[ind,:,:])
1450 smooth=3
1451 else:
1452 data = dataOut.data_pre[0]
1453 noise = dataOut.noise
1454 fwindow = None
1455 vers = 0 # old
1456 nIncohInt = None
1457 smooth=None
1458
1459 for ind in range(nChannel):
1460 data_param[ind,:,:] = self.__calculateMoments( data[ind,:,:] , absc , noise[ind], nicoh=nIncohInt, smooth=smooth, type1=type1, fwindow=fwindow, vers=vers)
1461 #print('snr:',data_param[:,0])
1462
1463 if proc_type == 1:
1464 dataOut.moments = data_param[:,1:,:]
1465 #dataOut.data_dop = data_param[:,0]
1466 dataOut.data_dop = data_param[:,2]
1467 dataOut.data_width = data_param[:,1]
1468 # dataOut.data_snr = data_param[:,2]
1469 dataOut.data_snr = data_param[:,0]
1470 dataOut.data_pow = data_param[:,6] # to compare with type0 proccessing
1471 dataOut.spcpar=numpy.stack((dataOut.data_dop,dataOut.data_width,dataOut.data_snr, data_param[:,3], data_param[:,4],data_param[:,5]),axis=2)
1472
1473 else:
1474 dataOut.moments = data_param[:,1:,:]
1475 dataOut.data_snr = data_param[:,0]
1476 dataOut.data_pow = data_param[:,1]
1477 dataOut.data_dop = data_param[:,2]
1478 dataOut.data_width = data_param[:,3]
1479 dataOut.spcpar=numpy.stack((dataOut.data_dop,dataOut.data_width,dataOut.data_snr, dataOut.data_pow),axis=2)
1480
1481 return dataOut
1482
1483 def __calculateMoments(self, oldspec, oldfreq, n0,
1429 def __calculateMoments(self, oldspec, oldfreq, n0,
1484 nicoh = None, graph = None, smooth = None, type1 = None, fwindow = None, snrth = None, dc = None, aliasing = None, oldfd = None, wwauto = None, vers= None):
1430 nicoh = None, graph = None, smooth = None, type1 = None, fwindow = None, \
1431 snrth = None, dc = None, aliasing = None, oldfd = None, wwauto = None, \
1432 vers= None, Hei= None, debug=False, dbg_hei=None, ymax=0.1, curr_ch=0, sel_ch=[0,1]):
1485
1433
1486 def __GAUSSWINFIT1(A, flagPDER=0):
1434 def __GAUSSWINFIT1(A, flagPDER=0):
1487 nonlocal truex, xvalid
1435 nonlocal truex, xvalid
@@ -1612,7 +1560,6 class SpectralMoments(Operation):
1612 #print('Failed to converge.')
1560 #print('Failed to converge.')
1613 chi2 = chisq # Return chi-squared (chi2 obsolete-still works)
1561 chi2 = chisq # Return chi-squared (chi2 obsolete-still works)
1614 if done_early: Niter -= 1
1562 if done_early: Niter -= 1
1615 #save_tp(chisq,Niter,yfit)
1616 return yfit, a, converge, sigma, chisq, chi2 # return result
1563 return yfit, a, converge, sigma, chisq, chi2 # return result
1617 ten=numpy.array([10.0],dtype='f4')[0]
1564 ten=numpy.array([10.0],dtype='f4')[0]
1618 flambda *= ten # Assume fit got worse
1565 flambda *= ten # Assume fit got worse
@@ -1625,21 +1572,30 class SpectralMoments(Operation):
1625 if ((chisq1-chisq)/chisq1) <= tol: # Finished?
1572 if ((chisq1-chisq)/chisq1) <= tol: # Finished?
1626 chi2 = chisq # Return chi-squared (chi2 obsolete-still works)
1573 chi2 = chisq # Return chi-squared (chi2 obsolete-still works)
1627 if done_early: Niter -= 1
1574 if done_early: Niter -= 1
1628 #save_tp(chisq,Niter,yfit)
1629 return yfit, a, converge, sigma, chisq, chi2 # return result
1575 return yfit, a, converge, sigma, chisq, chi2 # return result
1630 converge = 0
1576 converge = 0
1631 chi2 = chisq
1577 chi2 = chisq
1632 #print('Failed to converge.')
1578 #print('Failed to converge.')
1633 #save_tp(chisq,Niter,yfit)
1634 return yfit, a, converge, sigma, chisq, chi2
1579 return yfit, a, converge, sigma, chisq, chi2
1635
1580
1581
1582 def spectral_cut(Hei, ind, dbg_hei, freq, fd, snr, n1, w, ymax, spec, spec2, n0, max_spec, ss1, m, bb0, curr_ch, sel_ch):
1583 if Hei[ind] > dbg_hei[0] and Hei[ind] < dbg_hei[1] and (curr_ch in sel_ch):
1584 nsa=len(freq)
1585 aux='H=%iKm, dop: %4.1f, snr: %4.1f, noise: %4.1f, sw: %4.1f'%(Hei[ind],fd, 10*numpy.log10(snr),10*numpy.log10(n1), w)
1586 plt.subplots()
1587 plt.ylim(0,ymax)
1588 plt.plot(freq,spec,'b-',freq,spec2,'b--', freq,numpy.repeat(n1, nsa),'k--', freq,numpy.repeat(n0, nsa),'k-', freq,numpy.repeat(max_spec, nsa),'y.-', numpy.repeat(fd, nsa),numpy.linspace(0,ymax,nsa),'r--', numpy.repeat(freq[ss1], nsa),numpy.linspace(0,ymax,nsa),'g-.', numpy.repeat(freq[m + bb0], nsa),numpy.linspace(0,ymax,nsa),'g-.')
1589 plt.title(aux)
1590 plt.show()
1591
1592
1636 if (nicoh is None): nicoh = 1
1593 if (nicoh is None): nicoh = 1
1637 if (smooth is None): smooth = 0
1594 if (smooth is None): smooth = 0
1638 if (type1 is None): type1 = 0
1595 if (type1 is None): type1 = 0
1639 if (vers is None): vers = 0
1596 if (vers is None): vers = 0
1640 if (fwindow is None): fwindow = numpy.zeros(oldfreq.size) + 1
1597 if (fwindow is None): fwindow = numpy.zeros(oldfreq.size) + 1
1641 if (snrth is None): snrth = -20.0
1598 if (snrth is None): snrth = -20.0
1642 #if (snrth is None): snrth = -21.0 # abs test
1643 if (dc is None): dc = 0
1599 if (dc is None): dc = 0
1644 if (aliasing is None): aliasing = 0
1600 if (aliasing is None): aliasing = 0
1645 if (oldfd is None): oldfd = 0
1601 if (oldfd is None): oldfd = 0
@@ -1659,7 +1615,6 class SpectralMoments(Operation):
1659 vec_sigma_fd = numpy.empty(oldspec.shape[1])
1615 vec_sigma_fd = numpy.empty(oldspec.shape[1])
1660
1616
1661 for ind in range(oldspec.shape[1]):
1617 for ind in range(oldspec.shape[1]):
1662
1663 spec = oldspec[:,ind]
1618 spec = oldspec[:,ind]
1664 if (smooth == 0):
1619 if (smooth == 0):
1665 spec2 = spec
1620 spec2 = spec
@@ -1673,13 +1628,10 class SpectralMoments(Operation):
1673 if m > 2 and m < oldfreq.size - 3:
1628 if m > 2 and m < oldfreq.size - 3:
1674 newindex = m + numpy.array([-2,-1,0,1,2])
1629 newindex = m + numpy.array([-2,-1,0,1,2])
1675 newfreq = numpy.arange(20)/20.0*(numpy.max(freq[newindex])-numpy.min(freq[newindex]))+numpy.min(freq[newindex])
1630 newfreq = numpy.arange(20)/20.0*(numpy.max(freq[newindex])-numpy.min(freq[newindex]))+numpy.min(freq[newindex])
1676 #peakspec = SPLINE(,)
1677 tck = interpolate.splrep(freq[newindex], spec2[newindex])
1631 tck = interpolate.splrep(freq[newindex], spec2[newindex])
1678 peakspec = interpolate.splev(newfreq, tck)
1632 peakspec = interpolate.splev(newfreq, tck)
1679 # max_spec = MAX(peakspec,)
1680 max_spec = numpy.max(peakspec)
1633 max_spec = numpy.max(peakspec)
1681 mnew = numpy.argmax(peakspec)
1634 mnew = numpy.argmax(peakspec)
1682 #fp = newfreq(mnew)
1683 fp = newfreq[mnew]
1635 fp = newfreq[mnew]
1684 else:
1636 else:
1685 fp = freq[m]
1637 fp = freq[m]
@@ -1760,6 +1712,10 class SpectralMoments(Operation):
1760 except:
1712 except:
1761 w = float("NaN")
1713 w = float("NaN")
1762 snr = power/(n0*fwindow.sum())
1714 snr = power/(n0*fwindow.sum())
1715
1716 if debug:
1717 spectral_cut(Hei, ind, dbg_hei, freq, fd, snr, n1, w, ymax, spec, spec2, n0, max_spec, ss1, m, bb0, curr_ch, sel_ch)
1718
1763 if snr < 1.e-20: snr = 1.e-20
1719 if snr < 1.e-20: snr = 1.e-20
1764
1720
1765 # Here start gaussean adjustment
1721 # Here start gaussean adjustment
@@ -1814,8 +1770,7 class SpectralMoments(Operation):
1814 vec_power[ind] = power # to compare with type 0 proccessing
1770 vec_power[ind] = power # to compare with type 0 proccessing
1815
1771
1816 if vers==1:
1772 if vers==1:
1817 #return numpy.vstack((vec_fd, vec_w, vec_snr, vec_n1, vec_fp, vec_sigma_fd, vec_power))
1773 return numpy.vstack((vec_snr, vec_w, vec_fd, vec_n1, vec_fp, vec_sigma_fd, vec_power))
1818 return numpy.vstack((vec_snr, vec_w, vec_fd, vec_n1, vec_fp, vec_sigma_fd, vec_power)) # snr and fd exchanged to compare doppler of both types
1819 else:
1774 else:
1820 return numpy.vstack((vec_snr, vec_power, vec_fd, vec_w))
1775 return numpy.vstack((vec_snr, vec_power, vec_fd, vec_w))
1821
1776
@@ -1853,19 +1808,9 class SpectralMoments(Operation):
1853 Rutina para cΓ‘lculo de ruido por alturas(n0). Similar a IDL
1808 Rutina para cΓ‘lculo de ruido por alturas(n0). Similar a IDL
1854 '''
1809 '''
1855 num_pts = numpy.size(power)
1810 num_pts = numpy.size(power)
1856 #print('num_pts',num_pts)
1857 #print('power',power.shape)
1858 #print(power[256:267,0:2])
1859 fft_avg = fft_avg*1.0
1811 fft_avg = fft_avg*1.0
1860
1861 ind = numpy.argsort(power, axis=None, kind='stable')
1812 ind = numpy.argsort(power, axis=None, kind='stable')
1862 #ind = numpy.argsort(numpy.reshape(power,-1))
1863 #print(ind.shape)
1864 #print(ind[0:11])
1865 #print(numpy.reshape(power,-1)[ind[0:11]])
1866 ARR = numpy.reshape(power,-1)[ind]
1813 ARR = numpy.reshape(power,-1)[ind]
1867 #print('ARR',len(ARR))
1868 #print('ARR',ARR.shape)
1869 NUMS_MIN = num_pts//10
1814 NUMS_MIN = num_pts//10
1870 RTEST = (1.0+1.0/fft_avg)
1815 RTEST = (1.0+1.0/fft_avg)
1871 SUM = 0.0
1816 SUM = 0.0
@@ -1904,20 +1849,136 class SpectralMoments(Operation):
1904
1849
1905 if talk:
1850 if talk:
1906 print('noise =', noise)
1851 print('noise =', noise)
1907 return noise
1852 return noise, stdvnoise
1853
1854 def run(self, dataOut, proc_type=0, mode_fit=0, exp='150EEJ', debug=False, dbg_hei=None, ymax=1, sel_ch=[0,1]):
1855
1856 absc = dataOut.abscissaList[:-1]
1857 nChannel = dataOut.data_pre[0].shape[0]
1858 nHei = dataOut.data_pre[0].shape[2]
1859 Hei=dataOut.heightList
1860 data_param = numpy.zeros((nChannel, 4 + proc_type*3, nHei))
1861 nProfiles = dataOut.nProfiles
1862 nCohInt = dataOut.nCohInt
1863 nIncohInt = dataOut.nIncohInt
1864 M = numpy.power(numpy.array(1/(nProfiles * nCohInt) ,dtype='float32'),2)
1865 N = numpy.array(M / nIncohInt,dtype='float32')
1866
1867 if proc_type == 1:
1868 type1 = mode_fit
1869 fwindow = numpy.zeros(absc.size) + 1
1870 if exp == '150EEJ':
1871 b=64
1872 fwindow[0:absc.size//2 - b] = 0
1873 fwindow[absc.size//2 + b:] = 0
1874 vers = 1 # new
1875
1876 data = dataOut.data_pre[0] * N
1877
1878 noise = numpy.zeros(nChannel)
1879 stdvnoise = numpy.zeros(nChannel)
1880 for ind in range(nChannel):
1881 noise[ind], stdvnoise[ind] = self.__NoiseByChannel(nProfiles, nIncohInt, data[ind,:,:])
1882 smooth=3
1883 else:
1884 data = dataOut.data_pre[0]
1885 noise = dataOut.noise
1886 fwindow = None
1887 type1 = None
1888 vers = 0 # old
1889 nIncohInt = None
1890 smooth=None
1891
1892 for ind in range(nChannel):
1893 data_param[ind,:,:] = self.__calculateMoments(data[ind,:,:] , absc , noise[ind], nicoh=nIncohInt, smooth=smooth, type1=type1, fwindow=fwindow, vers=vers, Hei=Hei, debug=debug, dbg_hei=dbg_hei, ymax=ymax, curr_ch=ind, sel_ch=sel_ch)
1894 #data_param[ind,:,:] = self.__calculateMoments(data[ind,:,:] , absc , noise[ind], nicoh=nIncohInt, smooth=smooth, type1=type1, fwindow=fwindow, vers=vers, Hei=Hei, debug=debug)
1895 if exp == 'ESF_EW':
1896 data_param[ind,0,:]*=(noise[ind]/stdvnoise[ind])
1897 data_param[ind,3,:]*=(1.0/M)
1908
1898
1909 class JULIADriftsEstimation(Operation):
1899 if proc_type == 1:
1900 dataOut.moments = data_param[:,1:,:]
1901 dataOut.data_dop = data_param[:,2]
1902 dataOut.data_width = data_param[:,1]
1903 dataOut.data_snr = data_param[:,0]
1904 dataOut.data_pow = data_param[:,6] # to compare with type0 proccessing
1905 dataOut.spcpar=numpy.stack((dataOut.data_dop,dataOut.data_width,dataOut.data_snr, data_param[:,3], data_param[:,4],data_param[:,5]),axis=2)
1906
1907 if exp == 'ESF_EW':
1908 spc=dataOut.data_pre[0]* N
1909 cspc=dataOut.data_pre[1]* N
1910 nHei=dataOut.data_pre[1].shape[2]
1911 cross_pairs=dataOut.pairsList
1912 nDiffIncohInt = dataOut.nDiffIncohInt
1913 N2=numpy.array(1 / nDiffIncohInt,dtype='float32')
1914 diffcspectra=dataOut.data_diffcspc.copy()* N2 * M * M
1915 num_pairs=len(dataOut.pairsList)
1916
1917 if num_pairs >= 0:
1918 fbinv=numpy.where(absc != 0)[0]
1919 ccf=numpy.sum(cspc[:,fbinv,:], axis=1)
1920 jvpower=numpy.sum(spc[:,fbinv,:], axis=1)
1921 coh=ccf/numpy.sqrt(jvpower[cross_pairs[0][0],:]*jvpower[cross_pairs[0][1],:])
1922 dccf=numpy.sum(diffcspectra[:,fbinv,:], axis=1)
1923 dataOut.ccfpar = numpy.zeros((num_pairs,nHei,3))
1924 dataOut.ccfpar[:,:,0]=numpy.abs(coh)
1925 dataOut.ccfpar[:,:,1]=numpy.arctan(numpy.imag(coh)/numpy.real(coh))
1926 dataOut.ccfpar[:,:,2]=numpy.arctan(numpy.imag(dccf)/numpy.real(dccf))
1927 else:
1928 dataOut.moments = data_param[:,1:,:]
1929 dataOut.data_snr = data_param[:,0]
1930 dataOut.data_pow = data_param[:,1]
1931 dataOut.data_dop = data_param[:,2]
1932 dataOut.data_width = data_param[:,3]
1933 dataOut.spcpar=numpy.stack((dataOut.data_dop,dataOut.data_width,dataOut.data_snr, dataOut.data_pow),axis=2)
1934
1935 return dataOut
1936
1937
1938 class JULIA_DayVelocities(Operation):
1939 '''
1940 Function SpectralMoments()
1941
1942 From espectral parameters calculates:
1943
1944 1. Signal to noise level (SNL)
1945 2. Vertical velocity
1946 3. Zonal velocity
1947 4. Vertical velocity error
1948 5. Zonal velocity error.
1949
1950 Type of dataIn: SpectralMoments
1951
1952 Configuration Parameters:
1953
1954 zenith : Pairs of angles corresponding to the two beams related to the perpendicular to B from the center of the antenna.
1955 zenithCorrection : Adjustment angle for the zenith. Default 0.
1956 heights : Range to process 150kM echoes. By default [125,185].
1957 nchan : To process 2 or 1 channel. 2 by default.
1958 chan : If nchan = 1, chan indicates which of the 2 channels to process.
1959 clean : 2nd cleaning processing (Graphical). Default False
1960 driftstdv_th : Diferencia mΓ‘xima entre valores promedio consecutivos de vertical.
1961 zonalstdv_th : Diferencia mΓ‘xima entre valores promedio consecutivos de zonal.
1962
1963 Input:
1964
1965 Affected:
1966
1967 '''
1910
1968
1911 def __init__(self):
1969 def __init__(self):
1912 Operation.__init__(self)
1970 Operation.__init__(self)
1913
1971 self.old_drift=None
1972 self.old_zonal=None
1973 self.count_drift=0
1974 self.count_zonal=0
1975 self.oldTime_drift=None
1976 self.oldTime_zonal=None
1914
1977
1915 def newtotal(self, data):
1978 def newtotal(self, data):
1916 return numpy.nansum(data)
1979 return numpy.nansum(data)
1917
1980
1918 #def data_filter(self, parm, snrth=-19.5, swth=20, wErrth=500):
1919 def data_filter(self, parm, snrth=-20, swth=20, wErrth=500):
1981 def data_filter(self, parm, snrth=-20, swth=20, wErrth=500):
1920 #def data_filter(self, parm, snrth=-21, swth=20, wErrth=500): # abs test
1921
1982
1922 Sz0 = parm.shape # Sz0: h,p
1983 Sz0 = parm.shape # Sz0: h,p
1923 drift = parm[:,0]
1984 drift = parm[:,0]
@@ -1958,30 +2019,29 class JULIADriftsEstimation(Operation):
1958 cnovalid = len(novalid[0])
2019 cnovalid = len(novalid[0])
1959 if cnovalid > 0:
2020 if cnovalid > 0:
1960 mask[novalid] = numpy.nan
2021 mask[novalid] = numpy.nan
2022
1961 new_parm = numpy.zeros((Sz0[0],Sz0[1]))
2023 new_parm = numpy.zeros((Sz0[0],Sz0[1]))
1962 for h in range(Sz0[0]):
2024 for i in range(Sz0[1]):
1963 for p in range(Sz0[1]):
2025 new_parm[:,i] = parm[:,i] * mask
1964 if numpy.isnan(mask[h]):
1965 new_parm[h,p]=numpy.nan
1966 else:
1967 new_parm[h,p]=parm[h,p]
1968
2026
1969 return new_parm, th
2027 return new_parm, th
1970
2028
1971 def statistics150km(self, veloc , sigma , threshold , currTime=None, \
1972 amountdata=2, clearAll = None, timeFactor=None):
1973
2029
1974 step = threshold/2
2030 def statistics150km(self, veloc , sigma , threshold , old_veloc=None, count=0, \
2031 currTime=None, oldTime=None, amountdata=2, clearAll = None, timeFactor=1800, debug = False):
2032
2033 if oldTime == None:
2034 oldTime = currTime
2035
2036 step = (threshold/2)*(numpy.abs(currTime - oldTime)//timeFactor + 1)
1975 factor = 2
2037 factor = 2
1976 avg_threshold = 100
2038 avg_threshold = 100
1977
1978 # Calcula la mediana en todas las alturas por tiempo
2039 # Calcula la mediana en todas las alturas por tiempo
1979 val1=numpy.nanmedian(veloc)
2040 val1=numpy.nanmedian(veloc)
1980
2041
1981 # Calcula la media ponderada en todas las alturas por tiempo
2042 # Calcula la media ponderada en todas las alturas por tiempo
1982 val2 = self.newtotal(veloc/numpy.power(sigma,2))/self.newtotal(1/numpy.power(sigma,2))
2043 val2 = self.newtotal(veloc/numpy.power(sigma,2))/self.newtotal(1/numpy.power(sigma,2))
1983
2044
1984
1985 # Verifica la cercanΓ­a de los valores calculados de mediana y media, si son cercanos escoge la media ponderada
2045 # Verifica la cercanΓ­a de los valores calculados de mediana y media, si son cercanos escoge la media ponderada
1986 op1=numpy.abs(val2-val1)
2046 op1=numpy.abs(val2-val1)
1987 op2=threshold/factor
2047 op2=threshold/factor
@@ -1996,12 +2056,65 class JULIADriftsEstimation(Operation):
1996
2056
1997 # Se calcula nuevamente media ponderada, en base a estimado inicial de la media
2057 # Se calcula nuevamente media ponderada, en base a estimado inicial de la media
1998 # a fin de eliminar valores que estΓ‘n muy lejanos a dicho valor
2058 # a fin de eliminar valores que estΓ‘n muy lejanos a dicho valor
2059
2060 if debug:
2061 print('veloc_prof:', veloc_prof)
2062 print('veloc:',veloc)
2063 print('threshold:',threshold)
2064 print('factor:',factor)
2065 print('threshold/factor:',threshold/factor)
2066 print('numpy.abs(veloc-veloc_prof):', numpy.abs(veloc-veloc_prof))
2067 print('numpy.where(numpy.abs(veloc-veloc_prof) < threshold/factor)[0]:', numpy.where(numpy.abs(veloc-veloc_prof) < threshold/factor)[0])
2068
1999 junk = numpy.where(numpy.abs(veloc-veloc_prof) < threshold/factor)[0]
2069 junk = numpy.where(numpy.abs(veloc-veloc_prof) < threshold/factor)[0]
2070 if junk.size > 2:
2071 veloc_prof = self.newtotal(veloc[junk]/numpy.power(sigma[junk],2))/self.newtotal(1/numpy.power(sigma[junk],2))
2072 sigma_prof1 = numpy.sqrt(1/self.newtotal(1/numpy.power(sigma[junk],2)))
2073 sigma_prof2 = numpy.sqrt(self.newtotal(numpy.power(veloc[junk]-veloc_prof,2)/numpy.power(sigma[junk],2)))*sigma_prof1
2074 sigma_prof = numpy.sqrt(numpy.power(sigma_prof1,2)+numpy.power(sigma_prof2,2))
2075 sets = junk
2076
2077 # Compara con valor anterior para evitar presencia de "outliers"
2078 if debug:
2079 print('old_veloc:',old_veloc)
2080 print('step:', step)
2081
2082 if old_veloc == None:
2083 valid=numpy.isfinite(veloc_prof)
2084 else:
2085 valid=numpy.abs(veloc_prof-old_veloc) < step
2000
2086
2001 return junk
2087 if debug:
2088 print('valid:', valid)
2089
2090 if not valid:
2091 aver_veloc=numpy.nan
2092 aver_sigma=numpy.nan
2093 sets=numpy.array([-1])
2094 else:
2095 aver_veloc=veloc_prof
2096 aver_sigma=sigma_prof
2097 clearAll=0
2098 if old_veloc != None and count < 5:
2099 if numpy.abs(veloc_prof-old_veloc) > step:
2100 clearAll=1
2101 count=0
2102 old_veloc=None
2103 if numpy.isfinite(aver_veloc):
2104
2105 count+=1
2106 if old_veloc != None:
2107 old_veloc = (old_veloc + aver_veloc) * 0.5
2108 else:
2109 old_veloc=aver_veloc
2110 oldTime=currTime
2111 if debug:
2112 print('count:',count)
2113 print('sets:',sets)
2114 return sets, old_veloc, count, oldTime, aver_veloc, aver_sigma, clearAll
2002
2115
2003 def run(self, dataOut, zenith, zenithCorrection,heights=None, otype=0, nchan=2, chan=0):
2004
2116
2117 def run(self, dataOut, zenith, zenithCorrection=0.0, heights=[125, 185], nchan=2, chan=0, clean=False, driftstdv_th=100, zonalstdv_th=200):
2005
2118
2006 dataOut.lat=-11.95
2119 dataOut.lat=-11.95
2007 dataOut.lon=-76.87
2120 dataOut.lon=-76.87
@@ -2009,21 +2122,14 class JULIADriftsEstimation(Operation):
2009 nCh=dataOut.spcpar.shape[0]
2122 nCh=dataOut.spcpar.shape[0]
2010 nHei=dataOut.spcpar.shape[1]
2123 nHei=dataOut.spcpar.shape[1]
2011 nParam=dataOut.spcpar.shape[2]
2124 nParam=dataOut.spcpar.shape[2]
2012 # SelecciΓ³n de alturas
2013
2125
2014 if not heights:
2126 # SelecciΓ³n de alturas
2015 parm = numpy.zeros((nCh,nHei,nParam))
2016 parm[:] = dataOut.spcpar[:]
2017 else:
2018 hei=dataOut.heightList
2127 hei=dataOut.heightList
2019 hvalid=numpy.where([hei >= heights[0]][0] & [hei <= heights[1]][0])[0]
2128 hvalid=numpy.where([hei >= heights[0]][0] & [hei <= heights[1]][0])[0]
2020 nhvalid=len(hvalid)
2129 nhvalid=len(hvalid)
2021 dataOut.heightList = hei[hvalid]
2130 dataOut.heightList = hei[hvalid]
2022 parm = numpy.zeros((nCh,nhvalid,nParam))
2131 parm=numpy.empty((nCh,nhvalid,nParam)); parm[:]=numpy.nan
2023 parm[:] = dataOut.spcpar[:,hvalid,:]
2132 parm[:] = dataOut.spcpar[:,hvalid,:]
2024 print('parm:',parm.shape)
2025
2026
2027 # Primer filtrado: Umbral de SNR
2133 # Primer filtrado: Umbral de SNR
2028 for i in range(nCh):
2134 for i in range(nCh):
2029 parm[i,:,:] = self.data_filter(parm[i,:,:])[0]
2135 parm[i,:,:] = self.data_filter(parm[i,:,:])[0]
@@ -2071,18 +2177,21 class JULIADriftsEstimation(Operation):
2071 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))
2177 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))
2072
2178
2073 # 150Km statistics to clean data
2179 # 150Km statistics to clean data
2074
2075 clean_drift = drift.copy()
2180 clean_drift = drift.copy()
2076
2077 clean_drift[:] = numpy.nan
2181 clean_drift[:] = numpy.nan
2078 if nchan == 2:
2182 if nchan == 2:
2079 clean_zonal = zonal.copy()
2183 clean_zonal = zonal.copy()
2080 clean_zonal[:] = numpy.nan
2184 clean_zonal[:] = numpy.nan
2081
2185
2082 # Drifts
2186 # Vertical
2083 driftstdv_th = 20*2
2187 sets1, self.old_drift, self.count_drift, self.oldTime_drift, aver_veloc, aver_sigma, clearAll = self.statistics150km(drift, w_w_err, driftstdv_th, \
2084
2188 old_veloc=self.old_drift, count=self.count_drift, currTime=dataOut.utctime, \
2085 sets1 = self.statistics150km(drift, w_w_err, driftstdv_th)
2189 oldTime=self.oldTime_drift, timeFactor=120)
2190 if clearAll == 1:
2191 mean_zonal = numpy.nan
2192 sigma_zonal = numpy.nan
2193 mean_drift = aver_veloc
2194 sigma_drift = aver_sigma
2086
2195
2087 if sets1.size != 1:
2196 if sets1.size != 1:
2088 clean_drift[sets1] = drift[sets1]
2197 clean_drift[sets1] = drift[sets1]
@@ -2093,9 +2202,14 class JULIADriftsEstimation(Operation):
2093
2202
2094 # Zonal
2203 # Zonal
2095 if nchan == 2:
2204 if nchan == 2:
2096 zonalstdv_th = 30*2
2205 sets2, self.old_zonal, self.count_zonal, self.oldTime_zonal, aver_veloc, aver_sigma, clearAll = self.statistics150km(zonal, w_e_err, zonalstdv_th, \
2097 sets2 = self.statistics150km(zonal, w_e_err, zonalstdv_th)
2206 old_veloc=self.old_zonal, count=self.count_zonal, currTime=dataOut.utctime, \
2098
2207 oldTime=self.oldTime_zonal, timeFactor=600)
2208 if clearAll == 1:
2209 mean_zonal = numpy.nan
2210 sigma_zonal = numpy.nan
2211 mean_zonal = aver_veloc
2212 sigma_zonal = aver_sigma
2099 if sets2.size != 1:
2213 if sets2.size != 1:
2100 clean_zonal[sets2] = zonal[sets2]
2214 clean_zonal[sets2] = zonal[sets2]
2101
2215
@@ -2103,22 +2217,105 class JULIADriftsEstimation(Operation):
2103 if cnovalid > 0: zonal[novalid] = numpy.nan
2217 if cnovalid > 0: zonal[novalid] = numpy.nan
2104 if cnovalid > 0: snr[novalid] = numpy.nan
2218 if cnovalid > 0: snr[novalid] = numpy.nan
2105
2219
2220 n_avg_par=4
2221 avg_par=numpy.empty((n_avg_par,)); avg_par[:] = numpy.nan
2222 avg_par[0,]=mean_drift
2223 avg_par[1,]=mean_zonal
2224 avg_par[2,]=sigma_drift
2225 avg_par[3,]=sigma_zonal
2226
2227 set1 = 1.0
2228 navg = set1
2229 nci = dataOut.nCohInt
2230 # ----------------------------------
2231 ipp = 252.0
2232 nincoh = dataOut.nIncohInt
2233 nptsfft = dataOut.nProfiles
2234 hardcoded=False # if True, similar to IDL processing
2235 if hardcoded:
2236 ipp=200.1
2237 nincoh=22
2238 nptsfft=128
2239 # ----------------------------------
2240 nipp = ipp * nci
2241 height = dataOut.heightList
2242 nHei = len(height)
2243 kd = 213.6
2244 nint = nptsfft * nincoh
2245 drift1D = drift.copy()
2246 if nchan == 2:
2247 zonal1D=zonal.copy()
2248 snr1D = snr.copy()
2249 snr1D = 10*numpy.power(10, 0.1*snr1D)
2250 noise1D = noise.copy()
2251 noise0 = numpy.nanmedian(noise1D)
2252 noise = noise0 + noise0
2253 sw1D = sw.copy()
2254 pow0 = snr1D * noise0 + noise0
2255 acf0 = snr1D * noise0 * numpy.exp((-drift1D*nipp*numpy.pi/(1.5e5*1.5))*1j) * (1-0.5*numpy.power(sw1D*nipp*numpy.pi/(1.5e5*1.5),2))
2256 acf0 /= pow0
2257 acf1 = acf0
2258 dt= nint * nipp /1.5e5
2106
2259
2107 if otype == 0:
2260 if nchan == 2:
2108 winds = numpy.vstack((snr, drift, zonal, noise, sw, w_w_err, w_e_err)) # to process statistics drifts
2261 dccf = pow0 * pow0 * numpy.exp((zonal1D*kd*dt/(height*1e3))*(1j))
2109 elif otype == 2:
2262 else:
2110 winds = numpy.vstack((snr, drift)) # one channel good signal: 2 RTI's
2263 dccf = numpy.empty(nHei); dccf[:]=numpy.nan # complex?
2111 elif otype == 3:
2264 dccf /= pow0 * pow0
2112 winds = numpy.vstack((snr, drift, zonal)) # to generic plot: 3 RTI's
2265 sno=(pow0+pow0-noise)/noise
2113 elif otype == 4:
2266
2114 winds = numpy.vstack((snrCH0, drift, snrCH1, zonal)) # to generic plot: 4 RTI's
2267 # First parameter: Signal to noise ratio and its error
2268 sno = numpy.log10(sno)
2269 sno10 = 10 * sno
2270 dsno = 1.0/numpy.sqrt(nint*navg)*(1+1/sno10)
2271
2272 # Second parameter: Vertical Drifts
2273 s=numpy.sqrt(numpy.abs(acf0)*numpy.abs(acf1))
2274 sp = s*(1.0 + 1.0/sno10)
2275 vzo = -numpy.arctan2(numpy.imag(acf0+acf1),numpy.real(acf0+acf1))* \
2276 1.5e5*1.5/(nipp*numpy.pi)
2277 dvzo = numpy.sqrt(1-sp*sp)*0.338*1.5e5/(numpy.sqrt(nint*navg)*sp*nipp)
2278
2279 # Third parameter: Zonal Drifts
2280 dt = nint*nipp/1.5e5
2281 ss = numpy.sqrt(numpy.abs(dccf))
2282 vxo = numpy.arctan2(numpy.imag(dccf),numpy.real(dccf))*height*1e3/(kd*dt)
2283 dvxo = numpy.sqrt(1.0-ss*ss)*height*1e3/(numpy.sqrt(nint*navg)*ss*kd*dt)
2284
2285 npar = 5
2286 par = numpy.empty((npar, nHei)); par[:] = numpy.nan
2287
2288 par[0,:] = sno
2289 par[1,:] = vzo
2290 par[2,:] = vxo
2291 par[3,:] = dvzo
2292 par[4,:] = dvxo
2293
2294 # Segundo filtrado:
2295 # RemociΓ³n por altura: Menos de dos datos finitos no son considerados como eco 150Km.
2296 clean_par=numpy.empty((npar,nHei)); clean_par[:]=numpy.nan
2297 if clean:
2298
2299 for p in range(npar):
2300 ih=0
2301 while ih < nHei-1:
2302 j=ih
2303 if numpy.isfinite(snr1D[ih]):
2304 while numpy.isfinite(snr1D[j]):
2305 j+=1
2306 if j >= nHei:
2307 break
2308 if j > ih + 1:
2309 for k in range(ih,j):
2310 clean_par[p][k] = par[p][k]
2311 ih = j - 1
2312 ih+=1
2313 else:
2314 clean_par[:] = par[:]
2115
2315
2116 snr1 = numpy.vstack((snrCH0, snrCH1))
2316 winds = numpy.vstack((clean_par[0,:], clean_par[1,:], clean_par[2,:], clean_par[3,:], clean_par[4,:]))
2117 print('winds:',winds.shape)
2118 print('snrCH0:',snrCH0.shape)
2119 dataOut.data_output = winds
2317 dataOut.data_output = winds
2120 dataOut.data_snr = snr1
2318 dataOut.avg_output = avg_par
2121
2122 dataOut.utctimeInit = dataOut.utctime
2319 dataOut.utctimeInit = dataOut.utctime
2123 dataOut.outputInterval = dataOut.timeInterval
2320 dataOut.outputInterval = dataOut.timeInterval
2124
2321
@@ -2126,6 +2323,351 class JULIADriftsEstimation(Operation):
2126
2323
2127 return dataOut
2324 return dataOut
2128
2325
2326
2327 class JULIA_NightVelocities(Operation):
2328 '''
2329 Function SpreadFVelocities()
2330
2331 Calculates SNL and drifts
2332
2333 Type of dataIn: Parameters
2334
2335 Configuration Parameters:
2336
2337 mymode : (0) Interferometry,
2338 (1) Doppler beam swinging.
2339 myproc : (0) JULIA_V,
2340 (1) JULIA_EW.
2341 myantenna : (0) 1/4 antenna,
2342 (1) 1/2 antenna.
2343 jset : Number of Incoherent integrations.
2344
2345
2346 Input:
2347 channelList : simple channel list to select e.g. [2,3,7]
2348 self.dataOut.data_pre : Spectral data
2349 self.dataOut.abscissaList : List of frequencies
2350 self.dataOut.noise : Noise level per channel
2351
2352 Affected:
2353 self.dataOut.moments : Parameters per channel
2354 self.dataOut.data_snr : SNR per channel
2355
2356 '''
2357 def __init__(self):
2358 Operation.__init__(self)
2359
2360 def newtotal(self, data):
2361 return numpy.nansum(data)
2362
2363 def data_filter(self, parm, snrth=-17, swth=20, dopth=500.0, debug=False):
2364
2365 Sz0 = parm.shape # Sz0: h,p
2366 drift = parm[:,0]
2367 sw = 2*parm[:,1]
2368 snr = 10*numpy.log10(parm[:,2])
2369 Sz = drift.shape # Sz: h
2370 mask = numpy.ones((Sz[0]))
2371 th=0
2372 valid=numpy.where(numpy.isfinite(snr))
2373 cvalid = len(valid[0])
2374 if cvalid >= 1:
2375 # CΓ‘lculo del ruido promedio de snr para el i-Γ©simo grupo de alturas
2376 nbins = int(numpy.max(snr)-numpy.min(snr))+1 # bin size = 1, similar to IDL
2377 h = numpy.histogram(snr,bins=nbins)
2378 hist = h[0]
2379 values = numpy.round_(h[1])
2380 moda = values[numpy.where(hist == numpy.max(hist))]
2381 indNoise = numpy.where(numpy.abs(snr - numpy.min(moda)) < 3)[0]
2382
2383 noise = snr[indNoise]
2384 noise_mean = numpy.sum(noise)/len(noise)
2385 # CΓ‘lculo de media de snr
2386 med = numpy.median(snr)
2387 # Establece el umbral de snr
2388 if noise_mean > med + 3:
2389 th = med
2390 else:
2391 th = noise_mean + 3
2392 # Establece mΓ‘scara
2393 novalid = numpy.where(snr <= th)[0]
2394 mask[novalid] = numpy.nan
2395 # Elimina datos que no sobrepasen el umbral: PARAMETRO
2396 novalid = numpy.where(snr <= snrth)
2397 cnovalid = len(novalid[0])
2398 if cnovalid > 0:
2399 mask[novalid] = numpy.nan
2400 novalid = numpy.where(numpy.isnan(snr))
2401 cnovalid = len(novalid[0])
2402 if cnovalid > 0:
2403 mask[novalid] = numpy.nan
2404 # umbral de velocidad
2405 if dopth != None:
2406 novalid = numpy.where(numpy.logical_or(drift< dopth*(-1), drift > dopth))
2407 cnovalid = len(novalid[0])
2408 if cnovalid > 0:
2409 mask[novalid] = numpy.nan
2410 if debug:
2411 print('Descartados:%i de %i:' %(cnovalid, len(drift)))
2412 print('Porcentaje:%3.1f' %(100.0*cnovalid/len(drift)))
2413
2414 new_parm = numpy.zeros((Sz0[0],Sz0[1]))
2415 for i in range(Sz0[1]):
2416 new_parm[:,i] = parm[:,i] * mask
2417
2418 return new_parm, mask
2419
2420
2421 def run(self, dataOut, zenith, zenithCorrection, mymode=1, dbs_sel=0, myproc=0, myantenna=0, jset=None, clean=False):
2422
2423
2424 dataOut.lat=-11.95
2425 dataOut.lon=-76.87
2426 mode=mymode
2427 proc=myproc
2428 antenna=myantenna
2429 nci=dataOut.nCohInt
2430 nptsfft=dataOut.nProfiles
2431 navg= 3 if jset is None else jset
2432 nint=dataOut.nIncohInt//navg
2433 navg1=dataOut.nProfiles * nint * navg
2434 tau1=dataOut.ippSeconds
2435 nipp=dataOut.radarControllerHeaderObj.ipp
2436 jlambda=6
2437 kd=213.6
2438 hei=dataOut.heightList.copy()
2439
2440 nCh=dataOut.spcpar.shape[0]
2441 nHei=dataOut.spcpar.shape[1]
2442 nParam=dataOut.spcpar.shape[2]
2443
2444 parm = numpy.zeros((nCh,nHei,nParam))
2445 parm[:] = dataOut.spcpar[:]
2446 mask=numpy.ones(nHei)
2447 mask0=mask.copy()
2448 # Primer filtrado: Umbral de SNR
2449 for i in range(nCh):
2450 parm[i,:,:], mask = self.data_filter(parm[i,:,:], snrth = 0.1) # umbral 0.1 filtra seΓ±al que no corresponde a ESF, para interferometrΓ­a usar -17dB
2451 mask0 *= mask
2452
2453 ccf_results=numpy.transpose(dataOut.ccfpar,(2,1,0))
2454
2455 for i in range(3):
2456 ccf_results[i,:,0] *= mask0
2457
2458 zenith = numpy.array(zenith)
2459 zenith -= zenithCorrection
2460 zenith *= numpy.pi/180
2461 alpha = zenith[0]
2462 beta = zenith[1]
2463
2464 w_w = parm[0,:,0]
2465 w_e = parm[1,:,0]
2466
2467 if mode==1:
2468 # Vertical and zonal calculation
2469 sinB_A = numpy.sin(beta)*numpy.cos(alpha) - numpy.sin(alpha)* numpy.cos(beta)
2470 w = -(w_w * numpy.sin(beta) - w_e * numpy.sin(alpha))/ sinB_A
2471 u = (w_w * numpy.cos(beta) - w_e * numpy.cos(alpha))/ sinB_A
2472
2473 #Noise
2474 n0 = parm[0,:,3]
2475 n1 = parm[1,:,3]
2476 jn0_1 = numpy.nanmedian(n0)
2477 jn0_2 = numpy.nanmean(n0)
2478 jn1_1 = numpy.nanmedian(n1)
2479 jn1_2 = numpy.nanmean(n1)
2480 noise0 = jn0_2 if numpy.abs(jn0_1-jn0_2)/(jn0_1+jn0_2) <= 0.1 else jn0_1
2481 noise1 = jn1_2 if numpy.abs(jn1_1-jn1_2)/(jn1_1+jn1_2) <= 0.1 else jn1_1
2482
2483 noise = noise0 + noise0 if mode == 1 else noise0 + noise1
2484
2485 #Power
2486 apow1 = (parm[0,:,2]/numpy.sqrt(nint))*noise0 + n0
2487 apow2 = (parm[1,:,2]/numpy.sqrt(nint))*noise1 + n1
2488
2489 #SNR SNR=Detectability/ SQRT(nint) or (Pow-Noise)/Noise
2490 s_n0 = (apow1 - noise0)/noise0
2491 s_n1 = (apow2 - noise1)/noise1
2492
2493 swCH0 = parm[0,:,1]
2494 swCH1 = parm[1,:,1]
2495
2496 if mode == 1:
2497 aacf1=(1-numpy.square(tau1)*numpy.square(4*numpy.pi/jlambda*swCH0)/2)* \
2498 numpy.exp(-4*numpy.pi/jlambda*w*tau1*1j)* \
2499 apow1
2500 aacf2=(1-numpy.square(tau1)*numpy.square(4*numpy.pi/jlambda*swCH1)/2)* \
2501 numpy.exp(-4*numpy.pi/jlambda*w*tau1*1j)* \
2502 apow2
2503 dccf_0=numpy.zeros(nHei, dtype=complex)
2504
2505 else:
2506 aacf1=(1-numpy.square(tau1)*numpy.square(4*numpy.pi/jlambda*swCH0)/2)* \
2507 numpy.exp(4*numpy.pi/jlambda*w_w*tau1*1j)* \
2508 apow1
2509 aacf2=(1-numpy.square(tau1)*numpy.square(4*numpy.pi/jlambda*swCH1)/2)* \
2510 numpy.exp(4*numpy.pi/jlambda*w_e*tau1*1j)* \
2511 apow2
2512 dccf_0=numpy.power(ccf_results[0,:,0],2)*apow1*apow2* \
2513 numpy.exp( \
2514 ( \
2515 (1+1*(antenna==1))* \
2516 (-1+2*(proc == 1))* \
2517 ccf_results[2,:,0] \
2518 )*1j)
2519
2520 nsamp=len(hei)
2521 pow0 = numpy.empty(nsamp); pow0[:] = numpy.nan
2522 pow1 = numpy.empty(nsamp); pow1[:] = numpy.nan
2523 acf0 = numpy.empty(nsamp, dtype=complex); acf0[:] = numpy.nan
2524 acf1 = numpy.empty(nsamp, dtype=complex); acf1[:] = numpy.nan
2525 dccf = numpy.empty(nsamp, dtype=complex); dccf[:] = numpy.nan
2526 dop0 = numpy.empty(nsamp); dop0[:] = numpy.nan
2527 dop1 = numpy.empty(nsamp); dop1[:] = numpy.nan
2528 p_w = numpy.empty(nsamp); p_w[:] = numpy.nan
2529 p_u = numpy.empty(nsamp); p_u[:] = numpy.nan
2530
2531 if mode == 0 or (mode == 1 and dbs_sel == 0):
2532 ih=0
2533 while ih < nsamp-10:
2534 j=ih
2535 if numpy.isfinite(s_n0[ih]) and numpy.isfinite(s_n1[ih]):
2536 while numpy.isfinite(s_n0[j]) and numpy.isfinite(s_n1[j]):
2537 j+=1
2538 if j > ih + 2:
2539 for k in range(ih,j):
2540 pow0[k] = apow1[k]
2541 pow1[k] = apow2[k]
2542 acf0[k] = aacf1[k]
2543 acf1[k] = aacf2[k]
2544 dccf[k] = dccf_0[k]
2545 ih = j - 1
2546 ih+=1
2547 else:
2548 ih=0
2549 while ih < nsamp-10:
2550 j=ih
2551 if numpy.isfinite(s_n0[ih]):
2552 while numpy.isfinite(s_n0[j]) and j < nsamp-10:
2553 j+=1
2554 #if j > ih + 6:
2555 if j > ih + 2:
2556 #if j > ih + 3:
2557 for k in range(ih,j):
2558 pow0[k] = apow1[k]
2559 #acf0[k] = aacf1[k]
2560 #dccf[k] = dccf_0[k]
2561 p_w[k] = w[k]
2562 dop0[k] = w_w[k]
2563 ih = j - 1
2564 ih+=1
2565 ih=0
2566 while ih < nsamp-10:
2567 j=ih
2568 if numpy.isfinite(s_n1[ih]):
2569 while numpy.isfinite(s_n1[j]) and j < nsamp-10:
2570 j+=1
2571 #if j > ih + 6:
2572 if j > ih + 2:
2573 #if j > ih + 3:
2574 for k in range(ih,j):
2575 pow1[k] = apow2[k]
2576 #acf1[k] = aacf2[k]
2577 p_u[k] = u[k]
2578 dop1[k] = w_e[k]
2579 ih = j - 1
2580 ih+=1
2581
2582 acf0 = numpy.zeros(nsamp, dtype=complex)
2583 acf1 = numpy.zeros(nsamp, dtype=complex)
2584 dccf = numpy.zeros(nsamp, dtype=complex)
2585
2586 acf0 /= pow0
2587 acf1 /= pow1
2588 dccf /= pow0 * pow1
2589
2590 if mode == 0 or (mode == 1 and dbs_sel == 0):
2591 sno=(pow0+pow1-noise)/noise
2592 # First parameter: Signal to noise ratio and its error
2593 sno=numpy.log10(sno)
2594 dsno=1.0/numpy.sqrt(nint*navg)*(1+1/sno)
2595 # Second parameter: Vertical Drifts
2596 s=numpy.sqrt(numpy.abs(acf0)*numpy.abs(acf1))
2597 ind=numpy.where(numpy.abs(s)>=1.0)
2598 if numpy.size(ind)>0:
2599 s[ind]=numpy.sqrt(0.9999)
2600 sp=s*(1.0 + 1.0/sno)
2601 vzo=-numpy.arctan2(numpy.imag(acf0+acf1),numpy.real(acf0+acf1))* \
2602 1.5e5*1.5/(nipp*numpy.pi)
2603 dvzo=numpy.sqrt(1-sp*sp)*0.338*1.5e5/(numpy.sqrt(nint*navg)*sp*nipp)
2604 ind=numpy.where(dvzo<=0.1)
2605 if numpy.size(ind)>0:
2606 dvzo[ind]=0.1
2607 # Third parameter: Zonal Drifts
2608 dt=nint*nipp/1.5e5
2609 ss=numpy.sqrt(numpy.abs(dccf))
2610 ind=numpy.where(ss>=1.0)
2611 if numpy.size(ind)>0:
2612 ss[ind]=numpy.sqrt(0.99999)
2613 ind=numpy.where(ss<=0.1)
2614 if numpy.size(ind)>0:
2615 ss[ind]=numpy.sqrt(0.1)
2616 vxo=numpy.arctan2(numpy.imag(dccf),numpy.real(dccf))*hei*1e3/(kd*dt)
2617 dvxo=numpy.sqrt(1.0-ss*ss)*hei*1e3/(numpy.sqrt(nint*navg)*ss*kd*dt)
2618 ind=numpy.where(dvxo<=0.1)
2619 if numpy.size(ind)>0:
2620 dvxo[ind]=0.1
2621 else:
2622 sno0=(pow0-noise0)/noise0
2623 sno1=(pow1-noise1)/noise1
2624
2625 # First parameter: Signal to noise ratio and its error
2626 sno0=numpy.log10(sno0)
2627 dsno0=1.0/numpy.sqrt(nint*navg)*(1+1/sno0)
2628 sno1=numpy.log10(sno1)
2629 dsno1=1.0/numpy.sqrt(nint*navg)*(1+1/sno1)
2630
2631 npar=6
2632 par = numpy.empty((npar, nHei)); par[:] = numpy.nan
2633
2634 if mode == 0:
2635 par[0,:] = sno
2636 par[1,:] = vxo
2637 par[2,:] = dvxo
2638 par[3,:] = vzo
2639 par[4,:] = dvzo
2640
2641 elif mode == 1 and dbs_sel == 0:
2642 par[0,:] = sno
2643 par[1,:] = vzo
2644 else:
2645 par[0,:] = sno0
2646 par[1,:] = sno1
2647 par[2,:] = dop0
2648 par[3,:] = dop1
2649 #par[4,:] = p_w
2650 #par[5,:] = p_u
2651
2652 if mode == 0:
2653 winds = numpy.vstack((par[0,:], par[1,:], par[2,:], par[3,:], par[4,:]))
2654 elif mode == 1 and dbs_sel == 0:
2655 winds = numpy.vstack((par[0,:], par[1,:]))
2656 else:
2657 winds = numpy.vstack((par[0,:], par[1,:], par[2,:], par[3,:]))
2658
2659 dataOut.data_output = winds
2660 dataOut.data_snr = par[0,:]
2661
2662 dataOut.utctimeInit = dataOut.utctime
2663 dataOut.outputInterval = dataOut.timeInterval
2664
2665 aux1= numpy.all(numpy.isnan(dataOut.data_output[0])) # NAN vectors are not written
2666 aux2= numpy.all(numpy.isnan(dataOut.data_output[1])) # NAN vectors are not written
2667 dataOut.flagNoData = aux1 or aux2
2668
2669 return dataOut
2670
2129 class SALags(Operation):
2671 class SALags(Operation):
2130 '''
2672 '''
2131 Function GetMoments()
2673 Function GetMoments()
General Comments 0
You need to be logged in to leave comments. Login now