forked from blaylockbk/Herbie
-
Notifications
You must be signed in to change notification settings - Fork 0
Expand file tree
/
Copy pathcrs.py
More file actions
195 lines (158 loc) · 6.53 KB
/
Copy pathcrs.py
File metadata and controls
195 lines (158 loc) · 6.53 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
"""CF convention coordinate reference system (CRS).
Check out the notebook `check_pygrib_vs_herbie_crs_extraction.ipynb`
to test how Herbie and pygrib are extracting the CRS information.
"""
from typing import Any
import xarray as xr
from pyproj import CRS
def get_cf_crs(
ds: "xr.Dataset", variable: str | None = None, _return_projparams=False
) -> dict[str, Any]:
"""
Extract the CF coordinate reference system (CRS) from a cfgrib xarray dataset.
Note:
I originally used pygrib to do this, but it was hard to maintain an
additional grib package dependency. I had issues with pygrib after
it was updated to support Numpy version 2, so thought it would be
best to code this in Herbie. This may be incomplete.
"""
# Assume the first variable in the Dataset has the same grid crs
# as all other variables in the Dataset.
if variable is None:
variable = next(iter(ds.data_vars))
da = ds[variable]
# Shape of the Earth reference system
# https://codes.ecmwf.int/grib/format/grib2/ctables/3/2/
shapeOfTheEarth = da.GRIB_shapeOfTheEarth
if shapeOfTheEarth == 0:
# Earth assumed spherical with radius = 6 367 470.0 m
a = 6_367_470
b = 6_367_470
elif shapeOfTheEarth == 1 and ds.attrs["model"] == "graphcast":
# Earth assumed spherical with radius specified (in m) by data producer
# TODO: Why is model='graphcast' using this value?
a = 4326.0
b = 4326.0
elif shapeOfTheEarth == 1 and ds.attrs["model"] in {
"urma",
"rtma",
"nbm",
"nbmqmd",
}:
# Earth assumed spherical with radius specified (in m) by data producer
# TODO: Why is urma, rtma, and nbm using this value?
a = 6_371_200.0
b = 6_371_200.0
elif shapeOfTheEarth == 6: # i.e., valid for HRRR model
# Earth assumed spherical with radius of 6,371,229.0 m
a = 6_371_229
b = 6_371_229
else:
raise NotImplementedError(
f"{shapeOfTheEarth=} for {ds.attrs['model']=} is not implemented. "
"Please open an issue or pull request on GitHub."
)
# Grid type definition
# https://codes.ecmwf.int/grib/format/grib2/ctables/3/1/
if da.GRIB_gridType == "lambert":
projparams = {"proj": "lcc"}
projparams["a"] = a
projparams["b"] = b
projparams["lon_0"] = da.GRIB_LoVInDegrees
projparams["lat_0"] = da.GRIB_LaDInDegrees
projparams["lat_1"] = da.GRIB_Latin1InDegrees
projparams["lat_2"] = da.GRIB_Latin2InDegrees
elif da.GRIB_gridType == "regular_ll":
projparams = {"proj": "longlat"}
projparams["a"] = a
projparams["b"] = b
elif da.GRIB_gridType == "rotated_ll":
projparams = {"proj": "ob_tran"}
projparams["o_proj"] = "longlat"
projparams["a"] = a
projparams["b"] = b
projparams["lon_0"] = da.GRIB_longitudeOfSouthernPoleInDegrees
projparams["o_lon_p"] = 0.0
projparams["o_lat_p"] = -da.GRIB_latitudeOfSouthernPoleInDegrees
elif da.GRIB_gridType == "regular_gg":
projparams = {"proj": "longlat"}
projparams["a"] = a
projparams["b"] = b
elif da.GRIB_gridType == "polar_stereographic":
projparams = {"proj": "stere"}
projparams["a"] = a
projparams["b"] = b
projparams["lat_ts"] = da.GRIB_LaDInDegrees
projparams["lat_0"] = 90
projparams["lon_0"] = da.GRIB_orientationOfTheGridInDegrees
else:
raise NotImplementedError(f"gridType {da.GRIB_gridType} is not implemented.")
if _return_projparams:
return projparams
else:
return CRS(projparams).to_cf()
"""
Look at how pygrib parses with this...
with pygrib.open(str(ds.local_grib)) as grb:
msg = grb.message(1)
print(msg.projparams)
Also, look for clues by dumping all keys with:
grib_dump -j <filename.grib2> > filedump.json
-----------------------------------------
Model: HRRR
pygrib {'a': 6371229, 'b': 6371229, 'lat_0': 38.5, 'lat_1': 38.5, 'lat_2': 38.5, 'lon_0': 262.5, 'proj': 'lcc'}
Herbie {'a': 6371229, 'b': 6371229, 'lat_0': 38.5, 'lat_1': 38.5, 'lat_2': 38.5, 'lon_0': 262.5, 'proj': 'lcc'}
equal= True
Model: HRRRAK
pygrib {'a': 6371229, 'b': 6371229, 'lat_0': 90.0, 'lat_ts': 60.0, 'lon_0': 225.0, 'proj': 'stere'}
Herbie {'a': 6371229, 'b': 6371229, 'lat_0': 90, 'lat_ts': 60.0, 'lon_0': 225.0, 'proj': 'stere'}
equal= True
Model: GFS
pygrib {'a': 6371229, 'b': 6371229, 'proj': 'longlat'}
Herbie {'a': 6371229, 'b': 6371229, 'proj': 'longlat'}
equal= True
Model: GRAPHCAST
pygrib {'a': 4326.0, 'b': 4326.0, 'proj': 'longlat'}
Herbie {'a': 4326.0, 'b': 4326.0, 'proj': 'longlat'}
equal= True
Model: CFS
pygrib {'a': 6371229, 'b': 6371229, 'proj': 'longlat'}
Herbie {'a': 6371229, 'b': 6371229, 'proj': 'longlat'}
equal= True
Model: IFS
pygrib {'a': 6371229, 'b': 6371229, 'proj': 'longlat'}
Herbie {'a': 6371229, 'b': 6371229, 'proj': 'longlat'}
equal= True
Model: AIFS
pygrib {'a': 6371229, 'b': 6371229, 'proj': 'longlat'}
Herbie {'a': 6371229, 'b': 6371229, 'proj': 'longlat'}
equal= True
Model: GEFS
pygrib {'a': 6371229, 'b': 6371229, 'proj': 'longlat'}
Herbie {'a': 6371229, 'b': 6371229, 'proj': 'longlat'}
equal= True
Model: GEFS
pygrib {'a': 6371229, 'b': 6371229, 'proj': 'longlat'}
Herbie {'a': 6371229, 'b': 6371229, 'proj': 'longlat'}
equal= True
Model: HAFSA
pygrib {'a': 6371229, 'b': 6371229, 'proj': 'longlat'}
Herbie {'a': 6371229, 'b': 6371229, 'proj': 'longlat'}
equal= True
Model: HREF
pygrib {'a': 6371229, 'b': 6371229, 'lat_0': 25.0, 'lat_1': 25.0, 'lat_2': 25.0, 'lon_0': 265.0, 'proj': 'lcc'}
Herbie {'a': 6371229, 'b': 6371229, 'lat_0': 25.0, 'lat_1': 25.0, 'lat_2': 25.0, 'lon_0': 265.0, 'proj': 'lcc'}
equal= True
Model: NAM
pygrib {'a': 6371229, 'b': 6371229, 'lat_0': 38.5, 'lat_1': 38.5, 'lat_2': 38.5, 'lon_0': 262.5, 'proj': 'lcc'}
Herbie {'a': 6371229, 'b': 6371229, 'lat_0': 38.5, 'lat_1': 38.5, 'lat_2': 38.5, 'lon_0': 262.5, 'proj': 'lcc'}
equal= True
Model: URMA
pygrib {'a': 6371200.0, 'b': 6371200.0, 'lat_0': 25.0, 'lat_1': 25.0, 'lat_2': 25.0, 'lon_0': 265.0, 'proj': 'lcc'}
Herbie {'a': 6371200.0, 'b': 6371200.0, 'lat_0': 25.0, 'lat_1': 25.0, 'lat_2': 25.0, 'lon_0': 265.0, 'proj': 'lcc'}
equal= True
Model: RTMA
pygrib {'a': 6371200.0, 'b': 6371200.0, 'lat_0': 25.0, 'lat_1': 25.0, 'lat_2': 25.0, 'lon_0': 265.0, 'proj': 'lcc'}
Herbie {'a': 6371200.0, 'b': 6371200.0, 'lat_0': 25.0, 'lat_1': 25.0, 'lat_2': 25.0, 'lon_0': 265.0, 'proj': 'lcc'}
equal= True
"""