Issue with wind_speed calculating speed with reanalysis data

[graffmet11]

What went wrong?

Using 20th Century Reanalysis Data (V3) to make some simple charts for analysis activities for my students.

Code below produces 500mb chart with height lines and wind barbs, but when I attempt to compute windspeed from u and v components, the script fails. I have a similar issue with computing vorticity. If you comment out the line below "#this makes code fail", you'll see the chart works.

I feel like I'm missing something very simple and any help is appreciated.

Operating System

MacOS

Version

Latest

Python Version

3.12.2

Code to Reproduce

#Adapted from Script by Dr. Kevin Goebbert

from datetime import datetime, timezone, timedelta
from datetime import date
from pathlib import Path
from metpy.plots import declarative
from metpy.units import units
import xarray as xr
import matplotlib.pyplot as plt
import metpy.plots as mpplots
import metpy.calc as mpcalc
import numpy as np
import pandas as pd


#Initial Settings
#----------------------------------------------------------------------
# Set Level for Analysis Surfaceand Plotted Parameters
Chart_Title_Variables = 'HEIGHTS_WINDS'
Model_Type = '20thC_ReanV3'
field = 'hgt'
field2 = 'uwnd'
field3 = 'vwnd'
field4 = 'air'
level = 500
#area = 'mx'

maxlat = 52
minlat = 22
maxlon = -66
minlon = -127
dataminlon = (minlon-2)
datamaxlon = (maxlon+2)
dataminlat = (minlat-2)
datamaxlat = (maxlat+2)
area = (minlon, maxlon, minlat, maxlat)

path = str('http://psl.noaa.gov/thredds/dodsC/Datasets/20thC_ReanV3/prsSI/'+(field)+'.1847.nc')
Directory = str('/Users/Will/Desktop/AFIT Courses/MET 550/DATA/CHARTS/SFC/')
File_Time = ''
FileName =  str(level)+str('mb_')+str(Chart_Title_Variables)+str('_')+str(Model_Type)+str(File_Time[:-3])+str('.tiff')
FileTitle = (f'20CRV3 {level}mb Geopotential Heights, Temps (C), and Winds (Kt) ')
timepos = 2045 #time position in file

#Set DataSource and extract available times from NetCDF File for First Parameter
ds = xr.open_dataset(path).metpy.parse_cf()
dr = xr.open_dataset(path).metpy.parse_cf()
ds = ds['hgt'].metpy.time.data.astype('datetime64[ms]').astype('O')
time = pd.to_datetime(ds)
timekey = time[timepos] #Set Time Here as Timekey (Number Equals timestep value from NetCDF file)
File_Time = timekey.strftime('%m-%d-%Y%:%H')
#print(timekey)

#Set DataSource and extract available times from NetCDF File for Second Parameter
path2 = str('http://psl.noaa.gov/thredds/dodsC/Datasets/20thC_ReanV3/prsSI/'+(field2)+'.1847.nc')
dr2 = xr.open_dataset(path2).metpy.parse_cf()
#Set DataSource and extract available times from NetCDF File for Third Parameter
path3 = str('http://psl.noaa.gov/thredds/dodsC/Datasets/20thC_ReanV3/prsSI/'+(field3)+'.1847.nc')
path4 = str('http://psl.noaa.gov/thredds/dodsC/Datasets/20thC_ReanV3/prsSI/'+(field4)+'.1847.nc')
dr3 = xr.open_dataset(path3).metpy.parse_cf()
dr4 = xr.open_dataset(path4).metpy.parse_cf()
dr2['vwnd']= dr3['vwnd']

#Subset of plotting area for final image
#####STUCK HERE!
dr = dr.sel(lat=slice(dataminlat, datamaxlat), lon=slice(360+(dataminlon), 360+(datamaxlon)))
dr2 = dr2.sel(lat=slice(dataminlat, datamaxlat), lon=slice(360+(dataminlon), 360+(datamaxlon)))
dr3 = dr3.sel(lat=slice(dataminlat, datamaxlat), lon=slice(360+(dataminlon), 360+(datamaxlon)))
dr4 = dr4.sel(lat=slice(dataminlat, datamaxlat), lon=slice(360+(dataminlon), 360+(datamaxlon)))
#Countour and Chart Settings
#----------------------------------------------------------------------

