标准化降水指数（SPI）计算程序

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
//SPI_Calculator.h
#include<iostream>
#include<vector>
#include<utility>
#include<fstream>
#include<cmath>


using namespace std;

#ifndef SPICALCULATOR_H
#define SPICALCULATOR_H



class SPI_Calculator{
public:
	/**
	path is the file that provides the data of rain fall
	file format:
	year	amount of precipitation
	year	amount of precipitation
	.
	.
	.
	*/
	SPI_Calculator(char* path);

	/**
	parameter is the year to get the SPI
	return SPI
	*/
	double getSPI(int year);

	double f_probabilityDensity(double x);
	double f_gamma(double z, unsigned Ntermsx);
	double integral(double a,double b);

	void printYears();
	void printSPI();

	void test();

private:

	bool readData(char* path);
	void calculateParameter();
	
	vector<pair<int,double>> dataTable;
	double dataCounter,zeroCounter;

	double xAver; //average number of samples
	double c0,c1,c2,d1,d2,d3; //const parameter
	double para_A, para_Gamma, para_Beta;

};


#endif

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
#include"SPI_Calculator.h"

SPI_Calculator::SPI_Calculator(char* path){
	c0 = 2.515517;
	c1 = 0.802853;
	c2 = 0.010328;
	d1 = 1.432788;
	d2 = 0.189269;
	d3 = 0.001308;
	
	dataCounter = 0;
	zeroCounter = 0;
	if(!readData(path)){
		cout<<"Cannot read data successfully"<<endl;
		system("pause");
		exit(0);
	}

	calculateParameter();
}



bool SPI_Calculator::readData(char* path){
	fstream dataIn(path,ios::in);
	if(!dataIn.is_open()){
		cout<<"Cannot open file: "<<path<<endl;
		return false;
	}

	pair<int,double> tempPair;
	while(!dataIn.eof()){
		dataIn>>tempPair.first;
		dataIn>>tempPair.second;
		dataTable.push_back(tempPair);
		dataCounter ++;
		if(tempPair.second==0)
			zeroCounter ++;
		if(tempPair.second < 0){
			cout<<"Illigal perciption "<< tempPair.second<<" in year"<<tempPair.first<<endl;
			dataIn.close();
			return false;
		}
	}
	dataIn.close();
	return true;
}


void SPI_Calculator::printYears(){
	cout<<"nGet the following index: "<<endl;
	for(int i = 0; i < dataCounter; i ++)
		cout<<dataTable[i].first<<(((i+1)%5)?"t":"n");
	cout<<endl;
	
}

void SPI_Calculator::printSPI(){
	cout<<"nSPI: "<<endl;
	for(int i = 0; i < dataCounter; i ++){
		cout<<dataTable[i].first<< " SPI: "<<getSPI(dataTable[i].first)<<endl;;
	
	}
	cout<<endl;
	
}


void SPI_Calculator::calculateParameter(){
	double sum = 0;
	for(int i = 0; i < dataCounter; i ++)
		sum += dataTable[i].second;
	xAver = sum / (dataCounter-zeroCounter);

	para_A = log(xAver);
	double temp  = 0.0;
	for(int i = 0; i < dataCounter; i ++)
		if((dataTable[i].second!=0))
			temp += log(dataTable[i].second);
	

	para_A -= (temp / (dataCounter-zeroCounter));
	para_Gamma = (1 + sqrt(1+4*para_A/3))/ (4* para_A);
	para_Beta = xAver / para_Gamma;

}

double SPI_Calculator::getSPI(int year){
	double perciption;
	int i = 0;
	for(; i < dataCounter; i ++)
		if(year == dataTable[i].first){
			perciption = dataTable[i].second;
			break;
		}
	if(i == dataCounter){
		cout<<"Cannot find the year"<<endl;
		return -999;
	}

	double F = 0.0;
	if(perciption == 0){
		F=(zeroCounter / dataCounter);
	}
	else
		F = integral(0.0000001, perciption);
	
	double S = 0.0;

	if(F > 0.5){
		F = 1-F;
		S = 1;
	}
	else
		S = -1;
	double t = sqrt(log(1/(F*F)));


	double Z = S*(t - (c0+c1*t+c2*t*t)/(1+d1*t+d2*t*t+d3*t*t*t));
	return Z;

	
}

double SPI_Calculator::f_probabilityDensity(double x){
	double E = 2.71828;
	double t_a = 1/(pow(para_Beta,para_Gamma)* f_gamma(para_Gamma,1000));
	double t_b = pow(x, para_Gamma-1);
	double t_c = pow(E, -x/para_Beta);
	return (t_a *t_b *t_c);
}

double SPI_Calculator::f_gamma(double z, unsigned Nterms){
<span style="white-space:pre">	</span>double g = 0.577216;// Euler-Mascheroni constant
    double retVal = z*exp(g*z);

    // apply products
    for(unsigned n=1; n<=Nterms; ++n)
        retVal *= (1 + z/n)*exp(-z/n);

    retVal = 1.0/retVal;// invert

    return retVal;

}


double SPI_Calculator::integral(double a,double b)  {
	int N = 100000*para_A*para_A;
	double s = 0,h;
	s=(f_probabilityDensity(a) + f_probabilityDensity(b))/2.0;
	h=(b-a)/N;
	for(int i=1; i<N; i++)
		s += f_probabilityDensity(a+i*h);

	return (s*h);
}