-
Notifications
You must be signed in to change notification settings - Fork 279
/
aux_factory.py
1897 lines (1595 loc) · 69.1 KB
/
aux_factory.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
# Copyright Iris contributors
#
# This file is part of Iris and is released under the BSD license.
# See LICENSE in the root of the repository for full licensing details.
"""Definitions of derived coordinates."""
from abc import ABCMeta, abstractmethod
import cf_units
import dask.array as da
import numpy as np
from iris.common import CFVariableMixin, CoordMetadata, metadata_manager_factory
import iris.coords
from iris.exceptions import IrisIgnoringBoundsWarning, warn_once_at_level
class AuxCoordFactory(CFVariableMixin, metaclass=ABCMeta):
"""Represents a "factory" which can manufacture additional auxiliary coordinate.
Represents a "factory" which can manufacture an additional auxiliary
coordinate on demand, by combining the values of other coordinates.
Each concrete subclass represents a specific formula for deriving
values from other coordinates.
The `standard_name`, `long_name`, `var_name`, `units`, `attributes` and
`coord_system` of the factory are used to set the corresponding
properties of the resulting auxiliary coordinates.
"""
def __init__(self):
# Configure the metadata manager.
if not hasattr(self, "_metadata_manager"):
self._metadata_manager = metadata_manager_factory(CoordMetadata)
#: Descriptive name of the coordinate made by the factory
self.long_name = None
#: netCDF variable name for the coordinate made by the factory
self.var_name = None
self.coord_system = None
# See the climatological property getter.
self._metadata_manager.climatological = False
@property
def coord_system(self):
"""The coordinate-system (if any) of the coordinate made by the factory."""
return self._metadata_manager.coord_system
@coord_system.setter
def coord_system(self, value):
self._metadata_manager.coord_system = value
@property
def climatological(self):
"""Always returns False, as a factory itself can never have points/bounds.
Always returns False, as a factory itself can never have points/bounds
and therefore can never be climatological by definition.
"""
return self._metadata_manager.climatological
@property
@abstractmethod
def dependencies(self):
"""Return a dict mapping from constructor argument.
Return a dictionary mapping from constructor argument names to
the corresponding coordinates.
"""
@abstractmethod
def make_coord(self, coord_dims_func):
"""Return a new :class:`iris.coords.AuxCoord` as defined by this factory.
Parameters
----------
coord_dims_func :
A callable which can return the list of dimensions relevant
to a given coordinate.
See :meth:`iris.cube.Cube.coord_dims()`.
"""
def update(self, old_coord, new_coord=None):
"""Notify the factory of the removal/replacement of a coordinate.
Notify the factory of the removal/replacement of a coordinate
which might be a dependency.
Parameters
----------
old_coord :
The coordinate to be removed/replaced.
new_coord : optional
If None, any dependency using old_coord is removed, otherwise
any dependency using old_coord is updated to use new_coord.
"""
new_dependencies = self.dependencies
for name, coord in self.dependencies.items():
if old_coord is coord:
new_dependencies[name] = new_coord
try:
self._check_dependencies(**new_dependencies)
except ValueError as e:
msg = "Failed to update dependencies. " + str(e)
raise ValueError(msg)
else:
setattr(self, name, new_coord)
break
def __repr__(self):
def arg_text(item):
key, coord = item
return "{}={}".format(key, str(coord and repr(coord.name())))
items = sorted(self.dependencies.items(), key=lambda item: item[0])
args = map(arg_text, items)
return "<{}({})>".format(type(self).__name__, ", ".join(args))
def derived_dims(self, coord_dims_func):
"""Return the cube dimensions for the derived coordinate.
Parameters
----------
coord_dims_func :
A callable which can return the list of dimensions relevant to a given
coordinate.
See :meth:`iris.cube.Cube.coord_dims()`.
Returns
-------
A sorted list of cube dimension numbers.
"""
# Which dimensions are relevant?
# e.g. If sigma -> [1] and orog -> [2, 3] then result = [1, 2, 3]
derived_dims = set()
for coord in self.dependencies.values():
if coord:
derived_dims.update(coord_dims_func(coord))
# Apply a fixed order so we know how to map dependency dims to
# our own dims (and so the Cube can map them to Cube dims).
derived_dims = tuple(sorted(derived_dims))
return derived_dims
def updated(self, new_coord_mapping):
"""Create a new instance of this factory.
Create a new instance of this factory where the dependencies
are replaced according to the given mapping.
Parameters
----------
new_coord_mapping :
A dictionary mapping from the object IDs potentially used
by this factory, to the coordinate objects that should be
used instead.
"""
new_dependencies = {}
for key, coord in self.dependencies.items():
if coord:
coord = new_coord_mapping[id(coord)]
new_dependencies[key] = coord
return type(self)(**new_dependencies)
def xml_element(self, doc):
"""Return a DOM element describing this coordinate factory."""
element = doc.createElement("coordFactory")
for key, coord in self.dependencies.items():
element.setAttribute(key, coord._xml_id())
element.appendChild(self.make_coord().xml_element(doc))
return element
def _dependency_dims(self, coord_dims_func):
dependency_dims = {}
for key, coord in self.dependencies.items():
if coord:
dependency_dims[key] = coord_dims_func(coord)
return dependency_dims
@staticmethod
def _nd_bounds(coord, dims, ndim):
"""Return a lazy bounds array for a dependency coordinate, 'coord'.
The result is aligned to the first 'ndim' cube dimensions, and
expanded to the full ('ndim'+1)-dimensional shape.
The value of 'ndim' must be >= the highest cube dimension of the
dependency coordinate.
The extra final result dimension ('ndim'-th) is the bounds dimension.
Example::
coord.shape == (70,)
coord.nbounds = 2
dims == [3]
ndim == 5
results in::
nd_bounds.shape == (1, 1, 1, 70, 1, 2)
"""
# Transpose to be consistent with the Cube.
sorted_pairs = sorted(enumerate(dims), key=lambda pair: pair[1])
transpose_order = [pair[0] for pair in sorted_pairs] + [len(dims)]
bounds = coord.lazy_bounds()
if dims and transpose_order != list(range(len(dims))):
bounds = bounds.transpose(transpose_order)
# Figure out the n-dimensional shape.
nd_shape = [1] * ndim + [coord.nbounds]
for dim, size in zip(dims, coord.shape):
nd_shape[dim] = size
bounds = bounds.reshape(nd_shape)
return bounds
@staticmethod
def _nd_points(coord, dims, ndim):
"""Return a lazy points array for a dependency coordinate, 'coord'.
The result is aligned to the first 'ndim' cube dimensions, and
expanded to the full 'ndim'-dimensional shape.
The value of 'ndim' must be >= the highest cube dimension of the
dependency coordinate.
Examples
--------
::
coord.shape == (4, 3)
dims == [3, 2]
ndim == 5
results in::
nd_points.shape == (1, 1, 3, 4, 1)
"""
# Transpose to be consistent with the Cube.
sorted_pairs = sorted(enumerate(dims), key=lambda pair: pair[1])
transpose_order = [pair[0] for pair in sorted_pairs]
points = coord.lazy_points()
if dims and transpose_order != list(range(len(dims))):
points = points.transpose(transpose_order)
# Expand dimensionality to be consistent with the Cube.
if dims:
keys = [None] * ndim
for dim, size in zip(dims, coord.shape):
keys[dim] = slice(None)
points = points[tuple(keys)]
else:
# Scalar coordinates have one dimensional points despite
# mapping to zero dimensions, so we only need to add N-1
# new dimensions.
keys = (None,) * (ndim - 1)
points = points[keys]
return points
def _remap(self, dependency_dims, derived_dims):
"""Return a mapping from dependency names to coordinate points arrays.
For dependencies that are present, the values are all expanded and
aligned to the same dimensions, which is the full set of all the
dependency dimensions.
These non-missing values are all lazy arrays.
Missing dependencies, however, are assigned a scalar value of 0.0.
"""
if derived_dims:
ndim = max(derived_dims) + 1
else:
ndim = 1
nd_points_by_key = {}
for key, coord in self.dependencies.items():
if coord:
# Get the points as consistent with the Cube.
nd_points = self._nd_points(coord, dependency_dims[key], ndim)
# Restrict to just the dimensions relevant to the
# derived coord. NB. These are always in Cube-order, so
# no transpose is needed.
if derived_dims:
keys = tuple(
slice(None) if dim in derived_dims else 0 for dim in range(ndim)
)
nd_points = nd_points[keys]
else:
# If no coord, treat value as zero.
# Use a float16 to provide `shape` attribute and avoid
# promoting other arguments to a higher precision.
nd_points = np.float16(0)
nd_points_by_key[key] = nd_points
return nd_points_by_key
def _remap_with_bounds(self, dependency_dims, derived_dims):
"""Return a mapping from dependency names to coordinate bounds arrays.
For dependencies that are present, the values are all expanded and
aligned to the same dimensions, which is the full set of all the
dependency dimensions, plus an extra bounds dimension.
These non-missing values are all lazy arrays.
Missing dependencies, however, are assigned a scalar value of 0.0.
Where a dependency coordinate has no bounds, then the associated value
is taken from its points array, but reshaped to have an extra bounds
dimension of length 1.
"""
if derived_dims:
ndim = max(derived_dims) + 1
else:
ndim = 1
nd_values_by_key = {}
for key, coord in self.dependencies.items():
if coord:
# Get the bounds or points as consistent with the Cube.
if coord.nbounds:
nd_values = self._nd_bounds(coord, dependency_dims[key], ndim)
else:
nd_values = self._nd_points(coord, dependency_dims[key], ndim)
# Restrict to just the dimensions relevant to the
# derived coord. NB. These are always in Cube-order, so
# no transpose is needed.
shape = []
for dim in derived_dims:
shape.append(nd_values.shape[dim])
# Ensure the array always has at least one dimension to be
# compatible with normal coordinates.
if not derived_dims:
shape.append(1)
# Add on the N-bounds dimension
if coord.nbounds:
shape.append(nd_values.shape[-1])
else:
# NB. For a non-bounded coordinate we still need an
# extra dimension to make the shape compatible, so
# we just add an extra 1.
shape.append(1)
nd_values = nd_values.reshape(shape)
else:
# If no coord, treat value as zero.
# Use a float16 to provide `shape` attribute and avoid
# promoting other arguments to a higher precision.
nd_values = np.float16(0)
nd_values_by_key[key] = nd_values
return nd_values_by_key
class AtmosphereSigmaFactory(AuxCoordFactory):
"""Define an atmosphere sigma coordinate factory with the following formula.
.. math::
p = ptop + sigma * (ps - ptop)
"""
def __init__(self, pressure_at_top=None, sigma=None, surface_air_pressure=None):
"""Create an atmosphere sigma coordinate factory with a formula.
.. math::
p(n, k, j, i) = pressure_at_top + sigma(k) *
(surface_air_pressure(n, j, i) - pressure_at_top)
"""
# Configure the metadata manager.
self._metadata_manager = metadata_manager_factory(CoordMetadata)
super().__init__()
# Check that provided coordinates meet necessary conditions.
self._check_dependencies(pressure_at_top, sigma, surface_air_pressure)
# Initialize instance attributes
self.units = pressure_at_top.units
self.pressure_at_top = pressure_at_top
self.sigma = sigma
self.surface_air_pressure = surface_air_pressure
self.standard_name = "air_pressure"
self.attributes = {}
@staticmethod
def _check_dependencies(pressure_at_top, sigma, surface_air_pressure):
"""Check for sufficient coordinates."""
if any(
[
pressure_at_top is None,
sigma is None,
surface_air_pressure is None,
]
):
raise ValueError(
"Unable to construct atmosphere sigma coordinate factory due "
"to insufficient source coordinates"
)
# Check dimensions
if pressure_at_top.shape not in ((), (1,)):
raise ValueError(
f"Expected scalar 'pressure_at_top' coordinate, got shape "
f"{pressure_at_top.shape}"
)
# Check bounds
if sigma.nbounds not in (0, 2):
raise ValueError(
f"Invalid 'sigma' coordinate: must have either 0 or 2 bounds, "
f"got {sigma.nbounds:d}"
)
for coord in (pressure_at_top, surface_air_pressure):
if coord.nbounds:
msg = (
f"Coordinate '{coord.name()}' has bounds. These will "
"be disregarded"
)
warn_once_at_level(
msg, category=IrisIgnoringBoundsWarning, stacklevel=2
)
# Check units
if sigma.units.is_unknown():
# Be graceful, and promote unknown to dimensionless units.
sigma.units = cf_units.Unit("1")
if not sigma.units.is_dimensionless():
raise ValueError(
f"Invalid units: 'sigma' must be dimensionless, got " f"'{sigma.units}'"
)
if pressure_at_top.units != surface_air_pressure.units:
raise ValueError(
f"Incompatible units: 'pressure_at_top' and "
f"'surface_air_pressure' must have the same units, got "
f"'{pressure_at_top.units}' and "
f"'{surface_air_pressure.units}'"
)
if not pressure_at_top.units.is_convertible("Pa"):
raise ValueError(
"Invalid units: 'pressure_at_top' and 'surface_air_pressure' "
"must have units of pressure"
)
@property
def dependencies(self):
"""Return dependencies."""
dependencies = {
"pressure_at_top": self.pressure_at_top,
"sigma": self.sigma,
"surface_air_pressure": self.surface_air_pressure,
}
return dependencies
@staticmethod
def _derive(pressure_at_top, sigma, surface_air_pressure):
"""Derive coordinate."""
return pressure_at_top + sigma * (surface_air_pressure - pressure_at_top)
def make_coord(self, coord_dims_func):
"""Return a new :class:`iris.coords.AuxCoord` as defined by this factory.
Parameters
----------
coord_dims_func :
A callable which can return the list of dimensions relevant
to a given coordinate.
See :meth:`iris.cube.Cube.coord_dims()`.
"""
# Which dimensions are relevant?
derived_dims = self.derived_dims(coord_dims_func)
dependency_dims = self._dependency_dims(coord_dims_func)
# Build the points array
nd_points_by_key = self._remap(dependency_dims, derived_dims)
points = self._derive(
nd_points_by_key["pressure_at_top"],
nd_points_by_key["sigma"],
nd_points_by_key["surface_air_pressure"],
)
# Bounds
bounds = None
if self.sigma.nbounds:
nd_values_by_key = self._remap_with_bounds(dependency_dims, derived_dims)
pressure_at_top = nd_values_by_key["pressure_at_top"]
sigma = nd_values_by_key["sigma"]
surface_air_pressure = nd_values_by_key["surface_air_pressure"]
ok_bound_shapes = [(), (1,), (2,)]
if sigma.shape[-1:] not in ok_bound_shapes:
raise ValueError("Invalid sigma coordinate bounds")
if pressure_at_top.shape[-1:] not in [(), (1,)]:
warn_once_at_level(
"Pressure at top coordinate has bounds. These are being "
"disregarded",
category=IrisIgnoringBoundsWarning,
)
pressure_at_top_pts = nd_points_by_key["pressure_at_top"]
bds_shape = list(pressure_at_top_pts.shape) + [1]
pressure_at_top = pressure_at_top_pts.reshape(bds_shape)
if surface_air_pressure.shape[-1:] not in [(), (1,)]:
warn_once_at_level(
"Surface pressure coordinate has bounds. These are being "
"disregarded",
category=IrisIgnoringBoundsWarning,
)
surface_air_pressure_pts = nd_points_by_key["surface_air_pressure"]
bds_shape = list(surface_air_pressure_pts.shape) + [1]
surface_air_pressure = surface_air_pressure_pts.reshape(bds_shape)
bounds = self._derive(pressure_at_top, sigma, surface_air_pressure)
# Create coordinate
return iris.coords.AuxCoord(
points,
standard_name=self.standard_name,
long_name=self.long_name,
var_name=self.var_name,
units=self.units,
bounds=bounds,
attributes=self.attributes,
coord_system=self.coord_system,
)
class HybridHeightFactory(AuxCoordFactory):
"""Defines a hybrid-height coordinate factory."""
def __init__(self, delta=None, sigma=None, orography=None):
"""Create a hybrid-height coordinate factory with the following formula.
.. math::
z = a + b * orog
At least one of `delta` or `orography` must be provided.
Parameters
----------
delta : Coord, optional
The coordinate providing the `a` term.
sigma : Coord, optional
The coordinate providing the `b` term.
orography : Coord, optional
The coordinate providing the `orog` term.
"""
# Configure the metadata manager.
self._metadata_manager = metadata_manager_factory(CoordMetadata)
super().__init__()
if delta and delta.nbounds not in (0, 2):
raise ValueError(
"Invalid delta coordinate: must have either 0 or 2 bounds."
)
if sigma and sigma.nbounds not in (0, 2):
raise ValueError(
"Invalid sigma coordinate: must have either 0 or 2 bounds."
)
if orography and orography.nbounds:
msg = (
"Orography coordinate {!r} has bounds."
" These will be disregarded.".format(orography.name())
)
warn_once_at_level(msg, category=IrisIgnoringBoundsWarning, stacklevel=2)
self.delta = delta
self.sigma = sigma
self.orography = orography
self.standard_name = "altitude"
if delta is None and orography is None:
emsg = "Unable to determine units: no delta or orography available."
raise ValueError(emsg)
if delta and orography and delta.units != orography.units:
emsg = "Incompatible units: delta and orography must have the same units."
raise ValueError(emsg)
self.units = (delta and delta.units) or orography.units
if not self.units.is_convertible("m"):
emsg = (
"Invalid units: delta and/or orography must be expressed "
"in length units."
)
raise ValueError(emsg)
self.attributes = {"positive": "up"}
@property
def dependencies(self):
"""Return a dict mapping from constructor arg names to coordinates.
Return a dictionary mapping from constructor argument names to
the corresponding coordinates.
"""
return {
"delta": self.delta,
"sigma": self.sigma,
"orography": self.orography,
}
def _derive(self, delta, sigma, orography):
return delta + sigma * orography
def make_coord(self, coord_dims_func):
"""Return a new :class:`iris.coords.AuxCoord` as defined by this factory.
Parameters
----------
coord_dims_func :
A callable which can return the list of dimensions relevant
to a given coordinate.
See :meth:`iris.cube.Cube.coord_dims()`.
"""
# Which dimensions are relevant?
derived_dims = self.derived_dims(coord_dims_func)
dependency_dims = self._dependency_dims(coord_dims_func)
# Build the points array.
nd_points_by_key = self._remap(dependency_dims, derived_dims)
points = self._derive(
nd_points_by_key["delta"],
nd_points_by_key["sigma"],
nd_points_by_key["orography"],
)
bounds = None
if (self.delta and self.delta.nbounds) or (self.sigma and self.sigma.nbounds):
# Build the bounds array.
nd_values_by_key = self._remap_with_bounds(dependency_dims, derived_dims)
delta = nd_values_by_key["delta"]
sigma = nd_values_by_key["sigma"]
orography = nd_values_by_key["orography"]
ok_bound_shapes = [(), (1,), (2,)]
if delta.shape[-1:] not in ok_bound_shapes:
raise ValueError("Invalid delta coordinate bounds.")
if sigma.shape[-1:] not in ok_bound_shapes:
raise ValueError("Invalid sigma coordinate bounds.")
if orography.shape[-1:] not in [(), (1,)]:
warn_once_at_level(
"Orography coordinate has bounds. These are being disregarded.",
category=IrisIgnoringBoundsWarning,
stacklevel=2,
)
orography_pts = nd_points_by_key["orography"]
bds_shape = list(orography_pts.shape) + [1]
orography = orography_pts.reshape(bds_shape)
bounds = self._derive(delta, sigma, orography)
hybrid_height = iris.coords.AuxCoord(
points,
standard_name=self.standard_name,
long_name=self.long_name,
var_name=self.var_name,
units=self.units,
bounds=bounds,
attributes=self.attributes,
coord_system=self.coord_system,
)
return hybrid_height
def update(self, old_coord, new_coord=None):
"""Notify the factory of the removal/replacement of a coordinate.
Notify the factory of the removal/replacement of a coordinate
which might be a dependency.
Parameters
----------
old_coord :
The coordinate to be removed/replaced.
new_coord : optional
If None, any dependency using old_coord is removed, otherwise
any dependency using old_coord is updated to use new_coord.
"""
if self.delta is old_coord:
if new_coord and new_coord.nbounds not in (0, 2):
raise ValueError(
"Invalid delta coordinate: must have either 0 or 2 bounds."
)
self.delta = new_coord
elif self.sigma is old_coord:
if new_coord and new_coord.nbounds not in (0, 2):
raise ValueError(
"Invalid sigma coordinate: must have either 0 or 2 bounds."
)
self.sigma = new_coord
elif self.orography is old_coord:
if new_coord and new_coord.nbounds:
msg = (
"Orography coordinate {!r} has bounds."
" These will be disregarded.".format(new_coord.name())
)
warn_once_at_level(
msg, category=IrisIgnoringBoundsWarning, stacklevel=2
)
self.orography = new_coord
class HybridPressureFactory(AuxCoordFactory):
"""Define a hybrid-pressure coordinate factory."""
def __init__(self, delta=None, sigma=None, surface_air_pressure=None):
"""Create a hybrid-height coordinate factory with the following formula.
.. math::
p = ap + b * ps
At least one of `delta` or `surface_air_pressure` must be provided.
Parameters
----------
delta : Coord, optional
The coordinate providing the `ap` term.
sigma : Coord, optional
The coordinate providing the `b` term.
surface_air_pressure : Coord, optional
The coordinate providing the `ps` term.
"""
# Configure the metadata manager.
self._metadata_manager = metadata_manager_factory(CoordMetadata)
super().__init__()
# Check that provided coords meet necessary conditions.
self._check_dependencies(delta, sigma, surface_air_pressure)
self.units = (delta and delta.units) or surface_air_pressure.units
self.delta = delta
self.sigma = sigma
self.surface_air_pressure = surface_air_pressure
self.standard_name = "air_pressure"
self.attributes = {}
@staticmethod
def _check_dependencies(delta, sigma, surface_air_pressure):
# Check for sufficient coordinates.
if delta is None and (sigma is None or surface_air_pressure is None):
msg = (
"Unable to construct hybrid pressure coordinate factory "
"due to insufficient source coordinates."
)
raise ValueError(msg)
# Check bounds.
if delta and delta.nbounds not in (0, 2):
raise ValueError(
"Invalid delta coordinate: must have either 0 or 2 bounds."
)
if sigma and sigma.nbounds not in (0, 2):
raise ValueError(
"Invalid sigma coordinate: must have either 0 or 2 bounds."
)
if surface_air_pressure and surface_air_pressure.nbounds:
msg = (
"Surface pressure coordinate {!r} has bounds. These will"
" be disregarded.".format(surface_air_pressure.name())
)
warn_once_at_level(msg, category=IrisIgnoringBoundsWarning, stacklevel=2)
# Check units.
if sigma is not None and sigma.units.is_unknown():
# Be graceful, and promote unknown to dimensionless units.
sigma.units = cf_units.Unit("1")
if sigma is not None and not sigma.units.is_dimensionless():
raise ValueError("Invalid units: sigma must be dimensionless.")
if (
delta is not None
and surface_air_pressure is not None
and delta.units != surface_air_pressure.units
):
msg = (
"Incompatible units: delta and "
"surface_air_pressure must have the same units."
)
raise ValueError(msg)
if delta is not None:
units = delta.units
else:
units = surface_air_pressure.units
if not units.is_convertible("Pa"):
msg = (
"Invalid units: delta and "
"surface_air_pressure must have units of pressure."
)
raise ValueError(msg)
@property
def dependencies(self):
"""Return a dict mapping from constructor arg names to coordinates.
Returns a dictionary mapping from constructor argument names to
the corresponding coordinates.
"""
return {
"delta": self.delta,
"sigma": self.sigma,
"surface_air_pressure": self.surface_air_pressure,
}
def _derive(self, delta, sigma, surface_air_pressure):
return delta + sigma * surface_air_pressure
def make_coord(self, coord_dims_func):
"""Return a new :class:`iris.coords.AuxCoord` as defined by this factory.
Parameters
----------
coord_dims_func :
A callable which can return the list of dimensions relevant
to a given coordinate.
See :meth:`iris.cube.Cube.coord_dims()`.
"""
# Which dimensions are relevant?
derived_dims = self.derived_dims(coord_dims_func)
dependency_dims = self._dependency_dims(coord_dims_func)
# Build the points array.
nd_points_by_key = self._remap(dependency_dims, derived_dims)
points = self._derive(
nd_points_by_key["delta"],
nd_points_by_key["sigma"],
nd_points_by_key["surface_air_pressure"],
)
bounds = None
if (self.delta and self.delta.nbounds) or (self.sigma and self.sigma.nbounds):
# Build the bounds array.
nd_values_by_key = self._remap_with_bounds(dependency_dims, derived_dims)
delta = nd_values_by_key["delta"]
sigma = nd_values_by_key["sigma"]
surface_air_pressure = nd_values_by_key["surface_air_pressure"]
ok_bound_shapes = [(), (1,), (2,)]
if delta.shape[-1:] not in ok_bound_shapes:
raise ValueError("Invalid delta coordinate bounds.")
if sigma.shape[-1:] not in ok_bound_shapes:
raise ValueError("Invalid sigma coordinate bounds.")
if surface_air_pressure.shape[-1:] not in [(), (1,)]:
warn_once_at_level(
"Surface pressure coordinate has bounds. "
"These are being disregarded.",
category=IrisIgnoringBoundsWarning,
)
surface_air_pressure_pts = nd_points_by_key["surface_air_pressure"]
bds_shape = list(surface_air_pressure_pts.shape) + [1]
surface_air_pressure = surface_air_pressure_pts.reshape(bds_shape)
bounds = self._derive(delta, sigma, surface_air_pressure)
hybrid_pressure = iris.coords.AuxCoord(
points,
standard_name=self.standard_name,
long_name=self.long_name,
var_name=self.var_name,
units=self.units,
bounds=bounds,
attributes=self.attributes,
coord_system=self.coord_system,
)
return hybrid_pressure
class OceanSigmaZFactory(AuxCoordFactory):
"""Defines an ocean sigma over z coordinate factory."""
def __init__(
self,
sigma=None,
eta=None,
depth=None,
depth_c=None,
nsigma=None,
zlev=None,
):
"""Create an ocean sigma over z coordinate factory with the following formula.
if k < nsigma:
.. math::
z(n, k, j, i) = eta(n, j, i) + sigma(k) *
(min(depth_c, depth(j, i)) + eta(n, j, i))
if k >= nsigma:
.. math::
z(n, k, j, i) = zlev(k)
The `zlev` and 'nsigma' coordinates must be provided, and at least
either `eta`, or 'sigma' and `depth` and `depth_c` coordinates.
"""
# Configure the metadata manager.
self._metadata_manager = metadata_manager_factory(CoordMetadata)
super().__init__()
# Check that provided coordinates meet necessary conditions.
self._check_dependencies(sigma, eta, depth, depth_c, nsigma, zlev)
self.units = zlev.units
self.sigma = sigma
self.eta = eta
self.depth = depth
self.depth_c = depth_c
self.nsigma = nsigma
self.zlev = zlev
self.standard_name = "sea_surface_height_above_reference_ellipsoid"
self.attributes = {"positive": "up"}
@staticmethod
def _check_dependencies(sigma, eta, depth, depth_c, nsigma, zlev):
# Check for sufficient factory coordinates.
if zlev is None:
raise ValueError("Unable to determine units: no zlev coordinate available.")
if nsigma is None:
raise ValueError("Missing nsigma coordinate.")
if eta is None and (sigma is None or depth_c is None or depth is None):
msg = (
"Unable to construct ocean sigma over z coordinate "
"factory due to insufficient source coordinates."
)
raise ValueError(msg)
# Check bounds and shape.
for coord, term in ((sigma, "sigma"), (zlev, "zlev")):
if coord is not None and coord.nbounds not in (0, 2):
msg = (
"Invalid {} coordinate {!r}: must have either "
"0 or 2 bounds.".format(term, coord.name())
)
raise ValueError(msg)
if sigma and sigma.nbounds != zlev.nbounds:
msg = (
"The sigma coordinate {!r} and zlev coordinate {!r} "
"must be equally bounded.".format(sigma.name(), zlev.name())
)
raise ValueError(msg)
coords = (
(eta, "eta"),
(depth, "depth"),
(depth_c, "depth_c"),
(nsigma, "nsigma"),
)
for coord, term in coords:
if coord is not None and coord.nbounds:
msg = (
"The {} coordinate {!r} has bounds. "
"These are being disregarded.".format(term, coord.name())
)
warn_once_at_level(
msg, category=IrisIgnoringBoundsWarning, stacklevel=2
)
for coord, term in ((depth_c, "depth_c"), (nsigma, "nsigma")):
if coord is not None and coord.shape != (1,):
msg = "Expected scalar {} coordinate {!r}: got shape {!r}.".format(
term, coord.name(), coord.shape
)
raise ValueError(msg)
# Check units.
if not zlev.units.is_convertible("m"):
msg = (
"Invalid units: zlev coordinate {!r} "
"must have units of distance.".format(zlev.name())
)
raise ValueError(msg)
if sigma is not None and sigma.units.is_unknown():
# Be graceful, and promote unknown to dimensionless units.
sigma.units = cf_units.Unit("1")
if sigma is not None and not sigma.units.is_dimensionless():
msg = "Invalid units: sigma coordinate {!r} must be dimensionless.".format(
sigma.name()
)