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

之前遇到的一个客户需求就是计算标准化降水指数（SPI）。

主要参考的论文为：洪兴骏，等 “标准化降水指数SPI分布函数的适用性研究”。原文可以很容易搜索到。

论文将计算方法其实写的比较清楚了，本文主要提供其C++实现。

头文件

//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

实现文件：

#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);
}

其实数学的东西实现起来没有特别复杂，两个难点在于积分（概率密度）以及 gamma 函数，在代码里也写清楚了。

使用的时候只需要按代码注释中说明的格式传入每年的降水量就行了。调用 getSPI (int year) 函数可以获取当年的 SPI。

注意，本程序基于合理的数据进行检测。不合理的数据（如每年的降水量完全相同，被视为不合理的数据）将引起意料之外的结果。

最后

以上就是无辜老鼠为你收集整理的标准化降水指数（SPI）计算程序的全部内容，希望文章能够帮你解决标准化降水指数（SPI）计算程序所遇到的程序开发问题。

如果觉得靠谱客网站的内容还不错，欢迎将靠谱客网站推荐给程序员好友。