InterpTools/fast_average.py at master · UNESCO-IHE/InterpTools · GitHub

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
# -*- coding: utf-8 -*-
"""
==========================
Fast Average Precipitation
==========================
Implemented by Juan Chacon @ UNESCO-IHE
Integrated Water Systems and Governance Department
Hydroinformatics Laboratory

Average precipitation from multiple stations

* Pre requisites
    you will need the following libraries, not coming alongside with the\
    Anaconda ditribution (recommended)

* Functions
    * Run: Calculate average precipitation from several stations, as shown in\
    Lindström et al [1997]

* Use policy
    * you should include the respective citation to the authors
    * if you find this tool usefull, you will give the main author a beer next\
    time you see him :)

* References
    * Lindström et al., “Development and Test of the Distributed HBV-96\
    Hydrological Model.”
"""
import numpy


def Run(data, covariance=None):
    '''
    Calculates spatial average from several stations. Data consists in\
    a matrix element containing all the recordas from the stations. If\
    covariance matrix can be provided it can be passed as an optional argument.

    Paramters
    ----------
        **data** -- Matrix of recordings from stations. data has to be\
        oriented column wise

        **covariance** -- Matrix of covariance between stations. If it is not\
        provided, covariance data will be calculated out of the given dataset.

    Results
    -------
        **spatial_average** -- Average precipitation based on punctual \
        readings, for all\
        the precipitation events in the serie.
    '''
    pre_average_guess = numpy.average(data, 1)
    station_weights_updated = weights(data, covariance)
    spatial_average = []
    for i in xrange(len(data)):
        spatial_average.append(pre_average_guess[i] +
                               numpy.dot(station_weights_updated, data[i, :] -
                               pre_average_guess[i]))

    spatial_average = numpy.clip(spatial_average, 0, max(spatial_average))

    return spatial_average


def weights(data, covariance=None):
    '''
    Calculates spatial average from several stations. Data consists in\
    a matrix element containing all the recordas from the stations. If\
    covariance matrix can be provided it can be passed as an optional argument.

    Paramters
    ----------
        **data** -- Matrix of recordings from stations. data has to be\
        oriented column wise

        **covariance** -- Matrix of covariance between stations. If it is not\
        provided, covariance data will be calculated out of the given dataset.

    Results
    -------
        **spatial_average** -- Average precipitation based on punctual \
        readings, for all\
        the precipitation events in the serie.
    '''
    if covariance is None:
        covariance = numpy.cov(numpy.transpose(data))

    data = numpy.array(data)
    if numpy.linalg.det(covariance) == 0:
        print 'Singular covariance matrix... cannot make it'
        return 9999*numpy.ones(len(data))

    pre_average_guess = numpy.average(data, 1)
    station_weights = []
    for i in xrange(len(covariance)):
        station_weights.append(numpy.cov(data[:, i], pre_average_guess)[0][1])
    station_weights_updated = numpy.dot(numpy.linalg.inv(covariance),
                                        station_weights)

    return station_weights_updated