/******************************************************************************
 * Project:  Open GIS conference 2012/11
 * Purpose:  master class examples
 * Author:   Bishop (aka Baryshnikov Dmitriy), polimax@mail.ru
 ******************************************************************************
*   Copyright (C) 2012 Bishop
*
*    This program is free software: you can redistribute it and/or modify
*    it under the terms of the GNU General Public License as published by
*    the Free Software Foundation, either version 3 of the License, or
*    (at your option) any later version.
*
*    This program is distributed in the hope that it will be useful,
*    but WITHOUT ANY WARRANTY; without even the implied warranty of
*    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
*    GNU General Public License for more details.
*
*    You should have received a copy of the GNU General Public License
*    along with this program.  If not, see <http://www.gnu.org/licenses/>.
 ****************************************************************************/

#define VECTOR
#define RASTER

#define APPROX
#define APRROX_MAXERROR 0.125

#define MULTI

#include "cpl_conv.h" // for CPLMalloc()
#include "gdal_priv.h"
#include "ogrsf_frmts.h"
#include "ogr_api.h"
#include "gdalwarper.h"

//производный класс для перепроецирования фигуры обрезки
class CutlineTransformer : public OGRCoordinateTransformation
{
public:

    void         *hSrcImageTransformer;

    virtual OGRSpatialReference *GetSourceCS() { return NULL; }
    virtual OGRSpatialReference *GetTargetCS() { return NULL; }

    virtual int Transform( int nCount, 
                           double *x, double *y, double *z = NULL ) {
        int nResult;

        int *pabSuccess = (int *) CPLCalloc(sizeof(int),nCount);
        nResult = TransformEx( nCount, x, y, z, pabSuccess );
        CPLFree( pabSuccess );

        return nResult;
    }

    virtual int TransformEx( int nCount, 
                             double *x, double *y, double *z = NULL,
                             int *pabSuccess = NULL ) {
        return GDALGenImgProjTransform( hSrcImageTransformer, TRUE, 
                                        nCount, x, y, z, pabSuccess );
    }
};

