Интерполяция ветров GFS из изобари c в координаты высоты с использованием Metpy

Question

Мне было поручено составлять графики ветра на разных уровнях атмосферы для поддержки авиации. Несмотря на то, что мне удалось сделать несколько хороших графиков, используя данные модели GFS (см. Код ниже), мне действительно нужно приблизительное приближение высоты, используя координаты давления, доступные из GFS. Я использую ветер в 300 гПа, 700 гПа и 925 гПа, чтобы приблизить ветры в 30 000 футов, 9 000 футов и 3000 футов. Мой вопрос действительно для тех, кто является метпи-гуру ... есть ли способ, которым я могу интерполировать эти ветра на высоту поверхности? Было бы неплохо получить фактический ветер на этих уровнях высоты! Спасибо за любой свет, которым кто-то может поделиться на эту тему!

import cartopy.crs as ccrs
import cartopy.feature as cfeature
import matplotlib.pyplot as plt
from matplotlib.colors import ListedColormap
import numpy as np
from netCDF4 import num2date
from datetime import datetime, timedelta
from siphon.catalog import TDSCatalog
from siphon.ncss import NCSS
from PIL import Image
from matplotlib import cm


# For the vertical levels we want to grab with our queries
# Levels need to be in Pa not hPa
Levels = [30000,70000,92500]
# Time deltas for days
Deltas = [1,2,3]
#Deltas = [1]
# Levels in hPa for the file names
LevelDict = {30000:'300', 70000:'700', 92500:'925'}
# The path to where our banners are stored 
impath = 'C:\\Users\\shell\\Documents\\Python Scripts\\Banners\\'
# Final images saved here
imoutpath = 'C:\\Users\\shell\\Documents\\Python Scripts\\TVImages\\'

# Quick function for finding out which variable is the time variable in the
# netCDF files
def find_time_var(var, time_basename='time'):
    for coord_name in var.coordinates.split():
        if coord_name.startswith(time_basename):
            return coord_name
    raise ValueError('No time variable found for ' + var.name)

# Function to grab data at different levels from Siphon 
def grabData(level):

    query.var = set()
    query.variables('u-component_of_wind_isobaric', 'v-component_of_wind_isobaric')
    query.vertical_level(level)
    data = ncss.get_data(query)
    u_wind_var = data.variables['u-component_of_wind_isobaric']
    v_wind_var = data.variables['v-component_of_wind_isobaric']
    time_var = data.variables[find_time_var(u_wind_var)]
    lat_var = data.variables['lat']
    lon_var = data.variables['lon']

    return u_wind_var, v_wind_var, time_var, lat_var, lon_var

# Construct a TDSCatalog instance pointing to the gfs dataset
best_gfs = TDSCatalog('http://thredds-jetstream.unidata.ucar.edu/thredds/catalog/grib/'
                      'NCEP/GFS/Global_0p5deg/catalog.xml')

# Pull out the dataset you want to use and look at the access URLs
best_ds = list(best_gfs.datasets.values())[1]
#print(best_ds.access_urls)

# Create NCSS object to access the NetcdfSubset
ncss = NCSS(best_ds.access_urls['NetcdfSubset'])
print(best_ds.access_urls['NetcdfSubset'])

# Looping through the forecast times

