#include <gdal/gdal_priv.h>
#include <gdal/cpl_conv.h>
#include <gdal/ogr_spatialref.h>
#include <stdio.h>
#include <math.h>
#include <fstream>
#include <sys/stat.h>
#include <vector>
#include <algorithm>
#include <time.h>
#include <string>
#include <string.h>
#include <stdlib.h>
#include <stdint.h>
#include <iomanip>
#include <iostream>
#include <sstream>
using namespace std;

#define FILENUMBER 504
#define YEAR 14

int g_nHeight, g_nWidth, g_lImageSize;

char g_strCount[14][100];
char g_strBand[7][14][100]; //band number 7: band 3 4 5 7, minimum ndvi, bare percent and nbr
char g_strPath[7][200];

char g_strOut_feature[18][100]; // feature number for each band 18
char g_strOut_band[7][100];
char g_strPro_band[7][100];

char temp[200],temp1[200];

void Set_string()
{//set path name
	if (true)
	{
	strcpy(g_strPath[0],"/gfs2/d18/Google/V2_30m/01_tiles/annual_composite/"); //band3
	strcpy(g_strPath[1],"/gfs2/d18/Google/V2_30m/01_tiles/annual_composite/");//band4
	strcpy(g_strPath[2],"/gfs2/d18/Google/V2_30m/01_tiles/annual_composite/");//band5
	strcpy(g_strPath[3],"/gfs2/d18/Google/V2_30m/01_tiles/annual_composite/");//band7
	strcpy(g_strPath[4],"/gfs2/d18/Google/V2_30m/01_tiles/annual_min_ndvi/"); //minimum ndvi
	strcpy(g_strPath[5],"/gfs2/d18/Google/V2_30m/01_tiles/annual_bare/"); //bare percentage
	strcpy(g_strPath[6],"/gfs2/d13/data/nbr/"); // NBR
	}
	
	//for count
	if (true)
	{
	strcpy(g_strCount[0],"annual_count_1999-");
	strcpy(g_strCount[1],"annual_count_2000-");
	strcpy(g_strCount[2],"annual_count_2001-");
	strcpy(g_strCount[3],"annual_count_2002-");
	strcpy(g_strCount[4],"annual_count_2003-");
	strcpy(g_strCount[5],"annual_count_2004-");
	strcpy(g_strCount[6],"annual_count_2005-");
	strcpy(g_strCount[7],"annual_count_2006-");
	strcpy(g_strCount[8],"annual_count_2007-");
	strcpy(g_strCount[9],"annual_count_2008-");
	strcpy(g_strCount[10],"annual_count_2009-");
	strcpy(g_strCount[11],"annual_count_2010-");
	strcpy(g_strCount[12],"annual_count_2011-");
	strcpy(g_strCount[13],"annual_count_2012-");
	}
	
	//for g_strBand
	short I=0;
	if (true)
	{
	strcpy(g_strBand[I][0],	"1999/annual_1999_30-");
	strcpy(g_strBand[I][1],	"2000/annual_2000_30-");
	strcpy(g_strBand[I][2],	"2001/annual_2001_30-");
	strcpy(g_strBand[I][3],	"2002/annual_2002_30-");
	strcpy(g_strBand[I][4],	"2003/annual_2003_30-");
	strcpy(g_strBand[I][5],	"2004/annual_2004_30-");
	strcpy(g_strBand[I][6],	"2005/annual_2005_30-");
	strcpy(g_strBand[I][7],	"2006/annual_2006_30-");
	strcpy(g_strBand[I][8],	"2007/annual_2007_30-");
	strcpy(g_strBand[I][9],	"2008/annual_2008_30-");
	strcpy(g_strBand[I][10], "2009/annual_2009_30-");
	strcpy(g_strBand[I][11], "2010/annual_2010_30-");
	strcpy(g_strBand[I][12], "2011/annual_2011_30-");
	strcpy(g_strBand[I][13], "2012/annual_2012_30-");
	I=1;

	strcpy(g_strBand[I][0],	"1999/annual_1999_40-");
	strcpy(g_strBand[I][1],	"2000/annual_2000_40-");
	strcpy(g_strBand[I][2],	"2001/annual_2001_40-");
	strcpy(g_strBand[I][3],	"2002/annual_2002_40-");
	strcpy(g_strBand[I][4],	"2003/annual_2003_40-");
	strcpy(g_strBand[I][5],	"2004/annual_2004_40-");
	strcpy(g_strBand[I][6],	"2005/annual_2005_40-");
	strcpy(g_strBand[I][7],	"2006/annual_2006_40-");
	strcpy(g_strBand[I][8],	"2007/annual_2007_40-");
	strcpy(g_strBand[I][9],	"2008/annual_2008_40-");
	strcpy(g_strBand[I][10],	"2009/annual_2009_40-");
	strcpy(g_strBand[I][11],	"2010/annual_2010_40-");
	strcpy(g_strBand[I][12],	"2011/annual_2011_40-");
	strcpy(g_strBand[I][13],	"2012/annual_2012_40-");

	I=2;
	strcpy(g_strBand[I][0],	"1999/annual_1999_50-");
	strcpy(g_strBand[I][1],	"2000/annual_2000_50-");
	strcpy(g_strBand[I][2],	"2001/annual_2001_50-");
	strcpy(g_strBand[I][3],	"2002/annual_2002_50-");
	strcpy(g_strBand[I][4],	"2003/annual_2003_50-");
	strcpy(g_strBand[I][5],	"2004/annual_2004_50-");
	strcpy(g_strBand[I][6],	"2005/annual_2005_50-");
	strcpy(g_strBand[I][7],	"2006/annual_2006_50-");
	strcpy(g_strBand[I][8],	"2007/annual_2007_50-");
	strcpy(g_strBand[I][9],	"2008/annual_2008_50-");
	strcpy(g_strBand[I][10],	"2009/annual_2009_50-");
	strcpy(g_strBand[I][11],	"2010/annual_2010_50-");
	strcpy(g_strBand[I][12],	"2011/annual_2011_50-");
	strcpy(g_strBand[I][13],	"2012/annual_2012_50-");


	I=3;

	strcpy(g_strBand[I][0],	"1999/annual_1999_70-");
	strcpy(g_strBand[I][1],	"2000/annual_2000_70-");
	strcpy(g_strBand[I][2],	"2001/annual_2001_70-");
	strcpy(g_strBand[I][3],	"2002/annual_2002_70-");
	strcpy(g_strBand[I][4],	"2003/annual_2003_70-");
	strcpy(g_strBand[I][5],	"2004/annual_2004_70-");
	strcpy(g_strBand[I][6],	"2005/annual_2005_70-");
	strcpy(g_strBand[I][7],	"2006/annual_2006_70-");
	strcpy(g_strBand[I][8],	"2007/annual_2007_70-");
	strcpy(g_strBand[I][9],	"2008/annual_2008_70-");
	strcpy(g_strBand[I][10],	"2009/annual_2009_70-");
	strcpy(g_strBand[I][11],	"2010/annual_2010_70-");
	strcpy(g_strBand[I][12],	"2011/annual_2011_70-");
	strcpy(g_strBand[I][13],	"2012/annual_2012_70-");
	
	I=4;
	strcpy(g_strBand[I][0],"annual_min_ndvi_1999-");
	strcpy(g_strBand[I][1],"annual_min_ndvi_2000-");
	strcpy(g_strBand[I][2],"annual_min_ndvi_2001-");
	strcpy(g_strBand[I][3],"annual_min_ndvi_2002-");
	strcpy(g_strBand[I][4],"annual_min_ndvi_2003-");
	strcpy(g_strBand[I][5],"annual_min_ndvi_2004-");
	strcpy(g_strBand[I][6],"annual_min_ndvi_2005-");
	strcpy(g_strBand[I][7],"annual_min_ndvi_2006-");
	strcpy(g_strBand[I][8],"annual_min_ndvi_2007-");
	strcpy(g_strBand[I][9],"annual_min_ndvi_2008-");
	strcpy(g_strBand[I][10],"annual_min_ndvi_2009-");
	strcpy(g_strBand[I][11],"annual_min_ndvi_2010-");
	strcpy(g_strBand[I][12],"annual_min_ndvi_2011-");
	strcpy(g_strBand[I][13],"annual_min_ndvi_2012-");

	I=5;
	strcpy(g_strBand[I][0],"bare_1999-");
	strcpy(g_strBand[I][1],"bare_2000-");
	strcpy(g_strBand[I][2],"bare_2001-");
	strcpy(g_strBand[I][3],"bare_2002-");
	strcpy(g_strBand[I][4],"bare_2003-");
	strcpy(g_strBand[I][5],"bare_2004-");
	strcpy(g_strBand[I][6],"bare_2005-");
	strcpy(g_strBand[I][7],"bare_2006-");
	strcpy(g_strBand[I][8],"bare_2007-");
	strcpy(g_strBand[I][9],"bare_2008-");
	strcpy(g_strBand[I][10],"bare_2009-");
	strcpy(g_strBand[I][11],"bare_2010-");
	strcpy(g_strBand[I][12],"bare_2011-");
	strcpy(g_strBand[I][13],"bare_2012-");
	
	I=6;
	strcpy(g_strBand[I][0],"annual_1999-");
	strcpy(g_strBand[I][1],"annual_2000-");
	strcpy(g_strBand[I][2],"annual_2001-");
	strcpy(g_strBand[I][3],"annual_2002-");
	strcpy(g_strBand[I][4],"annual_2003-");
	strcpy(g_strBand[I][5],"annual_2004-");
	strcpy(g_strBand[I][6],"annual_2005-");
	strcpy(g_strBand[I][7],"annual_2006-");
	strcpy(g_strBand[I][8],"annual_2007-");
	strcpy(g_strBand[I][9],"annual_2008-");
	strcpy(g_strBand[I][10],"annual_2010-");
	strcpy(g_strBand[I][11],"annual_2010-");
	strcpy(g_strBand[I][12],"annual_2011-");
	strcpy(g_strBand[I][13],"annual_2012-");
	}
	if (true)
	{
	    strcpy(g_strPro_band[0],"b30"); //band3
		strcpy(g_strPro_band[1],"b40");//band4
		strcpy(g_strPro_band[2],"b50");//band5
		strcpy(g_strPro_band[3],"b70");//band7
		strcpy(g_strPro_band[4],"minndvi"); //minimum ndvi
		strcpy(g_strPro_band[5],"bare"); //bare percentage
		strcpy(g_strPro_band[6],"nbr"); // NBR
		
		strcpy(g_strOut_band[0],"b30-"); //band3
		strcpy(g_strOut_band[1],"b40-");//band4
		strcpy(g_strOut_band[2],"b50-");//band5
		strcpy(g_strOut_band[3],"b70-");//band7
		strcpy(g_strOut_band[4],"minndvi-"); //minimum ndvi
		strcpy(g_strOut_band[5],"bare-"); //bare percentage
		strcpy(g_strOut_band[6],"nbr-"); // NBR
		
	}
	
}	

