Unidata
gdxarray() returns unusable objects from "thinned GFS" files
Here's example 1:
from datetime import datetime
from metpy.io import GempakGrid
import metpy
print(metpy.__version__)
gem_file_name = '/ldmdata/gempak/model/gfs/21092206_thin.gem'
gem_file = GempakGrid(gem_file_name)
plot_time = datetime(2021, 9, 25, 6)
ht = gem_file.gdxarray(parameter='HGHT', date_time=plot_time, level=500)[0]
print(ht)
And the result:
1.1.0.post39+gd43d5eb02.d20210825
Traceback (most recent call last):
File "/home/decker/classes/met433/new_labs/lab3q3.py", line 13, in <module>
print(ht)
File "/home/decker/local/miniconda3/envs/devel21/lib/python3.9/site-packages/xarray/core/common.py", line 144, in __repr__
return formatting.array_repr(self)
File "/home/decker/local/miniconda3/envs/devel21/lib/python3.9/site-packages/xarray/core/formatting.py", line 529, in array_repr
summary.append(repr(arr.coords))
File "/home/decker/local/miniconda3/envs/devel21/lib/python3.9/site-packages/xarray/core/coordinates.py", line 79, in __repr__
return formatting.coords_repr(self)
File "/home/decker/local/miniconda3/envs/devel21/lib/python3.9/site-packages/xarray/core/formatting.py", line 422, in coords_repr
return _mapping_repr(
File "/home/decker/local/miniconda3/envs/devel21/lib/python3.9/site-packages/xarray/core/formatting.py", line 397, in _mapping_repr
summary += [summarizer(k, v, col_width) for k, v in mapping.items()]
File "/home/decker/local/miniconda3/envs/devel21/lib/python3.9/site-packages/xarray/core/formatting.py", line 397, in <listcomp>
summary += [summarizer(k, v, col_width) for k, v in mapping.items()]
File "/home/decker/local/miniconda3/envs/devel21/lib/python3.9/site-packages/xarray/core/formatting.py", line 334, in summarize_coord
return summarize_variable(name, var.variable, col_width, marker)
File "/home/decker/local/miniconda3/envs/devel21/lib/python3.9/site-packages/xarray/core/formatting.py", line 292, in summarize_variable
values_str = inline_variable_array_repr(var, values_width)
File "/home/decker/local/miniconda3/envs/devel21/lib/python3.9/site-packages/xarray/core/formatting.py", line 260, in inline_variable_array_repr
return format_array_flat(var, max_width)
File "/home/decker/local/miniconda3/envs/devel21/lib/python3.9/site-packages/xarray/core/formatting.py", line 178, in format_array_flat
relevant_front_items = format_items(
File "/home/decker/local/miniconda3/envs/devel21/lib/python3.9/site-packages/xarray/core/formatting.py", line 165, in format_items
formatted = [format_item(xi, timedelta_format) for xi in x]
File "/home/decker/local/miniconda3/envs/devel21/lib/python3.9/site-packages/xarray/core/formatting.py", line 165, in <listcomp>
formatted = [format_item(xi, timedelta_format) for xi in x]
File "/home/decker/local/miniconda3/envs/devel21/lib/python3.9/site-packages/xarray/core/formatting.py", line 148, in format_item
return str(x)
File "/home/decker/src/git_repos/metpy/src/metpy/plots/mapping.py", line 99, in __str__
return 'Projection: ' + self._attrs['grid_mapping_name']
KeyError: 'grid_mapping_name'
Example 2:
from datetime import datetime
from metpy.io import GempakGrid
from metpy.plots import ContourPlot, MapPanel, PanelContainer
import metpy
print(metpy.__version__)
gem_file_name = '/ldmdata/gempak/model/gfs/21092206_thin.gem'
gem_file = GempakGrid(gem_file_name)
plot_time = datetime(2021, 9, 25, 6)
ht = gem_file.gdxarray(parameter='HGHT', date_time=plot_time, level=500)[0]
hts = ContourPlot()
hts.data = ht
hts.contours = list(range(0, 6000, 60))
mp = MapPanel()
mp.plots = [hts]
pc = PanelContainer()
pc.size = (10,10)
pc.panels = [mp]
pc.show()
And the output:
1.1.0.post39+gd43d5eb02.d20210825
Traceback (most recent call last):
File "/home/decker/classes/met433/new_labs/lab3q3.py", line 26, 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 1137, in draw
self._build()
File "/home/decker/src/git_repos/metpy/src/metpy/plots/declarative.py", line 1302, in _build
kwargs = self.parent.plot_kwargs
File "/home/decker/src/git_repos/metpy/src/metpy/plots/declarative.py", line 814, in plot_kwargs
dataproj = self.plots[0].griddata.metpy.cartopy_crs
File "/home/decker/src/git_repos/metpy/src/metpy/xarray.py", line 242, in cartopy_crs
return self.crs.to_cartopy()
File "/home/decker/src/git_repos/metpy/src/metpy/plots/mapping.py", line 79, in to_cartopy
proj_name = self._attrs['grid_mapping_name']
KeyError: 'grid_mapping_name'
My expectation is that I get a printout or plot of the data.
Same error on my end. I can't manually specify the grid_mapping_name in the GempakGrid() call or after the fact since the output is a list. If I slice the list I get the error before the mepty.assign_crs() method acts.
GFS and JRA55 .gem files unusable. Is this a metpy bug or user error on my end (i.e. creating .gem with improper/lacking projection information).
It's entirely possible, even likely, that there are some untested corner cases.
@icbeckley Can you attach a sample file (assuming it's not too large) to this issue so we can reproduce locally? You might need to zip/gzip it in order to attach it.
Having trouble attaching, but I've opened the file for sharing on google drive. It is a 1 deg GFS analysis with a few added omega fields. gdinfo() reveals it is in CED.
This should work: https://drive.google.com/file/d/1BaOJRMR0kEODXmabjbUvWqbD7IZJVg7l/view?usp=sharing
Thanks in advance for your time!
So I finally dug into this one. The problem is that CED (Cylindrical Equidistant--also same for MCD) doesn't seem to have an equivalent representation in CF Metadata. Because of course.
Under the hood the projection is generated using PROJ, which we then ask to convert to CF, which then doesn't have a name.
@sgdecker AFAICT, CED has "native" coordinates that are definitively NOT lat/lon values, so I don't think so? And if so, what does "modified" CED mean in that case?
I'd say based on what PROJ lists for equidistant cylindrical the answer seems to be "no" as far as lat/lon is concerned. Probably possible to use PlateCarree in Cartopy, but that doesn't solve the CF metadata issue.
@dopplershift I haven't had a chance to see what numbers the gdxarray method is coming up with internally, but I suspect they are radians based on this GEMPAK function. It does look like GEMPAK in principle allows for a more general CED grid than the PROJ library (GEMPAK also allows for a rotation angle, in addition to the central latitude and central longitude), but in practice, grids like the GFS in question are not rotated, and have a central latitude of 0.
Regarding the MCD grid, that is GEMPAK notation for when +lat_ts (in PROJ-speak) is non-zero. I'm not aware of any grid that is natively MCD, although you can certainly plot any grid in MCD in GEMPAK.
My (admitedly untested) conclusion is that there is a PROJ specification that corresponds to the GFS grids, but a transformation from radians to degrees would need to happen to arrive at a grid that can be expressed in CF-speak.