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