-
Notifications
You must be signed in to change notification settings - Fork 19
/
regridoperator.py
811 lines (610 loc) · 24 KB
/
regridoperator.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
from copy import deepcopy
import numpy as np
from cfdm import Container
from ..decorators import _display_or_return
from ..functions import _DEPRECATION_ERROR_ATTRIBUTE, _DEPRECATION_ERROR_METHOD
from ..mixin_container import Container as mixin_Container
class RegridOperator(mixin_Container, Container):
"""A regridding operator between two grids.
Regridding is the process of interpolating from one grid
resolution to a different grid resolution.
The regridding operator stores the regridding weights; auxiliary
information, such as the grid shapes; the CF metadata for the
destination grid; and the source grid coordinates.
.. versionadded:: 3.10.0
"""
def __init__(
self,
weights=None,
row=None,
col=None,
coord_sys=None,
method=None,
src_shape=None,
dst_shape=None,
src_cyclic=None,
dst_cyclic=None,
src_mask=None,
dst_mask=None,
src_coords=None,
src_bounds=None,
start_index=0,
src_axes=None,
dst_axes=None,
dst=None,
weights_file=None,
src_mesh_location=None,
dst_mesh_location=None,
src_featureType=None,
dst_featureType=None,
dimensionality=None,
src_z=None,
dst_z=None,
ln_z=False,
):
"""**Initialisation**
:Parameters:
weights: array_like
The 1-d array of regridding weights for locations in
the 2-d dense weights matrix. The locations are defined
by the *row* and *col* parameters.
row, col: array_like, array_like
The 1-d arrays of the row and column indices of the
regridding weights in the dense weights matrix, which
has J rows and I columns, where J and I are the total
number of cells in the destination and source grids
respectively. See the *start_index* parameter.
coord_sys: `str`
The name of the coordinate system of the source and
destination grids. Either ``'spherical'`` or
``'Cartesian'``.
method: `str`
The name of the regridding method.
src_shape: sequence of `int`
The shape of the source grid.
dst_shape: sequence of `int`
The shape of the destination grid.
src_cyclic: `bool`
For spherical regridding, specifies whether or not the
source grid longitude axis is cyclic.
dst_cyclic: `bool`
For spherical regridding, specifies whether or not the
destination grid longitude axis is cyclic.
src_mask: `numpy.ndarray` or `None`, optional
If a `numpy.ndarray` with shape *src_shape* then this
is the source grid mask that was used during the
creation of the *weights*. If *src_mask* is a scalar
array with value `False`, then this is equivalent to a
source grid mask with shape *src_shape* entirely
populated with `False`.
If `None` (the default), then the weights are assumed
to have been created assuming no source grid masked
cells.
dst_mask: `numpy.ndarray` or `None`, optional
A destination grid mask to be applied to the weights
matrix, in addition to those destination grid cells
that have no non-zero weights. If `None` (the default)
then no additional destination grid cells are
masked. If a Boolean `numpy` array then it must have
shape *dst_shape*, and a value of `True` signifies a
masked destination grid cell.
start_index: `int`, optional
Specify whether the *row* and *col* parameters use 0-
or 1-based indexing. By default 0-based indexing is
used.
If *row* and *col* are to be read from a weights file
and their netCDF variables have ``start_index``
attributes, then these will be used in preference to
*start_index*.
parameters: Deprecated at version 3.14.0
Use keyword parameters instead.
dst: `Field` or `Domain`
The definition of the destination grid.
dst_axes: `dict` or sequence or `None`, optional
The destination grid axes to be regridded.
src_axes: `dict` or sequence or `None`, optional
The source grid axes to be regridded.
weights_file: `str` or `None`, optional
Path to a netCDF file that contained the regridding
weights. If `None`, the default, then the weights
were computed rather than read from a file.
.. versionadded:: 3.15.2
src_mesh_location: `str`, optional
The UGRID mesh element of the source grid
(e.g. ``'face'``).
.. versionadded:: 3.16.0
dst_mesh_location: `str`, optional
The UGRID mesh element of the destination grid
(e.g. ``'face'``).
.. versionadded:: 3.16.2
src_featureType: `str`, optional
The discrete sampling geometry (DSG) featureType of
the source grid (e.g. ``'trajectory'``).
.. versionadded:: 3.16.2
dst_featureType: `str`, optional
The DSG featureType of the destination grid
(e.g. ``'trajectory'``).
.. versionadded:: 3.16.2
src_z: optional
The identity of the source grid vertical coordinates
used to calculate the weights. If `None` then no
source grid vertical axis is identified.
.. versionadded:: 3.16.2
dst_z: optional
The identity of the destination grid vertical
coordinates used to calculate the weights. If `None`
then no destination grid vertical axis is identified.
.. versionadded:: 3.16.2
ln_z: `bool`, optional
Whether or not the weights were calculated with the
natural logarithm of vertical coordinates.
.. versionadded:: 3.16.2
dimensionality: `int`, optional
The number of physical regridding dimensions. This may
differ from the corresponding number of storage
dimensions in the source or destination grids, if
either has an unstructured mesh or a DSG featureType.
.. versionadded:: 3.16.2
"""
super().__init__()
if weights is None and weights_file is None:
# This to allow a no-arg init!
return
self._set_component("weights", weights, copy=False)
self._set_component("row", row, copy=False)
self._set_component("col", col, copy=False)
self._set_component("coord_sys", coord_sys, copy=False)
self._set_component("method", method, copy=False)
self._set_component("src_mask", src_mask, copy=False)
self._set_component("dst_mask", dst_mask, copy=False)
self._set_component("src_cyclic", bool(src_cyclic), copy=False)
self._set_component("dst_cyclic", bool(dst_cyclic), copy=False)
self._set_component("src_shape", tuple(src_shape), copy=False)
self._set_component("dst_shape", tuple(dst_shape), copy=False)
self._set_component("src_coords", tuple(src_coords), copy=False)
self._set_component("src_bounds", tuple(src_bounds), copy=False)
self._set_component("start_index", int(start_index), copy=False)
self._set_component("src_axes", src_axes, copy=False)
self._set_component("dst_axes", dst_axes, copy=False)
self._set_component("dst", dst, copy=False)
self._set_component("weights_file", weights_file, copy=False)
self._set_component("src_mesh_location", src_mesh_location, copy=False)
self._set_component("dst_mesh_location", dst_mesh_location, copy=False)
self._set_component("src_featureType", src_featureType, copy=False)
self._set_component("dst_featureType", dst_featureType, copy=False)
self._set_component("dimensionality", dimensionality, copy=False)
self._set_component("src_z", src_z, copy=False)
self._set_component("dst_z", dst_z, copy=False)
self._set_component("ln_z", bool(ln_z), copy=False)
def __repr__(self):
"""x.__repr__() <==> repr(x)"""
return (
f"<CF {self.__class__.__name__}: {self.coord_sys} {self.method}>"
)
@property
def col(self):
"""The 1-d array of the column indices of the regridding
weights.
.. versionadded:: 3.14.0
"""
return self._get_component("col")
@property
def coord_sys(self):
"""The name of the regridding coordinate system.
.. versionadded:: 3.14.0
"""
return self._get_component("coord_sys")
@property
def dimensionality(self):
"""The number of physical regridding dimensions.
.. versionadded:: 3.16.2
"""
return self._get_component("dimensionality")
@property
def dst(self):
"""The definition of the destination grid.
Either a `Field` or` Domain`.
.. versionadded:: 3.14.0
"""
return self._get_component("dst")
@property
def dst_axes(self):
"""The destination grid axes to be regridded.
If not required then this will be `None`.
.. versionadded:: 3.14.0
"""
return self._get_component("dst_axes")
@property
def dst_cyclic(self):
"""Whether or not the destination grid longitude axis is
cyclic."""
return self._get_component("dst_cyclic")
@property
def dst_featureType(self):
"""The DSG featureType of the destination grid.
.. versionadded:: 3.16.2
"""
return self._get_component("dst_featureType")
@property
def dst_mask(self):
"""A destination grid mask to be applied to the weights matrix.
If `None` then no additional destination grid cells are
masked.
If a Boolean `numpy` array then it is required that this mask
is applied to the weights matrix prior to use in a regridding
operation. The mask must have shape `!dst_shape`, and a value
of `True` signifies a masked destination grid cell.
.. versionadded:: 3.14.0
"""
return self._get_component("dst_mask")
@property
def dst_mesh_location(self):
"""The UGRID mesh element of the destination grid.
.. versionadded:: 3.16.0
"""
return self._get_component("dst_mesh_location")
@property
def dst_shape(self):
"""The shape of the destination grid.
.. versionadded:: 3.14.0
"""
return self._get_component("dst_shape")
@property
def dst_z(self):
"""The identity of the destination grid vertical coordinates.
.. versionadded:: 3.16.2
"""
return self._get_component("dst_z")
@property
def ln_z(self):
"""Whether or not vertical weights are based on ln(z).
.. versionadded:: 3.16.2
"""
return self._get_component("ln_z")
@property
def method(self):
"""The name of the regridding method.
.. versionadded:: 3.14.0
"""
return self._get_component("method")
@property
def name(self):
"""The name of the regridding method."""
_DEPRECATION_ERROR_ATTRIBUTE(
self,
"name",
version="3.14.0",
removed_at="5.0.0",
)
@property
def row(self):
"""The 1-d array of the row indices of the regridding weights.
.. versionadded:: 3.14.0
"""
return self._get_component("row")
@property
def src_axes(self):
"""The source grid axes to be regridded.
If not required then this will be `None`.
.. versionadded:: 3.14.0
"""
return self._get_component("src_axes")
@property
def src_bounds(self):
"""The bounds of the source grid cells.
.. versionadded:: 3.14.0
"""
return self._get_component("src_bounds")
@property
def src_coords(self):
"""The coordinates of the source grid cells.
.. versionadded:: 3.14.0
"""
return self._get_component("src_coords")
@property
def src_cyclic(self):
"""Whether or not the source grid longitude axis is cyclic.
.. versionadded:: 3.14.0
"""
return self._get_component("src_cyclic")
@property
def src_featureType(self):
"""The DSG featureType of the source grid.
.. versionadded:: 3.16.2
"""
return self._get_component("src_featureType")
@property
def src_mask(self):
"""The source grid mask that was applied during the weights
creation.
If `None` then the weights are assumed to have been created
assuming no source grid masked cells.
If a Boolean `numpy.ndarray` with shape `!src_shape` then this
is the source grid mask that was used during the creation of
the *weights*.
.. versionadded:: 3.14.0
"""
return self._get_component("src_mask")
@property
def src_mesh_location(self):
"""The UGRID mesh element of the source grid.
.. versionadded:: 3.16.0
"""
return self._get_component("src_mesh_location")
@property
def src_shape(self):
"""The shape of the source grid.
.. versionadded:: 3.14.0
"""
return self._get_component("src_shape")
@property
def src_z(self):
"""The identity of the source grid vertical coordinates.
.. versionadded:: 3.16.2
"""
return self._get_component("src_z")
@property
def start_index(self):
"""The start index of the row and column indices.
If `row` and `col` are to be read from a weights file and
their netCDF variables have ``start_index`` attributes, then
these will be used in preference to `start_index`.
.. versionadded:: 3.14.0
"""
return self._get_component("start_index")
@property
def weights(self):
"""The 1-d array of the regridding weights, or `None`.
If and only if it is a `scipy` sparse array that combines the
weights and the row and column indices (as opposed to a
`numpy` array of just the weights) then the `dst_mask`
attribute will have been updated to `True` for destination
grid points for which the weights are all zero.
.. versionadded:: 3.14.0
"""
return self._get_component("weights")
@property
def weights_file(self):
"""The file which contains the weights, or `None`.
.. versionadded:: 3.15.2
"""
return self._get_component("weights_file")
def copy(self):
"""Return a deep copy.
:Returns:
`RegridOperator`
The deep copy.
"""
row = self.row
if row is not None:
row = row.copy()
col = self.col
if col is not None:
col = col.copy()
return type(self)(
self.weights.copy(),
row,
col,
method=self.method,
src_shape=self.src_shape,
dst_shape=self.dst_shape,
src_cyclic=self.src_cyclic,
dst_cyclic=self.dst_cyclic,
src_mask=deepcopy(self.src_mask),
dst_mask=deepcopy(self.dst_mask),
src_coords=deepcopy(self.src_coords),
src_bounds=deepcopy(self.src_bounds),
coord_sys=self.coord_sys,
start_index=self.start_index,
src_axes=self.src_axes,
dst_axes=self.dst_axes,
dst=self.dst.copy(),
weights_file=self.weights_file,
src_mesh_location=self.src_mesh_location,
dst_mesh_location=self.dst_mesh_location,
src_featureType=self.src_featureType,
dst_featureType=self.dst_featureType,
src_z=self.src_z,
dst_z=self.dst_z,
ln_z=self.ln_z,
)
@_display_or_return
def dump(self, display=True):
"""A full description of the regrid operator.
Returns a description of all properties, including
the weights and their row and column indices.
.. versionadded:: 3.14.0
:Parameters:
display: `bool`, optional
If False then return the description as a string. By
default the description is printed.
:Returns:
{{returns dump}}
"""
_title = repr(self)[4:-1]
line = f"{''.ljust(len(_title), '-')}"
string = [line, _title, line]
for attr in (
"coord_sys",
"method",
"dimensionality",
"src_shape",
"dst_shape",
"src_cyclic",
"dst_cyclic",
"src_mask",
"dst_mask",
"src_coords",
"src_bounds",
"start_index",
"src_axes",
"dst_axes",
"src_mesh_location",
"dst_mesh_location",
"src_featureType",
"dst_featureType",
"src_z",
"dst_z",
"ln_z",
"dst",
"weights",
"row",
"col",
"weights_file",
):
string.append(f"{attr}: {getattr(self, attr)!r}")
return "\n".join(string)
def get_parameter(self, parameter, *default):
"""Return a regrid operation parameter.
Deprecated at version 3.14.0.
:Parameters:
parameter: `str`
The name of the parameter.
default: optional
Return the value of the *default* parameter if the
parameter has not been set.
If set to an `Exception` instance then it will be
raised instead.
.. versionadded:: 3.14.0
:Returns:
The value of the named parameter or the default value, if
set.
**Examples**
>>> r.get_parameter('dst_axes')
['domainaxis1', 'domainaxis0']
>>> r.get_parameter('x')
Traceback
...
ValueError: RegridOperator has no 'x' parameter
>>> r.get_parameter('x', 'missing')
'missing'
"""
_DEPRECATION_ERROR_METHOD(
self,
"get_parameter",
message="Use attributes directly.",
version="3.14.0",
removed_at="5.0.0",
)
def parameters(self):
"""Get the CF metadata parameters for the destination grid.
Deprecated at version 3.14.0.
Any parameter names and values are allowed, and it is assumed
that these are well defined during the creation and
subsequent use of a `RegridOperator` instance.
:Returns:
`dict`
The parameters.
**Examples**
>>> r.parameters()
{'dst': <CF Domain: {latitude(5), longitude(8), time(1)}>,
'dst_axes': ['domainaxis1', 'domainaxis2'],
'src_axes': None}
"""
_DEPRECATION_ERROR_METHOD(
self,
"parameters",
version="3.14.0",
removed_at="5.0.0",
)
def tosparse(self):
"""Convert the weights to `scipy` sparse array format in-place.
The `weights` attribute is set to a Compressed Sparse Row
(CSR) array (i.e. a `scipy.sparse._arrays.csr_array` instance)
that combines the weights and the row and column indices, and
the `row` and `col` attributes are set to `None`.
The `dst_mask` attribute is also updated to `True` for
destination grid points for which the weights are all zero.
A CSR array is used as the most efficient sparse array type
given that we expect no changes to the sparsity structure, and
any further modification of the weights to account for missing
values in the source grid will always involve row-slicing.
If the weights are already in a sparse array format then no
action is taken.
.. versionadded:: 3.13.0
:Returns:
`None`
"""
weights = self.weights
row = self.row
col = self.col
if weights is not None and row is None and col is None:
# Weights are already in sparse array format
return
from math import prod
from scipy.sparse import csr_array
start_index = self.start_index
col_start_index = None
row_start_index = None
if weights is None:
weights_file = self.weights_file
if weights_file is not None:
# Read the weights from the weights file
from netCDF4 import Dataset
from ..data.array.netcdfarray import _lock
_lock.acquire()
nc = Dataset(weights_file, "r")
weights = nc.variables["S"][...]
row = nc.variables["row"][...]
col = nc.variables["col"][...]
try:
col_start_index = nc.variables["col"].start_index
except AttributeError:
col_start_index = 1
try:
row_start_index = nc.variables["row"].start_index
except AttributeError:
row_start_index = 1
nc.close()
_lock.release()
else:
raise ValueError(
"Conversion to sparse array format requires at least "
"one of the 'weights' or 'weights_file' attributes to "
"be set"
)
# Convert to sparse array format
if col_start_index:
col = col - col_start_index
elif start_index:
col = col - start_index
if row_start_index:
row = row - row_start_index
elif start_index:
row = row - start_index
src_size = prod(self.src_shape)
dst_size = prod(self.dst_shape)
weights = csr_array((weights, (row, col)), shape=[dst_size, src_size])
self._set_component("weights", weights, copy=False)
self._set_component("row", None, copy=False)
self._set_component("col", None, copy=False)
del row, col
# Set the destination grid mask to True where the weights for
# destination grid points are all zero
dst_mask = self.dst_mask
if dst_mask is not None:
if dst_mask.dtype != bool or dst_mask.shape != self.dst_shape:
raise ValueError(
f"The {self.__class__.__name__}.dst_mask attribute must "
"be None or a Boolean numpy array with shape "
f"{self.dst_shape}. Got: dtype={dst_mask.dtype}, "
f"shape={dst_mask.shape}"
)
dst_mask = np.array(dst_mask).reshape((dst_size,))
else:
dst_mask = np.zeros((dst_size,), dtype=bool)
# Performance note:
#
# It is much more efficient to access 'weights.indptr' and
# 'weights.data' directly, rather than iterating over rows of
# 'weights' and using 'weights.getrow'.
indptr = weights.indptr.tolist()
data = weights.data
for j, (i0, i1) in enumerate(zip(indptr[:-1], indptr[1:])):
if not data[i0:i1].size:
dst_mask[j] = True
if not dst_mask.any():
dst_mask = None
else:
dst_mask = dst_mask.reshape(self.dst_shape)
self._set_component("dst_mask", dst_mask, copy=False)