-
Notifications
You must be signed in to change notification settings - Fork 78
/
Copy pathcrs.py
172 lines (140 loc) · 5.88 KB
/
crs.py
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
"""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
from pyproj import CRS
import xarray as xr
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"]:
# Earth assumed spherical with radius specified (in m) by data producer
# TODO: Why is urma and rtma using this value?
a = 6371200.0
b = 6371200.0
elif shapeOfTheEarth == 6:
# Earth assumed spherical with radius of 6,371,229.0 m
a = 6_371_229
b = 6_371_229
# 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 == "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
"""