ValueError: could not broadcast input array from shape (40,250,250) into shape (250,250)

General
, ,
1 
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
import numpy as np
import warnings
import wradlib as wrl

warnings.filterwarnings("ignore")

# Assuming grid_files is defined as in your code
grid_files = glob.glob("D:\\radar_data\\radinout\\gridded_20_jun_2022\\gridded_20_jun_2022grid_ZDR_*nc")

radar_lon= 75.8175
radar_radius_km=280.0
cir1=plt.Circle((radar_lon,radar_lat),radar_radius_km/111.0,edgecolor='black',facecolor='none')
cir2=plt.Circle((radar_lon,radar_lat),radar_radius_km/111.0,edgecolor='black',facecolor='none')




rain_rate_stratiform = np.zeros((len(grid_files), 250, 250))
rain_rate_convective= np.zeros((len(grid_files), 250, 250))



for i, file in enumerate(grid_files):
    grid_data = pyart.io.read_grid(file)
    ref_data = grid_data.fields['REF']['data']

    xg_files = grid_data.to_xarray()
    Z = wrl.trafo.idecibel(xg_files.REF[0, 0, :, :].values)
    #rain_rate = wrl.zr.z_to_r(Z, a=200, b=1.5)
    #total_rain_rate[i, :, :] = rain_rate

    grid_data = pyart.io.read_grid(grid_file)
    xg_files = grid_data.to_xarray()
    ref_data = grid_data.fields['REF']['data']
    feature_dict = pyart.retrieve.feature_detection(
        grid_data,
        dx=2000,
        dy=2000,
        level_m=500,
        field="REF",
        always_core_thres=30,
        bkg_rad_km=20,
        use_cosine=True,
        max_diff=5,
        zero_diff_cos_val=55,
        weak_echo_thres=10,
        max_rad_km=4,
    )

    # Reshape the feature_detection field data
    feature_data = feature_dict["feature_detection"]["data"][0, :, :]
    feature_data_reshaped = np.broadcast_to(feature_data, (grid_data.nz,) + feature_data.shape)

    # Mask the feature data
    feature_masked = np.ma.masked_equal(feature_data_reshaped, 0)
    feature_masked = np.ma.masked_equal(feature_masked, 3)

    # Add the masked data to the grid
    grid_data.add_field(
        "feature_detection", {"data": feature_masked}, replace_existing=True
    )

    # Masked array for stratiform and convective reflectivity
    stratiform_reflectivity = ma.masked_where(feature_masked != 1, ref_data)
    convective_reflectivity = ma.masked_where(feature_masked != 2, ref_data)

    # Convert reflectivity to rainfall rate using wradlib
    Zs_stratiform = wrl.trafo.idecibel(stratiform_reflectivity)
    Zs_convective = wrl.trafo.idecibel(convective_reflectivity)

    rain_rate_stratiform = wrl.zr.z_to_r(Zs_stratiform, a=251, b=1.48)
    rain_rate_convective = wrl.zr.z_to_r(Zs_convective, a=178, b=1.51)

    # Apply mask to rainfall rates
    rain_rate_stratiform[i, :, :] = np.ma.masked_greater(rain_rate_stratiform, 20)
    rain_rate_convective[i, :, :] = np.ma.masked_greater(rain_rate_convective, 20)

    # Define the Plate Carree projection
    proj = ccrs.PlateCarree()

    # Create the figure and axes instances
    fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(15, 5), subplot_kw={'projection': proj})

    # Plot the first contourf
    im1 = ax1.contourf(xg_files['lon'].values, xg_files['lat'].values, rain_rate_stratiform, 
                   levels=np.linspace(0, 20, 10), cmap='jet', origin="lower")
    # Plot the marker
    ax1.plot(75.8175, 26.8207, marker="x", color="black", markersize=15, markeredgewidth=2)
    # Add colorbar
    cbar1 = plt.colorbar(im1, ax=ax1, label='Rainfall rate in mm/h', ticks=np.arange(0, 21, 2))
    cbar1.ax.set_yticklabels(np.arange(0, 21, 2))
    # Add title
    ax1.set_title("Stratiform Rainfall")
    ax1.add_patch(cir1)
    # Add gridlines
    gl1 = ax1.gridlines(crs=proj, draw_labels=True, color='none')
    gl1.xformatter = LONGITUDE_FORMATTER
    gl1.yformatter = LATITUDE_FORMATTER
    # Add basemap
    ax1.add_feature(cfeature.COASTLINE)
    ax1.add_feature(cfeature.BORDERS, linewidth=1)
    ax1.add_feature(cfeature.STATES, linewidth=1)


    # Plot the second contourf
    im2 = ax2.contourf(xg_files['lon'].values, xg_files['lat'].values, rain_rate_convective, 
                       levels=np.linspace(0, 20, 10), cmap='jet', origin="lower")
    # Plot the marker
    ax2.plot(75.8175, 26.8207, marker="x", color="black", markersize=15, markeredgewidth=2)
    # Add colorbar
    cbar2 = plt.colorbar(im2, ax=ax2, label='Rainfall rate in mm/h', ticks=np.arange(0, 21, 2))
    cbar2.ax.set_yticklabels(np.arange(0, 21, 2))
    ax2.add_patch(cir2)
    # Add title
    ax2.set_title("Convective Rainfall")

    # Add gridlines
    gl2 = ax2.gridlines(crs=proj, draw_labels=True, color='none')
    gl2.xformatter = LONGITUDE_FORMATTER
    gl2.yformatter = LATITUDE_FORMATTER
    # Add basemap
    ax2.add_feature(cfeature.COASTLINE)
    ax2.add_feature(cfeature.BORDERS, linewidth=1)
    ax2.add_feature(cfeature.STATES, linewidth=1)

    # Show the plot
    plt.show()