# Set Contour Plot Parameters Field 1
contour = declarative.ContourPlot()
contour.data = dr
contour.time = timekey
contour.field = field
contour.level = level * units.hPa
contour.linecolor = 'midnightblue'
contour.linestyle = '-'
contour.linewidth = 1
contour.contours = list(range(0, 20000, 3))
contour.plot_units = 'dam'
contour.smooth_contour = 3
contour.clabels = True

#Set Countour Plot Parameters Field 4
contour2 = declarative.ContourPlot()
contour2.data = dr4
contour2.time = timekey
contour2.field = field4
contour2.level = level * units.hPa
contour2.linecolor = 'tab:red'
contour2.plot_units = 'degC'
contour2.linestyle = 'dashed'
contour2.linewidth = 1
contour2.contours = range(-100, 100, 1)
contour2.smooth_contour = 3
contour2.clabels = True

# Plot wind barbs
barb = mpplots.BarbPlot()
barb.data = dr2
barb.time = timekey
barb.field = (field2, field3)
barb.level = level * units.hPa
barb.skip = (3, 3)
barb.plot_units = 'knot'

# Panel for plot with Map features
panel = declarative.MapPanel()
panel.layout = (1, 1, 1)
panel.area = area
#panel.area = 'mx'
panel.projection = 'mer'
panel.layers = ['coastline', 'borders', 'states']
panel.layers_edgecolor = ['grey','grey','grey']
panel.title = FileTitle + str(' ') + str(timekey) + str(' UTC')
panel.plots = [contour, contour2, barb]


#this makes code fail
dr2_wind_speed['wind_speed'] = mpcalc.wind_speed(dr2['uwnd'], dr2['vwnd'])


# Bringing it all together
pc = declarative.PanelContainer()
pc.size = (15, 14)
pc.panels = [panel]

#Save the plot and Fix Name saved
title = 'Surface_Obs{date}.tiff'
pc.save('Surface_Obs.tiff', bbox_inches='tight', dpi= 210)


os.rename('Surface_Obs.tiff',(str(Directory)+str(FileName))+str(timekey)+str('UTC')+str('.tiff'))

pc.show()

Errors, Traceback, and Logs

syntax error, unexpected WORD_WORD, expecting SCAN_ATTR or SCAN_DATASET or SCAN_ERROR
context: Request^ Too Large: 755.4624 Mbytes, max=500.0
syntax error, unexpected WORD_WORD, expecting SCAN_ATTR or SCAN_DATASET or SCAN_ERROR
context: Request^ Too Large: 755.4624 Mbytes, max=500.0
syntax error, unexpected WORD_WORD, expecting SCAN_ATTR or SCAN_DATASET or SCAN_ERROR
context: Request^ Too Large: 755.4624 Mbytes, max=500.0
syntax error, unexpected WORD_WORD, expecting SCAN_ATTR or SCAN_DATASET or SCAN_ERROR
context: Request^ Too Large: 755.4624 Mbytes, max=500.0
syntax error, unexpected WORD_WORD, expecting SCAN_ATTR or SCAN_DATASET or SCAN_ERROR
context: Request^ Too Large: 755.4624 Mbytes, max=500.0
syntax error, unexpected WORD_WORD, expecting SCAN_ATTR or SCAN_DATASET or SCAN_ERROR
context: Request^ Too Large: 755.4624 Mbytes, max=500.0
syntax error, unexpected WORD_WORD, expecting SCAN_ATTR or SCAN_DATASET or SCAN_ERROR
context: Request^ Too Large: 755.4624 Mbytes, max=500.0
---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
Cell In[103], line 129
    126 panel.title = FileTitle + str(' ') + str(timekey) + str(' UTC')
    127 panel.plots = [contour, contour2, barb]
--> 129 dr2_wind_speed['wind_speed'] = mpcalc.wind_speed(dr2['uwnd'], dr2['vwnd'])
    132 # Bringing it all together
    133 pc = declarative.PanelContainer()

File /opt/anaconda3/lib/python3.12/site-packages/metpy/xarray.py:1323, in preprocess_and_wrap.<locals>.decorator.<locals>.wrapper(*args, **kwargs)
   1320     match = tuple(match)
   1322 # Cast all DataArrays to Pint Quantities
-> 1323 _mutate_arguments(bound_args, xr.DataArray, lambda arg, _: arg.metpy.unit_array)
   1325 # Optionally cast all Quantities to their magnitudes
   1326 if to_magnitude:

