Unidata
gdxarray() should assign units to x and y coordinates
Is there a way to construct the proper Dataset with the current MetPy functionality so my attempt to plot wind barbs works?
Here's my attempt:
from datetime import datetime
import xarray as xr
import metpy
from metpy.units import units
from metpy.io import GempakGrid
from metpy.plots.declarative import ContourPlot, FilledContourPlot, BarbPlot, MapPanel, PanelContainer
print(metpy.__version__)
MODEL = '/ldmdata/gempak/model/'
gem_file_name = MODEL + 'nam/21092012_nam211.gem'
gem_file = GempakGrid(gem_file_name)
plot_time = datetime(2021, 9, 20, 12)
u500 = gem_file.gdxarray(parameter='UREL', date_time=plot_time, level=500)[0]
v500 = gem_file.gdxarray(parameter='VREL', date_time=plot_time, level=500)[0]
u500 = u500 * units('m/s')
v500 = v500 * units('m/s')
wind = xr.merge([u500, v500])
barbs = BarbPlot()
barbs.data = wind
barbs.time = plot_time
barbs.field = ['urel', 'vrel']
barbs.earth_relative = False # I get a plot w/o this line, but obviously the winds aren't rotated correctly
barbs.skip = (3, 3)
barbs.plot_units = 'knot'
panel = MapPanel()
panel.area = [-120, -74, 22, 55]
panel.projection = 'lcc'
panel.layers = ['states', 'coastline', 'borders']
panel.title = f'500-mb Winds at {plot_time}'
panel.plots = [barbs]
pc = PanelContainer()
pc.size = (15, 15)
pc.panels = [panel]
pc.show()
The result:
1.1.0.post39+gd43d5eb02.d20210825
Traceback (most recent call last):
File "/home/decker/classes/met433/new_labs/barb_issue.py", line 42, in <module>
pc.show()
File "/home/decker/src/git_repos/metpy/src/metpy/plots/declarative.py", line 589, in show
self.draw()
File "/home/decker/src/git_repos/metpy/src/metpy/plots/declarative.py", line 576, in draw
panel.draw()
File "/home/decker/src/git_repos/metpy/src/metpy/plots/declarative.py", line 846, in draw
p.draw()
File "/home/decker/src/git_repos/metpy/src/metpy/plots/declarative.py", line 1480, in draw
self._build()
File "/home/decker/src/git_repos/metpy/src/metpy/plots/declarative.py", line 1504, in _build
x_like, y_like, u, v = self.plotdata
File "/home/decker/src/git_repos/metpy/src/metpy/plots/declarative.py", line 1465, in plotdata
if 'degree' in x.units:
File "/home/decker/local/miniconda3/envs/devel21/lib/python3.9/site-packages/xarray/core/common.py", line 239, in __getattr__
raise AttributeError(
AttributeError: 'DataArray' object has no attribute 'units'
There should be a way of doing this either through the metpy.dequantify(), which will take units attached to a variable and make the unit an attribute within the xarray object.
Try the following for the u500 and v500...
u500 = (u500 * units('m/s')).metpy.dequantify()
v500 = (v500 * units('m/s')).metpy.dequantify()
Here is the Units portion of the xarray tutorial that might also help.
https://unidata.github.io/MetPy/latest/tutorials/xarray_tutorial.html#units
Thanks. I had tried metpy.quantify(), which didn't appear to do anything (presumably because urel/vrel are not CF-kosher), and metpy.quantify('m/s'), which crashed IIRC, before settling on the unit multiplication approach. It didn't occur to me to try dequantify().
Unfortunately, that change had no effect!
Because the problem only occurs with barbs.earth_relative = False, I was thinking somehow the units for longitude were screwed up (since rotating the winds properly on this grid needs longitude info), and that's the missing attribute, but I haven't been able to confirm that yet, and if it is the problem, I haven't figured out how to fix it.
Ah, yes, it could be that the x/y (or lon/lat) values don't have the units needed.
You could try attaching the units to the coordinate variables. Assuming they are 'x' and 'y' and the units are meters...
wind['x'] = wind['x'] * units('meters')
wind['y'] = wind['y'] * units('meters')
@kgoebber thanks; I think we are getting somewhere! I was fixated on latitude/longitude, but the error traceback does mention x after all. Here is a revised program:
from datetime import datetime
import xarray as xr
import metpy
from metpy.units import units
from metpy.io import GempakGrid
from metpy.plots.declarative import ContourPlot, FilledContourPlot, BarbPlot, MapPanel, PanelContainer
print(metpy.__version__)
MODEL = '/ldmdata/gempak/model/'
gem_file_name = MODEL + 'nam/21092012_nam211.gem'
gem_file = GempakGrid(gem_file_name)
plot_time = datetime(2021, 9, 20, 12)
u500 = gem_file.gdxarray(parameter='UREL', date_time=plot_time, level=500)[0]
v500 = gem_file.gdxarray(parameter='VREL', date_time=plot_time, level=500)[0]
u500 = u500 * units('m/s')
v500 = v500 * units('m/s')
wind = xr.merge([u500, v500]).metpy.assign_latitude_longitude()
print(wind.longitude)
print(wind.x)
wind['x'] = wind['x'] * units('meters')
wind['y'] = wind['y'] * units('meters')
print(wind.x)
barbs = BarbPlot()
barbs.data = wind
barbs.time = plot_time
barbs.field = ['urel', 'vrel']
barbs.earth_relative = False
barbs.skip = (3, 3)
barbs.plot_units = 'knot'
panel = MapPanel()
panel.area = [-120, -74, 22, 55]
panel.projection = 'lcc'
panel.layers = ['states', 'coastline', 'borders']
panel.title = f'500-mb Winds at {plot_time}'
panel.plots = [barbs]
pc = PanelContainer()
pc.size = (15, 15)
pc.panels = [panel]
pc.show()
And the output:
1.1.0.post39+gd43d5eb02.d20210825
<xarray.DataArray 'longitude' (y: 65, x: 93)>
array([[-133.45899757, -132.75739652, -132.05371082, ..., -66.54178831,
-65.81553154, -65.09097213],
[-133.66429679, -132.95935346, -132.25229366, ..., -66.386047 ,
-65.65606298, -64.92780342],
[-133.87173891, -133.16342482, -132.4529618 , ..., -66.22861748,
-65.49486987, -64.76287439],
...,
[-151.97501727, -151.00076902, -150.02034207, ..., -52.20690183,
-51.15677614, -50.11207458],
[-152.41212945, -151.43223072, -150.44602669, ..., -51.86040261,
-50.80299078, -49.75112582],
[-152.85558563, -151.86999524, -150.87796976, ..., -51.50844159,
-50.44365796, -49.38454751]])
Coordinates:
* x (x) float32 -4.226e+06 -4.145e+06 ... 3.17e+06 3.251e+06
* y (y) float32 2.035e+06 2.117e+06 2.198e+06 ... 7.155e+06 7.237e+06
metpy_crs object Projection: lambert_conformal_conic
latitude (y, x) float64 12.19 12.39 12.58 12.77 ... 57.68 57.49 57.29
longitude (y, x) float64 -133.5 -132.8 -132.1 ... -51.51 -50.44 -49.38
Attributes:
units: degrees_east
standard_name: longitude
<xarray.DataArray 'x' (x: 93)>
array([-4.226088e+06, -4.144817e+06, -4.063546e+06, -3.982275e+06,
-3.901004e+06, -3.819733e+06, -3.738462e+06, -3.657191e+06,
-3.575920e+06, -3.494649e+06, -3.413378e+06, -3.332107e+06,
-3.250836e+06, -3.169565e+06, -3.088294e+06, -3.007023e+06,
-2.925752e+06, -2.844481e+06, -2.763210e+06, -2.681939e+06,
-2.600668e+06, -2.519397e+06, -2.438126e+06, -2.356855e+06,
-2.275584e+06, -2.194313e+06, -2.113042e+06, -2.031771e+06,
-1.950500e+06, -1.869229e+06, -1.787958e+06, -1.706687e+06,
-1.625416e+06, -1.544145e+06, -1.462874e+06, -1.381603e+06,
-1.300332e+06, -1.219061e+06, -1.137790e+06, -1.056519e+06,
-9.752478e+05, -8.939768e+05, -8.127058e+05, -7.314348e+05,
-6.501638e+05, -5.688928e+05, -4.876218e+05, -4.063508e+05,
-3.250798e+05, -2.438088e+05, -1.625378e+05, -8.126677e+04,
4.234092e+00, 8.127523e+04, 1.625462e+05, 2.438172e+05,
3.250882e+05, 4.063592e+05, 4.876302e+05, 5.689012e+05,
6.501722e+05, 7.314432e+05, 8.127142e+05, 8.939852e+05,
9.752562e+05, 1.056527e+06, 1.137798e+06, 1.219069e+06,
1.300340e+06, 1.381611e+06, 1.462882e+06, 1.544153e+06,
1.625424e+06, 1.706695e+06, 1.787966e+06, 1.869237e+06,
1.950508e+06, 2.031779e+06, 2.113050e+06, 2.194321e+06,
2.275592e+06, 2.356863e+06, 2.438134e+06, 2.519405e+06,
2.600676e+06, 2.681947e+06, 2.763218e+06, 2.844489e+06,
2.925760e+06, 3.007031e+06, 3.088302e+06, 3.169573e+06,
3.250844e+06], dtype=float32)
Coordinates:
* x (x) float32 -4.226e+06 -4.145e+06 ... 3.17e+06 3.251e+06
metpy_crs object Projection: lambert_conformal_conic
Attributes:
_metpy_axis: x
<xarray.DataArray 'x' (x: 93)>
array([-4.226088e+06, -4.144817e+06, -4.063546e+06, -3.982275e+06,
-3.901004e+06, -3.819733e+06, -3.738462e+06, -3.657191e+06,
-3.575920e+06, -3.494649e+06, -3.413378e+06, -3.332107e+06,
-3.250836e+06, -3.169565e+06, -3.088294e+06, -3.007023e+06,
-2.925752e+06, -2.844481e+06, -2.763210e+06, -2.681939e+06,
-2.600668e+06, -2.519397e+06, -2.438126e+06, -2.356855e+06,
-2.275584e+06, -2.194313e+06, -2.113042e+06, -2.031771e+06,
-1.950500e+06, -1.869229e+06, -1.787958e+06, -1.706687e+06,
-1.625416e+06, -1.544145e+06, -1.462874e+06, -1.381603e+06,
-1.300332e+06, -1.219061e+06, -1.137790e+06, -1.056519e+06,
-9.752478e+05, -8.939768e+05, -8.127058e+05, -7.314348e+05,
-6.501638e+05, -5.688928e+05, -4.876218e+05, -4.063508e+05,
-3.250798e+05, -2.438088e+05, -1.625378e+05, -8.126677e+04,
4.234092e+00, 8.127523e+04, 1.625462e+05, 2.438172e+05,
3.250882e+05, 4.063592e+05, 4.876302e+05, 5.689012e+05,
6.501722e+05, 7.314432e+05, 8.127142e+05, 8.939852e+05,
9.752562e+05, 1.056527e+06, 1.137798e+06, 1.219069e+06,
1.300340e+06, 1.381611e+06, 1.462882e+06, 1.544153e+06,
1.625424e+06, 1.706695e+06, 1.787966e+06, 1.869237e+06,
1.950508e+06, 2.031779e+06, 2.113050e+06, 2.194321e+06,
2.275592e+06, 2.356863e+06, 2.438134e+06, 2.519405e+06,
2.600676e+06, 2.681947e+06, 2.763218e+06, 2.844489e+06,
2.925760e+06, 3.007031e+06, 3.088302e+06, 3.169573e+06,
3.250844e+06], dtype=float32)
Coordinates:
* x (x) float32 -4.226e+06 -4.145e+06 ... 3.17e+06 3.251e+06
metpy_crs object Projection: lambert_conformal_conic
Traceback (most recent call last):
File "/home/decker/classes/met433/new_labs/barb_issue.py", line 47, in <module>
pc.show()
File "/home/decker/src/git_repos/metpy/src/metpy/plots/declarative.py", line 589, in show
self.draw()
File "/home/decker/src/git_repos/metpy/src/metpy/plots/declarative.py", line 576, in draw
panel.draw()
File "/home/decker/src/git_repos/metpy/src/metpy/plots/declarative.py", line 846, in draw
p.draw()
File "/home/decker/src/git_repos/metpy/src/metpy/plots/declarative.py", line 1480, in draw
self._build()
File "/home/decker/src/git_repos/metpy/src/metpy/plots/declarative.py", line 1504, in _build
x_like, y_like, u, v = self.plotdata
File "/home/decker/src/git_repos/metpy/src/metpy/plots/declarative.py", line 1465, in plotdata
if 'degree' in x.units:
File "/home/decker/local/miniconda3/envs/devel21/lib/python3.9/site-packages/xarray/core/common.py", line 239, in __getattr__
raise AttributeError(
AttributeError: 'DataArray' object has no attribute 'units'
My observations:
- For some reason, trying to attach units to
xandyvia multiplication doesn't seem to do anything. -
assign_latitude_longitude()does seem to bring in proper units for latitude and longitude.
When I get a chance, I'll look at the code in assign_latitude_longitude() to see how it is attaching units. Perhaps the same approach will successfully add meters to the x and y coordinates, and I'll get a plot.
My observations:
- For some reason, trying to attach units to
xandyvia multiplication doesn't seem to do anything.assign_latitude_longitude()does seem to bring in proper units for latitude and longitude.When I get a chance, I'll look at the code in
assign_latitude_longitude()to see how it is attaching units. Perhaps the same approach will successfully add meters to thexandycoordinates, and I'll get a plot.
Dimension coordinates (like x and y) do not support quantity arrays yet (xref https://github.com/pydata/xarray/projects/1), so multiplying by units doesn't do anything. Instead, you need to assign the units attribute:
wind['x'].attrs['units'] = 'meters'
wind['y'].attrs['units'] = 'meters'
You'll find this pattern (units-on-attribute for coordinates, units-in-quantity-array for data) often in MetPy's accessor code. Hopefully the in-progress changes upstream will make this easier/more intuitive.
Aha, that's the incantation I was looking for! I now have a plot. Thanks, everyone!
So the actual bug/feature request is that gdxarray should do this automatically.