this showing error:

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[178], line 79
     76 rain_rate_convective = wrl.zr.z_to_r(Zs_convective, a=178, b=1.51)
     78 # Apply mask to rainfall rates
---> 79 rain_rate_stratiform[i, :, :] = np.ma.masked_greater(rain_rate_stratiform, 20)
     80 rain_rate_convective[i, :, :] = np.ma.masked_greater(rain_rate_convective, 20)
     82 # Define the Plate Carree projection

File ~\PYTHON\envs\pcf\lib\contextlib.py:79, in ContextDecorator.__call__.<locals>.inner(*args, **kwds)
     76 @wraps(func)
     77 def inner(*args, **kwds):
     78     with self._recreate_cm():
---> 79         return func(*args, **kwds)

File ~\PYTHON\envs\pcf\lib\site-packages\numpy\ma\core.py:3397, in MaskedArray.__setitem__(self, indx, value)
   3395         _data[indx.data] = dval
   3396     else:
-> 3397         _data[indx] = dval
   3398         _mask[indx] = mval
   3399 elif hasattr(indx, 'dtype') and (indx.dtype == MaskType):

ValueError: could not broadcast input array from shape (40,250,250) into shape (250,250)
2 
rain_rate_stratiform = np.ma.masked_greater(rain_rate_stratiform, 20)

The error message is already very clear:
ValueError: could not broadcast input array from shape (40,250,250) into shape (250,250)

np.ma.masked_greater(rain_rate_stratiform, 20) is of shape (40, 250, 250) and rain_rate_stratiform[i, :, :] is of shape (250, 250) as there is the i-subset. Removing this subset will fix this error.

I know that sometimes Python error messages are hard to decipher. Please always try to get behind the error message. In many cases you will immediately (OK, you might struggle a while :grimacing:) see the issue solution.

Best,
Kai

1 Like
3

Thanks i tried that but the loop/code is still running without end and outputs(plots)

4

My above solution does fix the error, but it might not do the expected.

I’ve added missing markdown to your OP to correctly display the code.

To better get an idea what’s going on it would be great to get information about the shapes of the individual arrays within the loop.

Looking at this now it seems the issue starts here:

What’s the shape of ref_data, feature_data feature_data_reshaped, feature_masked?

