How can i add new DSR function

[yogeshkkolte]

I have calculated new axis ratio for our site and want to use that for calculate radar parameters so to test performance of new equation i replace new equation in DSR.py module in the def bc(D_eq) function. But when i run command "dsd.calculate_radar_parameters" then python kernel is dying def bc(D_eq): '''Beard and Chuang Drop Shape relationship model.

Implementation of the Beard and Chuang Drop Shape model given in [1]
. This gives the ratio of the major to minor axis.

Parameters
-----------
D_eq: float
    Volume Equivalent Drop Diameter

Returns
-------
axis_ratio: float
    The ratio of the major to minor axis.

See Also
--------
tb : Thurai and Bringi DSR
pb : Pruppacher and Beard DSR

References
----------
..[1]  Beard, Kenneth V., Catherine Chuang, 1987: A New Model for the
Equilibrium Shape of Raindrops. J. Atmos. Sci., 44, 1509-1524.
'''

return 1.0048 + 5.7e-04 * np.power(D_eq, 1) - \

2.628e-02 * np.power(D_eq, 2) + 3.682e-03 * np.power(D_eq, 3) - \

1.677e-04 * np.power(D_eq, 4)

Monsoon drop shape relationship for 0 < Deq <8

return 0.1479 + 2.413_(Deq) - 2.194(D_eq*2) + 0.7531(Deq**3) - 0.08305(Deq*4)

Monsoon drop shape relationship for 1 < Deq <8

return 0.9946 + 0.02297 * D_eq - 0.02051 * D_eq ** 2 + 0.002795 * D_eq ** 3 -\
    0.0001849 * D_eq ** 4

[jleinonen]

Can you formulate the code properly for viewing? The first monsoon relationship doesn't look like valid Python. Also, which frequency are you using? The scattering module (which I wrote around some pre-existing numerical code) may crash for some higher frequencies / aspect ratios.

[josephhardinee]

With the current formatting it is really hard to read what you are trying to do. It appears as if there are 3 return statements in one function? Can you provide an example of what you are doing in a self contained file so we can run it and see if we can figure out what is going on?

[yogeshkkolte]

I want to use this equation in DSR module
b/a = 0.9857 + 4.03e-02 * np.power(D_eq, 1) - 3.25e-02np.power(D_eq, 2) + \ 6.2e-03 * np.power(D_eq, 3) - 5.1e-04 \ np.power(D_eq, 4)
I am calculating radar parameters for four different frequency, so in have changed wavelengths in tmatrix_aux.py in pytmatrix module . previously i was getting radar parameters by changing wavelengths for Beard and Chuang equation.

typical wavelengths [mm] at different bands

wl_S = 111.0

wl_S = 100.0

wl_C = 53.5

wl_C = 50.0

wl_X = 33.3

wl_X = 31.46 wl_Ku = 22.0

wl_Ka = 8.43

wl_Ka = 8.50

Is it necessary to make change in "def dsr_bc(D_eq)" axis ratio function present in tmatrix_aux.py ? I tried that also but unexpectedly kernel is dying. I am not getting any error massage.

[josephhardinee]

Yogesh, In the current incarnation or pydisdrometer, it is limited to the wavelengths that were built into pytmatrix(Which we use for our scattering). In a newer branch release Nick and I are currently working on, this restriction will be relaxed.

To pass in your own DSR, you would define the relationship as you have above, say

def yogesh_dsr(D):
    return  0.9857 + 4.03e-02 * np.power(D_eq, 1) - 3.25e-02*np.power(D_eq, 2) + \
    6.2e-03 * np.power(D_eq, 3) - 5.1e-04 * np.power(D_eq, 4)

Then when issuing the call to calculate_radar_parameters, you can pass in dsr_func=yogesh_dsr So for instance, something like

dsd.calculate_radar_parameters(wavelength=tmatrix_aux.wl_X, dsr_func =yogesh_dsr)

Hopefully that helps. Until the next release though, I would stick with the built in radar frequencies. It does not make a particularly large difference for the most part if you deviate by a GHz here or there, especially not for practical purposes.

Also just as an FYI, I plotted that DSR real quick to compare it to the Pruppacher and Beard, and Thurai and Bringi relations and it differs by quite a bit, so I would be careful with interpreting the results and make sure that is what you intended.

[yogeshkkolte]

Hello Sir, The coed is running for default wavelengths present in tmatrix_aux and those wave length i have been added in tmatrix_aux for all Drop Shape relationship models present in DSR.py except yogesh_dsr. I dint understand why this is happening. I am having 12Gb RAM in computer.I thought that calculations need large apace so i run program on HPC(High-performance computing) system still getting same error.

dsd.change_scattering_temperature("20C") dsd.calculate_radar_parameters(wavelength=tmatrix_aux.wl_X, dsr_func = DSR.yogesh_dsr) Zhs = dsd.fields['Zh']['data'] Ais = dsd.fields['Ai']['data'] please help me.

