计算锋生

MeteoInfo

参照MetPy在MeteoInfoLab的meteolib Jython包里增加了锋生函数 frontogenesis，需要温度、风场U, V分量以及网格间距作为输入。示例脚本程序如下：

1. fn = 'D:/Temp/nc/GFS_20101026_1200.nc'
   
2. f = addfile(fn)
   
3. 
   
4. # Set subset slice for the geographic extent of data
   
5. lon_slice = slice(400, 650)
   
6. lat_slice = slice(50, 150, -1)
   
7. 
   
8. # Grab lat/lon values (GFS will be 1D)
   
9. lats = f['lat'][lat_slice]
   
10. lons = f['lon'][lon_slice]
    
11. 
    
12. level = '85000'
    
13. hght_850 = f['Geopotential_height_isobaric'][0, level, lat_slice, lon_slice]
    
14. tmpk_850 = f['Temperature_isobaric'][0, level, lat_slice, lon_slice]
    
15. uwnd_850 = f['u-component_of_wind_isobaric'][0, level, lat_slice, lon_slice]
    
16. vwnd_850 = f['v-component_of_wind_isobaric'][0, level, lat_slice, lon_slice]
    
17. 
    
18. # Convert temperatures to degree Celsius for plotting purposes
    
19. tmpc_850 = tmpk_850 - 273.15
    
20. 
    
21. # Calculate potential temperature for frontogenesis calculation
    
22. thta_850 = meteolib.potential_temperature(850, tmpk_850)
    
23. 
    
24. dx, dy = meteolib.lat_lon_grid_deltas(lons, lats)
    
25. fronto_850 = meteolib.frontogenesis(thta_850, uwnd_850, vwnd_850, dx, dy)
    
26. 
    
27. # A conversion factor to get frontogensis units of K per 100 km per 3 h
    
28. convert_to_per_100km_3h = 1000*100*3600*3
    
29. 
    
30. #plot
    
31. proj = projinfo(proj='lcc', lon_0=-100, lat_0=35, lat_1=30, lat_2=60)
    
32. ax = axesm(projinfo=proj, tickfontsize=12)
    
33. geoshow('us_states', edgecolor=(0,102,204))
    
34. geoshow('country', edgecolor=(0,102,204))
    
35. 
    
36. # Plot 850-hPa Frontogenesis
    
37. clevs_fronto = np.arange(-8, 8.5, 0.5)
    
38. cf = contourf(lons, lats, fronto_850*convert_to_per_100km_3h, clevs_fronto,
    
39.                  cmap='BlWhRe', smooth=False)
    
40. colorbar(cf, orientation='horizontal', shrink=0.8, aspect=50, fontsize=12,
    
41.     label='Frontogenesis K / 100 km / 3 h')
    
42. 
    
43. # Plot 850-hPa Temperature in Celsius
    
44. clevs_tmpc = np.arange(-40, 41, 2)
    
45. csf = contour(lons, lats, tmpc_850, clevs_tmpc, colors='gray',
    
46.                  linestyle='--')
    
47. clabel(csf, fmt='%d', fontsize=12, drawshadow=False)
    
48. 
    
49. # Plot 850-hPa Geopotential Heights
    
50. clevs_850_hght = np.arange(0, 8000, 30)
    
51. cs = contour(lons, lats, hght_850, clevs_850_hght, colors='black')
    
52. clabel(cs, fmt='%d', fontsize=12, drawshadow=False)
    
53. 
    
54. # Plot 850-hPa Wind Barbs only plotting every fifth barb
    
55. wind_slice = (slice(None, None, 5), slice(None, None, 5))
    
56. barbs(lons[wind_slice[0]], lats[wind_slice[1]],
    
57.          uwnd_850[wind_slice], vwnd_850[wind_slice],
    
58.          color='black')
    
59. 
    
60. axis([-128, -72, 16, 55])
    
61. 
    
62. # Plot some titles to tell people what is on the map
    
63. left_title(r'GFS 850-hPa Geopotential Heights (m), Temp (C), and Winds', fontsize=12, bold=False)
    
64. right_title('Valid Time: {}'.format(f.gettime(0)), fontsize=12, bold=False)

