Source code for disaggregation

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
#*******************************************************************************
# @Author: Anne Philipp (University of Vienna)
#
# @Date: March 2018
#
# @Change History:
#
#    November 2015 - Leopold Haimberger (University of Vienna):
#        - migration of the methods dapoly and darain from Fortran
#          (flex_extract_v6 and earlier) to Python
#
#    April 2018 - Anne Philipp (University of Vienna):
#        - applied PEP8 style guide
#        - added structured documentation
#        - outsourced the disaggregation functions dapoly and darain
#          to a new module named disaggregation
#        - added the new disaggregation method for precipitation
#
#    June 2020 - Anne Philipp (University of Vienna):
#        - reformulated formular for dapoly
#
# @License:
#    (C) Copyright 2014-2020.
#    Anne Philipp, Leopold Haimberger
#
#    SPDX-License-Identifier: CC-BY-4.0
#
#    This work is licensed under the Creative Commons Attribution 4.0
#    International License. To view a copy of this license, visit
#    http://creativecommons.org/licenses/by/4.0/ or send a letter to
#    Creative Commons, PO Box 1866, Mountain View, CA 94042, USA.
#
# @Methods:
#    - dapoly
#    - darain
#    - IA3
#*******************************************************************************
'''Disaggregation of deaccumulated flux data from an ECMWF model FG field.

Initially the flux data to be concerned are:
    - large-scale precipitation
    - convective precipitation
    - surface sensible heat flux
    - surface solar radiation
    - u stress
    - v stress

Different versions of disaggregation is provided for rainfall
data (darain, modified linear) and the surface fluxes and
stress data (dapoly, cubic polynomial).
'''

# ------------------------------------------------------------------------------
# MODULES
# ------------------------------------------------------------------------------

