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

General
,
1

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()

2

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.

1 Like
3

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.

4

Hey

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

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

# 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.

2 Likes