Tmean_Veg.py 8.84 KB
Newer Older
Klaus Zimmermann's avatar
Klaus Zimmermann committed
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
import glob
import numpy as np

import iris
import iris.coords
import iris.coord_categorisation
from iris.analysis import Aggregator
from iris.util import rolling_window
from cf_units import Unit

MISSVAL = 1.0e20


def veg_bounds(indata, axis, condition, start_end):
    """
    Args:
    * indata (array):
        thresholded (indicator) data.
    * axis (int):
        number of the array dimension mapping the time sequences.
        (Can also be negative, e.g. '-1' means last dimension)
    * condition (int):
        minimum number of consecutive days to start the vegetation period.
    * start_end (string):
        "start" | "end" indicating which bound
    """

    se = start_end.lower()
    if se.startswith('s'):
        data = indata
        s = True
    elif se.startswith('e'):
        data = np.flip(indata, axis=axis)
        s = False
    else:
        raise RuntimeError("Argument start_end not understood <{}>."
                           "".format(start_end))
    if axis < 0:
        # just cope with negative axis numbers
        axis += indata.ndim

    year_length = indata.shape[-1]
    # Make an array with data values "windowed" along the time axis.
    hit_windows = rolling_window(data, window=condition, axis=axis)
    # Find the windows "full of True-s" (along the added 'window axis').
    full_windows = np.all(hit_windows, axis=axis+1)
    # Find first occurence fulfilling the condition (along the time axis).
    doy = np.argmax(full_windows, axis=axis)
    # mask gridcells where doy == 0 as a result of the condition not being met
    ind = map(np.squeeze,
              np.split(np.mgrid[:doy.shape[0], :doy.shape[1]], 2, 0)) + [doy]
    mask = (~full_windows[ind]) | (doy > np.ceil(year_length / 2))
    masked_doy = np.ma.masked_array(np.float32(doy), mask)
    # convert to day of year in the interval 1..year_length
    # for end calculations count from end of year
    if s:
        masked_doy = masked_doy + 1
    else:
        masked_doy = year_length - masked_doy
    return masked_doy


def veg_gdd(data, start, end):
    doy = np.arange(data.shape[0]) + 1
    dayix = doy[:, None, None]
    L = (dayix < start[None, :, :]) | (dayix > end[None, :, :])
    data.data[L] = 0.0
    data.data = np.maximum(data.data - 5, 0.0)
    gdd = data.collapsed('time', iris.analysis.SUM)
    return gdd, data


