# Iam having problem in extracting(20 to 40 dbz) and potting reflectivity to rainrate

**URL:** <https://openradar.discourse.group/t/iam-having-problem-in-extracting-20-to-40-dbz-and-potting-reflectivity-to-rainrate/468>\
**Category:** General\
**Tags:** wradlib, py-art\
**Created:** [May 7, 2024, 5:24pm UTC](https://openradar.discourse.group/t/iam-having-problem-in-extracting-20-to-40-dbz-and-potting-reflectivity-to-rainrate/468 "2024-05-07T17:24:33Z")\
**Posts on this page:** 4\
**Page:** 1

<div class="post-metadata">

**Author:** ![ankithva](https://yyz2.discourse-cdn.com/free1/user_avatar/openradar.discourse.group/ankithva/32/317_2.png) [@ankithva](https://openradar.discourse.group/u/ankithva)\
**Post date:** [May 7, 2024, 5:24pm UTC](https://openradar.discourse.group/t/iam-having-problem-in-extracting-20-to-40-dbz-and-potting-reflectivity-to-rainrate/468/1 "2024-05-07T17:24:33Z")

</div>

I need to plot rainfall(IMD) (20dbz to 40 dbz) so when iam trying to plot it values greater than 40 dbz is also coming ,i will mention code bellow:

import numpy as np  
import matplotlib.pyplot as plt  
import cartopy.crs as ccrs  
import cartopy.feature as cfeature  
import glob  
import pyart  
import numpy.ma as ma  
import wradlib as wrl  
import warnings  
warnings.filterwarnings(“ignore”)

grid\_files = glob.glob(“D:\radar\_data\radinout\gridded\_23\_9\_2022\*nc”)

total\_rain\_rate = np.zeros((len(grid\_files), 250, 250))  
radar\_lat = 26.8206  
radar\_lon = 75.8186  
radar\_radius\_km = 255

for i, file in enumerate(grid\_files):  
grid\_data = pyart.io.read\_grid(file)  
ref\_data = grid\_data.fields[‘REF’][‘data’]  
# Filter reflectivity between 20 dBZ and 40 dBZ  
ref\_values\_between\_20\_and\_40 = ma.masked\_inside(ref\_data, 20, 40)

```
grid_data.fields['REF']['data'] = ref_values_between_20_and_40
grid_data.fields['REF']['long_name'] = 'Reflectivity'
grid_data.fields['REF']['units'] = 'dBZ'

xg_files = grid_data.to_xarray()
Z = wrl.trafo.idecibel(xg_files.REF[0, 0, :, :].values)
# Assuming Z-R relation parameters a=200 and b=1.6, as in the original code
rain_rate = wrl.zr.z_to_r(Z, a=200, b=1.6)
total_rain_rate[i, :, :] = rain_rate

```

missing\_mask = np.isnan(total\_rain\_rate)  
valid\_data = np.ma.masked\_array(total\_rain\_rate, mask=missing\_mask)  
daily\_aggregated\_rain = np.sum(valid\_data, axis=0)

# Define the Plate Carree projection

proj = ccrs.PlateCarree()

# Create the axes instance

plt.figure(figsize=(8, 8))  
ax = plt.axes(projection=proj)

# Plot the contourf

im = ax.contourf(xg\_files[‘lon’].values, xg\_files[‘lat’].values, daily\_aggregated\_rain, cmap=‘jet’, origin=“lower”, extend=‘max’)

# Add colorbar

cbar = plt.colorbar(im, ax=ax, label=‘Accumulated Rainfall (mm)’)

# Set labels and title

plt.title(‘IMD JAIPUR RADAR - Accumulated Rainfall (mm) / 23-9-2022’)

# Add gridlines

gl = ax.gridlines(crs=proj, draw\_labels=True, color=‘none’)  
gl.xformatter = LONGITUDE\_FORMATTER  
gl.yformatter = LATITUDE\_FORMATTER

# Add basemap

ax.add\_feature(cfeature.COASTLINE)  
ax.add\_feature(cfeature.BORDERS, linewidth=1)  
ax.add\_feature(cfeature.STATES, linewidth=1)

# Plot the marker

ax.plot(75.8186, 26.8206, ‘k+’, color=“black”, markersize=15, markeredgewidth=2, label=‘RADAR LOCATION’)  
ax.legend(loc=‘upper right’, fontsize=‘x-small’)  
cir1 = plt.Circle((radar\_lon, radar\_lat), radar\_radius\_km / 111.0, edgecolor=‘black’, facecolor=‘none’, ls=‘–’,  
label=‘RADAR RANGE(250km)’)  
ax.add\_patch(cir1)

# Show the plot

plt.show()

---

<div class="post-metadata">

**Author:** ![mgrover1](https://yyz2.discourse-cdn.com/free1/user_avatar/openradar.discourse.group/mgrover1/32/7_2.png) [@mgrover1](https://openradar.discourse.group/u/mgrover1)\
**Post date:** [May 10, 2024, 1:45pm UTC](https://openradar.discourse.group/t/iam-having-problem-in-extracting-20-to-40-dbz-and-potting-reflectivity-to-rainrate/468/2 "2024-05-10T13:45:51Z")

</div>

Is there a question? Could you please detail more about what you are asking for? We can respond to feedback about specific issues, but cannot fully debug scripts.

---

<div class="post-metadata">

**Author:** ![ankithva](https://yyz2.discourse-cdn.com/free1/user_avatar/openradar.discourse.group/ankithva/32/317_2.png) [@ankithva](https://openradar.discourse.group/u/ankithva)\
**Post date:** [May 12, 2024, 4:35am UTC](https://openradar.discourse.group/t/iam-having-problem-in-extracting-20-to-40-dbz-and-potting-reflectivity-to-rainrate/468/3 "2024-05-12T04:35:02Z")

</div>

Iam not asking for scripts ,The initial problem is mentioned at very beginning of the post and the code is only for a reference. Sorry for the misunderstanding.

---

<div class="post-metadata">

**Author:** ![jorahu](https://yyz2.discourse-cdn.com/free1/user_avatar/openradar.discourse.group/jorahu/32/19_2.png) [@jorahu](https://openradar.discourse.group/u/jorahu)\
**Post date:** [May 12, 2024, 11:54am UTC](https://openradar.discourse.group/t/iam-having-problem-in-extracting-20-to-40-dbz-and-potting-reflectivity-to-rainrate/468/4 "2024-05-12T11:54:56Z")

</div>

Hey

Seems like you need to debug your code line by line.

The following line in your code should do exactly what you want.

```auto
# Filter reflectivity between 20 dBZ and 40 dBZ
ref_values_between_20_and_40 = ma.masked_inside(ref_data, 20, 40)

```

Why is does not can depend on multiple problems. Maybe you have overwritten some variables before plotting (an eeeeeeasy mistake when using Jupyter or a similar IDE). Print out or plot the variable ref\_values\_between\_20\_and\_40 right after the definition with, for example, plt.imshow(ref\_values\_between\_20\_and\_40).

There is no way for others to debug your code by simply looking at it. Even more so when the functions used seem appropriate.

If the ma.masked\_inside doesn’t do the expected, try numpy.where(), for example.

Also, if the issue is strictly with filtering your data, there is no need to attach code about plotting. Always include the minimum amount of code needed to illustrate your issue, this tremendously increases your chances of getting informative answers.