for delta in Deltas:
    # Create lat/lon box and the time(s) for location you want to get data for
    now = datetime.utcnow()
    fcst = now + timedelta(days = delta)
    timestamp = datetime.strftime(fcst, '%A')
    query = ncss.query()
    query.lonlat_box(north=78, south=45, east=-90, west=-220).time(fcst)
    query.accept('netcdf4')


    # Now looping through the levels to create our plots

    for level in Levels:
        u_wind_var, v_wind_var, time_var, lat_var, lon_var = grabData(level)
        # Get actual data values and remove any size 1 dimensions
        lat = lat_var[:].squeeze()
        lon = lon_var[:].squeeze()
        u_wind = u_wind_var[:].squeeze()
        v_wind = v_wind_var[:].squeeze()
        #converting to knots
        u_windkt= u_wind*1.94384
        v_windkt= v_wind*1.94384
        wspd = np.sqrt(np.power(u_windkt,2)+np.power(v_windkt,2))

        # Convert number of hours since the reference time into an actual date
        time = num2date(time_var[:].squeeze(), time_var.units)
        print (time)
        # Combine 1D latitude and longitudes into a 2D grid of locations
        lon_2d, lat_2d = np.meshgrid(lon, lat)

        # Create new figure
        #fig = plt.figure(figsize = (18,12))
        fig = plt.figure()
        fig.set_size_inches(26.67,15)
        datacrs = ccrs.PlateCarree()
        plotcrs = ccrs.LambertConformal(central_longitude=-150,
                                       central_latitude=55,
                                       standard_parallels=(30, 60))

        # Add the map and set the extent
        ax = plt.axes(projection=plotcrs)
        ext = ax.set_extent([-195., -115., 50., 72.],datacrs)
        ext2 = ax.set_aspect('auto')
        ax.background_patch.set_fill(False)

        # Add state boundaries to plot
        ax.add_feature(cfeature.STATES, edgecolor='black', linewidth=2)

        # Add geopolitical boundaries for map reference
        ax.add_feature(cfeature.COASTLINE.with_scale('50m'))
        ax.add_feature(cfeature.OCEAN.with_scale('50m'))
        ax.add_feature(cfeature.LAND.with_scale('50m'),facecolor = '#cc9666', linewidth = 4)

        if level == 30000:
            spdrng_sped = np.arange(30, 190, 2)
            windlvl = 'Jet_Stream'
        elif level == 70000:
            spdrng_sped = np.arange(20, 100, 1)
            windlvl = '9000_Winds_Aloft'
        elif level == 92500:
            spdrng_sped = np.arange(20, 80, 1)
            windlvl = '3000_Winds_Aloft'
        else:
            pass


        top = cm.get_cmap('Greens')
        middle = cm.get_cmap('YlOrRd')
        bottom = cm.get_cmap('BuPu_r')
        newcolors = np.vstack((top(np.linspace(0, 1, 128)),
                       middle(np.linspace(0, 1, 128))))
        newcolors2 = np.vstack((newcolors,bottom(np.linspace(0,1,128))))

        cmap = ListedColormap(newcolors2)
        cf = ax.contourf(lon_2d, lat_2d, wspd, spdrng_sped, cmap=cmap,
                         transform=datacrs, extend = 'max', alpha=0.75)

        cbar = plt.colorbar(cf, orientation='horizontal', pad=0, aspect=50,
                            drawedges = 'true')
        cbar.ax.tick_params(labelsize=16)
        wslice = slice(1, None, 4)

        ax.quiver(lon_2d[wslice, wslice], lat_2d[wslice, wslice],
                       u_windkt[wslice, wslice], v_windkt[wslice, wslice], width=0.0015,
                       headlength=4, headwidth=3, angles='xy', color='black', transform = datacrs)

        plt.savefig(imoutpath+'TV_UpperAir'+LevelDict[level]+'_'+timestamp+'.png',bbox_inches= 'tight')

        # Now we use Pillow to overlay the banner with the appropriate day
        background = Image.open(imoutpath+'TV_UpperAir'+LevelDict[level]+'_'+timestamp+'.png')
        im = Image.open(impath+'Banner_'+windlvl+'_'+timestamp+'.png')

        # resize the image
        size = background.size
        im = im.resize(size,Image.ANTIALIAS)

        background.paste(im, (17, 8), im)
        background.save(imoutpath+'TV_UpperAir'+LevelDict[level]+'_'+timestamp+'.png','PNG')

Answer 1

Спасибо за вопрос! Мой подход заключается в том, чтобы для каждого отдельного столбца интерполировать координату давления высоты геопотенциала на выходе GFS на предоставленные вами высоты, чтобы оценить давление каждого уровня высоты для каждого столбца. Затем я могу использовать это давление для интерполяции GFS-вывода u, v на. GPH-выход GFS и ветер имеют немного разные координаты давления, поэтому я интерполировал дважды. Я выполнил интерполяцию, используя MetPy interpolate.log_interpolate_1d, который выполняет линейную интерполяцию в журнале входов. Вот код, который я использовал!

