[FIXED] How to properly calculate temperature advection with metpy, error with units

Issue

I'm kinda new with Metpy. I've been trying to calculate the temperature advection with Metpy but it's been unsuccessful. Since I'm new with this package, I don't understand why needs to have units to work properly. When I calculate temperature advection I end with some weird lines on my maps and I don't know why. I think it's because of the units or something but I'm not sure. I attach my script below:

import numpy as np
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import matplotlib.gridspec as gridspec
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
import scipy.ndimage as ndimage
import metpy.calc as mpcalc
from metpy.units import units
import matplotlib.patches as mpatches
from metpy.units import units
from netCDF4 import Dataset
from netCDF4 import num2date
from siphon.catalog import TDSCatalog
from pint import UnitRegistry

crs = ccrs.PlateCarree() #ccrs.Mercator()
def plot_background(ax):
    ax.set_extent([-80., -15., -10., -60.],ccrs.PlateCarree())
    ax.add_feature(cfeature.COASTLINE.with_scale('50m'), linewidth=0.5)
    ax.add_feature(cfeature.BORDERS, linewidth=0.5)
    gl = ax.gridlines(ccrs.PlateCarree(), draw_labels=True,
                  linewidth=2, color='gray', alpha=0.5, linestyle='--')
    gl.xlabels_top = False
    gl.xlabels_bottom = True
    gl.ylabels_left = True
    gl.ylabels_right = False
    gl.xlines = False; gl.ylines= False
    gl.xformatter = LONGITUDE_FORMATTER
    gl.yformatter = LATITUDE_FORMATTER
    gl.xlabel_style = {'size': 12, 'color': 'black'}
    gl.ylabel_style = {'size': 12, 'color': 'black'}
    return ax


def define_box_region(lon1,lon2,lat1,lat2,):
    lat_corners = np.array([lat1,  lat1, lat2, lat2]); 
    lon_corners = np.array([ lon1, lon2, lon2, lon1]) 
    poly_corners = np.zeros((len(lat_corners), 2))
    poly_corners[:,0] = lon_corners; poly_corners[:,1] = lat_corners
    poly = mpatches.Polygon(poly_corners, closed=True, ec='m', fill=False, lw=1, fc="yellow", transform=ccrs.PlateCarree())
    return poly
infilesRegEraDJF=list(open('DJFRegEra.txt', 'r'))


print(infilesRegEraDJF[1][:-1])
d1=infilesRegEraDJF[1][:-1]
data=Dataset(d1)
tahist1=data.variables['ta850'][:,:]
tahist1 = np.squeeze(tahist1, axis=None) #gsp =ERA
tahist1=np.transpose(tahist1)
# tahist1=tahist1-273.15
temp850=tahist1 *units.degK
# t850=tahist1.to('degC')
# t850=tahist1*units.degC
lats_data= data.variables['lat'][:]
lats=np.transpose(lats_data)
lons_data= data.variables['lon'][:]
lons=np.transpose(lons_data)


print(infilesRegEraDJF[2][:-1])
d1=infilesRegEraDJF[2][:-1]
data=Dataset(d1)
uahist1=data.variables['ua850'][:,:]
uahist1 = np.squeeze(uahist1, axis=None) #gsp =ERA
uahist1=np.transpose(uahist1)
u850=uahist1*units('m/s')

print(infilesRegEraDJF[3][:-1])
d1=infilesRegEraDJF[3][:-1]
data=Dataset(d1)
vahist1=data.variables['va850'][:,:]
vahist1 = np.squeeze(vahist1, axis=None) #gsp =ERA
vahist1=np.transpose(vahist1)
v850=vahist1*units('m/s')

dx1, dy1 = mpcalc.lat_lon_grid_deltas(lons, lats)

advRegEra = mpcalc.advection(temp850, [u850, v850],
                       (dx1, dy1), dim_order='yx')

advregera2=advRegEra.to(units('delta_degC/sec'))


values_adv=[-0.5,-0.4,-0.3,-0.2,-0.1,0,0.1,0.2,0.3,0.4,0.5]


fig, axarr = plt.subplots(nrows=3, ncols=2, figsize=(9.5, 11), constrained_layout=True,
                              subplot_kw={'projection': crs})
axlist = axarr.flatten()
for ax in axlist:
 plot_background(ax)

cf1=axlist[1].contourf(lons, lats,advRegEra,values_adv,cmap='RdBu_r',extend='both',transform=ccrs.PlateCarree())
cf3 = axlist[1].quiver(lons[::17,::17], lats[::17,::17],uahist1[::17,::17],vahist1[::17,::17],scale=200,headwidth=5, headlength=5,transform=ccrs.PlateCarree())
axlist[1].set_title('RegERA', fontsize=14)
axlist[1].quiverkey(cf3,X=0.9, Y=1.05,U=10,label='10 m/s', labelpos='E')
poly=define_box_region(-52,-38,-35,-22.5); axlist[1].add_patch(poly)  #Box SBR
poly=define_box_region(-67.5,-52,-37.5,-25); axlist[1].add_patch(poly)  #Box LPB
poly=define_box_region(-72.5,-57.5,-37.5,-55); axlist[1].add_patch(poly)  #Box ARG


cax=fig.add_axes([0.95,0.30,0.025,0.40])
ax3=fig.colorbar(cf1,ticks=values_adv,cax=cax,orientation='vertical',label='(°C/hr)')
cax.tick_params(labelsize=14)
fig.set_constrained_layout_pads(w_pad=0.1, h_pad=0.1, hspace=0.1, wspace=0.1)

When I run the line where I convert the adv (where the units should be K/s) into advregera2=advRegEra.to(units('delta_degC/sec')) I get the following error

DimensionalityError: Cannot convert from '1 / meter' (1 / [length]) to 'delta_degree_Celsius / second' ([temperature] / [time])

I tried also to put the units from the beginning on the left side (tahist1=units.K*tahis1) and this works but when I plot the map I get some weird lines.

Can someone kindly help me? It's my first time using metpy so it's been difficult for me to understand how the function of advection temperature works :(. I'm using spyder 3.7 on Linux OpenSUSE 15.1 leap

Solution

advection is definitely one of the trickier functions to use in MetPy. Since you're using netcdf4-python to open the files, you definitely want to multiply with the units on the left, like:

u850 = units('m/s') * uahist1

That should fix your unit issues. Regarding the weird lines, I'm guessing they look like the image in this similar GitHub issue? If so, that's caused by a problem where your plot wrapping around from +180 to -180 longitude, but the data start at 0 and end at 360 longitude (at least that's what I'm guessing). Assuming that's the issue, since you're not actually needing to cross those cuts, you should be able to remove those points by only plotting sub slices of your data around the region of interest. Something like:

lon_subset = slice(200, 300)
cf1=axlist[1].contourf(lons[lon_subset], lats, advRegEra[:, lon_subset],
    values_adv, cmap='RdBu_r',extend='both',transform=ccrs.PlateCarree())

Replace 200 and 300 above with indices that are correct.

Answered By - DopplerShift