source: flex_extract.git/source/python/mods/disaggregation.py @ 79729d5

ctbtodev
Last change on this file since 79729d5 was 79729d5, checked in by Anne Philipp <anne.philipp@…>, 5 years ago

switched from python2 to python3

  • Property mode set to 100644
File size: 16.1 KB
RevLine 
[79729d5]1#!/usr/bin/env python3
[efdb01a]2# -*- coding: utf-8 -*-
[991df6a]3#*******************************************************************************
4# @Author: Anne Philipp (University of Vienna)
5#
6# @Date: March 2018
7#
8# @Change History:
[54a8a01]9#
[991df6a]10#    November 2015 - Leopold Haimberger (University of Vienna):
11#        - migration of the methods dapoly and darain from Fortran
12#          (flex_extract_v6 and earlier) to Python
13#
14#    April 2018 - Anne Philipp (University of Vienna):
15#        - applied PEP8 style guide
16#        - added structured documentation
17#        - outsourced the disaggregation functions dapoly and darain
[ff99eae]18#          to a new module named disaggregation
[6f951ca]19#        - added the new disaggregation method for precipitation
[991df6a]20#
21# @License:
[6f951ca]22#    (C) Copyright 2014-2019.
23#    Anne Philipp, Leopold Haimberger
[991df6a]24#
[6f951ca]25#    This work is licensed under the Creative Commons Attribution 4.0
26#    International License. To view a copy of this license, visit
27#    http://creativecommons.org/licenses/by/4.0/ or send a letter to
28#    Creative Commons, PO Box 1866, Mountain View, CA 94042, USA.
[991df6a]29#
[6f951ca]30# @Methods:
[991df6a]31#    - dapoly
32#    - darain
[ae88f7d]33#    - IA3
[991df6a]34#*******************************************************************************
[6f951ca]35'''Disaggregation of deaccumulated flux data from an ECMWF model FG field.
36
37Initially the flux data to be concerned are:
38    - large-scale precipitation
39    - convective precipitation
40    - surface sensible heat flux
41    - surface solar radiation
42    - u stress
43    - v stress
44
45Different versions of disaggregation is provided for rainfall
46data (darain, modified linear) and the surface fluxes and
47stress data (dapoly, cubic polynomial).
48'''
[efdb01a]49
50# ------------------------------------------------------------------------------
51# MODULES
52# ------------------------------------------------------------------------------
53
54# ------------------------------------------------------------------------------
55# FUNCTIONS
56# ------------------------------------------------------------------------------
57def dapoly(alist):
[708c667]58    """Cubic polynomial interpolation of deaccumulated fluxes.
59
60    Interpolation of deaccumulated fluxes of an ECMWF model FG field
61    using a cubic polynomial solution which conserves the integrals
62    of the fluxes within each timespan.
63    Disaggregation is done for 4 accumluated timespans which
64    generates a new, disaggregated value which is output at the
65    central point of the 4 accumulation timespans.
66    This new point is used for linear interpolation of the complete
67    timeseries afterwards.
68
69    Parameters
70    ----------
[6f951ca]71    alist : list of array of float
[708c667]72        List of 4 timespans as 2-dimensional, horizontal fields.
73        E.g. [[array_t1], [array_t2], [array_t3], [array_t4]]
74
75    Return
76    ------
[6f951ca]77    nfield : array of float
[708c667]78        Interpolated flux at central point of accumulation timespan.
79
80    Note
81    ----
82    March 2000    : P. JAMES
83        Original author
84
85    June 2003     : A. BECK
86        Adaptations
87
88    November 2015 : Leopold Haimberger (University of Vienna)
89        Migration from Fortran to Python
90
91    """
92
[efdb01a]93    pya = (alist[3] - alist[0] + 3. * (alist[1] - alist[2])) / 6.
94    pyb = (alist[2] + alist[0]) / 2. - alist[1] - 9. * pya / 2.
95    pyc = alist[1] - alist[0] - 7. * pya / 2. - 2. * pyb
96    pyd = alist[0] - pya / 4. - pyb / 3. - pyc / 2.
97    nfield = 8. * pya + 4. * pyb + 2. * pyc + pyd
98
99    return nfield
100
101
102def darain(alist):
[708c667]103    """Linear interpolation of deaccumulated fluxes.
104
105    Interpolation of deaccumulated fluxes of an ECMWF model FG rainfall
106    field using a modified linear solution which conserves the integrals
107    of the fluxes within each timespan.
108    Disaggregation is done for 4 accumluated timespans which generates
109    a new, disaggregated value which is output at the central point
110    of the 4 accumulation timespans. This new point is used for linear
111    interpolation of the complete timeseries afterwards.
112
113    Parameters
114    ----------
[6f951ca]115    alist : list of array of float
[708c667]116        List of 4 timespans as 2-dimensional, horizontal fields.
117        E.g. [[array_t1], [array_t2], [array_t3], [array_t4]]
118
119    Return
120    ------
[6f951ca]121    nfield : array of float
[708c667]122        Interpolated flux at central point of accumulation timespan.
123
124    Note
125    ----
126    March 2000    : P. JAMES
127        Original author
128
129    June 2003     : A. BECK
130        Adaptations
131
132    November 2015 : Leopold Haimberger (University of Vienna)
133        Migration from Fortran to Python
134    """
135
[efdb01a]136    xa = alist[0]
137    xb = alist[1]
138    xc = alist[2]
139    xd = alist[3]
140    xa[xa < 0.] = 0.
141    xb[xb < 0.] = 0.
142    xc[xc < 0.] = 0.
143    xd[xd < 0.] = 0.
144
145    xac = 0.5 * xb
146    mask = xa + xc > 0.
147    xac[mask] = xb[mask] * xc[mask] / (xa[mask] + xc[mask])
148    xbd = 0.5 * xc
149    mask = xb + xd > 0.
150    xbd[mask] = xb[mask] * xc[mask] / (xb[mask] + xd[mask])
151    nfield = xac + xbd
152
153    return nfield
[ae88f7d]154
155def IA3(g):
[708c667]156    """ Interpolation with a non-negative geometric mean based algorithm.
157
158    The original grid is reconstructed by adding two sampling points in each
159    data series interval. This subgrid is used to keep all information during
160    the interpolation within the associated interval. Additionally, an advanced
161    monotonicity filter is applied to improve the monotonicity properties of
162    the series.
163
164    Note
165    ----
[6f951ca]166    (C) Copyright 2017-2019
[708c667]167    Sabine Hittmeir, Anne Philipp, Petra Seibert
168
169    This work is licensed under the Creative Commons Attribution 4.0
170    International License. To view a copy of this license, visit
171    http://creativecommons.org/licenses/by/4.0/ or send a letter to
172    Creative Commons, PO Box 1866, Mountain View, CA 94042, USA.
173
174    Parameters
175    ----------
[6f951ca]176    g : list of float
[708c667]177        Complete data series that will be interpolated having
178        the dimension of the original raw series.
179
180    Return
181    ------
[6f951ca]182    f : list of float
[708c667]183        The interpolated data series with additional subgrid points.
184        Its dimension is equal to the length of the input data series
185        times three.
186
187
188    References
189    ----------
190    For more information see article:
191    Hittmeir, S.; Philipp, A.; Seibert, P. (2017): A conservative
192    interpolation scheme for extensive quantities with application to the
193    Lagrangian particle dispersion model FLEXPART.,
194    Geoscientific Model Development
[ae88f7d]195    """
196
197    #######################  variable description #############################
198    #                                                                         #
199    # i      - index variable for looping over the data series                #
200    # g      - input data series                                              #
201    # f      - interpolated and filtered data series with additional          #
202    #          grid points                                                    #
203    # fi     - function value at position i, f_i                              #
204    # fi1    - first  sub-grid function value f_i^1                           #
205    # fi2    - second sub-grid function value f_i^2                           #
206    # fip1   - next function value at position i+1, f_(i+1)                   #
207    # dt     - time step                                                      #
208    # fmon   - monotonicity filter                                            #
209    #                                                                         #
210    ###########################################################################
211
212
213    import numpy as np
214
215    # time step
[0e2f93e]216    dt = 1.0
[ae88f7d]217
218    ############### Non-negative Geometric Mean Based Algorithm ###############
219
220    # for the left boundary the following boundary condition is valid:
221    # the value at t=0 of the interpolation algorithm coincides with the
222    # first data value according to the persistence hypothesis
[0e2f93e]223    f = [g[0]]
[ae88f7d]224
225    # compute two first sub-grid intervals without monotonicity check
226    # go through the data series and extend each interval by two sub-grid
227    # points and interpolate the corresponding data values
228    # except for the last interval due to boundary conditions
[0e2f93e]229    for i in range(0, 2):
[ae88f7d]230
231        # as a requirement:
232        # if there is a zero data value such that g[i]=0, then the whole
233        # interval in f has to be zero to such that f[i+1]=f[i+2]=f[i+3]=0
234        # according to Eq. (6)
[0e2f93e]235        if g[i] == 0.:
236            f.extend([0., 0., 0.])
[ae88f7d]237
238        # otherwise the sub-grid values are calculated and added to the list
239        else:
240            # temporal save of last value in interpolated list
241            # since it is the left boundary and hence the new (fi) value
242            fi = f[-1]
243
244            # the value at the end of the interval (fip1) is prescribed by the
245            # geometric mean, restricted such that non-negativity is guaranteed
246            # according to Eq. (25)
247            fip1=min( 3.*g[i] , 3.*g[i+1] , np.sqrt(g[i+1]*g[i]) )
248
249            # the function value at the first sub-grid point (fi1) is determined
250            # according to the equal area condition with Eq. (19)
[0e2f93e]251            fi1 = 3./2.*g[i]-5./12.*fip1-1./12.*fi
[ae88f7d]252
253            # the function value at the second sub-grid point (fi2) is determined
254            # according Eq. (18)
[0e2f93e]255            fi2 = fi1+1./3.*(fip1-fi)
[ae88f7d]256
257            # add next interval of interpolated (sub-)grid values
258            f.append(fi1)
259            f.append(fi2)
260            f.append(fip1)
261
262    # compute rest of the data series intervals
263    # go through the data series and extend each interval by two sub-grid
264    # points and interpolate the corresponding data values
265    # except for the last interval due to boundary conditions
[0e2f93e]266    for i in range(2, len(g)-1):
[ae88f7d]267
268        # as a requirement:
269        # if there is a zero data value such that g[i]=0, then the whole
270        # interval in f has to be zero to such that f[i+1]=f[i+2]=f[i+3]=0
271        # according to Eq. (6)
[0e2f93e]272        if g[i] == 0.:
[ae88f7d]273            # apply monotonicity filter for interval before
274            # check if there is "M" or "W" shape
[0e2f93e]275            if     np.sign(f[-5]-f[-6]) * np.sign(f[-4]-f[-5]) == -1 \
276               and np.sign(f[-4]-f[-5]) * np.sign(f[-3]-f[-4]) == -1 \
277               and np.sign(f[-3]-f[-4]) * np.sign(f[-2]-f[-3]) == -1:
[ae88f7d]278
279                # the monotonicity filter corrects the value at (fim1) by
280                # substituting (fim1) with (fmon), see Eq. (27), (28) and (29)
281                fmon = min(3.*g[i-2], \
282                           3.*g[i-1], \
283                           np.sqrt(max(0,(18./13.*g[i-2] - 5./13.*f[-7]) *
284                                         (18./13.*g[i-1] - 5./13.*f[-1]))))
285
286                # recomputation of the sub-grid interval values while the
287                # interval boundaries (fi) and (fip2) remains unchanged
288                # see Eq. (18) and (19)
[0e2f93e]289                f[-4] = fmon
290                f[-6] = 3./2.*g[i-2]-5./12.*fmon-1./12.*f[-7]
291                f[-5] = f[-6]+(fmon-f[-7])/3.
292                f[-3] = 3./2.*g[i-1]-5./12.*f[-1]-1./12.*fmon
293                f[-2] = f[-3]+(f[-1]-fmon)/3.
[ae88f7d]294
295            f.extend([0.,0.,0.])
296
297        # otherwise the sub-grid values are calculated and added to the list
298        else:
299            # temporal save of last value in interpolated list
300            # since it is the left boundary and hence the new (fi) value
301            fi = f[-1]
302
303            # the value at the end of the interval (fip1) is prescribed by the
304            # geometric mean, restricted such that non-negativity is guaranteed
305            # according to Eq. (25)
[0e2f93e]306            fip1 = min( 3.*g[i] , 3.*g[i+1] , np.sqrt(g[i+1]*g[i]) )
[ae88f7d]307
308            # the function value at the first sub-grid point (fi1) is determined
309            # according to the equal area condition with Eq. (19)
[0e2f93e]310            fi1 = 3./2.*g[i]-5./12.*fip1-1./12.*fi
[ae88f7d]311
312            # the function value at the second sub-grid point (fi2) is determined
313            # according Eq. (18)
[0e2f93e]314            fi2 = fi1+1./3.*(fip1-fi)
[ae88f7d]315
316            # apply monotonicity filter for interval before
317            # check if there is "M" or "W" shape
[0e2f93e]318            if     np.sign(f[-5]-f[-6]) * np.sign(f[-4]-f[-5]) == -1 \
319               and np.sign(f[-4]-f[-5]) * np.sign(f[-3]-f[-4]) == -1 \
320               and np.sign(f[-3]-f[-4]) * np.sign(f[-2]-f[-3]) == -1:
[ae88f7d]321
322                # the monotonicity filter corrects the value at (fim1) by
323                # substituting (fim1) with fmon, see Eq. (27), (28) and (29)
324                fmon = min(3.*g[i-2], \
325                           3.*g[i-1], \
326                           np.sqrt(max(0,(18./13.*g[i-2] - 5./13.*f[-7]) *
327                                         (18./13.*g[i-1] - 5./13.*f[-1]))))
328
329                # recomputation of the sub-grid interval values while the
330                # interval boundaries (fi) and (fip2) remains unchanged
331                # see Eq. (18) and (19)
[0e2f93e]332                f[-4] = fmon
333                f[-6] = 3./2.*g[i-2]-5./12.*fmon-1./12.*f[-7]
334                f[-5] = f[-6]+(fmon-f[-7])/3.
335                f[-3] = 3./2.*g[i-1]-5./12.*f[-1]-1./12.*fmon
336                f[-2] = f[-3]+(f[-1]-fmon)/3.
[ae88f7d]337
338            # add next interval of interpolated (sub-)grid values
339            f.append(fi1)
340            f.append(fi2)
341            f.append(fip1)
342
343    # separate treatment of the final interval
344
345    # as a requirement:
346    # if there is a zero data value such that g[i]=0, then the whole
347    # interval in f has to be zero to such that f[i+1]=f[i+2]=f[i+3]=0
348    # according to Eq. (6)
[0e2f93e]349    if g[-1] == 0.:
[ae88f7d]350        # apply monotonicity filter for interval before
351        # check if there is "M" or "W" shape
[0e2f93e]352        if     np.sign(f[-5]-f[-6]) * np.sign(f[-4]-f[-5]) == -1 \
353           and np.sign(f[-4]-f[-5]) * np.sign(f[-3]-f[-4]) == -1 \
354           and np.sign(f[-3]-f[-4]) * np.sign(f[-2]-f[-3]) == -1:
[ae88f7d]355
356            # the monotonicity filter corrects the value at (fim1) by
357            # substituting (fim1) with (fmon), see Eq. (27), (28) and (29)
358            fmon = min(3.*g[-3], \
359                       3.*g[-2], \
360                       np.sqrt(max(0,(18./13.*g[-3] - 5./13.*f[-7]) *
361                                     (18./13.*g[-2] - 5./13.*f[-1]))))
362
363            # recomputation of the sub-grid interval values while the
364            # interval boundaries (fi) and (fip2) remains unchanged
365            # see Eq. (18) and (19)
[0e2f93e]366            f[-4] = fmon
367            f[-6] = 3./2.*g[-3]-5./12.*fmon-1./12.*f[-7]
368            f[-5] = f[-6]+(fmon-f[-7])/3.
369            f[-3] = 3./2.*g[-2]-5./12.*f[-1]-1./12.*fmon
370            f[-2] = f[-3]+(f[-1]-fmon)/3.
[ae88f7d]371
372        f.extend([0.,0.,0.])
373
374    # otherwise the sub-grid values are calculated and added to the list
375    # using the persistence hypothesis as boundary condition
376    else:
377        # temporal save of last value in interpolated list
378        # since it is the left boundary and hence the new (fi) value
379        fi = f[-1]
380        # since last interval in series, last value is also fip1
381        fip1 = g[-1]
382        # the function value at the first sub-grid point (fi1) is determined
383        # according to the equal area condition with Eq. (19)
384        fi1 = 3./2.*g[-1]-5./12.*fip1-1./12.*fi
385        # the function value at the second sub-grid point (fi2) is determined
386        # according Eq. (18)
387        fi2 = fi1+dt/3.*(fip1-fi)
388
389        # apply monotonicity filter for interval before
390        # check if there is "M" or "W" shape
[0e2f93e]391        if     np.sign(f[-5]-f[-6]) * np.sign(f[-4]-f[-5]) == -1 \
392           and np.sign(f[-4]-f[-5]) * np.sign(f[-3]-f[-4]) == -1 \
393           and np.sign(f[-3]-f[-4]) * np.sign(f[-2]-f[-3]) == -1:
[ae88f7d]394
395            # the monotonicity filter corrects the value at (fim1) by
396            # substituting (fim1) with (fmon), see Eq. (27), (28) and (29)
397            fmon = min(3.*g[-3], \
398                       3.*g[-2], \
399                       np.sqrt(max(0,(18./13.*g[-3] - 5./13.*f[-7]) *
400                                     (18./13.*g[-2] - 5./13.*f[-1]))))
401
402            # recomputation of the sub-grid interval values while the
403            # interval boundaries (fi) and (fip2) remains unchanged
404            # see Eq. (18) and (19)
[0e2f93e]405            f[-4] = fmon
406            f[-6] = 3./2.*g[-3]-5./12.*fmon-1./12.*f[-7]
407            f[-5] = f[-6]+(fmon-f[-7])/3.
408            f[-3] = 3./2.*g[-2]-5./12.*f[-1]-1./12.*fmon
409            f[-2] = f[-3]+(f[-1]-fmon)/3.
[ae88f7d]410
411        # add next interval of interpolated (sub-)grid values
412        f.append(fi1)
413        f.append(fi2)
414        f.append(fip1)
415
416    return f
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG