使用metpy绘制探空图

bsb0

最近在使用metpy的SkewT功能绘制探空图，但是三条曲线始终存在着问题，包括官网上的图，在高空300hPa上气温还出奇地高，网上搜到还有高空比地面温度还高的，想求助一下是什么原因，代码如下，

1. import pandas as pd
   
2. import metpy.calc as mpcalc
   
3. from metpy.units import units
   
4. import numpy as np
   
5. from scipy import interpolate
   
6. import matplotlib.pyplot as plt
   
7. 
   
8. df = pd.read_csv('d:/微波辐射计数据/ZP2021-05-28_00-00-00LV2.csv',skiprows=(1,2))#读取dataframe，跳过1,2行
   
9. df1 = df.drop(['Date/Time'],axis=1).astype('float')#删除时间列，并转换为float格式
   
10. time = df['Date/Time']#读取时间并去重
    
11. IsDuplicated = time.duplicated()
    
12. time = time.drop_duplicates()
    
13. h1 = df.columns[11:94]#读取高度列表并转换为float，单位为米
    
14. h2 = h1.str.replace(r'km','0')
    
15. h3 = h2.astype('float')*1000
    
16. #方便计算，转换为array
    
17. ds_T = df1.loc[0].values
    
18. ds_ρ = df1.loc[1].values
    
19. ds_rh = df1.loc[2].values
    
20. T = ds_T[10:93]+273.5
    
21. ρ = ds_ρ[10:93]/1000
    
22. P0 = ds_T[4]#地面气压
    
23. #计算虚温
    
24. Tv = (1+0.378*(ρ*287.05*T/P0))*T
    
25. x = 9.8/(-287.05*Tv)*h3
    
26. #压高公式，将高度换算为气压
    
27. P1 = P0*np.exp(x)
    
28. P = units.Quantity(P1,"hPa")
    
29. e = units.Quantity(ρ*287.05*T/100,"hPa")#计算水汽压，单位为J/m*3等于pa
    
30. e = mpcalc.vapor_pressure(P,e/P)
    
31. t = units.Quantity(ds_T[10:93],"degC")
    
32. rh = units.Quantity(ds_rh[10:93],"%")
    
33. Td = mpcalc.dewpoint_from_relative_humidity(t,rh)
    
34. nd = np.array([P,t,Td])
    
35. #插值计算
    
36. x = nd[0]
    
37. y = nd[1]
    
38. z = nd[2]
    
39. P_new = np.arange(max(x),min(x),-0.1)
    
40. f1 = interpolate.interp1d(x, y, kind = "slinear")
    
41. f2 = interpolate.interp1d(x, z, kind = "slinear")
    
42. t_new = f1(P_new)
    
43. Td_new = f2(P_new)
    
44. nd = np.array([P_new,t_new,Td_new])
    
45. P = units.Quantity(nd[0],"hPa")
    
46. t = units.Quantity(nd[1],"degC")
    
47. Td = units.Quantity(nd[2],"degC")
    
48. fig = plt.figure(figsize=(18, 9))
    
49. plt.rcParams['font.sans-serif']=['SimHei'] #添加中文字体
    
50. plt.rcParams['axes.unicode_minus']=False #显示负号
    
51. # ax = fig.add_subplot(1,1,1)
    
52. skew = SkewT(fig, rotation=45)
    
53. skew.plot(P, t, 'r')
    
54. skew.plot(P, Td, 'g')
    
55. skew.ax.set_ylim(1000, 100)
    
56. skew.ax.set_xlim(-40, 60)
    
57. # Add the relevant special lines
    
58. skew.plot_dry_adiabats()
    
59. skew.plot_moist_adiabats()
    
60. skew.plot_mixing_lines()
    
61. plt.title('%s'%(time[0]),loc = 'right',fontsize = 18)
    
62. plt.suptitle('温度对数压力图',fontsize = 25)

复制代码

图一是楼主所绘制，图二是metpy官网。

edwardli

这个是斜温图，Skew-T Diagram。特点就是：温度坐标是左下到右上倾斜的。对于你的图，地面是21℃左右，500hPa是-8℃左右，300hPa是-28℃左右。
对于下方的官网示例图，分别是22.5℃，-15℃，-43℃左右。
不是图的问题，是“识图”的问题。

werewolf

edwardli 发表于 2021-6-17 18:52
这个是斜温图，Skew-T Diagram。特点就是：温度坐标是左下到右上倾斜的。对于你的图，地面是21℃左右，500h ...

哈哈哈哈，孝死

bsb0

