Changeset ae88f7d in flex_extract.git


Ignore:
Timestamp:
Oct 8, 2018, 3:34:53 PM (8 months ago)
Author:
Anne Philipp <anne.philipp@…>
Branches:
dev
Children:
0e576fc
Parents:
5bad6ec
Message:

added function IA3

File:
1 edited

Legend:

Unmodified
Added
Removed
  • source/python/mods/disaggregation.py

    r25b14be rae88f7d  
    4040#    - dapoly
    4141#    - darain
     42#    - IA3
    4243#
    4344#*******************************************************************************
     
    140141
    141142    return nfield
     143
     144def IA3(g):
     145    """
     146
     147    ***************************************************************************
     148    * Copyright 2017                                                          *
     149    * Sabine Hittmeir, Anne Philipp, Petra Seibert                            *
     150    *                                                                         *
     151    * This work is licensed under the Creative Commons Attribution 4.0        *
     152    * International License. To view a copy of this license, visit            *
     153    * http://creativecommons.org/licenses/by/4.0/ or send a letter to         *
     154    * Creative Commons, PO Box 1866, Mountain View, CA 94042, USA.            *
     155    ***************************************************************************
     156
     157    @Description
     158       The given data series will be interpolated with a non-negative geometric
     159       mean based algorithm. The original grid is reconstructed by adding two
     160       sampling points in each data series interval. This subgrid is used to
     161       keep all information during the interpolation within the associated
     162       interval. Additionally, an advanced monotonicity filter is applied to
     163       improve the monotonicity properties of the series.
     164
     165       For more information see article:
     166       Hittmeir, S.; Philipp, A.; Seibert, P. (2017): A conservative
     167       interpolation scheme for extensive quantities with application to the
     168       Lagrangian particle dispersion model FLEXPART.,
     169       Geoscientific Model Development
     170
     171    @Input
     172       g: list of float values
     173          A list of float values which represents the complete data series that
     174          will be interpolated having the dimension of the original raw series.
     175
     176    @Return
     177       f: list of float values
     178          The interpolated data series with additional subgrid points.
     179          Its dimension is equal to the length of the input data series
     180          times three.
     181    """
     182
     183    #######################  variable description #############################
     184    #                                                                         #
     185    # i      - index variable for looping over the data series                #
     186    # g      - input data series                                              #
     187    # f      - interpolated and filtered data series with additional          #
     188    #          grid points                                                    #
     189    # fi     - function value at position i, f_i                              #
     190    # fi1    - first  sub-grid function value f_i^1                           #
     191    # fi2    - second sub-grid function value f_i^2                           #
     192    # fip1   - next function value at position i+1, f_(i+1)                   #
     193    # dt     - time step                                                      #
     194    # fmon   - monotonicity filter                                            #
     195    #                                                                         #
     196    ###########################################################################
     197
     198
     199    import numpy as np
     200
     201    # time step
     202    dt=1.0
     203
     204    ############### Non-negative Geometric Mean Based Algorithm ###############
     205
     206    # for the left boundary the following boundary condition is valid:
     207    # the value at t=0 of the interpolation algorithm coincides with the
     208    # first data value according to the persistence hypothesis
     209    f=[g[0]]
     210
     211    # compute two first sub-grid intervals without monotonicity check
     212    # go through the data series and extend each interval by two sub-grid
     213    # points and interpolate the corresponding data values
     214    # except for the last interval due to boundary conditions
     215    for i in range(0,2):
     216
     217        # as a requirement:
     218        # if there is a zero data value such that g[i]=0, then the whole
     219        # interval in f has to be zero to such that f[i+1]=f[i+2]=f[i+3]=0
     220        # according to Eq. (6)
     221        if g[i]==0.:
     222            f.extend([0.,0.,0.])
     223
     224        # otherwise the sub-grid values are calculated and added to the list
     225        else:
     226            # temporal save of last value in interpolated list
     227            # since it is the left boundary and hence the new (fi) value
     228            fi = f[-1]
     229
     230            # the value at the end of the interval (fip1) is prescribed by the
     231            # geometric mean, restricted such that non-negativity is guaranteed
     232            # according to Eq. (25)
     233            fip1=min( 3.*g[i] , 3.*g[i+1] , np.sqrt(g[i+1]*g[i]) )
     234
     235            # the function value at the first sub-grid point (fi1) is determined
     236            # according to the equal area condition with Eq. (19)
     237            fi1=3./2.*g[i]-5./12.*fip1-1./12.*fi
     238
     239            # the function value at the second sub-grid point (fi2) is determined
     240            # according Eq. (18)
     241            fi2=fi1+1./3.*(fip1-fi)
     242
     243            # add next interval of interpolated (sub-)grid values
     244            f.append(fi1)
     245            f.append(fi2)
     246            f.append(fip1)
     247
     248    # compute rest of the data series intervals
     249    # go through the data series and extend each interval by two sub-grid
     250    # points and interpolate the corresponding data values
     251    # except for the last interval due to boundary conditions
     252    for i in range(2,len(g)-1):
     253
     254        # as a requirement:
     255        # if there is a zero data value such that g[i]=0, then the whole
     256        # interval in f has to be zero to such that f[i+1]=f[i+2]=f[i+3]=0
     257        # according to Eq. (6)
     258        if g[i]==0.:
     259            # apply monotonicity filter for interval before
     260            # check if there is "M" or "W" shape
     261            if     np.sign(f[-5]-f[-6]) * np.sign(f[-4]-f[-5])==-1 \
     262               and np.sign(f[-4]-f[-5]) * np.sign(f[-3]-f[-4])==-1 \
     263               and np.sign(f[-3]-f[-4]) * np.sign(f[-2]-f[-3])==-1:
     264
     265                # the monotonicity filter corrects the value at (fim1) by
     266                # substituting (fim1) with (fmon), see Eq. (27), (28) and (29)
     267                fmon = min(3.*g[i-2], \
     268                           3.*g[i-1], \
     269                           np.sqrt(max(0,(18./13.*g[i-2] - 5./13.*f[-7]) *
     270                                         (18./13.*g[i-1] - 5./13.*f[-1]))))
     271
     272                # recomputation of the sub-grid interval values while the
     273                # interval boundaries (fi) and (fip2) remains unchanged
     274                # see Eq. (18) and (19)
     275                f[-4]=fmon
     276                f[-6]=3./2.*g[i-2]-5./12.*fmon-1./12.*f[-7]
     277                f[-5]=f[-6]+(fmon-f[-7])/3.
     278                f[-3]=3./2.*g[i-1]-5./12.*f[-1]-1./12.*fmon
     279                f[-2]=f[-3]+(f[-1]-fmon)/3.
     280
     281            f.extend([0.,0.,0.])
     282
     283        # otherwise the sub-grid values are calculated and added to the list
     284        else:
     285            # temporal save of last value in interpolated list
     286            # since it is the left boundary and hence the new (fi) value
     287            fi = f[-1]
     288
     289            # the value at the end of the interval (fip1) is prescribed by the
     290            # geometric mean, restricted such that non-negativity is guaranteed
     291            # according to Eq. (25)
     292            fip1=min( 3.*g[i] , 3.*g[i+1] , np.sqrt(g[i+1]*g[i]) )
     293
     294            # the function value at the first sub-grid point (fi1) is determined
     295            # according to the equal area condition with Eq. (19)
     296            fi1=3./2.*g[i]-5./12.*fip1-1./12.*fi
     297
     298            # the function value at the second sub-grid point (fi2) is determined
     299            # according Eq. (18)
     300            fi2=fi1+1./3.*(fip1-fi)
     301
     302            # apply monotonicity filter for interval before
     303            # check if there is "M" or "W" shape
     304            if     np.sign(f[-5]-f[-6]) * np.sign(f[-4]-f[-5])==-1 \
     305               and np.sign(f[-4]-f[-5]) * np.sign(f[-3]-f[-4])==-1 \
     306               and np.sign(f[-3]-f[-4]) * np.sign(f[-2]-f[-3])==-1:
     307
     308                # the monotonicity filter corrects the value at (fim1) by
     309                # substituting (fim1) with fmon, see Eq. (27), (28) and (29)
     310                fmon = min(3.*g[i-2], \
     311                           3.*g[i-1], \
     312                           np.sqrt(max(0,(18./13.*g[i-2] - 5./13.*f[-7]) *
     313                                         (18./13.*g[i-1] - 5./13.*f[-1]))))
     314
     315                # recomputation of the sub-grid interval values while the
     316                # interval boundaries (fi) and (fip2) remains unchanged
     317                # see Eq. (18) and (19)
     318                f[-4]=fmon
     319                f[-6]=3./2.*g[i-2]-5./12.*fmon-1./12.*f[-7]
     320                f[-5]=f[-6]+(fmon-f[-7])/3.
     321                f[-3]=3./2.*g[i-1]-5./12.*f[-1]-1./12.*fmon
     322                f[-2]=f[-3]+(f[-1]-fmon)/3.
     323
     324            # add next interval of interpolated (sub-)grid values
     325            f.append(fi1)
     326            f.append(fi2)
     327            f.append(fip1)
     328
     329    # separate treatment of the final interval
     330
     331    # as a requirement:
     332    # if there is a zero data value such that g[i]=0, then the whole
     333    # interval in f has to be zero to such that f[i+1]=f[i+2]=f[i+3]=0
     334    # according to Eq. (6)
     335    if g[-1]==0.:
     336        # apply monotonicity filter for interval before
     337        # check if there is "M" or "W" shape
     338        if     np.sign(f[-5]-f[-6]) * np.sign(f[-4]-f[-5])==-1 \
     339           and np.sign(f[-4]-f[-5]) * np.sign(f[-3]-f[-4])==-1 \
     340           and np.sign(f[-3]-f[-4]) * np.sign(f[-2]-f[-3])==-1:
     341
     342            # the monotonicity filter corrects the value at (fim1) by
     343            # substituting (fim1) with (fmon), see Eq. (27), (28) and (29)
     344            fmon = min(3.*g[-3], \
     345                       3.*g[-2], \
     346                       np.sqrt(max(0,(18./13.*g[-3] - 5./13.*f[-7]) *
     347                                     (18./13.*g[-2] - 5./13.*f[-1]))))
     348
     349            # recomputation of the sub-grid interval values while the
     350            # interval boundaries (fi) and (fip2) remains unchanged
     351            # see Eq. (18) and (19)
     352            f[-4]=fmon
     353            f[-6]=3./2.*g[-3]-5./12.*fmon-1./12.*f[-7]
     354            f[-5]=f[-6]+(fmon-f[-7])/3.
     355            f[-3]=3./2.*g[-2]-5./12.*f[-1]-1./12.*fmon
     356            f[-2]=f[-3]+(f[-1]-fmon)/3.
     357
     358        f.extend([0.,0.,0.])
     359
     360    # otherwise the sub-grid values are calculated and added to the list
     361    # using the persistence hypothesis as boundary condition
     362    else:
     363        # temporal save of last value in interpolated list
     364        # since it is the left boundary and hence the new (fi) value
     365        fi = f[-1]
     366        # since last interval in series, last value is also fip1
     367        fip1 = g[-1]
     368        # the function value at the first sub-grid point (fi1) is determined
     369        # according to the equal area condition with Eq. (19)
     370        fi1 = 3./2.*g[-1]-5./12.*fip1-1./12.*fi
     371        # the function value at the second sub-grid point (fi2) is determined
     372        # according Eq. (18)
     373        fi2 = fi1+dt/3.*(fip1-fi)
     374
     375        # apply monotonicity filter for interval before
     376        # check if there is "M" or "W" shape
     377        if     np.sign(f[-5]-f[-6]) * np.sign(f[-4]-f[-5])==-1 \
     378           and np.sign(f[-4]-f[-5]) * np.sign(f[-3]-f[-4])==-1 \
     379           and np.sign(f[-3]-f[-4]) * np.sign(f[-2]-f[-3])==-1:
     380
     381            # the monotonicity filter corrects the value at (fim1) by
     382            # substituting (fim1) with (fmon), see Eq. (27), (28) and (29)
     383            fmon = min(3.*g[-3], \
     384                       3.*g[-2], \
     385                       np.sqrt(max(0,(18./13.*g[-3] - 5./13.*f[-7]) *
     386                                     (18./13.*g[-2] - 5./13.*f[-1]))))
     387
     388            # recomputation of the sub-grid interval values while the
     389            # interval boundaries (fi) and (fip2) remains unchanged
     390            # see Eq. (18) and (19)
     391            f[-4]=fmon
     392            f[-6]=3./2.*g[-3]-5./12.*fmon-1./12.*f[-7]
     393            f[-5]=f[-6]+(fmon-f[-7])/3.
     394            f[-3]=3./2.*g[-2]-5./12.*f[-1]-1./12.*fmon
     395            f[-2]=f[-3]+(f[-1]-fmon)/3.
     396
     397        # add next interval of interpolated (sub-)grid values
     398        f.append(fi1)
     399        f.append(fi2)
     400        f.append(fip1)
     401
     402    return f
Note: See TracChangeset for help on using the changeset viewer.
hosted by ZAMG