int main(int argc, char **argv)
{
    //настроим переменные окружения GDAL
    CPLSetConfigOption( "GDAL_FILENAME_IS_UTF8", "NO" );
    CPLSetConfigOption( "GDAL_DATA", "D:/work/gis-lab/conf/files/gdal-1.9.2/gdal_data" );
    CPLSetConfigOption( "CPL_DEBUG", "ON" );
    CPLSetConfigOption( "CPL_TIMESTAMP", "ON" );
    CPLSetConfigOption( "CPL_LOG_ERRORS", "ON" );

    //создадим выходную СК
    OGRSpatialReference oDstSpaRef(SRS_WKT_WGS84);
    OGRErr eErr = oDstSpaRef.SetUTM(37, 1);
    
    //создадим обрезающую геометрию
    double dfPM = oDstSpaRef.GetProjParm("central_meridian");
	OGRLinearRing ring;
    double dfWidth = 35;
	ring.addPoint(dfPM - dfWidth, -89);
	ring.addPoint(dfPM - dfWidth, 89);
	ring.addPoint(dfPM + dfWidth, 89);
	ring.addPoint(dfPM + dfWidth, -89);
    ring.closeRings();

    OGRPolygon* pCutPoly = new OGRPolygon();
    pCutPoly->addRing(&ring);
    pCutPoly->flattenTo2D();
    pCutPoly->assignSpatialReference(new OGRSpatialReference(SRS_WKT_WGS84));
    pCutPoly->segmentize(1);

    OGRGeometry* pCutGeom = static_cast<OGRGeometry*>(pCutPoly);

    //определим выходной охват
    OGRGeometry* pCutGeomDst = pCutGeom->clone();
    pCutGeomDst->transformTo(&oDstSpaRef);
    OGREnvelope DstEnv;
    pCutGeomDst->getEnvelope(&DstEnv);

#ifdef RASTER
    //загрузка драйверов поддержки растровых форматов
    GDALAllRegister();
    
    //открытие рабочего набора
    GDALDataset  *poSrcDataset = (GDALDataset *) GDALOpen( "/vsizip/D:/work/gis-lab/conf/files/data/TrueMarble.4km.10800x5400.tif.zip/TrueMarble.4km.10800x5400.tif", GA_ReadOnly );
    if( NULL == poSrcDataset )
    {
        const char* err = CPLGetLastErrorMsg();
        printf( "Open failed.%s\n", err );
        exit( 1 );
    }

    //получаем СК из файла
    OGRSpatialReference SrcSpaRef( poSrcDataset->GetProjectionRef() );
    char* pszSpaRefDef(NULL);
    if( SrcSpaRef.exportToPrettyWkt(&pszSpaRefDef) == OGRERR_NONE)
    {
        printf(pszSpaRefDef);
        printf("\n");
        CPLFree( pszSpaRefDef );
    }
    
    char **papszTO = NULL;

    if( oDstSpaRef.exportToWkt(&pszSpaRefDef) == OGRERR_NONE)
    {
        papszTO = CSLSetNameValue( papszTO, "DST_SRS", pszSpaRefDef );
        CPLFree( pszSpaRefDef );
    } 

    //создание выходного растра
    GDALDriver* poDriver = poSrcDataset->GetDriver();
    GDALRasterBand *poGDALRasterBand = poSrcDataset->GetRasterBand(1);

    if(NULL == poGDALRasterBand)
    {
        printf( "Band get failed\n" );
        exit( 1 );
    }
    GDALDataType eDT = poGDALRasterBand->GetRasterDataType();
    
    CPLSetConfigOption( "CHECK_WITH_INVERT_PROJ", "TRUE" );
    void* hTransformArg = GDALCreateGenImgProjTransformer2( poSrcDataset, NULL, papszTO );
    GDALTransformerInfo* psInfo = (GDALTransformerInfo*)hTransformArg;
    
    double adfThisGeoTransform[6];
    double adfExtent[4];
    int nThisPixels, nThisLines;
    
    //оцениваемы размеры выходного файла
    if( GDALSuggestedWarpOutput2( poSrcDataset, psInfo->pfnTransform, hTransformArg, adfThisGeoTransform, &nThisPixels, &nThisLines, adfExtent, 0 ) != CE_None )
    {
        printf( "Suggest Output failed\n" );
        exit( 1 );
    }

    //корректируем размер согласно обрезающей фигуры
    adfThisGeoTransform[0] = DstEnv.MinX;
    adfThisGeoTransform[3] = DstEnv.MaxY;

    int nPixels = (int) ((DstEnv.MaxX - DstEnv.MinX) / adfThisGeoTransform[1] + 0.5);
    int nLines = (int) ((DstEnv.MaxY - DstEnv.MinY) / -adfThisGeoTransform[5] + 0.5);

    GDALDataset  *poDstDataset = poDriver->Create("D:/work/gis-lab/conf/files/data/test.tif", nPixels, nLines, poSrcDataset->GetRasterCount(), eDT, NULL );
    if( NULL == poDstDataset )
    {
        printf( "Create Output failed\n" );
        exit( 1 );
    }
    
    GDALSetGenImgProjTransformerDstGeoTransform( hTransformArg, adfThisGeoTransform);
    
    char* dstSrs;
    oDstSpaRef.exportToWkt(&dstSrs);
    poDstDataset->SetProjection( dstSrs );
    poDstDataset->SetGeoTransform( adfThisGeoTransform );

#ifdef APPROX
    hTransformArg = GDALCreateApproxTransformer( GDALGenImgProjTransform,  hTransformArg, APRROX_MAXERROR);
    GDALTransformerFunc pfnTransformer = GDALApproxTransform;
#else //APPROX
    GDALTransformerFunc pfnTransformer = GDALGenImgProjTransform;
#endif //APPROX

    //подготовка обрезающей фигуры
    papszTO = CSLSetNameValue( papszTO, "INSERT_CENTER_LONG", "FALSE" );

    CutlineTransformer oTransformer;
    oTransformer.hSrcImageTransformer = GDALCreateGenImgProjTransformer2( poSrcDataset, NULL, papszTO );

    pCutGeomDst->transform((OGRCoordinateTransformation*)&oTransformer);
    char* pszWKT = NULL;
    if(pCutGeomDst->exportToWkt(&pszWKT) != OGRERR_NONE )
    {
        const char* err = CPLGetLastErrorMsg();
        printf( "Open failed.%s\n", err );
        exit( 1 );
    }


    //создание и настройка Warp
    GDALWarpOptions *psWO = GDALCreateWarpOptions();
    psWO->papszWarpOptions = CSLSetNameValue(psWO->papszWarpOptions, "INIT_DEST", "0"); 

    psWO->papszWarpOptions = CSLSetNameValue( psWO->papszWarpOptions, "CUTLINE", pszWKT );
    CPLFree( pszWKT );

    psWO->eWorkingDataType = eDT;
    psWO->eResampleAlg = GRA_Bilinear;

    psWO->hSrcDS = poSrcDataset;
    psWO->hDstDS = poDstDataset;

    psWO->pfnTransformer = pfnTransformer;
    psWO->pTransformerArg = hTransformArg;

    psWO->pfnProgress = GDALTermProgress;
    psWO->nBandCount = poSrcDataset->GetRasterCount();

    psWO->panSrcBands = (int *) CPLMalloc(psWO->nBandCount*sizeof(int));
    psWO->panDstBands = (int *) CPLMalloc(psWO->nBandCount*sizeof(int));

    for(int i = 0; i < psWO->nBandCount; ++i )
    {
        psWO->panSrcBands[i] = i+1;
        psWO->panDstBands[i] = i+1;
    }

    //выполнение Warp
    GDALWarpOperation oWO;
    if( oWO.Initialize( psWO ) == CE_None )
    {
#ifdef MULTI
        if( oWO.ChunkAndWarpMulti( 0, 0, poDstDataset->GetRasterXSize(), poDstDataset->GetRasterYSize() ) != CE_None)
#else //MULTI
        if( oWO.ChunkAndWarpImage( 0, 0, poDstDataset->GetRasterXSize(), poDstDataset->GetRasterYSize() ) != CE_None)
#endif //MULTI
        {
            const char* err = CPLGetLastErrorMsg();
            printf( "Warp failed.%s\n", err );
            exit( 1 );
        }
    }

    GDALDestroyWarpOptions( psWO );
    CSLDestroy( papszTO );

    //выгрузка драйверов поддержки растровых форматов
    GDALDestroyDriverManager();
#endif //RASTER

#ifdef VECTOR
    //загрузка драйверов поддержки векторных форматов
    OGRRegisterAll();

    //выведем список драйверов которые поддерживает GDAL 
	for(size_t i = 0; i < OGRSFDriverRegistrar::GetRegistrar()->GetDriverCount(); ++i)
	{
		OGRSFDriver *pDrv = OGRSFDriverRegistrar::GetRegistrar()->GetDriver(i);
		const char* pszDrvName = pDrv->GetName();
        printf(pszDrvName);
        printf("\n");
	}


    //откроем тестовый векторный файл
    OGRDataSource *poDS = OGRSFDriverRegistrar::Open( "/vsizip/D:/work/gis-lab/conf/files/data/ne_50m_admin_0_countries_lakes.zip/ne_50m_admin_0_countries_lakes.shp", FALSE );
    if( NULL == poDS )
    {
        const char* err = CPLGetLastErrorMsg();
        printf( "Open failed.%s\n", err );
        exit( 1 );
    }

    OGRLayer *pSrcLayer = poDS->GetLayer(0);
    if( NULL == pSrcLayer )
    {
        const char* err = CPLGetLastErrorMsg();
        printf( "Open failed.%s\n", err );
        exit( 1 );
    }

    //выведем текущую проекцию
    OGRSpatialReference* pSrcSpaRef = pSrcLayer->GetSpatialRef();
    if( NULL == pSrcSpaRef )
    {
        const char* err = CPLGetLastErrorMsg();
        printf( "Open failed.%s\n", err );
        exit( 1 );
    }

    char* pSpaRefDef(NULL);
    if( pSrcSpaRef->exportToPrettyWkt(&pSpaRefDef) == OGRERR_NONE)
    {
        printf(pSpaRefDef);
        CPLFree( pSpaRefDef );
    }

    //создание выходного файла

    OGRSFDriver *poVecDriver = poDS->GetDriver();
    if( NULL == poVecDriver )
    {
        printf( "Driver not available.\n");
        exit( 1 );
    }

    //создаем выходной OGRDataSource
    OGRDataSource *poDstDS = poVecDriver->CreateDataSource( "D:/work/gis-lab/conf/files/data/test.shp", NULL );
    if( NULL == poDstDS )
    {
        printf( "Creation of output file failed.\n" );
        exit( 1 );
    }
    
    //создаем выходной OGRLayer
    OGRLayer *poDstLayer = poDstDS->CreateLayer( "test", oDstSpaRef.Clone(), pSrcLayer->GetGeomType(), NULL );
    if( poDstLayer == NULL )
    {
        printf( "Layer creation failed.\n" );
        exit( 1 );
    }

    //добавляем поля
    OGRFeatureDefn* poFields = pSrcLayer->GetLayerDefn();
    for(size_t i = 0; i < poFields->GetFieldCount(); ++i)
    {
        OGRFieldDefn *pField = poFields->GetFieldDefn(i);
	    if( poDstLayer->CreateField( pField ) != OGRERR_NONE )
	    {
            printf("Error create the output layer '%s'! OGR error: %s", CPLGetLastErrorMsg());
            exit( 1 );
        }
    }

    //обрезаем и перепроецируем
    OGRCoordinateTransformation *poCT = OGRCreateCoordinateTransformation( pSrcSpaRef, &oDstSpaRef);
    if(poCT)
    { 
        OGRFeature* pSrcFeature;
        while((pSrcFeature = pSrcLayer->GetNextFeature()) != NULL)
        {
            OGRGeometry *pSrcGeom = pSrcFeature->GetGeometryRef();
            if(pSrcGeom)
            {
                OGRGeometry* pDstGeom = pSrcGeom->Intersection(pCutGeom);
                if(pDstGeom && pDstGeom->transform(poCT) == OGRERR_NONE)
                {
                    OGRFeature *poDstFeature = OGRFeature::CreateFeature( poDstLayer->GetLayerDefn() );                     
                    OGRFeatureDefn* pSrcFeatureDefn = pSrcFeature->GetDefnRef();
                    for(size_t j = 0; j < pSrcFeatureDefn->GetFieldCount(); ++j)
                    {
                        poDstFeature->SetField(j, pSrcFeature->GetRawFieldRef(j));
                    }
                    poDstFeature->SetGeometryDirectly(pDstGeom);

                    if( poDstLayer->CreateFeature( poDstFeature ) != OGRERR_NONE )
                    {
                        printf( "Failed to create feature in shapefile.\n" );
                        exit( 1 );
                    }

                    OGRFeature::DestroyFeature( poDstFeature );
                }
            }
        }

        OCTDestroyCoordinateTransformation(poCT);
    }

    OGRDataSource::DestroyDataSource( poDS );

    OGRDataSource::DestroyDataSource( poDstDS );

    //выгрузка драйверов поддержки векторных форматов
	OGRCleanupAll();
#endif //VECTOR

    return 0;
}