I’m assuming something is going wrong when doing the above quoted assignments as you are overwriting (at least partly, already defined arrays.

HTH,
Kai

5 
ref_data.shape=(40, 250, 250)
feature_data.shape=(250, 250)
feature_data_reshaped.shape=(40, 250, 250)
feature_masked.shape=(40, 250, 250)
6

So, the question is, do you want to have a 3D output (z,y,x)? Then your rain_rate_stratiform would need to be 4D (files, z, y, x). If you need (files, y, x), then you need to go to 2D here:

Does that make sense?

Kai

7

all i need to plot is the rain_rate_stratiform and rain_rate_convective, so how can i reshape the arrays to plot the result?

8

iam getting the output for a single time step(with out loop) like this:

import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
import numpy as np
import warnings
import wradlib as wrl

warnings.filterwarnings("ignore")

grid_file = "D:\\radar_data\\radinout\\gridded_20_jun_2022\\gridded_20_jun_2022grid_ZDR_JPR220620000253-IMD-B.nc"
radar_lat= 26.8207
radar_lon= 75.8175
radar_radius_km=280.0
cir1=plt.Circle((radar_lon,radar_lat),radar_radius_km/111.0,edgecolor='black',facecolor='none')
cir2=plt.Circle((radar_lon,radar_lat),radar_radius_km/111.0,edgecolor='black',facecolor='none')

stotal_rain_rate = np.zeros((len(grid_files), 250, 250))
ctotal_rain_rate = np.zeros((len(grid_files), 250, 250))
grid_data = pyart.io.read_grid(grid_file)
xg_files = grid_data.to_xarray()
ref_data = grid_data.fields['REF']['data']
feature_dict = pyart.retrieve.feature_detection(
    grid_data,
    dx=2000,
    dy=2000,
    level_m=500,
    field="REF",
    always_core_thres=30,
    bkg_rad_km=20,
    use_cosine=True,
    max_diff=5,
    zero_diff_cos_val=55,
    weak_echo_thres=10,
    max_rad_km=4,
)

# Reshape the feature_detection field data
feature_data = feature_dict["feature_detection"]["data"][0, :, :]
feature_data_reshaped = np.broadcast_to(feature_data, (grid_data.nz,) + feature_data.shape)

# Mask the feature data
feature_masked = np.ma.masked_equal(feature_data_reshaped, 0)
feature_masked = np.ma.masked_equal(feature_masked, 3)

# Add the masked data to the grid
grid_data.add_field(
    "feature_detection", {"data": feature_masked}, replace_existing=True
)

# Masked array for stratiform and convective reflectivity
stratiform_reflectivity = ma.masked_where(feature_masked != 1, ref_data)
convective_reflectivity = ma.masked_where(feature_masked != 2, ref_data)

# Convert reflectivity to rainfall rate using wradlib
Zs_stratiform = wrl.trafo.idecibel(stratiform_reflectivity)
Zs_convective = wrl.trafo.idecibel(convective_reflectivity)

rain_rate_stratiform = wrl.zr.z_to_r(Zs_stratiform, a=251, b=1.48)
rain_rate_convective = wrl.zr.z_to_r(Zs_convective, a=178, b=1.51)

# Define the Plate Carree projection
proj = ccrs.PlateCarree()

# Create the figure and axes instances
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(15, 5), subplot_kw={'projection': proj})

# Plot the first contourf
im1 = ax1.contourf(xg_files['lon'].values, xg_files['lat'].values, rain_rate_stratiform[1], 
                   levels=np.linspace(0, 100, 10), cmap='jet', origin="lower")
# Plot the marker
ax1.plot(75.8175, 26.8207, marker="x", color="black", markersize=15, markeredgewidth=2)
# Add colorbar
cbar1 = plt.colorbar(im1, ax=ax1, label='Rainfall rate in mm/h', ticks=np.arange(0, 100, 10))
cbar1.ax.set_yticklabels(np.arange(0, 100, 10))
# Add title
ax1.set_title("Stratiform Rainfall")
ax1.add_patch(cir1)
# Add gridlines
gl1 = ax1.gridlines(crs=proj, draw_labels=True, color='none')
gl1.xformatter = LONGITUDE_FORMATTER
gl1.yformatter = LATITUDE_FORMATTER
# Add basemap
ax1.add_feature(cfeature.COASTLINE)
ax1.add_feature(cfeature.BORDERS, linewidth=1)
ax1.add_feature(cfeature.STATES, linewidth=1)