edwardli 发表于 2021-6-17 18:52
这个是斜温图，Skew-T Diagram。特点就是：温度坐标是左下到右上倾斜的。对于你的图，地面是21℃左右，500h ...

谢谢指出我这么低级的错误，赶紧去补课了

独孤酒见

楼主及各位高人，求助问题。
我用网上的脚本画探空图。
........
p = df['pressure'].values * units.hPa                # 单位：hPa
T = df['temperature'].values * units.degC            # 单位：℃
Td = df['dewpoint'].values * units.degC              # 单位：℃
wind_speed = df['speed'].values * units.knots        # 单位：knot
wind_dir = df['direction'].values * units.degrees    # 单位：°
u, v = mpcalc.wind_components(wind_speed, wind_dir)  # 计算水平风速u和v

求u，v的时候出问题了。
提示：

TypeError                                 Traceback (most recent call last)
~\AppData\Local\Temp/ipykernel_746504/612715894.py in <module>
      4 wind_speed = df['speed'].values * units.knots        # 单位：knot
      5 wind_dir = df['direction'].values * units.degrees    # 单位：°
----> 6 u, v = mpcalc.wind_components(wind_speed, wind_dir)  # 计算水平风速u和v

C:\Anaconda\anaconda3\lib\site-packages\metpy\xarray.py in wrapper(*args, **kwargs)
   1228
   1229             # Evaluate inner calculation
-> 1230             result = func(*bound_args.args, **bound_args.kwargs)
   1231
   1232             # Wrap output based on match and match_unit

C:\Anaconda\anaconda3\lib\site-packages\metpy\units.py in wrapper(*args, **kwargs)
    286         def wrapper(*args, **kwargs):
    287             _check_units_inner_helper(func, sig, defaults, dims, *args, **kwargs)
--> 288             return func(*args, **kwargs)
    289
    290         return wrapper

C:\Anaconda\anaconda3\lib\site-packages\metpy\calc\basic.py in wind_components(speed, wind_direction)
    147
    148     """
--> 149     wind_direction = _check_radians(wind_direction, max_radians=4 * np.pi)
    150     u = -speed * np.sin(wind_direction)
    151     v = -speed * np.cos(wind_direction)

C:\Anaconda\anaconda3\lib\site-packages\metpy\calc\basic.py in _check_radians(value, max_radians)
   1269     """
   1270     with contextlib.suppress(AttributeError):
-> 1271         value = value.to('radians').m
   1272     if np.any(np.greater(np.abs(value), max_radians)):
   1273         warnings.warn('Input over {} radians. '

C:\Anaconda\anaconda3\lib\site-packages\pint\quantity.py in to(self, other, *contexts, **ctx_kwargs)
    722         other = to_units_container(other, self._REGISTRY)
    723
--> 724         magnitude = self._convert_magnitude_not_inplace(other, *contexts, **ctx_kwargs)
    725
    726         return self.__class__(magnitude, other)

C:\Anaconda\anaconda3\lib\site-packages\pint\quantity.py in _convert_magnitude_not_inplace(self, other, *contexts, **ctx_kwargs)
    671                 return self._REGISTRY.convert(self._magnitude, self._units, other)
    672
--> 673         return self._REGISTRY.convert(self._magnitude, self._units, other)
    674
    675     def _convert_magnitude(self, other, *contexts, **ctx_kwargs):

C:\Anaconda\anaconda3\lib\site-packages\pint\registry.py in convert(self, value, src, dst, inplace)
   1001             return value
   1002
-> 1003         return self._convert(value, src, dst, inplace)
   1004
   1005     def _convert(self, value, src, dst, inplace=False, check_dimensionality=True):

C:\Anaconda\anaconda3\lib\site-packages\pint\registry.py in _convert(self, value, src, dst, inplace)
   1915                 value, src = src._magnitude, src._units
   1916
-> 1917         return super()._convert(value, src, dst, inplace)
   1918
   1919     def _get_compatible_units(self, input_units, group_or_system):

C:\Anaconda\anaconda3\lib\site-packages\pint\registry.py in _convert(self, value, src, dst, inplace)
   1516
   1517         if not (src_offset_unit or dst_offset_unit):
-> 1518             return super()._convert(value, src, dst, inplace)
   1519
   1520         src_dim = self._get_dimensionality(src)

C:\Anaconda\anaconda3\lib\site-packages\pint\registry.py in _convert(self, value, src, dst, inplace, check_dimensionality)
   1050             value *= factor
   1051         else:
-> 1052             value = value * factor
   1053
   1054         return value

TypeError: can't multiply sequence by non-int of type 'float'