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