# Plot the second contourf
im2 = ax2.contourf(xg_files['lon'].values, xg_files['lat'].values, rain_rate_convective[1], 
                   levels=np.linspace(0, 100, 10), cmap='jet', origin="lower")
# Plot the marker
ax2.plot(75.8175, 26.8207, marker="x", color="black", markersize=15, markeredgewidth=2)
# Add colorbar
cbar2 = plt.colorbar(im2, ax=ax2, label='Rainfall rate in mm/h', ticks=np.arange(0, 100, 10))
cbar2.ax.set_yticklabels(np.arange(0, 100, 10))
ax2.add_patch(cir2)
# Add title
ax2.set_title("Convective Rainfall")

# Add gridlines
gl2 = ax2.gridlines(crs=proj, draw_labels=True, color='none')
gl2.xformatter = LONGITUDE_FORMATTER
gl2.yformatter = LATITUDE_FORMATTER
# Add basemap
ax2.add_feature(cfeature.COASTLINE)
ax2.add_feature(cfeature.BORDERS, linewidth=1)
ax2.add_feature(cfeature.STATES, linewidth=1)

# Show the plot
plt.show()

AND HOW TO APPLY IT IN A LOOP TO PLOT ACCUMULATED /RAIN RATE?

download
download1142×467 57.3 KB

9

@ankithva

If you add code, please wrap it like this (3 backticks):

   ```python
   # code goes here
   ```
1 Like
10 
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
import numpy as np
import warnings
import wradlib as wrl

warnings.filterwarnings("ignore")

grid_file = "D:\\radar_data\\radinout\\gridded_20_jun_2022\\gridded_20_jun_2022grid_ZDR_JPR220620000253-IMD-B.nc"
radar_lat= 26.8207
radar_lon= 75.8175
radar_radius_km=280.0
cir1=plt.Circle((radar_lon,radar_lat),radar_radius_km/111.0,edgecolor='black',facecolor='none')
cir2=plt.Circle((radar_lon,radar_lat),radar_radius_km/111.0,edgecolor='black',facecolor='none')

stotal_rain_rate = np.zeros((len(grid_files), 250, 250))
ctotal_rain_rate = np.zeros((len(grid_files), 250, 250))
grid_data = pyart.io.read_grid(grid_file)
xg_files = grid_data.to_xarray()
ref_data = grid_data.fields['REF']['data']
feature_dict = pyart.retrieve.feature_detection(
  grid_data,
  dx=2000,
  dy=2000,
  level_m=500,
  field="REF",
  always_core_thres=30,
  bkg_rad_km=20,
  use_cosine=True,
  max_diff=5,
  zero_diff_cos_val=55,
  weak_echo_thres=10,
  max_rad_km=4,
)

# Reshape the feature_detection field data
feature_data = feature_dict["feature_detection"]["data"][0, :, :]
feature_data_reshaped = np.broadcast_to(feature_data, (grid_data.nz,) + feature_data.shape)

# Mask the feature data
feature_masked = np.ma.masked_equal(feature_data_reshaped, 0)
feature_masked = np.ma.masked_equal(feature_masked, 3)

# Add the masked data to the grid
grid_data.add_field(
  "feature_detection", {"data": feature_masked}, replace_existing=True
)

# Masked array for stratiform and convective reflectivity
stratiform_reflectivity = ma.masked_where(feature_masked != 1, ref_data)
convective_reflectivity = ma.masked_where(feature_masked != 2, ref_data)

# Convert reflectivity to rainfall rate using wradlib
Zs_stratiform = wrl.trafo.idecibel(stratiform_reflectivity)
Zs_convective = wrl.trafo.idecibel(convective_reflectivity)

