%pylab inline

import pydisdrometer
from matplotlib.colors import LogNorm

filename = '/Users/hard505/Dropbox/Projects/NetworkRetrieval/data/Disdrometer/IFloodS_APU01_2013_0502_dsd.txt' #Parsivel 05, March 2nd

dsd = pydisdrometer.read_parsivel_nasa_gv(filename, campaign='ifloods')

"""
Explanation: Quick usage examples for PyDisdrometer
This notebook shows some (very) brief examples of current PyDisdrometer functionality and how to interact with it.
NASA GV Field Campaign Parsivel Format
Let's start with a NASA Ground Validation Parsivel formatted file.
End of explanation
"""


dsd.calculate_radar_parameters()  

"""
Explanation: So at this point we have the drop size distribution read in. NASA strips out rainfall information though.  Let's do some T-Matrix scattering though. This should take a little bit(up to a minute or so on my laptop depending on the file).
End of explanation
"""


plt.figure(figsize=(12,12))

plt.subplot(2,2,1)
plt.plot(dsd.time/60.0, dsd.fields['Zh']['data'])
plt.xlabel('Time(hrs)')
plt.ylabel('Reflectivity(dBZ)')
plt.xlim(5,24)

plt.subplot(2,2,2)
plt.plot(dsd.time/60.0, dsd.fields['Zdr']['data'])
plt.xlabel('Time(hrs)')
plt.ylabel('Differential Reflectivity(dB)')
plt.xlim(5,24)

plt.subplot(2,2,3)
plt.plot(dsd.time/60.0, dsd.fields['Kdp']['data'])
plt.xlabel('Time(hrs)')
plt.ylabel('Specific Differential Phase(deg/km)')
plt.xlim(5,24)

plt.subplot(2,2,4)
plt.pcolor(dsd.time/60.0, dsd.diameter, np.log10(dsd.Nd.T), vmin=0, vmax=6)
plt.xlabel('Time(hrs)')
plt.ylabel('Diameter')
plt.colorbar()
plt.ylim(0,5) #Zoom in some
plt.xlim(5,24)

plt.tight_layout()

plt.show()

"""
Explanation: This assumes BC shape relationship, X band, 10C. You can pass in a new wavelength to change that.
 Let's plot some of the parameters, and then try to do something with the data.
End of explanation
"""


dsd.calculate_dsd_parameterization()

plt.figure(figsize=(12,12))

plt.subplot(2,2,1)
plt.plot(dsd.time/60.0, dsd.fields['D0']['data'])
plt.xlabel('Time(hrs)')
plt.ylabel('$D_0$')
plt.xlim(5,24)

plt.subplot(2,2,2)
plt.plot(dsd.time/60.0, np.log10(dsd.fields['Nw']['data']))
plt.xlabel('Time(hrs)')
plt.ylabel('$log_{10}(N_w)$')
plt.xlim(5,24)

plt.subplot(2,2,3)
plt.plot(dsd.time/60.0, dsd.fields['Dmax']['data'])
plt.xlabel('Time(hrs)')
plt.ylabel('Shape parameters')
plt.xlim(5,24)

plt.subplot(2,2,4)
plt.pcolor(dsd.time/60.0, dsd.diameter, np.log10(dsd.Nd.T), vmin=0, vmax=6)
plt.xlabel('Time(hrs)')
plt.ylabel('Dmax')
plt.colorbar()
plt.ylim(0,5) #Zoom in some
plt.xlim(5,24)

plt.tight_layout()

plt.show()

"""
Explanation: Next let's estimate some microphysical parameterizations.
End of explanation
"""


plt.figure(figsize=(12,12))

plt.subplot(2,2,1)
plt.hist2d(dsd.fields['D0']['data'], np.log10(dsd.fields['Nw']['data']), 
           range=((0,6),(0,6)), vmin=1, bins=50, norm=LogNorm())
plt.ylabel('log10(Nw)')
plt.xlabel('$D_0$')

plt.subplot(2,2,2)
plt.hist2d(dsd.fields['Zh']['data'], dsd.fields['Zdr']['data'], 
           range=((-10,40),(-1,1)), vmin=1, bins=50, norm=LogNorm())
plt.xlabel('$Z_h$')
plt.ylabel('$Z_{dr}$')

plt.subplot(2,2,3)
plt.hist2d(dsd.fields['Dm']['data'], dsd.fields['Zdr']['data'],
        range=((0,2),(-1,1)), vmin=1, bins=50, norm=LogNorm())
plt.xlabel('$D_m$')
plt.ylabel('$Z_{dr}$')

plt.subplot(2,2,4)
plt.hist2d(dsd.fields['Dmax']['data'], dsd.fields['Zh']['data'],
        range=((0,8),(-10,40)), vmin=1, bins=50, norm=LogNorm())
plt.ylabel('$Z_h$')
plt.xlabel('$Dmax$')

plt.tight_layout()

plt.show()

"""
Explanation: At hour twenty we see some drops that are probably instrument inaccuracies. Let's look at the overall behavior of the parameters with respect to each other.
End of explanation
"""


plot(dsd.time/60.0 , 0<sum(dsd.Nd[:,dsd.diameter>4], axis=1))
plt.xlim(5,24)

"""
Explanation: Let's take what might be another common use case. Finding areas with certain size drops. In this case everything above 4mm. There are a few ways to do this but let's just create an indicator function for now.
End of explanation
"""


time_hrs = dsd.time/60.0
print(time_hrs[0<sum(dsd.Nd[:,dsd.diameter>4], axis=1)])

"""
Explanation: Or maybe we just want a listing of the time steps where there are drops greater than 4mm
End of explanation
"""