File /opt/anaconda3/lib/python3.12/site-packages/metpy/units.py:225, in _mutate_arguments(bound_args, check_type, mutate_arg)
    223 for arg_name, arg_val in bound_args.arguments.items():
    224     if isinstance(arg_val, check_type):
--> 225         bound_args.arguments[arg_name] = mutate_arg(arg_val, arg_name)
    227 if isinstance(bound_args.arguments.get('args'), tuple):
    228     bound_args.arguments['args'] = tuple(
    229         mutate_arg(arg_val, '(unnamed)') if isinstance(arg_val, check_type) else arg_val
    230         for arg_val in bound_args.arguments['args'])

File /opt/anaconda3/lib/python3.12/site-packages/metpy/xarray.py:1323, in preprocess_and_wrap.<locals>.decorator.<locals>.wrapper.<locals>.<lambda>(arg, _)
   1320     match = tuple(match)
   1322 # Cast all DataArrays to Pint Quantities
-> 1323 _mutate_arguments(bound_args, xr.DataArray, lambda arg, _: arg.metpy.unit_array)
   1325 # Optionally cast all Quantities to their magnitudes
   1326 if to_magnitude:

File /opt/anaconda3/lib/python3.12/site-packages/metpy/xarray.py:162, in MetPyDataArrayAccessor.unit_array(self)
    152 @property
    153 def unit_array(self):
    154     """Return the data values of this DataArray as a `pint.Quantity`.
    155 
    156     Notes
   (...)
    160     large-sized remote datasets before subsetting!
    161     """
--> 162     if is_quantity(self._data_array.data):
    163         return self._data_array.data
    164     else:

File /opt/anaconda3/lib/python3.12/site-packages/xarray/core/dataarray.py:795, in DataArray.data(self)
    783 @property
    784 def data(self) -> Any:
    785     """
    786     The DataArray's data as an array. The underlying array type
    787     (e.g. dask, sparse, pint) is preserved.
   (...)
    793     DataArray.values
    794     """
--> 795     return self.variable.data

File /opt/anaconda3/lib/python3.12/site-packages/xarray/core/variable.py:474, in Variable.data(self)
    472     return self._data
    473 elif isinstance(self._data, indexing.ExplicitlyIndexed):
--> 474     return self._data.get_duck_array()
    475 else:
    476     return self.values

File /opt/anaconda3/lib/python3.12/site-packages/xarray/core/indexing.py:840, in MemoryCachedArray.get_duck_array(self)
    839 def get_duck_array(self):
--> 840     self._ensure_cached()
    841     return self.array.get_duck_array()

File /opt/anaconda3/lib/python3.12/site-packages/xarray/core/indexing.py:837, in MemoryCachedArray._ensure_cached(self)
    836 def _ensure_cached(self):
--> 837     self.array = as_indexable(self.array.get_duck_array())

File /opt/anaconda3/lib/python3.12/site-packages/xarray/core/indexing.py:794, in CopyOnWriteArray.get_duck_array(self)
    793 def get_duck_array(self):
--> 794     return self.array.get_duck_array()

File /opt/anaconda3/lib/python3.12/site-packages/xarray/core/indexing.py:664, in LazilyIndexedArray.get_duck_array(self)
    659 # self.array[self.key] is now a numpy array when
    660 # self.array is a BackendArray subclass
    661 # and self.key is BasicIndexer((slice(None, None, None),))
    662 # so we need the explicit check for ExplicitlyIndexed
    663 if isinstance(array, ExplicitlyIndexed):
--> 664     array = array.get_duck_array()
    665 return _wrap_numpy_scalars(array)

File /opt/anaconda3/lib/python3.12/site-packages/xarray/coding/variables.py:81, in _ElementwiseFunctionArray.get_duck_array(self)
     80 def get_duck_array(self):
---> 81     return self.func(self.array.get_duck_array())

File /opt/anaconda3/lib/python3.12/site-packages/xarray/core/indexing.py:657, in LazilyIndexedArray.get_duck_array(self)
    653     array = apply_indexer(self.array, self.key)
    654 else:
    655     # If the array is not an ExplicitlyIndexedNDArrayMixin,
    656     # it may wrap a BackendArray so use its __getitem__
