calculate the MCS index

Question

I would like to optimize the code to calculate the MCS index for the Southern Hemisphere using the formula ( # SA − MCS index = [Shear (0–6 km) –20.01] ∕7.87 + [ gradiente_horizontal (775 hPa) –4.84 × 10–5] ∕5.65 × 10–5 + { – [ omega (800 hPa) + 0.27] ∕0.29] } + { – [ LI (four layers) + 2.17] ∕2.23}).

import matplotlib.pyplot as plt
import numpy as np
import xarray as xr
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import metpy
import metpy.calc as mpcalc
from metpy.units import units

# Abrir el archivo GRIB con los datos
ds = xr.open_dataset('cdas1_2022011700_.t00z.pgrbh00.grib2', engine='cfgrib', backend_kwargs={
                     'filter_by_keys': {'typeOfLevel': 'isobaricInhPa'}})

ds1 = xr.open_dataset('cdas1_2022011700_.t00z.pgrbh00.grib2', engine='cfgrib', backend_kwargs={
                     'filter_by_keys': {'typeOfLevel': 'pressureFromGroundLayer'}})
ds = ds.metpy.parse_cf()
ds1 = ds1.metpy.parse_cf()

# Filtrar los datos para obtener solo el nivel de 1000 hPa
u_1000hPa = ds.u.sel(isobaricInhPa=1000)
v_1000hPa = ds.v.sel(isobaricInhPa=1000)

# Filtrar los datos para obtener solo el nivel de 700 hPa
u_400hPa = ds.u.sel(isobaricInhPa=400)
v_400hPa = ds.v.sel(isobaricInhPa=400)

# Filtrar los datos para obtener solo el nivel de 750 hPa
u_775hPa = ds.u.sel(isobaricInhPa=775)
v_775hPa = ds.v.sel(isobaricInhPa=775)
t_775hPa = ds.t.sel(isobaricInhPa=775)

# Filtrar los datos para obtener solo el nivel de 800 hPa
w_800hPa = ds.w.sel(isobaricInhPa=800)

# Extraer los datos del LI (Lifted Index)
pli = ds1.pli
# pli = pli.where(pli > 0)

# Extraer los datos de agua precipitable
pwat = ds1.pwat

# Calcular la velocidad del viento a 650 hPa y 1000 hPa con metpy
mag400=metpy.calc.wind_speed(u_400hPa,v_400hPa)
mag1000=metpy.calc.wind_speed(u_1000hPa,v_1000hPa)

# Calcular el shear
shear=(mag400-mag1000)

# Calcular las derivadas parciales de la temperatura en las direcciones x e y (advección horizontal de temperatura)
# Inicialización de dTdX con ceros
dT_dx = np.zeros_like(t_775hPa.values)
dT_dy = np.zeros_like(t_775hPa.values)
# Usamos numpy.gradient para calcular las diferencias finitas para las derivadas espaciales
dT_dx, dT_dy = np.gradient(t_775hPa,axis=(0,1))

# Calcular el producto entre v y el gradiente horizontal de T en las direcciones x e y (Tadv=-(u750*dTdX+v750*dTdY))
gradiente_horizontal = -u_775hPa * dT_dx - v_775hPa * dT_dy

# # Calcular el MCS index
mcs=((shear.values-13.66)/5.5) + ((gradiente_horizontal-4.28e-5)/5.19e-5) + (-(w_800hPa+0.269)/0.286) + (-(pli+2.142)/2.175)

# Crear la figura y el subplot
fig = plt.figure(figsize=(16, 14))
ax = fig.add_subplot(1, 1, 1, projection=ccrs.PlateCarree())

# Agregar características geográficas
ax.add_feature(cfeature.COASTLINE)
ax.add_feature(cfeature.BORDERS)
ax.add_feature(cfeature.STATES)

# Agregar las líneas de latitud y longitud
ax.gridlines(draw_labels=True, linewidth=0.5, color='gray', alpha=0.5, linestyle='--')

# Definir algunas propiedades del mapa
central_lon, central_lat = -33, -55
extent = [-75, -45, -40, -16]
ax.set_extent(extent, crs=ccrs.PlateCarree())
ax.set_title(
    "Parcel MCS index", fontsize=16, fontweight='bold')
    
# Cambiar los valores en el rango (-0.5, 0.5) a NaN
masked_mcs_pwat = np.where((pwat < 27) & (pli > 0) & (mcs >= -0.5) & (mcs <= 0.5), np.nan, mcs)

# Dibujar el mapa de mcs
clevs_mcs = np.arange(-1.5, 5)  # levels=[-1.15, -0.12, 1.58, 2.74, 5],
contourf_mcs = ax.contourf(ds['longitude'], ds['latitude'], masked_mcs_pwat, levels = clevs_mcs,
                           cmap= 'rainbow' , transform=ccrs.PlateCarree())

#ploteo la barra de referencias
cbar_temp = plt.colorbar(
    contourf_mcs, ax=ax, orientation='horizontal',aspect=40, pad=0.02)
cbar_temp.set_label('Parcel MCS index ', size=14)

# Mostrar el mapa
plt.show()

Please edit the question to limit it to a specific problem with enough detail to identify an adequate answer. — Community, Aug 07 '23 at 18:57

score 0 · Answer 1 · answered Aug 07 '23 at 16:10

To optimize the code for calculating the MCS index, you can make use of MetPy's built-in functionality for calculations and utilize xarray's broadcasting capabilities. Additionally, you can vectorize some operations to improve the efficiency. Here's an optimized version of the code:

import matplotlib.pyplot as plt
import numpy as np
import xarray as xr
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import metpy.calc as mpcalc
from metpy.units import units

# Load the GRIB file data
ds = xr.open_dataset('cdas1_2022011700_.t00z.pgrbh00.grib2', engine='cfgrib', backend_kwargs={
                     'filter_by_keys': {'typeOfLevel': 'isobaricInhPa'}})

ds1 = xr.open_dataset('cdas1_2022011700_.t00z.pgrbh00.grib2', engine='cfgrib', backend_kwargs={
                     'filter_by_keys': {'typeOfLevel': 'pressureFromGroundLayer'}})
ds = ds.metpy.parse_cf()
ds1 = ds1.metpy.parse_cf()

# Filtrar los datos para obtener solo los niveles requeridos
u_400hPa = ds.u.sel(isobaricInhPa=400)
v_400hPa = ds.v.sel(isobaricInhPa=400)
u_775hPa = ds.u.sel(isobaricInhPa=775)
v_775hPa = ds.v.sel(isobaricInhPa=775)
t_775hPa = ds.t.sel(isobaricInhPa=775)
w_800hPa = ds.w.sel(isobaricInhPa=800)
pli = ds1.pli
pwat = ds1.pwat

# Calculate wind speed at 400 hPa and 1000 hPa
mag400 = mpcalc.wind_speed(u_400hPa, v_400hPa)
mag1000 = mpcalc.wind_speed(u_1000hPa, v_1000hPa)

# Calculate the shear
shear = mag400 - mag1000

# Calculate horizontal temperature advection (dT/dx and dT/dy)
dT_dx, dT_dy = np.gradient(t_775hPa.values, axis=(2, 1))
gradiente_horizontal = -u_775hPa * dT_dx - v_775hPa * dT_dy

# Calculate the MCS index
mcs = ((shear - 13.66) / 5.5) + ((gradiente_horizontal - 4.28e-5) / 5.19e-5) + (-(w_800hPa + 0.269) / 0.286) + (-(pli + 2.142) / 2.175)

# Create the figure and subplot
fig = plt.figure(figsize=(16, 14))
ax = fig.add_subplot(1, 1, 1, projection=ccrs.PlateCarree())

# Add geographic features
ax.add_feature(cfeature.COASTLINE)
ax.add_feature(cfeature.BORDERS)
ax.add_feature(cfeature.STATES)

# Add latitude and longitude gridlines
ax.gridlines(draw_labels=True, linewidth=0.5, color='gray', alpha=0.5, linestyle='--')

# Define some map properties
central_lon, central_lat = -33, -55
extent = [-75, -45, -40, -16]
ax.set_extent(extent, crs=ccrs.PlateCarree())
ax.set_title("Parcel MCS index", fontsize=16, fontweight='bold')

# Apply masking to the MCS index data
masked_mcs_pwat = np.where((pwat < 27) & (pli > 0) & (mcs >= -0.5) & (mcs <= 0.5), np.nan, mcs)

# Draw the MCS index map
clevs_mcs = np.arange(-1.5, 5)
contourf_mcs = ax.contourf(ds['longitude'], ds['latitude'], masked_mcs_pwat, levels=clevs_mcs, cmap='rainbow', transform=ccrs.PlateCarree())

# Add colorbar
cbar_temp = plt.colorbar(contourf_mcs, ax=ax, orientation='horizontal', aspect=40, pad=0.02)
cbar_temp.set_label('Parcel MCS index', size=14)

# Show the map
plt.show()

Hi @jspaeth . Thanks for your reply. I tried the code that you recommended in your answer. And I have an error in the calculation of the horizontal temperature advection: I would like a recommendation to improve the calculation of horizontal temperature advection. AxisError Traceback (most recent call last) CellIn[2], line 38 37 # Calculate horizontal temperature advection (dT/dx and dT/dy) ---> 38 `dT_dx, dT_dy = np.gradient(t_775hPa.values, axis=(2, 1))` AxisError: axis 2 is out of bounds for array of dimension 2 — leojim19, Aug 08 '23 at 02:26
I answered your other question about advection (https://stackoverflow.com/questions/76860347/calculate-the-partial-derivatives-of-temperature-horizontal-advection-of-temper) Here your error comes from the `axis` argument being incorrect when passing in a 2D array. I'd recommend incorporating the answer to the other question before proceeding further. — DopplerShift, Aug 20 '23 at 22:06
This open Pull Request on MetPy may also be of interest to calculating MCS motion: https://github.com/Unidata/MetPy/pull/3116 — DopplerShift, Aug 20 '23 at 22:06

calculate the MCS index

1 Answers1