fn = 'D:/Temp/nc/GFS_20101026_1200.nc'<br /> f = addfile(fn)<br /> <br /> # Set subset slice for the geographic extent of data<br /> lon_slice = slice(400, 650)<br /> lat_slice = slice(50, 150, -1)<br /> <br /> # Grab lat/lon values (GFS will be 1D)<br /> lats = f['lat'][lat_slice]<br /> lons = f['lon'][lon_slice]<br /> <br /> level = '85000'<br /> hght_850 = f['Geopotential_height_isobaric'][0, level, lat_slice, lon_slice]<br /> tmpk_850 = f['Temperature_isobaric'][0, level, lat_slice, lon_slice]<br /> uwnd_850 = f['u-component_of_wind_isobaric'][0, level, lat_slice, lon_slice]<br /> vwnd_850 = f['v-component_of_wind_isobaric'][0, level, lat_slice, lon_slice]<br /> <br /> # Convert temperatures to degree Celsius for plotting purposes<br /> tmpc_850 = tmpk_850 - 273.15<br /> <br /> # Calculate potential temperature for frontogenesis calculation<br /> thta_850 = meteolib.potential_temperature(850, tmpk_850)<br /> <br /> dx, dy = meteolib.lat_lon_grid_deltas(lons, lats)<br /> fronto_850 = meteolib.frontogenesis(thta_850, uwnd_850, vwnd_850, dx, dy)<br /> <br /> # A conversion factor to get frontogensis units of K per 100 km per 3 h<br /> convert_to_per_100km_3h = 1000*100*3600*3<br /> <br /> #plot<br /> proj = projinfo(proj='lcc', lon_0=-100, lat_0=35, lat_1=30, lat_2=60)<br /> ax = axesm(projinfo=proj, tickfontsize=12)<br /> geoshow('us_states', edgecolor=(0,102,204))<br /> geoshow('country', edgecolor=(0,102,204))<br /> <br /> # Plot 850-hPa Frontogenesis<br /> clevs_fronto = np.arange(-8, 8.5, 0.5)<br /> cf = contourf(lons, lats, fronto_850*convert_to_per_100km_3h, clevs_fronto,<br />                  cmap='BlWhRe', smooth=False)<br /> colorbar(cf, orientation='horizontal', shrink=0.8, aspect=50, fontsize=12,<br />     label='Frontogenesis K / 100 km / 3 h')<br /> <br /> # Plot 850-hPa Temperature in Celsius<br /> clevs_tmpc = np.arange(-40, 41, 2)<br /> csf = contour(lons, lats, tmpc_850, clevs_tmpc, colors='gray',<br />                  linestyle='--')<br /> clabel(csf, fmt='%d', fontsize=12, drawshadow=False)<br /> <br /> # Plot 850-hPa Geopotential Heights<br /> clevs_850_hght = np.arange(0, 8000, 30)<br /> cs = contour(lons, lats, hght_850, clevs_850_hght, colors='black')<br /> clabel(cs, fmt='%d', fontsize=12, drawshadow=False)<br /> <br /> # Plot 850-hPa Wind Barbs only plotting every fifth barb<br /> wind_slice = (slice(None, None, 5), slice(None, None, 5))<br /> barbs(lons[wind_slice[0]], lats[wind_slice[1]],<br />          uwnd_850[wind_slice], vwnd_850[wind_slice],<br />          color='black')<br /> <br /> axis([-128, -72, 16, 55])<br /> <br /> # Plot some titles to tell people what is on the map<br /> left_title(r'GFS 850-hPa Geopotential Heights (m), Temp (C), and Winds', fontsize=12, bold=False)<br /> right_title('Valid Time: {}'.format(f.gettime(0)), fontsize=12, bold=False)

子非鱼。

{:5_213:}{:5_213:}{:5_213:}王老师太厉害了！

lixiliang2008

王老师，参照MetPy在MeteoInfoLab的meteolib Jython包里增加了锋生函数 frontogenesis这个怎么操作？

MeteoInfo

lixiliang2008 发表于 2023-9-3 10:49
王老师，参照MetPy在MeteoInfoLab的meteolib Jython包里增加了锋生函数 frontogenesis这个怎么操作？

参考这里：https://github.com/meteoinfo/Met ... /kinematics.py#L428

lixiliang2008

王老师，按照您给的提示我做了处理。还是没搞明白。提示'module' object has no attribute 'lat_lon_grid_deltas'以及no attribute 'frontogenesis'。具体步骤是什么呢？麻烦您了

MeteoInfo

lixiliang2008 发表于 2023-9-3 16:20
王老师，按照您给的提示我做了处理。还是没搞明白。提示'module' object has no attribute 'lat_lon_grid_d ...

下载MeteoInfo最新版本试试

lixiliang2008

做出来了，多谢王老师指点。

司机是猪