Getting GFS forecast variables for sounding-plotting

[kylejgillett]

Hey Metpy-ers,

I would like to transform my observational sounding plotter script into a modeled plot (first tests with GFS data), and I am having trouble getting ahold of the data from the THREDDS server. I have followed along with MetPy Monday video # 164 on accessing the THREDDS server. I think I am on the right track (maybe?) but I am not sure how to declare the specific lat and lon I want for the data, I tested it with temp_raw.isel(time=0, lat=##.####, lon=-##.####) but that doesn't seem to work like I thought it should lol. I get this error:

"TypeError: invalid indexer array, does not have integer dtype: array(42.242516)"

I am likely thinking this is more simple than it really is, any ideas?

here if my code thus far:

--- define TDSCatalog as GFS 0p25deg Best---
`best_gfs = TDSCatalog('http://thredds.ucar.edu/thredds/catalog/grib/NCEP/GFS/Global_0p25deg/catalog.xml?dataset=grib/NCEP/GFS/Global_0p25deg/Best')
best_ds = best_gfs.datasets[0]
best_ds

---access the ds and define variables ---
ds = xr.open_dataset(best_ds.access_urls['OPENDAP'])
gfs_variables = ['Temperature_isobaric','Relative_humidity_isobaric',
'u-component_of_wind_isobaric', 'v-component_of_wind_isobaric',
'Geopotential_height_isobaric']

--- create temp, RH (to calc depoint), u, v, and geopotential height variables from the GFS ----
temp_raw = ds.metpy.parse_cf(gfs_variables[0])
rh_raw = ds.metpy.parse_cf(gfs_variables[1])
ucomp_raw = ds.metpy.parse_cf(gfs_variables[2])
vcomp_raw = ds.metpy.parse_cf(gfs_variables[3])
geohght_raw = ds.metpy.parse_cf(gfs_variables[4])

---- get temp_raw data for time of zero.... eventually also add lat= and lon= to get temp_raw for desired sounding location??? ----
temp_raw.isel(time=0)`

Thanks!

[DanielAdriaansen]

When using isel from Xarray, you have to provide integers which correspond to the index into the lat/lon array that you want (like you do for time, using time=0). If you want to use actual lat/lon values, you have to use plain sel:
https://docs.xarray.dev/en/latest/user-guide/indexing.html

You can string together your selection I think like this:
temp_raw.isel(time=0).sel(lat=##.####,lon=-##.####)

Note that I am not very familiar with lat/lon and using sel. If the exact lat/lon value of your sounding site isn't found within the lat/lon array, there may be ways you can tell Xarray how to select the "best" one by interpolating, or doing nearest, etc. For example: https://docs.xarray.dev/en/latest/user-guide/indexing.html#nearest-neighbor-lookups

[DanielAdriaansen]

It's just how the data are stored. If you look at your coordinate variable isobaric1, it's ordered differently than a human thinks about it, and the first pressure value (pressure[0]) is 1.0 hPa. It shouldn't matter that it's stored this way, until you get to plotting at which point you can reverse the data for the plotting step only. One way to do that is:

skew.plot(temp_rawK.isobaric1[::-1],temp_rawK.Temperature_isobaric[::-1])

Note the [::-1] is one way to reverse the data. Xarray may also have ways of doing this, for example: https://docs.xarray.dev/en/latest/generated/xarray.DataArray.transpose.html

For your longitude you can just add 360.0 to the sel step like:
temp_rawK = temp_raw.isel(time=0).sel(lat=##.####,lon=(-##.####+360.0))

[kylejgillett]

Ah prefect! Thanks!

[kylejgillett]

While I have you here, any idea how I can get pressure from this GFS catalog? Im am noticing that there isn't a pressure with altitude variable (or something similar). Any thoughts?

[DanielAdriaansen]

Pressure is the vertical coordinate of the data, so it is the same for every sounding (at every grid cell) and is a coordinate variable in the Xarray Dataset.

pres = ds['isobaric1']
print(pres)

Then you can use pres for every sounding.

Note it will not always be called "isobaric1", this is just the case for these data.

[dopplershift]

You should be able to get the pressure values consistently through the metpy xarray accessor as the vertical coordinates, so:

pres = temp_rawK.metpy.vertical