计算锋生的函数 frontogenesis

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

fn = 'D:/Temp/nc/GFS_20101026_1200.nc'
f = addfile(fn)

# Set subset slice for the geographic extent of data
lon_slice = slice(400, 650)
lat_slice = slice(50, 150, -1)

# Grab lat/lon values (GFS will be 1D)
lats = f['lat'][lat_slice]
lons = f['lon'][lon_slice]

level = '85000'
hght_850 = f['Geopotential_height_isobaric'][0, level, lat_slice, lon_slice]
tmpk_850 = f['Temperature_isobaric'][0, level, lat_slice, lon_slice]
uwnd_850 = f['u-component_of_wind_isobaric'][0, level, lat_slice, lon_slice]
vwnd_850 = f['v-component_of_wind_isobaric'][0, level, lat_slice, lon_slice]

# Convert temperatures to degree Celsius for plotting purposes
tmpc_850 = tmpk_850 - 273.15

# Calculate potential temperature for frontogenesis calculation
thta_850 = meteolib.potential_temperature(850, tmpk_850)

dx, dy = meteolib.lat_lon_grid_deltas(lons, lats)
fronto_850 = meteolib.frontogenesis(thta_850, uwnd_850, vwnd_850, dx, dy)

# A conversion factor to get frontogensis units of K per 100 km per 3 h
convert_to_per_100km_3h = 1000*100*3600*3

#plot
proj = projinfo(proj='lcc', lon_0=-100, lat_0=35, lat_1=30, lat_2=60)
ax = axesm(projinfo=proj, tickfontsize=12)
geoshow('us_states', edgecolor=(0,102,204))
geoshow('country', edgecolor=(0,102,204))

# Plot 850-hPa Frontogenesis
clevs_fronto = np.arange(-8, 8.5, 0.5)
cf = contourf(lons, lats, fronto_850*convert_to_per_100km_3h, clevs_fronto,
                 cmap='BlWhRe', smooth=False)
colorbar(cf, orientation='horizontal', shrink=0.8, aspect=50, fontsize=12,
    label='Frontogenesis K / 100 km / 3 h')

# Plot 850-hPa Temperature in Celsius
clevs_tmpc = np.arange(-40, 41, 2)
csf = contour(lons, lats, tmpc_850, clevs_tmpc, colors='gray',
                 linestyle='--')
clabel(csf, fmt='%d', fontsize=12, drawshadow=False)

# Plot 850-hPa Geopotential Heights
clevs_850_hght = np.arange(0, 8000, 30)
cs = contour(lons, lats, hght_850, clevs_850_hght, colors='black')
clabel(cs, fmt='%d', fontsize=12, drawshadow=False)

# Plot 850-hPa Wind Barbs only plotting every fifth barb
wind_slice = (slice(None, None, 5), slice(None, None, 5))
barbs(lons[wind_slice[0]], lats[wind_slice[1]],
         uwnd_850[wind_slice], vwnd_850[wind_slice],
         color='black')

axis([-128, -72, 16, 55])

# Plot some titles to tell people what is on the map
left_title(r'GFS 850-hPa Geopotential Heights (m), Temp (C), and Winds', fontsize=12, bold=False)
right_title('Valid Time: {}'.format(f.gettime(0)), fontsize=12, bold=False)