int main(void)
{
	char *Maskpath="/gfs2/d18/Google/V2_30m/01_tiles/datamask/data_mask-";
	char *Countpath="/gfs2/d18/Google/V2_30m/01_tiles/annual_count/"; 
	char *Outpath="/gfs2/d14/test/";
	Set_string();
	
	int I,J,K,L,M,N,P;
	int itemp;
	unsigned char m_value[14], m_year[14];
	unsigned char m_first[3], m_last[3];
	unsigned char m_byTemp,m_byMin;
	
	char str_flag[200]="end";
	char filename[200];
	
	char strCount[14][200];
	char strWork[14][200];
	char strMask[200];

	char strOut_profile[200];
	
	int band_sum[2][7][14], band_count[2][7][14];
	int band_ave[2][7][14];
	double Sx,Sy,Sxx,Sxy,Syy;
	float x1;
	
	GDALAllRegister();
	GDALDataset  *poDataset,*poDataset1,*poDataset2;
	GDALRasterBand *poBand;
	
	char *ppszOptions[] ={NULL} ;
	GDALDriver *poDriver = GetGDALDriverManager()->GetDriverByName("GTiff");
	GDALDriver *poDriver1 = GetGDALDriverManager()->GetDriverByName("ENVI");
		if( poDriver == NULL )
			exit(1);

	GDALDataset *m_pDatasetWriting; 	
	double plfTrans[6];
	CPLErr errNum;

	g_lImageSize =36000 *36000;
	
	unsigned char **m_pbyCount, **m_pbyBand,*m_pbyMask;

	m_pbyCount=new unsigned char *[14];
	m_pbyBand=new unsigned char *[14];
	
	
	for (I=0;I<14;I++)
	{
		m_pbyCount[I]	=new unsigned char [g_lImageSize];
		m_pbyBand[I]	=new unsigned char [g_lImageSize];
	}
			
	m_pbyMask=new unsigned char [g_lImageSize];
	unsigned char  m_pbyOut_avail=0;

	
				
		strcpy(temp,Maskpath);
		strcat(temp,"0000144000-0000324000");	
		strcat(temp,".tif");	
		strcpy(strMask,temp);


		for(I=0;I<14;I++)
		{
		strcpy(temp,Countpath);
		strcat(temp,g_strCount[I]);
		strcat(temp,"0000144000-0000324000");	
		strcat(temp,".tif");	
		strcpy(strCount[I],temp);		
		}
		
		strcpy(temp,Outpath);				
		strcat(temp,"1948");	
		strcat(temp,".txt");	
		strcpy(strOut_profile,temp);
		
		//output profile text
        ofstream out_profile;
		out_profile.open(strOut_profile, ios::app);
		//initialize result-carrying arrays
		for (M=0;M<2;M++) {for (I=0;I<7;I++) {for (J=0;J<14;J++) {band_sum[M][I][J]=0; band_count[M][I][J]=0; band_ave[M][I][J]=0;}}}
		
	poDataset = (GDALDataset *) GDALOpen( strMask, GA_ReadOnly );
    g_nHeight=poDataset->GetRasterYSize(); //get image height
	g_nWidth=poDataset->GetRasterXSize();  //get image width
	
	poBand=poDataset->GetRasterBand(1);
	poBand->RasterIO(GF_Read, 0, 0, g_nWidth, g_nHeight, m_pbyMask,
		g_nWidth, g_nHeight, GDT_Byte, 0, 0);
	
	// read count data for a certain file
		for (I=0;I<14;I++)
		{				
		poDataset1 = (GDALDataset *) GDALOpen( strCount[I], GA_ReadOnly );    				
		poBand= poDataset1->GetRasterBand(1);
		poBand->RasterIO(GF_Read, 0, 0, g_nWidth, g_nHeight, m_pbyCount[I],
							g_nWidth, g_nHeight, GDT_Byte, 0, 0);// read the NDVI tiles.
		GDALClose( (GDALDatasetH) poDataset1 );
		}
	//count valuable data in 14 years within land. 

	I=34938*36000+8565;
		if(*(m_pbyMask+I)>0) //land
		{
			for(J=0;J<14;J++)
			if(*(m_pbyCount[J]+I)>0)
					m_pbyOut_avail=m_pbyOut_avail+1;		
		}
		for(M=0;M<7;M++) //band
		{
			
			if(true)
			{//set file name
				for(I=0;I<14;I++) //year
				{
				strcpy(temp,g_strPath[M]);
				strcat(temp,g_strBand[M][I]);
				strcat(temp,"0000144000-0000324000");	
				strcat(temp,".tif");	
				strcpy(strWork[I],temp);
				}	
			
			}
		
		//*****************************************
		
			for (I=0;I<14;I++)
			{	//read the data from file			
			poDataset1 = (GDALDataset *) GDALOpen( strWork[I], GA_ReadOnly );    				
			poBand= poDataset1->GetRasterBand(1);
			poBand->RasterIO(GF_Read, 0, 0, g_nWidth, g_nHeight, m_pbyBand[I],
								g_nWidth, g_nHeight, GDT_Byte, 0, 0);// read the NDVI tiles.
			GDALClose( (GDALDatasetH) poDataset1 );
			}
			
		//*****************************************	//
		I=34938*36000+8565;
			if(m_pbyOut_avail>0)  //at least have data
			{
			// band profiles
			if (*(m_pbyMask+I)==1){
			   for (J=0;J<14;J++){
			     if(*(m_pbyCount[J]+I)>0){
				    band_sum[0][M][J] = *(m_pbyBand[J]+I); 
				 }}}
            }


	//write out profile
    I=0;
	  out_profile<<g_strPro_band[M];
	  for(J=0;J<14;J++){
             out_profile<<","<<band_sum[I][M][J];
		}
		out_profile<<endl;
	
	}//M all bands
	
	GDALClose( (GDALDatasetH) poDataset );	
	out_profile.close();	

		for (I=0;I<14;I++)
	{
		delete []m_pbyCount[I];
		delete []m_pbyBand[I];
	}
	
	delete []m_pbyCount;
	delete []m_pbyBand;		
	delete []m_pbyMask;
	return 0;
}