from datetime import datetime
import numpy as np
import metpy.calc as mpcalc
from metpy.units import units
from metpy.interpolate import log_interpolate_1d
from siphon.catalog import TDSCatalog


gfs_url = 'https://tds.scigw.unidata.ucar.edu/thredds/catalog/grib/NCEP/GFS/Global_0p5deg/catalog.xml'
cat = TDSCatalog(gfs_url)

now = datetime.utcnow()

# A shortcut to NCSS
ncss = cat.datasets['Best GFS Half Degree Forecast Time Series'].subset()

query = ncss.query()
query.var = set()
query.variables('u-component_of_wind_isobaric', 'v-component_of_wind_isobaric', 'Geopotential_height_isobaric')
query.lonlat_box(north=78, south=45, east=-90, west=-220)
query.time(now)
query.accept('netcdf4')

data = ncss.get_data(query)


# Reading in the u(isobaric), v(isobaric), isobaric vars and the GPH(isobaric6) and isobaric6 vars
# These are two slightly different vertical pressure coordinates.
# We will also assign units here, and this can allow us to go ahead and convert to knots
lat = units.Quantity(data.variables['lat'][:].squeeze(), units('degrees'))
lon = units.Quantity(data.variables['lon'][:].squeeze(), units('degrees'))
iso_wind = units.Quantity(data.variables['isobaric'][:].squeeze(), units('Pa'))
iso_gph = units.Quantity(data.variables['isobaric6'][:].squeeze(), units('Pa'))
u = units.Quantity(data.variables['u-component_of_wind_isobaric'][:].squeeze(), units('m/s')).to(units('knots'))
v = units.Quantity(data.variables['v-component_of_wind_isobaric'][:].squeeze(), units('m/s')).to(units('knots'))
gph = units.Quantity(data.variables['Geopotential_height_isobaric'][:].squeeze(), units('gpm'))


# Here we will select our altitudes to interpolate onto and convert them to geopotential meters 
altitudes = ([30000., 9000., 3000.] * units('ft')).to(units('gpm'))

# Now we will interpolate the pressure coordinate for model output geopotential height to 
# estimate the pressure level for our given altitudes at each grid point
pressures_of_alts = np.zeros((len(altitudes), len(lat), len(lon)))

for ilat in range(len(lat)):
    for ilon in range(len(lon)):
        pressures_of_alts[:, ilat, ilon] = log_interpolate_1d(altitudes,
                                                              gph[:, ilat, ilon],
                                                              iso_gph)

pressures_of_alts = pressures_of_alts * units('Pa')


# Similarly, we will use our interpolated pressures to interpolate
# our u and v winds across their given pressure coordinates.
# This will provide u, v at each of our interpolated pressure
# levels corresponding to our provided initial altitudes
u_at_levs = np.zeros((len(altitudes), len(lat), len(lon)))
v_at_levs = np.zeros((len(altitudes), len(lat), len(lon)))

for ilat in range(len(lat)):
    for ilon in range(len(lon)):
        u_at_levs[:, ilat, ilon], v_at_levs[:, ilat, ilon] = log_interpolate_1d(pressures_of_alts[:, ilat, ilon],
                                                                                iso_wind,
                                                                                u[:, ilat, ilon],
                                                                                v[:, ilat, ilon])

u_at_levs = u_at_levs * units('knots')
v_at_levs = v_at_levs * units('knots')


# We can use mpcalc to calculate a wind speed array from these
wspd = mpcalc.wind_speed(u_at_levs, v_at_levs)

Я смог взять свой вывод из этого и привести его к вашему графическому коду (с некоторым удалением единиц).

Ваш GFS на 300 гПа обмоток
Мои "30000-футовые" ветры GFS

Здесь - это то, как выглядят мои интерполированные поля давления на каждом предполагаемом уровне высоты.

Надеюсь, это поможет!

Answer 2

Я не уверен, что это то, что вы ищете (я очень плохо знаком с Metpy), но я использовал функцию metpy height_to_pressure_std (altitude). Он выражается в единицах hPa, которые затем я конвертирую в Паскали, а затем в значения без единиц измерения для использования в функции Siphon vertical_level (float).