学步园

//

#include "stdafx.h"
#include "fangshibo.h"

#include <vector>

#include<math.h>
///gdal头文件
#include "..//include//gdal.h"
#include "..//include//gdal_priv.h"
#include "..//include//ogr_srs_api.h"
#include "..//include//cpl_string.h"
#include "..//include//cpl_conv.h"
#pragma comment (lib,"..//lib//gdal_i.lib")
/////////////////////////////////////////////////////////////////////////////

#ifdef _DEBUG
#define new DEBUG_NEW
#undef THIS_FILE
static char THIS_FILE[] = __FILE__;
#endif

/////////////////////////////////////////////////////////////////////////////
// The one and only application object

CWinApp theApp;
using namespace std;

///////////////////////
void  MaxNDVI(vector<CString> files,CString maxNDVfile)
{    
//    GDALAllRegister();
    vector<GDALDataset *> datasets;    
        //打开文件
        for(int i=0;i<files.size();i++)
        {    
            CString filename=files[i];            
            GDALDataset *tmpDataset = (GDALDataset *) GDALOpen( files[i], GA_ReadOnly);
            if(tmpDataset!=NULL)
            {                
                datasets.push_back(tmpDataset);                
            }
            else
            {
                fprintf( stderr, "GDALOpen failed - %d/n%s/n", CPLGetLastErrorNo(), CPLGetLastErrorMsg() );
            }    
            tmpDataset=NULL;//以后再释放
        }
        ////读取数据咯，按行来读取
        //GDALDir
        if(datasets.size()==0)return;
        GDALDataset *tmpDataset=datasets[0];
        if(tmpDataset!=NULL)
        {    
            GDALDriver * driver=NULL;
            int index = maxNDVfile.ReverseFind('.');
            if (index < 0)
                return ;
            if (index == maxNDVfile.GetLength()-1)
                return ;
            CString suffix = maxNDVfile.Right(maxNDVfile.GetLength()-1-index);            
            suffix.MakeLower();
            if (suffix == "bmp")
                driver = GetGDALDriverManager()->GetDriverByName("BMP");
            else if (suffix == "jpg")
                driver = GetGDALDriverManager()->GetDriverByName("JPEG");
            else if (suffix == "tif")
                driver = GetGDALDriverManager()->GetDriverByName("GTiff");
            else if (suffix == "img")
                driver = GetGDALDriverManager()->GetDriverByName("HFA");
            else if (suffix == "bt")
                driver = GetGDALDriverManager()->GetDriverByName("BT");
            else if (suffix == "ecw")
                driver = GetGDALDriverManager()->GetDriverByName("ECW");
            else if (suffix == "fits")
                driver = GetGDALDriverManager()->GetDriverByName("FITS");
            else if (suffix == "gif")
                driver = GetGDALDriverManager()->GetDriverByName("GIF");
            else if (suffix == "hdf")
                driver = GetGDALDriverManager()->GetDriverByName("HDF4");
            else if (suffix == "hdr")
                driver = GetGDALDriverManager()->GetDriverByName("EHdr");
            
    
            int w=tmpDataset->GetRasterXSize();
            int h=tmpDataset->GetRasterYSize();            
            GDALDataset *maxNDV=driver->Create(maxNDVfile,w,h,1,GDT_Float32,NULL);        
            int xOff=0;
            int yOff=0;
            int width=w;
            int height=1;//一行一行地读取
            vector<float*> bufs;
            for(int i=0;i<datasets.size();i++)
            {
                float *buf=new float[width*height];
                bufs.push_back(buf);
                buf=NULL;
            }        
            float *newbuf=new float[width*height];
            GDALRasterBand * newpoband=maxNDV->GetRasterBand(1);        
            for(int j=0;j<h;j++)
            {                
                for(int i=0;i<datasets.size();i++)
                {
                    GDALDataset *tmpDt=datasets[i];
                    GDALRasterBand * poband=tmpDt->GetRasterBand(1);                
                    float *buf=bufs[i];
                    xOff=0;
                    yOff=j;
                    poband->RasterIO(GF_Read,xOff,yOff,width,height,buf,width,height,GDT_Float32,0,0);
                    buf=NULL;
                    tmpDt=NULL;//在后面再释放,因为很多个指针指向同一一段内存，不能随便释放内存哦
                }            
                //////////
                //在这可以开始运算咯
                for(int k=0;k<w;k++)//对于每个像元
                {
                    //下面这段代码以后要支持均值，极差，方差
                    for(int kk=0;kk<bufs.size();kk++)
                    {
                        float *buf=bufs[kk];
                        if(kk==0)//初始化
                        {
                            newbuf[k]=buf[k];
                        }
                        if(buf[k]>newbuf[k])//最大的,稍微改下就是最小值，
                        {
                            newbuf[k]=buf[k];
                        }                    
                    }
                }
                //////////写入数据
                newpoband->RasterIO(GF_Write,xOff,yOff,width,height,newbuf,width,height,GDT_Float32,0,0);
            }
            delete newbuf;
            newbuf=NULL;
            for( i=0;i<bufs.size();i++)
            {
                if(bufs[i]!=NULL)
                    delete bufs[i];            
            }
            bufs.resize(0);
            delete maxNDV;
            maxNDV=NULL;
            delete driver;
            driver=NULL;
        }
        tmpDataset=NULL;
        //释放掉内存
        for( i=0;i<datasets.size();i++)
        {
            if(datasets[i]!=NULL)
            {
                delete datasets[i];
                datasets[i]=NULL;
            }
        }
        datasets.resize(0);
}