|
| 1 | +# -*- coding: utf-8 -*- |
| 2 | +# ----------------------------------------------------------------------------- |
| 3 | +# (C) British Crown Copyright 2017-2021 Met Office. |
| 4 | +# All rights reserved. |
| 5 | +# |
| 6 | +# Redistribution and use in source and binary forms, with or without |
| 7 | +# modification, are permitted provided that the following conditions are met: |
| 8 | +# |
| 9 | +# * Redistributions of source code must retain the above copyright notice, this |
| 10 | +# list of conditions and the following disclaimer. |
| 11 | +# |
| 12 | +# * Redistributions in binary form must reproduce the above copyright notice, |
| 13 | +# this list of conditions and the following disclaimer in the documentation |
| 14 | +# and/or other materials provided with the distribution. |
| 15 | +# |
| 16 | +# * Neither the name of the copyright holder nor the names of its |
| 17 | +# contributors may be used to endorse or promote products derived from |
| 18 | +# this software without specific prior written permission. |
| 19 | +# |
| 20 | +# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" |
| 21 | +# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE |
| 22 | +# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE |
| 23 | +# ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE |
| 24 | +# LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR |
| 25 | +# CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF |
| 26 | +# SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS |
| 27 | +# INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN |
| 28 | +# CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) |
| 29 | +# ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE |
| 30 | +# POSSIBILITY OF SUCH DAMAGE. |
| 31 | +"""Module containing a plugin to calculate the modal weather code in a period.""" |
| 32 | + |
| 33 | +import numpy as np |
| 34 | +from iris.analysis import Aggregator |
| 35 | +from iris.cube import Cube, CubeList |
| 36 | +from numpy import ndarray |
| 37 | +from scipy import stats |
| 38 | + |
| 39 | +from improver import BasePlugin |
| 40 | +from improver.utilities.cube_manipulation import MergeCubes |
| 41 | + |
| 42 | +from .utilities import DAYNIGHT_CODES, GROUPED_CODES |
| 43 | + |
| 44 | +CODE_MAX = 100 |
| 45 | +UNSET_CODE_INDICATOR = -99 |
| 46 | + |
| 47 | + |
| 48 | +class ModalWeatherCode(BasePlugin): |
| 49 | + """Plugin that returns the modal code over the period spanned by the |
| 50 | + input data. In cases of a tie in the mode values, scipy returns the smaller |
| 51 | + value. The opposite is desirable in this case as the significance / |
| 52 | + importance of the weather codes generally increases with the value. To |
| 53 | + achieve this the codes are subtracted from an arbitrarily larger |
| 54 | + number prior to calculating the mode, and this operation reversed in the |
| 55 | + final output. |
| 56 | +
|
| 57 | + If there are many different codes for a single point over the time |
| 58 | + spanned by the input cubes it may be that the returned mode is not robust. |
| 59 | + Given the preference to return more significant codes explained above, |
| 60 | + a 12 hour period with 12 different codes, one of which is thunder, will |
| 61 | + return a thunder code to describe the whole period. This is likely not a |
| 62 | + good representation. In these cases grouping is used to try and select |
| 63 | + a suitable weather code (e.g. a rain shower if the codes include a mix of |
| 64 | + rain showers and dynamic rain) by providing a more robust mode. The lowest |
| 65 | + number (least significant) member of the group is returned as the code. |
| 66 | + Use of the least significant member reflects the lower certainty in the |
| 67 | + forecasts. |
| 68 | +
|
| 69 | + Where there are different weather codes available for night and day, the |
| 70 | + modal code returned is always a day code, regardless of the times |
| 71 | + covered by the input files. |
| 72 | + """ |
| 73 | + |
| 74 | + def __init__(self): |
| 75 | + """Create an aggregator instance for reuse""" |
| 76 | + self.aggregator_instance = Aggregator("mode", self.mode_aggregator) |
| 77 | + |
| 78 | + @staticmethod |
| 79 | + def _unify_day_and_night(cube: Cube): |
| 80 | + """Remove distinction between day and night codes so they can each |
| 81 | + contribute when calculating the modal code. The cube of weather |
| 82 | + codes is modified in place with all night codes made into their |
| 83 | + daytime equivalents. |
| 84 | +
|
| 85 | + Args: |
| 86 | + A cube of weather codes. |
| 87 | + """ |
| 88 | + night_codes = np.array(DAYNIGHT_CODES) - 1 |
| 89 | + for code in night_codes: |
| 90 | + cube.data[cube.data == code] += 1 |
| 91 | + |
| 92 | + @staticmethod |
| 93 | + def _group_codes(modal: Cube, cube: Cube): |
| 94 | + """In instances where the mode returned is not significant, i.e. the |
| 95 | + weather code chosen occurs infrequently in the period, the codes can be |
| 96 | + grouped to yield a more definitive period code. Given the uncertainty, |
| 97 | + the least significant weather type (lowest number in a group that is |
| 98 | + found in the data) is used to replace the other data values that belong |
| 99 | + to that group prior to recalculating the modal code. |
| 100 | +
|
| 101 | + The modal cube is modified in place. |
| 102 | +
|
| 103 | + Args: |
| 104 | + modal: |
| 105 | + The modal weather code cube which contains UNSET_CODE_INDICATOR |
| 106 | + values that need to be replaced with a more definitive period |
| 107 | + code. |
| 108 | + cube: |
| 109 | + The original input data. Data relating to unset points will be |
| 110 | + grouped and the mode recalculated.""" |
| 111 | + |
| 112 | + undecided_points = np.argwhere(modal.data == UNSET_CODE_INDICATOR) |
| 113 | + |
| 114 | + for point in undecided_points: |
| 115 | + data = cube.data[(..., *point)].copy() |
| 116 | + |
| 117 | + for _, codes in GROUPED_CODES.items(): |
| 118 | + default_code = sorted([code for code in data if code in codes]) |
| 119 | + if default_code: |
| 120 | + data[np.isin(data, codes)] = default_code[0] |
| 121 | + mode_result, counts = stats.mode(CODE_MAX - data) |
| 122 | + modal.data[tuple(point)] = CODE_MAX - mode_result |
| 123 | + |
| 124 | + @staticmethod |
| 125 | + def mode_aggregator(data: ndarray, axis: int) -> ndarray: |
| 126 | + """An aggregator for use with iris to calculate the mode along the |
| 127 | + specified axis. If the modal value selected comprises less than 10% |
| 128 | + of data along the dimension being collapsed, the value is set to the |
| 129 | + UNSET_CODE_INDICATOR to indicate that the uncertainty was too high to |
| 130 | + return a mode. |
| 131 | +
|
| 132 | + Args: |
| 133 | + data: |
| 134 | + The data for which a mode is to be calculated. |
| 135 | + axis: |
| 136 | + The axis / dimension over which to calculate the mode. |
| 137 | +
|
| 138 | + Returns: |
| 139 | + The data array collapsed over axis, containing the calculated modes. |
| 140 | + """ |
| 141 | + # Iris aggregators support indexing from the end of the array. |
| 142 | + if axis < 0: |
| 143 | + axis += data.ndim |
| 144 | + # Aggregation coordinate is moved to the -1 position in initialisation; |
| 145 | + # move this back to the leading coordinate |
| 146 | + data = np.moveaxis(data, [axis], [0]) |
| 147 | + minimum_significant_count = 0.1 * data.shape[0] |
| 148 | + mode_result, counts = stats.mode(CODE_MAX - data, axis=0) |
| 149 | + mode_result[counts < minimum_significant_count] = ( |
| 150 | + CODE_MAX - UNSET_CODE_INDICATOR |
| 151 | + ) |
| 152 | + return CODE_MAX - np.squeeze(mode_result) |
| 153 | + |
| 154 | + def process(self, cubes: CubeList) -> Cube: |
| 155 | + """Calculate the modal weather code, with handling for edge cases. |
| 156 | +
|
| 157 | + Args: |
| 158 | + cubes: |
| 159 | + A list of weather code cubes at different times. A modal |
| 160 | + code will be calculated over the time coordinate to return |
| 161 | + the most comon code, which is taken to be the best |
| 162 | + representation of the whole period. |
| 163 | +
|
| 164 | + Returns: |
| 165 | + A single weather code cube with time bounds that span those of |
| 166 | + the input weather code cubes. |
| 167 | + """ |
| 168 | + # Handle case in which a single time is provided. |
| 169 | + if len(cubes) == 1: |
| 170 | + return cubes[0] |
| 171 | + |
| 172 | + cube = MergeCubes()(cubes) |
| 173 | + self._unify_day_and_night(cube) |
| 174 | + |
| 175 | + result = cube.collapsed("time", self.aggregator_instance) |
| 176 | + result.coord("time").points = result.coord("time").bounds[0][-1] |
| 177 | + |
| 178 | + # Handle any unset points where it was hard to determine a suitable mode |
| 179 | + if (result.data == UNSET_CODE_INDICATOR).any(): |
| 180 | + self._group_codes(result, cube) |
| 181 | + |
| 182 | + return result |
0 commit comments