--> 657     array = self.array[self.key]
    659 # self.array[self.key] is now a numpy array when
    660 # self.array is a BackendArray subclass
    661 # and self.key is BasicIndexer((slice(None, None, None),))
    662 # so we need the explicit check for ExplicitlyIndexed
    663 if isinstance(array, ExplicitlyIndexed):

File /opt/anaconda3/lib/python3.12/site-packages/xarray/backends/netCDF4_.py:103, in NetCDF4ArrayWrapper.__getitem__(self, key)
    102 def __getitem__(self, key):
--> 103     return indexing.explicit_indexing_adapter(
    104         key, self.shape, indexing.IndexingSupport.OUTER, self._getitem
    105     )

File /opt/anaconda3/lib/python3.12/site-packages/xarray/core/indexing.py:1018, in explicit_indexing_adapter(key, shape, indexing_support, raw_indexing_method)
    996 """Support explicit indexing by delegating to a raw indexing method.
    997 
    998 Outer and/or vectorized indexers are supported by indexing a second time
   (...)
   1015 Indexing result, in the form of a duck numpy-array.
   1016 """
   1017 raw_key, numpy_indices = decompose_indexer(key, shape, indexing_support)
-> 1018 result = raw_indexing_method(raw_key.tuple)
   1019 if numpy_indices.tuple:
   1020     # index the loaded np.ndarray
   1021     indexable = NumpyIndexingAdapter(result)

File /opt/anaconda3/lib/python3.12/site-packages/xarray/backends/netCDF4_.py:116, in NetCDF4ArrayWrapper._getitem(self, key)
    114     with self.datastore.lock:
    115         original_array = self.get_array(needs_lock=False)
--> 116         array = getitem(original_array, key)
    117 except IndexError as err:
    118     # Catch IndexError in netCDF4 and return a more informative
    119     # error message.  This is most often called when an unsorted
    120     # indexer is used before the data is loaded from disk.
    121     msg = (
    122         "The indexing operation you are attempting to perform "
    123         "is not valid on netCDF4.Variable object. Try loading "
    124         "your data into memory first by calling .load()."
    125     )

File /opt/anaconda3/lib/python3.12/site-packages/xarray/backends/common.py:249, in robust_getitem(array, key, catch, max_retries, initial_delay)
    247 for n in range(max_retries + 1):
    248     try:
--> 249         return array[key]
    250     except catch:
    251         if n == max_retries:

File src/netCDF4/_netCDF4.pyx:4972, in netCDF4._netCDF4.Variable.__getitem__()

File src/netCDF4/_netCDF4.pyx:5930, in netCDF4._netCDF4.Variable._get()

File src/netCDF4/_netCDF4.pyx:2034, in netCDF4._netCDF4._ensure_nc_success()

RuntimeError: NetCDF: Authorization failure

[dopplershift]

If you look at the top of the errors, see: context: Request^ Too Large: 755.4624 Mbytes, max=500.0. That implies that too much data is being requested from the server, which limits individual requests to 500MB. These are likely caused by the calls to parse_cf() with no variable names, which requests all variables from the server. You can reduce the data downloaded by passing specific field names to the call to parse_cf(), such as:

dr = xr.open_dataset(path).metpy.parse_cf(field)
dr2 = xr.open_dataset(path2).metpy.parse_cf('uwnd', 'vwnd')

[graffmet11]

@dopplershift,

Thanks for the assist! I tried the parse_cf angle, but didn't have any luck. It put me on the right path though. Re-examining the file structure, I realized these were different than some of the over nc files that users probably plug in (such as today's GFS data). Each file has one weather parameter (like wind), but has coordinates of levels, times, lat, and lon. I ended up doing slices by level, lat, lon, and time. It looks like that cut down on the data I was pulling extensively. For anyone else who might have this issue, I'm pasting what that looks like below:

#Set Remaining File time
time = pd.to_datetime(ds)
timekey = time[timepos] #Set Time Here as Timekey (Number Equals timestep value from NetCDF file)
File_Time = timekey.strftime('%m-%d-%Y%:%H')
File_Time_Var = timekey.strftime('%Y-%m-%dT%H:%S:00')

#Open data sources
dr = xr.open_dataset(path).metpy.parse_cf().sel(level=slice(pres_level), time=slice(File_Time_Var, File_Time_Var), lat=slice(dataminlat, datamaxlat), lon=slice(360+(dataminlon), 360+(datamaxlon)))