calculate the partial derivatives of temperature (horizontal advection of temperature)

I would like to know which way is the most correct to calculate the partial derivatives of temperature in the x and y directions (horizontal advection of temperature)
The second code used the data matrix of temperature, zonal wind and meridional wind. Extract the data for temperature (T), zonal wind component (u) and meridional wind component (v)

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'}})
ds = ds.metpy.parse_cf()
# 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)
# 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)
# 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
# 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(
"horizontal advection of temperature", fontsize=16, fontweight='bold')
# Dibujar el mapa de gradiente_horizontal
contourf = ax.contourf(ds['longitude'], ds['latitude'], gradiente_horizontal,
cmap= 'rainbow' , transform=ccrs.PlateCarree())
#ploteo la barra de referencias
cbar_temp = plt.colorbar(
contourf, ax=ax, orientation='horizontal',aspect=40, pad=0.02)
cbar_temp.set_label('horizontal advection of temperature ', size=14)
# Mostrar el mapa
plt.show()

o

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'}})
ds = ds.metpy.parse_cf()
# calcular advección de temperatura para la matriz de datos de temperatura, viento zonal y viento meridional
# Extraer los datos de temperatura (T), componente zonal del viento (u) y componente meridional del viento (v)
T = ds.t
u = ds.u
v = ds.v
# Calcular las derivadas parciales de temperatura en las direcciones x e y
dT_dx, dT_dy = np.gradient(T, axis=(1,2))
# Calcular la advección de temperatura
advection_of_temperature = -u * dT_dx - v * dT_dy
# Ahora puedes acceder a los valores de la advección de temperatura:
advection_of_temperature = advection_of_temperature.sel(isobaricInhPa=775)
# 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(
"horizontal advection of temperature", fontsize=16, fontweight='bold')
# Dibujar el mapa de gradiente_horizontal
contourf = ax.contourf(ds['longitude'], ds['latitude'], advection_of_temperature,
cmap= 'rainbow' , transform=ccrs.PlateCarree())
#ploteo la barra de referencias
cbar_temp = plt.colorbar(
contourf, ax=ax, orientation='horizontal',aspect=40, pad=0.02)
cbar_temp.set_label('horizontal advection of temperature ', size=14)
# Mostrar el mapa
plt.show()

The use np.gradient should be fine regardless of which approach you choose, but there is another bug in your code. You assume the return from np.gradient is dt_dx, dt_dy. However, the original order of the data is pressure x latitude x longitude. This means you've incorrectly calculated advection above because you're mismatching the wind components and gradient components.

Also, you're using dT/d(lat) as dT/dy, which isn't unit-wise correct to multiply with wind and to get the right dimensionality of temperature/time.

Since you're already using MetPy, the easiest way to accomplish your goal is to use MetPy's advection function, which can take advantage of xarray coordinates and other metadata to do the "right thing" more or less automatically (using parse_cf() in advance isn't even really necessary):

import metpy.calc as mpcalc
mpcalc.advection(ds.t, u=ds.u, v=ds.v)