[josephhardinee]

By error you mean the python kernel is dying correct?(It's dropping you back to a shell in other words if doing this on the command line). As Jussi pointed out, this happens when drops that are very large relative to the wavelength are passed in, which is a deficiency of the original Fortran scattering code, or drops with extreme oblateness. What is the largest entry in dsd.diameter, and the output of DSR.yogesh_dsr for that diameter?

[yogeshkkolte]

yes it is dropping me out of ipython console. I fitted 4th order equation for 1mm <= D <7mm.

[josephhardinee]

Can you put the output of dsd.diameter in the comments box? At X-Band you should not be seeing pytmatrix crash for those drop sizes so I'm guessing maybe a large diameter snuck in or something like that.

[yogeshkkolte]

DSR.yogesh_dsr not at all working for any wavelength.

dsd.change_scattering_temperature("20C") dsd.calculate_radar_parameters(wavelength=tmatrix_aux.wl_X, dsr_func = DSR.bc)

dsd.diameter Out[2]: array([ 0.2, 0.4, 0.6, 0.8, 1. , 1.2, 1.4, 1.6, 1.8,

2. , 2.2, 2.4, 2.6, 2.8, 3. , 3.2, 3.4, 3.6, 3.8, 4. , 4.2, 4.4, 4.6, 4.8, 5. , 5.2, 5.4, 5.6, 5.8, 6. , 6.2, 6.4, 6.6, 6.8, 7. , 7.2, 7.4, 7.6, 7.8, 8. , 8.2, 8.4, 8.6, 8.8, 9. , 9.2, 9.4, 9.6, 9.8, 10. ])

[josephhardinee]

Can you upload a small self contained and runnable example that I can run(So a .py and a data file) that reproduces this. I think that will make this easier. Also can you post the output of

pip freeze | grep PyDisdrometer

My guess is that your DSR goes very oblate at 10mm, which is causing the generated scattering tables to have issues.

[yogeshkkolte]

PyDisdrometer==0.1.13.0

[josephhardinee]

Actually at 10mm your oblateness ratio is negative which is physically nonsensical and causing the code to crash. I would consider something like an if statement for whether the drop is greater than 7mm, and if so just returning the axis ratio of a 7mm drop as a first approximation.

[yogeshkkolte]

In our observations drop diameter is not grater than 7mm . less than 100 drops present in 6-7mm range out of 3 years monsoon observations. I am so sorry for by mistake i written wrong equation and successive conversations i did not recognized my mistakes.

correct equation : def dsr_yogesh(D_eq): return 0.9857 + 4.03e-02 * np.power(D_eq, 1) - 3.25e-02np.power(D_eq, 2) + \ 6.2e-03 * np.power(D_eq, 3) - 5.1e-04 \ np.power(D_eq, 4)

How can i limit T-Matrix simulation within 0 to 7mm? please forgive me.

[josephhardinee]

It's no problem. The easiest thing here is probably to just make sure your DSR does not go negative. The quickest way to do this is to just cap a maximum oblateness ratio. So something like

def dsr_yogesh(D_eq):
    ratio= 0.9857 + 4.03e-02 * np.power(D_eq, 1) - 3.25e-02*np.power(D_eq, 2) + \
    6.2e-03 * np.power(D_eq, 3) - 5.1e-04 * np.power(D_eq, 4) 
    ratio[ratio<0.577] =0.577 #dsr_yogesh(7)
    return ratio

As you don't have any drops larger than 7, the capping should not matter.

The second approach would be to delete the relevant entries from dsd.diameter and delete the relevant columns from Nd. This is a little more challenging. A third approach is to limit what goes into the scattering code. While the most elegant and correct, it requires some changes within PyDisdrometer. We'll throw this on a feature request list as this would help handling the 2DVD's which have sample locations up to much larger sizes.

[yogeshkkolte]

I was getting following error ''numpy.float64' object does not support item assignment'' so tried this: def dsr_yogesh(D_eq): if (D_eq <= 7.0): ratio= 0.9857 + 4.03e-02 * np.power(D_eq, 1) - 3.25e-02np.power(D_eq, 2) + \ 6.2e-03 * np.power(D_eq, 3) - 5.1e-04 \ np.power(D_eq, 4) else: ratio = 0.577 return ratio

[josephhardinee]

Is that working for you now?

[yogeshkkolte]

yes it is working. when i compare simulated reflectivity between dsr_yogesh and dsr_bc not seen any difference in reflectivity. plot for both function is linear.

[josephhardinee]

That seems surprising. Can you paste the code sample you are using to compare the two?

[yogeshkkolte]

import matplotlib.pylab as plt import pydisdrometer as pyd from pytmatrix import tmatrix_aux from pydisdrometer import DSR filename = '/home/user/anaconda/DSD/Dsd_dbl/DSD_CSV/V12281_1.csv' reader = pyd.aux_readers.HACPL_2DVD_reader.HACPL_2DVD_dsd_reader(filename,'hacpl') dsd = pyd.aux_readers.HACPL_2DVD_reader.read_2dvd_dsd_hacpl(filename,'hacpl') dsd.change_scattering_temperature("20C") dsd.calculate_radar_parameters(wavelength=tmatrix_aux.wl_X, dsr_func = DSR.dsr_yogesh) Zhsy = dsd.fields['Zh']['data'] Aisy = dsd.fields['Ai']['data']

dsd.change_scattering_temperature("20C") dsd.calculate_radar_parameters(wavelength=tmatrix_aux.wl_X, dsr_func = DSR.bc) Zhs = dsd.fields['Zh']['data'] Ais = dsd.fields['Ai']['data'] plt.plot(Zhs,Zhsy,'r-',label = 'yogesh - BC',linewidth = 2)

plt.legend(loc = 9)

[josephhardinee]

Try replacing the zhsy = dsd.fields['Zh]['data'].copy() and see if you get the same results.

Sent from my iPhone

On Aug 15, 2016, at 10:53 PM, yogesh kolte notifications@github.com wrote:

import matplotlib.pylab as plt import pydisdrometer as pyd from pytmatrix import tmatrix_aux from pydisdrometer import DSR filename = '/home/user/anaconda/DSD/Dsd_dbl/DSD_CSV/V12281_1.csv' reader = pyd.aux_readers.HACPL_2DVD_reader.HACPL_2DVD_dsd_reader(filename,'hacpl') dsd = pyd.aux_readers.HACPL_2DVD_reader.read_2dvd_dsd_hacpl(filename,'hacpl') dsd.change_scattering_temperature("20C") dsd.calculate_radar_parameters(wavelength=tmatrix_aux.wl_X, dsr_func = DSR.dsr_yogesh) Zhsy = dsd.fields['Zh']['data'] Aisy = dsd.fields['Ai']['data']

dsd.change_scattering_temperature("20C") dsd.calculate_radar_parameters(wavelength=tmatrix_aux.wl_X, dsr_func = DSR.bc) Zhs = dsd.fields['Zh']['data'] Ais = dsd.fields['Ai']['data'] plt.plot(Zhs,Zhsy,'r-',label = 'yogesh - BC',linewidth = 2)

plt.legend(loc = 9)

— You are receiving this because you commented. Reply to this email directly, view it on GitHub, or mute the thread.

[yogeshkkolte]

Getting same result: dsd.change_scattering_temperature("20C") dsd.calculate_radar_parameters(wavelength=tmatrix_aux.wl_X, dsr_func = DSR.dsr_yogesh) Aisy = dsd.fields['Ai']['data'] zhsy = dsd.fields['Zh']['data'].copy() dsd.change_scattering_temperature("20C") dsd.calculate_radar_parameters(wavelength=tmatrix_aux.wl_X, dsr_func = DSR.bc) Zhs = dsd.fields['Zh']['data'] Ais = dsd.fields['Ai']['data']

plt.plot(Zhs,zhsy,'r-',label = 'yogesh - BC',linewidth = 2)

plt.plot(Zhsy,'b-',label = 'yogesh',linewidth = 2)

plt.legend(loc = 9) plt.xlabel('BC-Zh',fontsize = 20) plt.ylabel('yogesh-Zh',fontsize = 20)

I verified DSR.py def dsr_yogesh(D_eq): if (D_eq <= 7.0): ratio= 0.9857 + 4.03e-02 * np.power(D_eq, 1) - 3.25e-02np.power(D_eq, 2) + \ 6.2e-03 * np.power(D_eq, 3) - 5.1e-04 \ np.power(D_eq, 4) else: ratio = 0.577 return ratio

[josephhardinee]

Can you try comparing Zdr. It should definitely be different. There was a bug a while back where the DSR was not getting changed properly and was defaulting to BC. This was fixed and I want to make sure that is not present in your version.

[yogeshkkolte]

please go through figure . please suggest any correction. could i share the csv file of input date on gmail.

[yogeshkkolte]

Sorry Joseph, Again I made mistake in observing figures.The dsr_yogesh under estimating ZDR Values. Here X axis is time index

[josephhardinee]

Okay that looks more like what I expected based on the DSR. It had much less oblateness in the middle size drops if I remember right(Which can have a strong effect on Zdr) and so Zdr would be suppressed. I'm going to close this issue as I believe everything is working fine at this point(with PyDisdrometer at least). If you run into any other issues with PyDisdrometer feel free to open another issue.

[yogeshkkolte]

Thank you for your valuable support.