def main():
    # Make an aggregator from the user function.
    VEG_BOUNDS = Aggregator('veg_bounds',
                            veg_bounds,
                            units_func=lambda units: 1)

    # Define the parameters of the index.
    threshold_temperature = 5.0
    consecutive_days = 4

    var = 'tas'
    indir = '/home/rossby/imports/obs/UERRA/SURFEX-MESCAN/input/day/'
    inroot = '_EUR-055_UERRA-HARMONIE_RegRean_v1d1-v1d2_SURFEX-MESCAN_v1_day_*.nc'  # noqa
    filepattern = '{}{}{}'.format(indir, var, inroot)
    filelist = sorted(glob.glob(filepattern))

    creation_date = ''
    tracking_id = ''
    vegStart = iris.cube.CubeList([])
    vegEnd = iris.cube.CubeList([])
    vegLen = iris.cube.CubeList([])
    vegGDD = iris.cube.CubeList([])

    for f in filelist:  # one file / year at the time ...
        datacube = iris.load_cube(f)
        creation_date += datacube.attributes.pop('creation_date', None) + ' '
        tracking_id += datacube.attributes.pop('tracking_id', None) + ' '
        datacube.convert_units(Unit('degC'))
        iris.coord_categorisation.add_day_of_year(datacube, 'time', name='doy')
        iris.coord_categorisation.add_year(datacube, 'time', name='year')
        gtThreshold = datacube.copy()
        gtThreshold.data = gtThreshold.data > threshold_temperature

        Vstart = gtThreshold.collapsed('time', VEG_BOUNDS,
                                       condition=consecutive_days,
                                       start_end='start')

        mask = np.ma.getmaskarray(Vstart.data)
        # iris.util.mask_cube(Vstart, mask)
        Vend = gtThreshold.collapsed('time', VEG_BOUNDS,
                                     condition=consecutive_days,
                                     start_end='end')
        mask = mask | np.ma.getmaskarray(Vend.data)
        Vend.data = np.ma.masked_array(Vend.data, mask)
        # iris.util.mask_cube(Vend, mask)
        Vlen = Vend - Vstart + 1
        Vlen.data[mask] = 0
        (Vgdd, testgdd) = veg_gdd(datacube.copy(), Vstart.data, Vend.data)
        Vgdd.data[mask] = 0
        Vstart.remove_coord('doy')
        Vend.remove_coord('doy')
        Vlen.remove_coord('doy')
        Vgdd.remove_coord('doy')
        Vstart.remove_coord('year')
        Vend.remove_coord('year')
        Vlen.remove_coord('year')
        Vgdd.remove_coord('year')
        vegStart.append(Vstart.copy())
        vegEnd.append(Vend.copy())
        vegLen.append(Vlen.copy())
        vegGDD.append(Vgdd.copy())

    dayVegStart5_ANN = vegStart.merge_cube()
    dayVegStart5_ANN.units = Unit('days')
    dayVegStart5_ANN._FillValue = MISSVAL
    dayVegStart5_ANN.valid_min = 1.
    dayVegStart5_ANN.valid_max = 183.
    dayVegStart5_ANN.standard_name = None
    dayVegStart5_ANN.long_name = u'Start of Vegetation Period (first 4 days above 5 degC)'  # noqa
    dayVegStart5_ANN.var_name = u'dayVegStart5_ANN'
    dayVegStart5_ANN.attributes = datacube.attributes
    dayVegStart5_ANN.attributes['creation_date'] = creation_date
    dayVegStart5_ANN.attributes['tracking_id'] = tracking_id

    dayVegEnd5_ANN = vegEnd.merge_cube()
    dayVegEnd5_ANN.units = Unit('days')
    dayVegEnd5_ANN._FillValue = MISSVAL
    dayVegEnd5_ANN.valid_min = 183.
    dayVegEnd5_ANN.valid_max = 366.
    dayVegEnd5_ANN.standard_name = None
    dayVegEnd5_ANN.long_name = u'End of Vegetation Period (last 4 days above 5 degC)'  # noqa
    dayVegEnd5_ANN.var_name = u'dayVegEnd5_ANN'
    dayVegEnd5_ANN.attributes = datacube.attributes
    dayVegEnd5_ANN.attributes['creation_date'] = creation_date
    dayVegEnd5_ANN.attributes['tracking_id'] = tracking_id

    dayVegLen5_ANN = vegLen.merge_cube()
    dayVegLen5_ANN.units = Unit('days')
    dayVegLen5_ANN._FillValue = MISSVAL
    dayVegLen5_ANN.valid_min = 0.
    dayVegLen5_ANN.valid_max = 366.
    dayVegLen5_ANN.standard_name = None
    dayVegLen5_ANN.long_name = u'Length of Vegetation Period'  # noqa
    dayVegLen5_ANN.var_name = u'dayVegLen5_ANN'
    dayVegLen5_ANN.attributes = datacube.attributes
    dayVegLen5_ANN.attributes['creation_date'] = creation_date
    dayVegLen5_ANN.attributes['tracking_id'] = tracking_id

    GDD5_ANN = vegGDD.merge_cube()
    GDD5_ANN.units = Unit('degC days')
    GDD5_ANN._FillValue = MISSVAL
    GDD5_ANN.valid_min = 0.
    GDD5_ANN.valid_max = 15000.
    GDD5_ANN.standard_name = None
    GDD5_ANN.long_name = u'Growing Degree Days (over 5 degC) during Vegetation Period'  # noqa
    GDD5_ANN.var_name = u'GDD5_ANN'
    GDD5_ANN.attributes = datacube.attributes
    GDD5_ANN.attributes['creation_date'] = creation_date
    GDD5_ANN.attributes['tracking_id'] = tracking_id

    outdir = '/home/sm_lbarr/PROJ/B4EST/'
    outroot = '_EUR-055_UERRA-HARMONIE_RegRean_v1d1-v1d2_SURFEX-MESCAN_v1_ann_1961-2015.nc'  # noqa

    iris.coord_categorisation.add_year(dayVegStart5_ANN, 'time', name='year')  # noqa
    iris.coord_categorisation.add_year(dayVegEnd5_ANN, 'time', name='year')
    iris.coord_categorisation.add_year(dayVegLen5_ANN, 'time', name='year')
    iris.coord_categorisation.add_year(GDD5_ANN, 'time', name='year')
    iris.save(dayVegStart5_ANN, '{}{}{}'.format(outdir, 'dayVegStart5_ANN', outroot), fill_value=MISSVAL)  # noqa
    iris.save(dayVegEnd5_ANN, '{}{}{}'.format(outdir, 'dayVegEnd5_ANN', outroot), fill_value=MISSVAL)  # noqa
    iris.save(dayVegLen5_ANN, '{}{}{}'.format(outdir, 'dayVegLen5_ANN', outroot), fill_value=MISSVAL)  # noqa
    iris.save(GDD5_ANN, '{}{}{}'.format(outdir, 'GDD5_ANN', outroot), fill_value=MISSVAL)  # noqa

    c = dayVegStart5_ANN.extract(iris.Constraint(
        coord_values={'year': lambda cell: 1960 < cell < 1991}))
    dayVegStart5_6190 = c.collapsed('year', iris.analysis.MEAN)
    c = dayVegEnd5_ANN.extract(iris.Constraint(
        coord_values={'year': lambda cell: 1960 < cell < 1991}))
    dayVegEnd5_6190 = c.collapsed('year', iris.analysis.MEAN)
    c = dayVegLen5_ANN.extract(iris.Constraint(
        coord_values={'year': lambda cell: 1960 < cell < 1991}))
    dayVegLen5_6190 = c.collapsed('year', iris.analysis.MEAN)
    c = GDD5_ANN.extract(iris.Constraint(
        coord_values={'year': lambda cell: 1960 < cell < 1991}))
    GDD5_6190 = c.collapsed('year', iris.analysis.MEAN)
    iris.save(dayVegStart5_6190, '{}{}{}'.format(outdir, 'dayVegStart5_6190', outroot), fill_value=MISSVAL)  # noqa
    iris.save(dayVegEnd5_6190, '{}{}{}'.format(outdir, 'dayVegEnd5_6190', outroot), fill_value=MISSVAL)  # noqa
    iris.save(dayVegLen5_6190, '{}{}{}'.format(outdir, 'dayVegLen5_6190', outroot), fill_value=MISSVAL)  # noqa
    iris.save(GDD5_6190, '{}{}{}'.format(outdir, 'GDD5_6190', outroot), fill_value=MISSVAL)  # noqa


if __name__ == "__main__":
    main()