Error ( TypeError: 'GeometryCollection' object is not subscriptable ) when creating the figure with global GFS data gfs.t1xz.pgrb2.0p25.f00x. With some files I don't get this error.
The error is in the next line of the code.
contourf_mcs = ax.contourf(ds['longitude'], ds['latitude'], masked_mcs_pwat, levels=[-1.15, -0.12, 1.58, 2.74], cmap=cmap, transform=ccrs.PlateCarree(),extend= 'both')
I use python3.11 in Anaconda environment, matplotlib 3.7.2 metpy 1.5.1 cartopy 0.22.0
Code fragment:
Pyhton
import pygrib
import xarray as xr
import metpy.calc as mpcalc
import matplotlib.pyplot as plt
import matplotlib as mpl
import numpy as np
import cartopy.crs as ccrs
import cartopy.feature as cfeature
def read_f00x_data (file):
indices = []
variable_mapping = {}
grbs = pygrib.open(file)
ds = xr.Dataset()
for idx in grbs:
# buscar los nombres de las variables
if idx.name == 'U component of wind' and idx.level == 500:
indices.append(idx.messagenumber)
variable_mapping[idx.messagenumber] = 'u500hPa'
if idx.name == 'V component of wind' and idx.level == 500:
indices.append(idx.messagenumber)
variable_mapping[idx.messagenumber] = 'v500hPa'
if idx.name == 'Temperature' and idx.level == 750:
indices.append(idx.messagenumber)
variable_mapping[idx.messagenumber] = 't750hPa'
if idx.name == 'U component of wind' and idx.level == 750:
indices.append(idx.messagenumber)
variable_mapping[idx.messagenumber] = 'u750hPa'
if idx.name == 'V component of wind' and idx.level == 750:
indices.append(idx.messagenumber)
variable_mapping[idx.messagenumber] = 'v750hPa'
if idx.name == 'U component of wind' and idx.level ==1000:
indices.append(idx.messagenumber)
variable_mapping[idx.messagenumber] = 'u1000hPa'
if idx.name == 'V component of wind' and idx.level == 1000:
indices.append(idx.messagenumber)
variable_mapping[idx.messagenumber] = 'v1000hPa'
if idx.name == 'Vertical velocity' and idx.level == 800:
indices.append(idx.messagenumber)
variable_mapping[idx.messagenumber] = 'w800hPa'
if idx.name == 'Precipitable water':
indices.append(idx.messagenumber)
variable_mapping[idx.messagenumber] = 'pwatcapa'
if idx.name == 'Best (4-layer) lifted index':
indices.append(idx.messagenumber)
variable_mapping[idx.messagenumber] = 'lftx4'
grbs.close()
grbs = pygrib.open(file)
ds = xr.Dataset()
for idx in indices:
grb = grbs.message(idx)
variable_short_name = variable_mapping.get(idx, grb.shortName)
unit_of_measure = grb.units
START_DATE = grb.dataDate
latitudes, longitudes = grb.latlons()
values = grb.values
data = xr.DataArray(values, coords={'latitude': latitudes[:, 0], 'longitude': longitudes[0, :]},
dims=['latitude', 'longitude'], attrs={'units': unit_of_measure, 'long_name': grb.name, 'short_name': variable_short_name, 'START_DATE': START_DATE})
ds[variable_short_name] = data
grbs.close()
ds.metpy.parse_cf()
return ds
def calculate_mcs_index(ds):
u_500hPa = ds.u500hPa
v_500hPa = ds.v500hPa
u_1000hPa = ds.u1000hPa
v_1000hPa = ds.v1000hPa
u_750hPa = ds.u750hPa
v_750hPa = ds.v750hPa
t_750hPa = ds.t750hPa
w_800hPa = ds.w800hPa
pli = ds.lftx4
pwat = ds.pwatcapa
mag500 = mpcalc.wind_speed(u_500hPa, v_500hPa)
mag1000 = mpcalc.wind_speed(u_1000hPa, v_1000hPa)
shear = mag500 - mag1000
gradiente_horizontal = mpcalc.advection(t_750hPa, u_750hPa, v_750hPa)
mcs = ((shear.values - 13.66) / 5.5) + ((gradiente_horizontal.values - 4.28e-5) / 5.19e-5) + (-(w_800hPa + 0.269) / 0.286) + (-(pli + 2.142) / 2.175)
return mcs
def create_mcs_map(ds, mcs, pwat, pli):
fig = plt.figure(figsize=(16, 14))
ax = fig.add_subplot(1, 1, 1, projection=ccrs.PlateCarree())
ax.add_feature(cfeature.COASTLINE)
ax.add_feature(cfeature.BORDERS)
# ax.add_feature(cfeature.STATES)
ax.gridlines(draw_labels=True, linewidth=0.5, color='gray', alpha=0.5, linestyle='--')
extent = [-75, -45, -40, -16]
ax.set_extent(extent, crs=ccrs.PlateCarree())
ax.set_title(f"Parcel MCS index {fecha_actual}.{archivo[20:]} {fecha_actual_title}.{hora_local}", fontsize=16, fontweight='bold')
masked_mcs_pwat = np.where((pwat < 27) & (pli > 0.9), np.nan, mcs)
cmap = (mpl.colors.ListedColormap(['green', 'yellow', 'orange'])
.with_extremes(over='red', under='white'))
contourf_mcs = ax.contourf(ds['longitude'], ds['latitude'], masked_mcs_pwat, levels=[-1.15, -0.12, 1.58, 2.74], cmap=cmap, transform=ccrs.PlateCarree(),extend='both')
bounds = [-1.15, -0.12, 1.58, 2.74]
norm = mpl.colors.BoundaryNorm(bounds, cmap.N)
cbar_temp = plt.colorbar(
mpl.cm.ScalarMappable(cmap=cmap, norm=norm),
ax=ax,
extend='both',
extendfrac='auto',
ticks=bounds,
spacing='uniform',
orientation='horizontal',
label='Parcel MCS index',
aspect=30,
pad=0.02,
shrink=0.8,
fraction=0.1,
drawedges=True
)
cbar_temp.set_ticklabels([-1.15, -0.12, 1.58, '2.74'])
cbar_temp.set_label('Parcel MCS index', fontsize=14)
ds = read_f00x_data(archivo)
mcs = calculate_mcs_index(ds)
create_mcs_map(ds, mcs, ds.pwatcapa, ds.lftx4)