rain_rate_stratiform = wrl.zr.z_to_r(Zs_stratiform, a=251, b=1.48)
rain_rate_convective = wrl.zr.z_to_r(Zs_convective, a=178, b=1.51)

# Define the Plate Carree projection
proj = ccrs.PlateCarree()

# Create the figure and axes instances
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(15, 5), subplot_kw={'projection': proj})

# Plot the first contourf
im1 = ax1.contourf(xg_files['lon'].values, xg_files['lat'].values, rain_rate_stratiform[1], 
                 levels=np.linspace(0, 100, 10), cmap='jet', origin="lower")
# Plot the marker
ax1.plot(75.8175, 26.8207, marker="x", color="black", markersize=15, markeredgewidth=2)
# Add colorbar
cbar1 = plt.colorbar(im1, ax=ax1, label='Rainfall rate in mm/h', ticks=np.arange(0, 100, 10))
cbar1.ax.set_yticklabels(np.arange(0, 100, 10))
# Add title
ax1.set_title("Stratiform Rainfall")
ax1.add_patch(cir1)
# Add gridlines
gl1 = ax1.gridlines(crs=proj, draw_labels=True, color='none')
gl1.xformatter = LONGITUDE_FORMATTER
gl1.yformatter = LATITUDE_FORMATTER
# Add basemap
ax1.add_feature(cfeature.COASTLINE)
ax1.add_feature(cfeature.BORDERS, linewidth=1)
ax1.add_feature(cfeature.STATES, linewidth=1)


# Plot the second contourf
im2 = ax2.contourf(xg_files['lon'].values, xg_files['lat'].values, rain_rate_convective[1], 
                 levels=np.linspace(0, 100, 10), cmap='jet', origin="lower")
# Plot the marker
ax2.plot(75.8175, 26.8207, marker="x", color="black", markersize=15, markeredgewidth=2)
# Add colorbar
cbar2 = plt.colorbar(im2, ax=ax2, label='Rainfall rate in mm/h', ticks=np.arange(0, 100, 10))
cbar2.ax.set_yticklabels(np.arange(0, 100, 10))
ax2.add_patch(cir2)
# Add title
ax2.set_title("Convective Rainfall")

# Add gridlines
gl2 = ax2.gridlines(crs=proj, draw_labels=True, color='none')
gl2.xformatter = LONGITUDE_FORMATTER
gl2.yformatter = LATITUDE_FORMATTER
# Add basemap
ax2.add_feature(cfeature.COASTLINE)
ax2.add_feature(cfeature.BORDERS, linewidth=1)
ax2.add_feature(cfeature.STATES, linewidth=1)

# Show the plot
plt.show()
11

any solution for this?

12

@ankithva

I think there are two problems you need to fix to get this to work. Firstly the arrays you are defining to hold your results before the loop are being overwritten within the loop as you’ve reused variable names.
Outside the loop you have:

rain_rate_stratiform = np.zeros((len(grid_files), 250, 250))
rain_rate_convective= np.zeros((len(grid_files), 250, 250))

and within the loop:

rain_rate_stratiform = wrl.zr.z_to_r(Zs_stratiform, a=251, b=1.48)
rain_rate_convective = wrl.zr.z_to_r(Zs_convective, a=178, b=1.51)

This leads to your results array being replaced when you first run the loop.

Also you are then trying to place a 3D array into a 2D slice of your results array. You either need to take a 2D slice of your individual timestamp rain rates (like you are plotting in your working example) to add to the results array or make your results array 4D to hold the entire 3D rain rates.

To do the first of those you could simply replace the results arrays as follows:

rain_rate_stratiform_timeseries = np.zeros((len(grid_files), 250, 250))
rain_rate_convective_timeseries = np.zeros((len(grid_files), 250, 250))

and then change the code within the loop to the following to store the same level as you are plotting in your example for a single file:

rain_rate_stratiform_timeseries[i, :, :] = np.ma.masked_greater(rain_rate_stratiform, 20)[1]
rain_rate_convective_timeseries[i, :, :] = np.ma.masked_greater(rain_rate_convective, 20)[1]

Hopefully that helps you get closer to the solution you are looking for.