# ------------------------------------------------------------------------------
# FUNCTIONS
# ------------------------------------------------------------------------------
[docs]def dapoly(alist): """Cubic polynomial interpolation of deaccumulated fluxes. Interpolation of deaccumulated fluxes of an ECMWF model FG field using a cubic polynomial solution which conserves the integrals of the fluxes within each timespan. Disaggregation is done for 4 accumluated timespans which generates a new, disaggregated value which is output at the central point of the 4 accumulation timespans. This new point is used for linear interpolation of the complete timeseries afterwards. Parameters ---------- alist : list of array of float List of 4 timespans as 2-dimensional, horizontal fields. E.g. [[array_t1], [array_t2], [array_t3], [array_t4]] Return ------ nfield : array of float Interpolated flux at central point of accumulation timespan. Note ---- March 2000 : P. JAMES Original author June 2003 : A. BECK Adaptations November 2015 : Leopold Haimberger (University of Vienna) Migration from Fortran to Python """ nfield = -1./12.*alist[0] + \ 7./12.*alist[1] + \ 7./12.*alist[2] - \ 1./12.*alist[3] return nfield
[docs]def darain(alist): """Linear interpolation of deaccumulated fluxes. Interpolation of deaccumulated fluxes of an ECMWF model FG rainfall field using a modified linear solution which conserves the integrals of the fluxes within each timespan. Disaggregation is done for 4 accumluated timespans which generates a new, disaggregated value which is output at the central point of the 4 accumulation timespans. This new point is used for linear interpolation of the complete timeseries afterwards. Parameters ---------- alist : list of array of float List of 4 timespans as 2-dimensional, horizontal fields. E.g. [[array_t1], [array_t2], [array_t3], [array_t4]] Return ------ nfield : array of float Interpolated flux at central point of accumulation timespan. Note ---- March 2000 : P. JAMES Original author June 2003 : A. BECK Adaptations November 2015 : Leopold Haimberger (University of Vienna) Migration from Fortran to Python """ xa = alist[0] xb = alist[1] xc = alist[2] xd = alist[3] xa[xa < 0.] = 0. xb[xb < 0.] = 0. xc[xc < 0.] = 0. xd[xd < 0.] = 0. xac = 0.5 * xb mask = xa + xc > 0. xac[mask] = xb[mask] * xc[mask] / (xa[mask] + xc[mask]) xbd = 0.5 * xc mask = xb + xd > 0. xbd[mask] = xb[mask] * xc[mask] / (xb[mask] + xd[mask]) nfield = xac + xbd return nfield
[docs]def IA3(g): """ Interpolation with a non-negative geometric mean based algorithm. The original grid is reconstructed by adding two sampling points in each data series interval. This subgrid is used to keep all information during the interpolation within the associated interval. Additionally, an advanced monotonicity filter is applied to improve the monotonicity properties of the series. Note ---- (C) Copyright 2017-2019 Sabine Hittmeir, Anne Philipp, Petra Seibert This work is licensed under the Creative Commons Attribution 4.0 International License. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/ or send a letter to Creative Commons, PO Box 1866, Mountain View, CA 94042, USA. Parameters ---------- g : list of float Complete data series that will be interpolated having the dimension of the original raw series. Return ------ f : list of float The interpolated data series with additional subgrid points. Its dimension is equal to the length of the input data series times three. References ---------- For more information see article: Hittmeir, S.; Philipp, A.; Seibert, P. (2017): A conservative interpolation scheme for extensive quantities with application to the Lagrangian particle dispersion model FLEXPART., Geoscientific Model Development """ ####################### variable description ############################# # # # i - index variable for looping over the data series # # g - input data series # # f - interpolated and filtered data series with additional # # grid points # # fi - function value at position i, f_i # # fi1 - first sub-grid function value f_i^1 # # fi2 - second sub-grid function value f_i^2 # # fip1 - next function value at position i+1, f_(i+1) # # dt - time step # # fmon - monotonicity filter # # # ########################################################################### import numpy as np # time step dt = 1.0 ############### Non-negative Geometric Mean Based Algorithm ############### # for the left boundary the following boundary condition is valid: # the value at t=0 of the interpolation algorithm coincides with the # first data value according to the persistence hypothesis f = [g[0]] # compute two first sub-grid intervals without monotonicity check # go through the data series and extend each interval by two sub-grid # points and interpolate the corresponding data values # except for the last interval due to boundary conditions for i in range(0, 2): # as a requirement: # if there is a zero data value such that g[i]=0, then the whole # interval in f has to be zero to such that f[i+1]=f[i+2]=f[i+3]=0 # according to Eq. (6) if g[i] == 0.: f.extend([0., 0., 0.]) # otherwise the sub-grid values are calculated and added to the list else: # temporal save of last value in interpolated list # since it is the left boundary and hence the new (fi) value fi = f[-1] # the value at the end of the interval (fip1) is prescribed by the # geometric mean, restricted such that non-negativity is guaranteed # according to Eq. (25) fip1 = min(3. * g[i], 3. * g[i + 1], np.sqrt(g[i + 1] * g[i])) # the function value at the first sub-grid point (fi1) is determined # according to the equal area condition with Eq. (19) fi1 = 3./2.*g[i]-5./12.*fip1-1./12.*fi # the function value at the second sub-grid point (fi2) is determined # according Eq. (18) fi2 = fi1+1./3.*(fip1-fi) # add next interval of interpolated (sub-)grid values f.append(fi1) f.append(fi2) f.append(fip1) # compute rest of the data series intervals # go through the data series and extend each interval by two sub-grid # points and interpolate the corresponding data values # except for the last interval due to boundary conditions for i in range(2, len(g)-1): # as a requirement: # if there is a zero data value such that g[i]=0, then the whole # interval in f has to be zero to such that f[i+1]=f[i+2]=f[i+3]=0 # according to Eq. (6) if g[i] == 0.: # apply monotonicity filter for interval before # check if there is "M" or "W" shape if np.sign(f[-5]-f[-6]) * np.sign(f[-4]-f[-5]) == -1 \ and np.sign(f[-4]-f[-5]) * np.sign(f[-3]-f[-4]) == -1 \ and np.sign(f[-3]-f[-4]) * np.sign(f[-2]-f[-3]) == -1: # the monotonicity filter corrects the value at (fim1) by # substituting (fim1) with (fmon), see Eq. (27), (28) and (29) fmon = min(3. * g[i - 2], 3. * g[i - 1], np.sqrt(max(0, (18. / 13. * g[i - 2] - 5. / 13. * f[-7]) * (18. / 13. * g[i - 1] - 5. / 13. * f[-1])))) # recomputation of the sub-grid interval values while the # interval boundaries (fi) and (fip2) remains unchanged # see Eq. (18) and (19) f[-4] = fmon f[-6] = 3./2.*g[i-2]-5./12.*fmon-1./12.*f[-7] f[-5] = f[-6]+(fmon-f[-7])/3. f[-3] = 3./2.*g[i-1]-5./12.*f[-1]-1./12.*fmon f[-2] = f[-3]+(f[-1]-fmon)/3. f.extend([0., 0., 0.]) # otherwise the sub-grid values are calculated and added to the list else: # temporal save of last value in interpolated list # since it is the left boundary and hence the new (fi) value fi = f[-1] # the value at the end of the interval (fip1) is prescribed by the # geometric mean, restricted such that non-negativity is guaranteed # according to Eq. (25) fip1 = min(3. * g[i], 3. * g[i + 1], np.sqrt(g[i + 1] * g[i])) # the function value at the first sub-grid point (fi1) is determined # according to the equal area condition with Eq. (19) fi1 = 3./2.*g[i]-5./12.*fip1-1./12.*fi # the function value at the second sub-grid point (fi2) is determined # according Eq. (18) fi2 = fi1+1./3.*(fip1-fi) # apply monotonicity filter for interval before # check if there is "M" or "W" shape if np.sign(f[-5]-f[-6]) * np.sign(f[-4]-f[-5]) == -1 \ and np.sign(f[-4]-f[-5]) * np.sign(f[-3]-f[-4]) == -1 \ and np.sign(f[-3]-f[-4]) * np.sign(f[-2]-f[-3]) == -1: # the monotonicity filter corrects the value at (fim1) by # substituting (fim1) with fmon, see Eq. (27), (28) and (29) fmon = min(3. * g[i - 2], 3. * g[i - 1], np.sqrt(max(0, (18. / 13. * g[i - 2] - 5. / 13. * f[-7]) * (18. / 13. * g[i - 1] - 5. / 13. * f[-1])))) # recomputation of the sub-grid interval values while the # interval boundaries (fi) and (fip2) remains unchanged # see Eq. (18) and (19) f[-4] = fmon f[-6] = 3./2.*g[i-2]-5./12.*fmon-1./12.*f[-7] f[-5] = f[-6]+(fmon-f[-7])/3. f[-3] = 3./2.*g[i-1]-5./12.*f[-1]-1./12.*fmon f[-2] = f[-3]+(f[-1]-fmon)/3. # add next interval of interpolated (sub-)grid values f.append(fi1) f.append(fi2) f.append(fip1) # separate treatment of the final interval # as a requirement: # if there is a zero data value such that g[i]=0, then the whole # interval in f has to be zero to such that f[i+1]=f[i+2]=f[i+3]=0 # according to Eq. (6) if g[-1] == 0.: # apply monotonicity filter for interval before # check if there is "M" or "W" shape if np.sign(f[-5]-f[-6]) * np.sign(f[-4]-f[-5]) == -1 \ and np.sign(f[-4]-f[-5]) * np.sign(f[-3]-f[-4]) == -1 \ and np.sign(f[-3]-f[-4]) * np.sign(f[-2]-f[-3]) == -1: # the monotonicity filter corrects the value at (fim1) by # substituting (fim1) with (fmon), see Eq. (27), (28) and (29) fmon = min(3. * g[-3], 3. * g[-2], np.sqrt(max(0, (18. / 13. * g[-3] - 5. / 13. * f[-7]) * (18. / 13. * g[-2] - 5. / 13. * f[-1])))) # recomputation of the sub-grid interval values while the # interval boundaries (fi) and (fip2) remains unchanged # see Eq. (18) and (19) f[-4] = fmon f[-6] = 3./2.*g[-3]-5./12.*fmon-1./12.*f[-7] f[-5] = f[-6]+(fmon-f[-7])/3. f[-3] = 3./2.*g[-2]-5./12.*f[-1]-1./12.*fmon f[-2] = f[-3]+(f[-1]-fmon)/3. f.extend([0., 0., 0.]) # otherwise the sub-grid values are calculated and added to the list # using the persistence hypothesis as boundary condition else: # temporal save of last value in interpolated list # since it is the left boundary and hence the new (fi) value fi = f[-1] # since last interval in series, last value is also fip1 fip1 = g[-1] # the function value at the first sub-grid point (fi1) is determined # according to the equal area condition with Eq. (19) fi1 = 3./2.*g[-1]-5./12.*fip1-1./12.*fi # the function value at the second sub-grid point (fi2) is determined # according Eq. (18) fi2 = fi1+dt/3.*(fip1-fi) # apply monotonicity filter for interval before # check if there is "M" or "W" shape if np.sign(f[-5]-f[-6]) * np.sign(f[-4]-f[-5]) == -1 \ and np.sign(f[-4]-f[-5]) * np.sign(f[-3]-f[-4]) == -1 \ and np.sign(f[-3]-f[-4]) * np.sign(f[-2]-f[-3]) == -1: # the monotonicity filter corrects the value at (fim1) by # substituting (fim1) with (fmon), see Eq. (27), (28) and (29) fmon = min(3. * g[-3], 3. * g[-2], np.sqrt(max(0, (18. / 13. * g[-3] - 5. / 13. * f[-7]) * (18. / 13. * g[-2] - 5. / 13. * f[-1])))) # recomputation of the sub-grid interval values while the # interval boundaries (fi) and (fip2) remains unchanged # see Eq. (18) and (19) f[-4] = fmon f[-6] = 3./2.*g[-3]-5./12.*fmon-1./12.*f[-7] f[-5] = f[-6]+(fmon-f[-7])/3. f[-3] = 3./2.*g[-2]-5./12.*f[-1]-1./12.*fmon f[-2] = f[-3]+(f[-1]-fmon)/3. # add next interval of interpolated (sub-)grid values f.append(fi1) f.append(fi2) f.append(fip1) return f