-
Notifications
You must be signed in to change notification settings - Fork 43
/
_save_rules.py
1430 lines (1171 loc) · 54.8 KB
/
_save_rules.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
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
# (C) British Crown Copyright 2010 - 2018, Met Office
#
# This file is part of iris-grib.
#
# iris-grib is free software: you can redistribute it and/or modify it under
# the terms of the GNU Lesser General Public License as published by the
# Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# iris-grib is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with iris-grib. If not, see <http://www.gnu.org/licenses/>.
"""
Grib save implementation.
:mod:`iris_grib._save_rules` is a private module with no public API.
It is invoked from :meth:`iris_grib.save_grib2`.
"""
from __future__ import (absolute_import, division, print_function)
from six.moves import (filter, input, map, range, zip) # noqa
import warnings
import cf_units
import gribapi
import numpy as np
import numpy.ma as ma
import cartopy.crs as ccrs
import iris
from iris.aux_factory import HybridHeightFactory, HybridPressureFactory
from iris.coord_systems import (GeogCS, RotatedGeogCS, Mercator,
TransverseMercator, LambertConformal)
import iris.exceptions
from ._iris_mercator_support import confirm_extended_mercator_supported
from . import grib_phenom_translation as gptx
from ._load_convert import (_STATISTIC_TYPE_NAMES, _TIME_RANGE_UNITS)
from iris.util import is_regular, regular_step
# Invert code tables from :mod:`iris_grib._load_convert`.
_STATISTIC_TYPE_NAMES = {val: key for key, val in
_STATISTIC_TYPE_NAMES.items()}
_TIME_RANGE_UNITS = {val: key for key, val in _TIME_RANGE_UNITS.items()}
def fixup_float32_as_int32(value):
"""
Workaround for use when the ECMWF GRIB API treats an IEEE 32-bit
floating-point value as a signed, 4-byte integer.
Returns the integer value which will result in the on-disk
representation corresponding to the IEEE 32-bit floating-point
value.
"""
value_as_float32 = np.array(value, dtype='f4')
value_as_uint32 = value_as_float32.view(dtype='u4')
if value_as_uint32 >= 0x80000000:
# Convert from two's-complement to sign-and-magnitude.
# NB. Because of the silly representation of negative
# integers in GRIB2, there is no value we can pass to
# grib_set that will result in the bit pattern 0x80000000.
# But since that bit pattern corresponds to a floating
# point value of negative-zero, we can safely treat it as
# positive-zero instead.
value_as_grib_int = 0x80000000 - int(value_as_uint32)
else:
value_as_grib_int = int(value_as_uint32)
return value_as_grib_int
def fixup_int32_as_uint32(value):
"""
Workaround for use when the ECMWF GRIB API treats a signed, 4-byte
integer value as an unsigned, 4-byte integer.
Returns the unsigned integer value which will result in the on-disk
representation corresponding to the signed, 4-byte integer value.
"""
value = int(value)
if -0x7fffffff <= value <= 0x7fffffff:
if value < 0:
# Convert from two's-complement to sign-and-magnitude.
value = 0x80000000 - value
else:
msg = '{} out of range -2147483647 to 2147483647.'.format(value)
raise ValueError(msg)
return value
def ensure_set_int32_value(grib, key, value):
"""
Ensure the workaround function :func:`fixup_int32_as_uint32` is applied as
necessary to problem keys.
"""
try:
gribapi.grib_set(grib, key, value)
except gribapi.GribInternalError:
value = fixup_int32_as_uint32(value)
gribapi.grib_set(grib, key, value)
###############################################################################
#
# Constants
#
###############################################################################
# Reference Flag Table 3.3
_RESOLUTION_AND_COMPONENTS_GRID_WINDS_BIT = 3 # NB "bit5", from MSB=1.
# Reference Regulation 92.1.6
_DEFAULT_DEGREES_UNITS = 1.0e-6
###############################################################################
#
# Identification Section 1
#
###############################################################################
def centre(cube, grib):
# TODO: read centre from cube
gribapi.grib_set_long(grib, "centre", 74) # UKMO
gribapi.grib_set_long(grib, "subCentre", 0) # exeter is not in the spec
def reference_time(cube, grib):
# Set the reference time.
# (analysis, forecast start, verify time, obs time, etc)
try:
fp_coord = cube.coord("forecast_period")
except iris.exceptions.CoordinateNotFoundError:
fp_coord = None
if fp_coord is not None:
rt, rt_meaning, _, _ = _non_missing_forecast_period(cube)
else:
rt, rt_meaning, _, _ = _missing_forecast_period(cube)
gribapi.grib_set_long(grib, "significanceOfReferenceTime", rt_meaning)
gribapi.grib_set_long(
grib, "dataDate", "%04d%02d%02d" % (rt.year, rt.month, rt.day))
gribapi.grib_set_long(
grib, "dataTime", "%02d%02d" % (rt.hour, rt.minute))
# TODO: Set the calendar, when we find out what happened to the proposal!
# http://tinyurl.com/oefqgv6
# I was sure it was approved for pre-operational use but it's not there.
def identification(cube, grib):
centre(cube, grib)
reference_time(cube, grib)
# operational product, operational test, research product, etc
# (missing for now)
gribapi.grib_set_long(grib, "productionStatusOfProcessedData", 255)
# Code table 1.4
# analysis, forecast, processed satellite, processed radar,
if cube.coords('realization'):
# assume realization will always have 1 and only 1 point
# as cubes saving to GRIB2 a 2D horizontal slices
if cube.coord('realization').points[0] != 0:
gribapi.grib_set_long(grib, "typeOfProcessedData", 4)
else:
gribapi.grib_set_long(grib, "typeOfProcessedData", 3)
else:
gribapi.grib_set_long(grib, "typeOfProcessedData", 2)
###############################################################################
#
# Grid Definition Section 3
#
###############################################################################
def shape_of_the_earth(cube, grib):
# assume latlon
cs = cube.coord(dimensions=[0]).coord_system
# Initially set shape_of_earth keys to missing (255 for byte, -1 for long).
gribapi.grib_set_long(grib, "scaleFactorOfRadiusOfSphericalEarth", 255)
gribapi.grib_set_long(grib, "scaledValueOfRadiusOfSphericalEarth", -1)
gribapi.grib_set_long(grib, "scaleFactorOfEarthMajorAxis", 255)
gribapi.grib_set_long(grib, "scaledValueOfEarthMajorAxis", -1)
gribapi.grib_set_long(grib, "scaleFactorOfEarthMinorAxis", 255)
gribapi.grib_set_long(grib, "scaledValueOfEarthMinorAxis", -1)
if isinstance(cs, GeogCS):
ellipsoid = cs
else:
ellipsoid = cs.ellipsoid
if ellipsoid is None:
msg = "Could not determine shape of the earth from coord system "\
"of horizontal grid."
raise iris.exceptions.TranslationError(msg)
# Spherical earth.
if ellipsoid.inverse_flattening == 0.0:
gribapi.grib_set_long(grib, "shapeOfTheEarth", 1)
gribapi.grib_set_long(grib, "scaleFactorOfRadiusOfSphericalEarth", 0)
gribapi.grib_set_long(grib, "scaledValueOfRadiusOfSphericalEarth",
ellipsoid.semi_major_axis)
gribapi.grib_set_long(grib, "scaleFactorOfEarthMajorAxis", 0)
gribapi.grib_set_long(grib, "scaledValueOfEarthMajorAxis", 0)
gribapi.grib_set_long(grib, "scaleFactorOfEarthMinorAxis", 0)
gribapi.grib_set_long(grib, "scaledValueOfEarthMinorAxis", 0)
# Oblate spheroid earth.
else:
gribapi.grib_set_long(grib, "shapeOfTheEarth", 7)
gribapi.grib_set_long(grib, "scaleFactorOfEarthMajorAxis", 0)
gribapi.grib_set_long(grib, "scaledValueOfEarthMajorAxis",
ellipsoid.semi_major_axis)
gribapi.grib_set_long(grib, "scaleFactorOfEarthMinorAxis", 0)
gribapi.grib_set_long(grib, "scaledValueOfEarthMinorAxis",
ellipsoid.semi_minor_axis)
def grid_dims(x_coord, y_coord, grib, x_str, y_str):
gribapi.grib_set_long(grib, x_str, x_coord.shape[0])
gribapi.grib_set_long(grib, y_str, y_coord.shape[0])
def latlon_first_last(x_coord, y_coord, grib):
if x_coord.has_bounds() or y_coord.has_bounds():
warnings.warn("Ignoring xy bounds")
# XXX Pending #1125
# gribapi.grib_set_double(grib, "latitudeOfFirstGridPointInDegrees",
# float(y_coord.points[0]))
# gribapi.grib_set_double(grib, "latitudeOfLastGridPointInDegrees",
# float(y_coord.points[-1]))
# gribapi.grib_set_double(grib, "longitudeOfFirstGridPointInDegrees",
# float(x_coord.points[0]))
# gribapi.grib_set_double(grib, "longitudeOfLastGridPointInDegrees",
# float(x_coord.points[-1]))
# WORKAROUND
gribapi.grib_set_long(grib, "latitudeOfFirstGridPoint",
int(y_coord.points[0]*1000000))
gribapi.grib_set_long(grib, "latitudeOfLastGridPoint",
int(y_coord.points[-1]*1000000))
gribapi.grib_set_long(grib, "longitudeOfFirstGridPoint",
int((x_coord.points[0] % 360)*1000000))
gribapi.grib_set_long(grib, "longitudeOfLastGridPoint",
int((x_coord.points[-1] % 360)*1000000))
def dx_dy(x_coord, y_coord, grib):
x_step = regular_step(x_coord)
y_step = regular_step(y_coord)
# Set x and y step. For degrees, this is encoded as an integer:
# 1 * 10^6 * floating point value.
# WMO Manual on Codes regulation 92.1.6
if x_coord.units == 'degrees':
gribapi.grib_set(grib, "iDirectionIncrement",
round(1e6 * float(abs(x_step))))
else:
raise ValueError('X coordinate must be in degrees, not {}'
'.'.format(x_coord.units))
if y_coord.units == 'degrees':
gribapi.grib_set(grib, "jDirectionIncrement",
round(1e6 * float(abs(y_step))))
else:
raise ValueError('Y coordinate must be in degrees, not {}'
'.'.format(y_coord.units))
def scanning_mode_flags(x_coord, y_coord, grib):
gribapi.grib_set_long(grib, "iScansPositively",
int(x_coord.points[1] - x_coord.points[0] > 0))
gribapi.grib_set_long(grib, "jScansPositively",
int(y_coord.points[1] - y_coord.points[0] > 0))
def horizontal_grid_common(cube, grib, xy=False):
nx_str, ny_str = ("Nx", "Ny") if xy else ("Ni", "Nj")
# Grib encoding of the sequences of X and Y points.
y_coord = cube.coord(dimensions=[0])
x_coord = cube.coord(dimensions=[1])
shape_of_the_earth(cube, grib)
grid_dims(x_coord, y_coord, grib, nx_str, ny_str)
scanning_mode_flags(x_coord, y_coord, grib)
def latlon_points_regular(cube, grib):
y_coord = cube.coord(dimensions=[0])
x_coord = cube.coord(dimensions=[1])
latlon_first_last(x_coord, y_coord, grib)
dx_dy(x_coord, y_coord, grib)
def latlon_points_irregular(cube, grib):
y_coord = cube.coord(dimensions=[0])
x_coord = cube.coord(dimensions=[1])
# Distinguish between true-north and grid-oriented vectors.
is_grid_wind = cube.name() in ('x_wind', 'y_wind', 'grid_eastward_wind',
'grid_northward_wind')
# Encode in bit "5" of 'resolutionAndComponentFlags' (other bits unused).
component_flags = 0
if is_grid_wind:
component_flags |= 2 ** _RESOLUTION_AND_COMPONENTS_GRID_WINDS_BIT
gribapi.grib_set(grib, 'resolutionAndComponentFlags', component_flags)
# Record the X and Y coordinate values.
# NOTE: there is currently a bug in the gribapi which means that the size
# of the longitudes array does not equal 'Nj', as it should.
# See : https://software.ecmwf.int/issues/browse/SUP-1096
# So, this only works at present if the x and y dimensions are **equal**.
lon_values = x_coord.points / _DEFAULT_DEGREES_UNITS
lat_values = y_coord.points / _DEFAULT_DEGREES_UNITS
gribapi.grib_set_array(grib, 'longitudes',
np.array(np.round(lon_values), dtype=np.int64))
gribapi.grib_set_array(grib, 'latitudes',
np.array(np.round(lat_values), dtype=np.int64))
def rotated_pole(cube, grib):
# Grib encoding of a rotated pole coordinate system.
cs = cube.coord(dimensions=[0]).coord_system
if cs.north_pole_grid_longitude != 0.0:
raise iris.exceptions.TranslationError(
'Grib save does not yet support Rotated-pole coordinates with '
'a rotated prime meridian.')
# XXX Pending #1125
# gribapi.grib_set_double(grib, "latitudeOfSouthernPoleInDegrees",
# float(cs.n_pole.latitude))
# gribapi.grib_set_double(grib, "longitudeOfSouthernPoleInDegrees",
# float(cs.n_pole.longitude))
# gribapi.grib_set_double(grib, "angleOfRotationInDegrees", 0)
# WORKAROUND
latitude = cs.grid_north_pole_latitude / _DEFAULT_DEGREES_UNITS
longitude = (((cs.grid_north_pole_longitude + 180) % 360) /
_DEFAULT_DEGREES_UNITS)
gribapi.grib_set(grib, "latitudeOfSouthernPole", - int(round(latitude)))
gribapi.grib_set(grib, "longitudeOfSouthernPole", int(round(longitude)))
gribapi.grib_set(grib, "angleOfRotation", 0)
def points_in_unit(coord, unit):
points = coord.units.convert(coord.points, unit)
points = np.around(points).astype(int)
return points
def step(points, atol):
diffs = points[1:] - points[:-1]
mean_diff = np.mean(diffs).astype(points.dtype)
if not np.allclose(diffs, mean_diff, atol=atol):
raise ValueError()
return int(mean_diff)
def grid_definition_template_0(cube, grib):
"""
Set keys within the provided grib message based on
Grid Definition Template 3.0.
Template 3.0 is used to represent "latitude/longitude (or equidistant
cylindrical, or Plate Carree)".
The coordinates are regularly spaced, true latitudes and longitudes.
"""
# Constant resolution, aka 'regular' true lat-lon grid.
gribapi.grib_set_long(grib, "gridDefinitionTemplateNumber", 0)
horizontal_grid_common(cube, grib)
latlon_points_regular(cube, grib)
def grid_definition_template_1(cube, grib):
"""
Set keys within the provided grib message based on
Grid Definition Template 3.1.
Template 3.1 is used to represent "rotated latitude/longitude (or
equidistant cylindrical, or Plate Carree)".
The coordinates are regularly spaced, rotated latitudes and longitudes.
"""
# Constant resolution, aka 'regular' rotated lat-lon grid.
gribapi.grib_set_long(grib, "gridDefinitionTemplateNumber", 1)
# Record details of the rotated coordinate system.
rotated_pole(cube, grib)
# Encode the lat/lon points.
horizontal_grid_common(cube, grib)
latlon_points_regular(cube, grib)
def grid_definition_template_4(cube, grib):
"""
Set keys within the provided grib message based on
Grid Definition Template 3.4.
Template 3.4 is used to represent "variable resolution latitude/longitude".
The coordinates are irregularly spaced latitudes and longitudes.
"""
# XXX: will we need to set `Ni` and `Nj`?
gribapi.grib_set(grib, "gridDefinitionTemplateNumber", 4)
horizontal_grid_common(cube, grib)
latlon_points_irregular(cube, grib)
def grid_definition_template_5(cube, grib):
"""
Set keys within the provided grib message based on
Grid Definition Template 3.5.
Template 3.5 is used to represent "variable resolution rotated
latitude/longitude".
The coordinates are irregularly spaced, rotated latitudes and longitudes.
"""
# NOTE: we must set Ni=Nj=1 before establishing the template.
# Without this, setting "gridDefinitionTemplateNumber" = 5 causes an
# immediate error.
# See: https://software.ecmwf.int/issues/browse/SUP-1095
# This is acceptable, as the subsequent call to 'horizontal_grid_common'
# will set these to the correct horizontal dimensions
# (by calling 'grid_dims').
gribapi.grib_set(grib, "Ni", 1)
gribapi.grib_set(grib, "Nj", 1)
gribapi.grib_set(grib, "gridDefinitionTemplateNumber", 5)
# Record details of the rotated coordinate system.
rotated_pole(cube, grib)
# Encode the lat/lon points.
horizontal_grid_common(cube, grib)
latlon_points_irregular(cube, grib)
def grid_definition_template_10(cube, grib):
"""
Set keys within the provided grib message based on
Grid Definition Template 3.10.
Template 3.10 is used to represent a Mercator grid.
"""
gribapi.grib_set(grib, "gridDefinitionTemplateNumber", 10)
# Retrieve some information from the cube.
y_coord = cube.coord(dimensions=[0])
x_coord = cube.coord(dimensions=[1])
cs = y_coord.coord_system
# Normalise the coordinate values to millimetres - the resolution
# used in the GRIB message.
y_mm = points_in_unit(y_coord, 'mm')
x_mm = points_in_unit(x_coord, 'mm')
# Encode the horizontal points.
# NB. Since we're already in millimetres, our tolerance for
# discrepancy in the differences is 1.
try:
x_step = step(x_mm, atol=1)
y_step = step(y_mm, atol=1)
except ValueError:
msg = 'Irregular coordinates not supported for Mercator.'
raise iris.exceptions.TranslationError(msg)
gribapi.grib_set(grib, 'Di', abs(x_step))
gribapi.grib_set(grib, 'Dj', abs(y_step))
horizontal_grid_common(cube, grib)
# Transform first and last points into geographic CS.
geog = cs.ellipsoid if cs.ellipsoid is not None else GeogCS(1)
first_x, first_y, = geog.as_cartopy_crs().transform_point(
x_coord.points[0],
y_coord.points[0],
cs.as_cartopy_crs())
last_x, last_y = geog.as_cartopy_crs().transform_point(
x_coord.points[-1],
y_coord.points[-1],
cs.as_cartopy_crs())
first_x = first_x % 360
last_x = last_x % 360
gribapi.grib_set(grib, "latitudeOfFirstGridPoint",
int(np.round(first_y / _DEFAULT_DEGREES_UNITS)))
gribapi.grib_set(grib, "longitudeOfFirstGridPoint",
int(np.round(first_x / _DEFAULT_DEGREES_UNITS)))
gribapi.grib_set(grib, "latitudeOfLastGridPoint",
int(np.round(last_y / _DEFAULT_DEGREES_UNITS)))
gribapi.grib_set(grib, "longitudeOfLastGridPoint",
int(np.round(last_x / _DEFAULT_DEGREES_UNITS)))
# Check and raise a more intelligible error, if the Iris version is too old
# to support the Mercator 'standard_parallel' property.
confirm_extended_mercator_supported()
# Encode the latitude at which the projection intersects the Earth.
gribapi.grib_set(grib, 'LaD',
cs.standard_parallel / _DEFAULT_DEGREES_UNITS)
# Encode resolution and component flags
gribapi.grib_set(grib, 'resolutionAndComponentFlags',
0x1 << _RESOLUTION_AND_COMPONENTS_GRID_WINDS_BIT)
def grid_definition_template_12(cube, grib):
"""
Set keys within the provided grib message based on
Grid Definition Template 3.12.
Template 3.12 is used to represent a Transverse Mercator grid.
"""
gribapi.grib_set(grib, "gridDefinitionTemplateNumber", 12)
# Retrieve some information from the cube.
y_coord = cube.coord(dimensions=[0])
x_coord = cube.coord(dimensions=[1])
cs = y_coord.coord_system
# Normalise the coordinate values to centimetres - the resolution
# used in the GRIB message.
y_cm = points_in_unit(y_coord, 'cm')
x_cm = points_in_unit(x_coord, 'cm')
# Set some keys specific to GDT12.
# Encode the horizontal points.
# NB. Since we're already in centimetres, our tolerance for
# discrepancy in the differences is 1.
try:
x_step = step(x_cm, atol=1)
y_step = step(y_cm, atol=1)
except ValueError:
msg = ('Irregular coordinates not supported for Transverse '
'Mercator.')
raise iris.exceptions.TranslationError(msg)
gribapi.grib_set(grib, 'Di', abs(x_step))
gribapi.grib_set(grib, 'Dj', abs(y_step))
horizontal_grid_common(cube, grib)
# GRIBAPI expects unsigned ints in X1, X2, Y1, Y2 but it should accept
# signed ints, so work around this.
# See https://software.ecmwf.int/issues/browse/SUP-1101
ensure_set_int32_value(grib, 'Y1', int(y_cm[0]))
ensure_set_int32_value(grib, 'Y2', int(y_cm[-1]))
ensure_set_int32_value(grib, 'X1', int(x_cm[0]))
ensure_set_int32_value(grib, 'X2', int(x_cm[-1]))
# Lat and lon of reference point are measured in millionths of a degree.
gribapi.grib_set(grib, "latitudeOfReferencePoint",
cs.latitude_of_projection_origin / _DEFAULT_DEGREES_UNITS)
gribapi.grib_set(grib, "longitudeOfReferencePoint",
cs.longitude_of_central_meridian / _DEFAULT_DEGREES_UNITS)
# Convert a value in metres into the closest integer number of
# centimetres.
def m_to_cm(value):
return int(round(value * 100))
# False easting and false northing are measured in units of (10^-2)m.
gribapi.grib_set(grib, 'XR', m_to_cm(cs.false_easting))
gribapi.grib_set(grib, 'YR', m_to_cm(cs.false_northing))
# GRIBAPI expects a signed int for scaleFactorAtReferencePoint
# but it should accept a float, so work around this.
# See https://software.ecmwf.int/issues/browse/SUP-1100
value = cs.scale_factor_at_central_meridian
key_type = gribapi.grib_get_native_type(grib,
"scaleFactorAtReferencePoint")
if key_type is not float:
value = fixup_float32_as_int32(value)
gribapi.grib_set(grib, "scaleFactorAtReferencePoint", value)
def grid_definition_template_30(cube, grib):
"""
Set keys within the provided grib message based on
Grid Definition Template 3.30.
Template 3.30 is used to represent a Lambert Conformal grid.
"""
gribapi.grib_set(grib, "gridDefinitionTemplateNumber", 30)
# Retrieve some information from the cube.
y_coord = cube.coord(dimensions=[0])
x_coord = cube.coord(dimensions=[1])
cs = y_coord.coord_system
# Normalise the coordinate values to millimetres - the resolution
# used in the GRIB message.
y_mm = points_in_unit(y_coord, 'mm')
x_mm = points_in_unit(x_coord, 'mm')
# Encode the horizontal points.
# NB. Since we're already in millimetres, our tolerance for
# discrepancy in the differences is 1.
try:
x_step = step(x_mm, atol=1)
y_step = step(y_mm, atol=1)
except ValueError:
msg = ('Irregular coordinates not supported for Lambert '
'Conformal.')
raise iris.exceptions.TranslationError(msg)
gribapi.grib_set(grib, 'Dx', abs(x_step))
gribapi.grib_set(grib, 'Dy', abs(y_step))
horizontal_grid_common(cube, grib, xy=True)
# Transform first point into geographic CS
geog = cs.ellipsoid if cs.ellipsoid is not None else GeogCS(1)
first_x, first_y = geog.as_cartopy_crs().transform_point(
x_coord.points[0],
y_coord.points[0],
cs.as_cartopy_crs())
first_x = first_x % 360
central_lon = cs.central_lon % 360
gribapi.grib_set(grib, "latitudeOfFirstGridPoint",
int(np.round(first_y / _DEFAULT_DEGREES_UNITS)))
gribapi.grib_set(grib, "longitudeOfFirstGridPoint",
int(np.round(first_x / _DEFAULT_DEGREES_UNITS)))
gribapi.grib_set(grib, "LaD", cs.central_lat / _DEFAULT_DEGREES_UNITS)
gribapi.grib_set(grib, "LoV", central_lon / _DEFAULT_DEGREES_UNITS)
latin1, latin2 = cs.secant_latitudes
gribapi.grib_set(grib, "Latin1", latin1 / _DEFAULT_DEGREES_UNITS)
gribapi.grib_set(grib, "Latin2", latin2 / _DEFAULT_DEGREES_UNITS)
gribapi.grib_set(grib, 'resolutionAndComponentFlags',
0x1 << _RESOLUTION_AND_COMPONENTS_GRID_WINDS_BIT)
# Which pole are the parallels closest to? That is the direction
# that the cone converges.
poliest_sec = latin1 if abs(latin1) > abs(latin2) else latin2
centre_flag = 0x0 if poliest_sec > 0 else 0x1
gribapi.grib_set(grib, 'projectionCentreFlag', centre_flag)
gribapi.grib_set(grib, "latitudeOfSouthernPole", 0)
gribapi.grib_set(grib, "longitudeOfSouthernPole", 0)
def grid_definition_section(cube, grib):
"""
Set keys within the grid definition section of the provided grib message,
based on the properties of the cube.
"""
x_coord = cube.coord(dimensions=[1])
y_coord = cube.coord(dimensions=[0])
cs = x_coord.coord_system # N.B. already checked same cs for x and y.
regular_x_and_y = is_regular(x_coord) and is_regular(y_coord)
if isinstance(cs, GeogCS):
if regular_x_and_y:
grid_definition_template_0(cube, grib)
else:
grid_definition_template_4(cube, grib)
elif isinstance(cs, RotatedGeogCS):
# Rotated coordinate system cases.
# Choose between GDT 3.1 and 3.5 according to coordinate regularity.
if regular_x_and_y:
grid_definition_template_1(cube, grib)
else:
grid_definition_template_5(cube, grib)
elif isinstance(cs, Mercator):
# Mercator coordinate system (template 3.10)
grid_definition_template_10(cube, grib)
elif isinstance(cs, TransverseMercator):
# Transverse Mercator coordinate system (template 3.12).
grid_definition_template_12(cube, grib)
elif isinstance(cs, LambertConformal):
# Lambert Conformal coordinate system (template 3.30).
grid_definition_template_30(cube, grib)
else:
name = cs.grid_mapping_name.replace('_', ' ').title()
emsg = 'Grib saving is not supported for coordinate system {!r}'
raise ValueError(emsg.format(name))
###############################################################################
#
# Product Definition Section 4
#
###############################################################################
def set_discipline_and_parameter(cube, grib):
# NOTE: for now, can match by *either* standard_name or long_name.
# This allows workarounds for data with no identified standard_name.
grib2_info = gptx.cf_phenom_to_grib2_info(cube.standard_name,
cube.long_name)
if grib2_info is not None:
gribapi.grib_set(grib, "discipline", grib2_info.discipline)
gribapi.grib_set(grib, "parameterCategory", grib2_info.category)
gribapi.grib_set(grib, "parameterNumber", grib2_info.number)
else:
gribapi.grib_set(grib, "discipline", 255)
gribapi.grib_set(grib, "parameterCategory", 255)
gribapi.grib_set(grib, "parameterNumber", 255)
warnings.warn('Unable to determine Grib2 parameter code for cube.\n'
'discipline, parameterCategory and parameterNumber '
'have been set to "missing".')
def _non_missing_forecast_period(cube):
# Calculate "model start time" to use as the reference time.
fp_coord = cube.coord("forecast_period")
# Convert fp and t to hours so we can subtract to calculate R.
cf_fp_hrs = fp_coord.units.convert(fp_coord.points[0], 'hours')
t_coord = cube.coord("time").copy()
hours_since = cf_units.Unit("hours since epoch",
calendar=t_coord.units.calendar)
t_coord.convert_units(hours_since)
rt_num = t_coord.points[0] - cf_fp_hrs
rt = hours_since.num2date(rt_num)
rt_meaning = 1 # "start of forecast"
# Forecast period
if fp_coord.units == cf_units.Unit("hours"):
grib_time_code = 1
elif fp_coord.units == cf_units.Unit("minutes"):
grib_time_code = 0
elif fp_coord.units == cf_units.Unit("seconds"):
grib_time_code = 13
else:
raise iris.exceptions.TranslationError(
"Unexpected units for 'forecast_period' : %s" % fp_coord.units)
if not t_coord.has_bounds():
fp = fp_coord.points[0]
else:
if not fp_coord.has_bounds():
raise iris.exceptions.TranslationError(
"bounds on 'time' coordinate requires bounds on"
" 'forecast_period'.")
fp = fp_coord.bounds[0][0]
if fp - int(fp):
warnings.warn("forecast_period encoding problem: "
"scaling required.")
fp = int(fp)
return rt, rt_meaning, fp, grib_time_code
def _missing_forecast_period(cube):
"""
Returns a reference time and significance code together with a forecast
period and corresponding units type code.
"""
t_coord = cube.coord("time")
if cube.coords('forecast_reference_time'):
# Make copies and convert them to common "hours since" units.
hours_since = cf_units.Unit('hours since epoch',
calendar=t_coord.units.calendar)
frt_coord = cube.coord('forecast_reference_time').copy()
frt_coord.convert_units(hours_since)
t_coord = t_coord.copy()
t_coord.convert_units(hours_since)
# Extract values.
t = t_coord.bounds[0, 0] if t_coord.has_bounds() else t_coord.points[0]
frt = frt_coord.points[0]
# Calculate GRIB parameters.
rt = frt_coord.units.num2date(frt)
rt_meaning = 1 # Forecast reference time.
fp = t - frt
integer_fp = int(fp)
if integer_fp != fp:
msg = 'Truncating floating point forecast period {} to ' \
'integer value {}'
warnings.warn(msg.format(fp, integer_fp))
fp = integer_fp
fp_meaning = 1 # Hours
else:
# With no forecast period or forecast reference time set assume a
# reference time significance of "Observation time" and set the
# forecast period to 0h.
t = t_coord.bounds[0, 0] if t_coord.has_bounds() else t_coord.points[0]
rt = t_coord.units.num2date(t)
rt_meaning = 3 # Observation time
fp = 0
fp_meaning = 1 # Hours
return rt, rt_meaning, fp, fp_meaning
def set_forecast_time(cube, grib):
"""
Set the forecast time keys based on the forecast_period coordinate. In
the absence of a forecast_period and forecast_reference_time,
the forecast time is set to zero.
"""
try:
fp_coord = cube.coord("forecast_period")
except iris.exceptions.CoordinateNotFoundError:
fp_coord = None
if fp_coord is not None:
_, _, fp, grib_time_code = _non_missing_forecast_period(cube)
else:
_, _, fp, grib_time_code = _missing_forecast_period(cube)
gribapi.grib_set(grib, "indicatorOfUnitOfTimeRange", grib_time_code)
gribapi.grib_set(grib, "forecastTime", fp)
def set_fixed_surfaces(cube, grib, full3d_cube=None):
# Look for something we can export
v_coord = grib_v_code = output_unit = None
# Detect factories for hybrid vertical coordinates.
hybrid_factories = [
factory for factory in cube.aux_factories
if isinstance(factory, (HybridHeightFactory, HybridPressureFactory))]
if not hybrid_factories:
hybrid_factory = None
elif len(hybrid_factories) > 1:
msg = 'Data contains >1 vertical coordinate factory : {}'
raise ValueError(msg.format(hybrid_factories))
else:
factory = hybrid_factories[0]
# Fetch the matching 'complete' factory from the *full* 3d cube, so we
# have all the level information.
hybrid_factory = full3d_cube.aux_factory(factory.name())
# Handle various different styles of vertical coordinate.
# hybrid height / pressure
if hybrid_factory is not None:
# N.B. in this case, there are additional operations, besides just
# encoding v_coord : see below at end ..
v_coord = cube.coord('model_level_number')
output_unit = cf_units.Unit("1")
if isinstance(hybrid_factory, HybridHeightFactory):
grib_v_code = 118
elif isinstance(hybrid_factory, HybridPressureFactory):
grib_v_code = 119
else:
msg = 'Unrecognised factory type : {}'
raise ValueError(msg.format(hybrid_factory))
# pressure
elif cube.coords("air_pressure") or cube.coords("pressure"):
grib_v_code = 100
output_unit = cf_units.Unit("Pa")
v_coord = (cube.coords("air_pressure") or cube.coords("pressure"))[0]
# altitude
elif cube.coords("altitude"):
grib_v_code = 102
output_unit = cf_units.Unit("m")
v_coord = cube.coord("altitude")
# height
elif cube.coords("height"):
grib_v_code = 103
output_unit = cf_units.Unit("m")
v_coord = cube.coord("height")
elif cube.coords("air_potential_temperature"):
grib_v_code = 107
output_unit = cf_units.Unit('K')
v_coord = cube.coord("air_potential_temperature")
# unknown / absent
else:
# check for *ANY* height coords at all...
v_coords = cube.coords(axis='z')
if v_coords:
# There are vertical coordinate(s), but we don't understand them...
v_coords_str = ' ,'.join(["'{}'".format(c.name())
for c in v_coords])
raise iris.exceptions.TranslationError(
'The vertical-axis coordinate(s) ({}) '
'are not recognised or handled.'.format(v_coords_str))
# What did we find?
if v_coord is None:
# No vertical coordinate: record as 'surface' level (levelType=1).
# NOTE: may *not* be truly correct, but seems to be common practice.
# Still under investigation :
# See https://github.com/SciTools/iris/issues/519
gribapi.grib_set(grib, "typeOfFirstFixedSurface", 1)
gribapi.grib_set(grib, "scaleFactorOfFirstFixedSurface", 0)
gribapi.grib_set(grib, "scaledValueOfFirstFixedSurface", 0)
# Set secondary surface = 'missing'.
gribapi.grib_set(grib, "typeOfSecondFixedSurface", -1)
gribapi.grib_set(grib, "scaleFactorOfSecondFixedSurface", 255)
gribapi.grib_set(grib, "scaledValueOfSecondFixedSurface", -1)
elif not v_coord.has_bounds():
# No second surface
output_v = v_coord.units.convert(v_coord.points[0], output_unit)
if output_v - abs(output_v):
warnings.warn("Vertical level encoding problem: scaling required.")
output_v = int(output_v)
gribapi.grib_set(grib, "typeOfFirstFixedSurface", grib_v_code)
gribapi.grib_set(grib, "scaleFactorOfFirstFixedSurface", 0)
gribapi.grib_set(grib, "scaledValueOfFirstFixedSurface", output_v)
gribapi.grib_set(grib, "typeOfSecondFixedSurface", -1)
gribapi.grib_set(grib, "scaleFactorOfSecondFixedSurface", 255)
gribapi.grib_set(grib, "scaledValueOfSecondFixedSurface", -1)
else:
# bounded : set lower+upper surfaces
output_v = v_coord.units.convert(v_coord.bounds[0], output_unit)
if output_v[0] - abs(output_v[0]) or output_v[1] - abs(output_v[1]):
warnings.warn("Vertical level encoding problem: scaling required.")
gribapi.grib_set(grib, "typeOfFirstFixedSurface", grib_v_code)
gribapi.grib_set(grib, "typeOfSecondFixedSurface", grib_v_code)
gribapi.grib_set(grib, "scaleFactorOfFirstFixedSurface", 0)
gribapi.grib_set(grib, "scaleFactorOfSecondFixedSurface", 0)
gribapi.grib_set(grib, "scaledValueOfFirstFixedSurface",
int(output_v[0]))
gribapi.grib_set(grib, "scaledValueOfSecondFixedSurface",
int(output_v[1]))
if hybrid_factory is not None:
# Need to record ALL the level coefficents in a 'PV' vector.
level_delta_coord = hybrid_factory.delta
sigma_coord = hybrid_factory.sigma
model_levels = full3d_cube.coord('model_level_number').points
# Just check these make some kind of sense (!)
if model_levels.dtype.kind not in 'iu':
msg = 'model_level_number is not an integer: dtype={}.'
raise ValueError(msg.format(model_levels.dtype))
if np.min(model_levels) < 1:
msg = 'model_level_number must be > 0: mininum value = {}.'
raise ValueError(msg.format(np.min(model_levels)))
# Need to save enough levels for indexes up to [max(model_levels)]
n_levels = np.max(model_levels)
max_valid_nlevels = 9999
if n_levels > max_valid_nlevels:
msg = ('model_level_number values are > {} : '
'maximum value = {}.')
raise ValueError(msg.format(max_valid_nlevels, n_levels))
# In sample data we have seen, there seems to be an extra missing data
# value after each set of n-levels coefficients.
# For now, we retain that encoding. This should not cause problems as
# values are indexed according to model-level-number,
# I.E. sigma, delta = PV[i], PV[NV/2+i]
n_coeffs = n_levels + 1
coeffs_array = np.zeros(n_coeffs * 2, dtype=np.float32)
for n_lev, height, sigma in zip(model_levels,
level_delta_coord.points,
sigma_coord.points):
# Record all the level coefficients coming from the 'full' cube.
# Note: if some model levels are missing, we must still have the
# coeffs at the correct index according to the model_level_number
# value, i.e. at [level - 1] and [NV // 2 + level - 1].
# Thus, we can *not* paste the values in a block: each one needs to
# go in the index corresponding to its 'model_level_number' value.
coeffs_array[n_lev - 1] = height
coeffs_array[n_coeffs + n_lev - 1] = sigma
pv_values = [float(el) for el in coeffs_array]
gribapi.grib_set(grib, "NV", n_coeffs * 2)
gribapi.grib_set_array(grib, "pv", pv_values)
def set_time_range(time_coord, grib):
"""
Set the time range keys in the specified message
based on the bounds of the provided time coordinate.
"""
if len(time_coord.points) != 1:
msg = 'Expected length one time coordinate, got {} points'
raise ValueError(msg.format(len(time_coord.points)))
if time_coord.nbounds != 2:
msg = 'Expected time coordinate with two bounds, got {} bounds'
raise ValueError(msg.format(time_coord.nbounds))
# Set type to hours and convert period to this unit.
gribapi.grib_set(grib, "indicatorOfUnitForTimeRange",
_TIME_RANGE_UNITS['hours'])
hours_since_units = cf_units.Unit('hours since epoch',
calendar=time_coord.units.calendar)
start_hours, end_hours = time_coord.units.convert(time_coord.bounds[0],