C++ VRTDataset::SetGeoTransform方法代码示例

static GDALDataset *
VRTCreateCopy( const char * pszFilename,
               GDALDataset *poSrcDS,
               int /* bStrict */,
               char ** /* papszOptions */,
               GDALProgressFunc /* pfnProgress */,
               void * /* pProgressData */ )
{
    CPLAssert( NULL != poSrcDS );

/* -------------------------------------------------------------------- */
/*      If the source dataset is a virtual dataset then just write      */
/*      it to disk as a special case to avoid extra layers of           */
/*      indirection.                                                    */
/* -------------------------------------------------------------------- */
    if( poSrcDS->GetDriver() != NULL &&
        EQUAL(poSrcDS->GetDriver()->GetDescription(),"VRT") )
    {

    /* -------------------------------------------------------------------- */
    /*      Convert tree to a single block of XML text.                     */
    /* -------------------------------------------------------------------- */
        char *pszVRTPath = CPLStrdup(CPLGetPath(pszFilename));
        reinterpret_cast<VRTDataset *>(
            poSrcDS )->UnsetPreservedRelativeFilenames();
        CPLXMLNode *psDSTree = reinterpret_cast<VRTDataset *>(
            poSrcDS )->SerializeToXML( pszVRTPath );

        char *pszXML = CPLSerializeXMLTree( psDSTree );

        CPLDestroyXMLNode( psDSTree );

        CPLFree( pszVRTPath );

    /* -------------------------------------------------------------------- */
    /*      Write to disk.                                                  */
    /* -------------------------------------------------------------------- */
        GDALDataset* pCopyDS = NULL;

        if( 0 != strlen( pszFilename ) )
        {
            VSILFILE *fpVRT = VSIFOpenL( pszFilename, "wb" );
            if( fpVRT == NULL )
            {
                CPLError(CE_Failure, CPLE_AppDefined,
                         "Cannot create %s", pszFilename);
                CPLFree( pszXML );
                return NULL;
            }

            bool bRet = VSIFWriteL( pszXML, strlen(pszXML), 1, fpVRT ) > 0;
            if( VSIFCloseL( fpVRT ) != 0 )
                bRet = false;

            if( bRet )
                pCopyDS = reinterpret_cast<GDALDataset *>(
                    GDALOpen( pszFilename, GA_Update ) );
        }
        else
        {
            /* No destination file is given, so pass serialized XML directly. */
            pCopyDS = reinterpret_cast<GDALDataset *>(
                GDALOpen( pszXML, GA_Update ) );
        }

        CPLFree( pszXML );

        return pCopyDS;
    }

/* -------------------------------------------------------------------- */
/*      Create the virtual dataset.                                     */
/* -------------------------------------------------------------------- */
    VRTDataset *poVRTDS = reinterpret_cast<VRTDataset *>(
        VRTDataset::Create( pszFilename,
                            poSrcDS->GetRasterXSize(),
                            poSrcDS->GetRasterYSize(),
                            0, GDT_Byte, NULL ) );
    if( poVRTDS == NULL )
        return NULL;

/* -------------------------------------------------------------------- */
/*      Do we have a geotransform?                                      */
/* -------------------------------------------------------------------- */
    double adfGeoTransform[6] = { 0.0 };

    if( poSrcDS->GetGeoTransform( adfGeoTransform ) == CE_None )
    {
        poVRTDS->SetGeoTransform( adfGeoTransform );
    }

/* -------------------------------------------------------------------- */
/*      Copy projection                                                 */
/* -------------------------------------------------------------------- */
    poVRTDS->SetProjection( poSrcDS->GetProjectionRef() );

/* -------------------------------------------------------------------- */
/*      Emit dataset level metadata.                                    */
/* -------------------------------------------------------------------- */
    poVRTDS->SetMetadata( poSrcDS->GetMetadata() );
//.........这里部分代码省略.........

//.........这里部分代码省略.........
/* ==================================================================== */
/*      Create a virtual dataset.                                       */
/* ==================================================================== */
    VRTDataset *poVDS;

/* -------------------------------------------------------------------- */
/*      Make a virtual clone.                                           */
/* -------------------------------------------------------------------- */
    poVDS = (VRTDataset*) VRTCreate(nOXSize, nOYSize);

    if (nGCPCount == 0)
    {
        if (pszOutputSRS != NULL)
        {
            poVDS->SetProjection(pszOutputSRS);
        }
        else
        {
            pszProjection = GDALGetProjectionRef(hDataset);
            if (pszProjection != NULL && strlen(pszProjection) > 0)
                poVDS->SetProjection(pszProjection);
        }
    }

    if (bGotBounds)
    {
        adfGeoTransform[0] = adfULLR[0];
        adfGeoTransform[1] = (adfULLR[2] - adfULLR[0]) / nOXSize;
        adfGeoTransform[2] = 0.0;
        adfGeoTransform[3] = adfULLR[1];
        adfGeoTransform[4] = 0.0;
        adfGeoTransform[5] = (adfULLR[3] - adfULLR[1]) / nOYSize;

        poVDS->SetGeoTransform(adfGeoTransform);
    }

    else if (GDALGetGeoTransform(hDataset, adfGeoTransform) == CE_None
             && nGCPCount == 0)
    {
        adfGeoTransform[0] += anSrcWin[0] * adfGeoTransform[1]
                              + anSrcWin[1] * adfGeoTransform[2];
        adfGeoTransform[3] += anSrcWin[0] * adfGeoTransform[4]
                              + anSrcWin[1] * adfGeoTransform[5];

        adfGeoTransform[1] *= anSrcWin[2] / (double) nOXSize;
        adfGeoTransform[2] *= anSrcWin[3] / (double) nOYSize;
        adfGeoTransform[4] *= anSrcWin[2] / (double) nOXSize;
        adfGeoTransform[5] *= anSrcWin[3] / (double) nOYSize;

        poVDS->SetGeoTransform(adfGeoTransform);
    }

    if (nGCPCount != 0)
    {
        const char *pszGCPProjection = pszOutputSRS;

        if (pszGCPProjection == NULL)
            pszGCPProjection = GDALGetGCPProjection(hDataset);

        if (pszGCPProjection == NULL)
            pszGCPProjection = "";

        poVDS->SetGCPs(nGCPCount, pasGCPs, pszGCPProjection);

        GDALDeinitGCPs(nGCPCount, pasGCPs);
        CPLFree(pasGCPs);

//.........这里部分代码省略.........

        exit( 0 );
    }

    if (padfScaleMin == NULL || padfScaleMax == NULL)
    {
        fprintf( stderr, "-equalize or -config filename command line options must be specified.\n");
        exit(1);
    }

/* ==================================================================== */
/*      Create a virtual dataset.                                       */
/* ==================================================================== */
    VRTDataset *poVDS;
    EnhanceCBInfo *pasEInfo = (EnhanceCBInfo *) 
        CPLCalloc(nBandCount, sizeof(EnhanceCBInfo));
        
/* -------------------------------------------------------------------- */
/*      Make a virtual clone.                                           */
/* -------------------------------------------------------------------- */
    poVDS = new VRTDataset( GDALGetRasterXSize(hDataset),
                            GDALGetRasterYSize(hDataset) );

    if( GDALGetGCPCount(hDataset) == 0 )
    {
        const char *pszProjection;
        double adfGeoTransform[6];

        pszProjection = GDALGetProjectionRef( hDataset );
        if( pszProjection != NULL && strlen(pszProjection) > 0 )
            poVDS->SetProjection( pszProjection );

        if( GDALGetGeoTransform( hDataset, adfGeoTransform ) == CE_None )
            poVDS->SetGeoTransform( adfGeoTransform );
    }
    else
    {
        poVDS->SetGCPs( GDALGetGCPCount(hDataset), 
                        GDALGetGCPs(hDataset),
                        GDALGetGCPProjection( hDataset ) );
    }
    
    poVDS->SetMetadata( ((GDALDataset*)hDataset)->GetMetadata() );

    for( iBand = 0; iBand < nBandCount; iBand++ )
    {
        VRTSourcedRasterBand   *poVRTBand;
        GDALRasterBand  *poSrcBand;
        GDALDataType    eBandType;

        poSrcBand = ((GDALDataset *) hDataset)->GetRasterBand(iBand+1);

/* -------------------------------------------------------------------- */
/*      Select output data type to match source.                        */
/* -------------------------------------------------------------------- */
        if( eOutputType == GDT_Unknown )
            eBandType = GDT_Byte;
        else
            eBandType = eOutputType;

/* -------------------------------------------------------------------- */
/*      Create this band.                                               */
/* -------------------------------------------------------------------- */
        poVDS->AddBand( eBandType, NULL );
        poVRTBand = (VRTSourcedRasterBand *) poVDS->GetRasterBand( iBand+1 );
            

bool CreateSubRaster( wxGISRasterDatasetSPtr pSrcRasterDataSet, OGREnvelope &Env, const OGRGeometry *pGeom, GDALDriver* pDriver, CPLString &szDstPath, GDALDataType eOutputType, int nBandCount, int *panBandList, double dfOutResX, double dfOutResY, bool bCopyNodata, bool bSkipSourceMetadata, char** papszOptions, ITrackCancel* pTrackCancel )
{
	GDALDataset* pDset = pSrcRasterDataSet->GetRaster();
	if(!pDset)
	{
		if(pTrackCancel)
			pTrackCancel->PutMessage(_("Get raster failed"), -1, enumGISMessageErr);
		return false;
	}

    double adfGeoTransform[6] = { 0, 0, 0, 0, 0, 0 };
	CPLErr err = pDset->GetGeoTransform(adfGeoTransform);
	if(err == CE_Fatal)
	{
		if(pTrackCancel)
			pTrackCancel->PutMessage(_("Get raster failed"), -1, enumGISMessageErr);
		return false;
	}
	if( adfGeoTransform[2] != 0.0 || adfGeoTransform[4] != 0.0 )
	{
		if(pTrackCancel)
			pTrackCancel->PutMessage(_("The geotransform is rotated. This configuration is not supported."), -1, enumGISMessageErr);
		return false;
    }
	int anSrcWin[4] = {0, 0, 0, 0};

    anSrcWin[0] = floor ((Env.MinX - adfGeoTransform[0]) / adfGeoTransform[1] + 0.001);
    anSrcWin[1] = floor ((Env.MaxY - adfGeoTransform[3]) / adfGeoTransform[5] + 0.001);
	anSrcWin[2] = ceil ((Env.MaxX - Env.MinX) / adfGeoTransform[1]);
	anSrcWin[3] = ceil ((Env.MinY - Env.MaxY) / adfGeoTransform[5]);
	if(pTrackCancel)
		pTrackCancel->PutMessage(wxString::Format(_("Computed source pixel window %d %d %d %d from geographic window."), anSrcWin[0], anSrcWin[1], anSrcWin[2], anSrcWin[3] ), -1, enumGISMessageInfo);

	if( anSrcWin[0] < 0 || anSrcWin[1] < 0 || anSrcWin[0] + anSrcWin[2] > pSrcRasterDataSet->GetWidth() || anSrcWin[1] + anSrcWin[3] > pSrcRasterDataSet->GetHeight() )
    {
		if(pTrackCancel)
			pTrackCancel->PutMessage(wxString::Format(_("Computed source pixel window falls outside raster size of %dx%d."), pSrcRasterDataSet->GetWidth(), pSrcRasterDataSet->GetHeight()), -1, enumGISMessageErr);
		return false;
    }

	int nOXSize = 0, nOYSize = 0;

    if(IsDoubleEquil(dfOutResX, -1) && IsDoubleEquil(dfOutResY, -1))
    {
        nOXSize = anSrcWin[2];
        nOYSize = anSrcWin[3];
    }
    else
    {
        nOXSize = ceil ((Env.MaxX - Env.MinX) / dfOutResX);
        nOYSize = ceil ((Env.MinY - Env.MaxY) / (adfGeoTransform[5] < 0 ? dfOutResY * -1 : dfOutResY));
    }

/* ==================================================================== */
/*      Create a virtual dataset.                                       */
/* ==================================================================== */
    VRTDataset *poVDS;
/* -------------------------------------------------------------------- */
/*      Make a virtual clone.                                           */
/* -------------------------------------------------------------------- */
    poVDS = (VRTDataset *) VRTCreate( nOXSize, nOYSize );

    if( pSrcRasterDataSet->GetSpatialReference() != NULL )
    {
		poVDS->SetProjection( pDset->GetProjectionRef() );
    }

	adfGeoTransform[0] += anSrcWin[0] * adfGeoTransform[1] + anSrcWin[1] * adfGeoTransform[2];
    adfGeoTransform[3] += anSrcWin[0] * adfGeoTransform[4] + anSrcWin[1] * adfGeoTransform[5];

    adfGeoTransform[1] *= anSrcWin[2] / (double) nOXSize;
    adfGeoTransform[2] *= anSrcWin[3] / (double) nOYSize;
    adfGeoTransform[4] *= anSrcWin[2] / (double) nOXSize;
    adfGeoTransform[5] *= anSrcWin[3] / (double) nOYSize;

    poVDS->SetGeoTransform( adfGeoTransform );

    int nGCPs = pDset->GetGCPCount();
    if( nGCPs > 0 )
    {
        GDAL_GCP *pasGCPs = GDALDuplicateGCPs( nGCPs, pDset->GetGCPs() );

        for(size_t i = 0; i < nGCPs; ++i )
        {
            pasGCPs[i].dfGCPPixel -= anSrcWin[0];
            pasGCPs[i].dfGCPLine  -= anSrcWin[1];
            pasGCPs[i].dfGCPPixel *= (nOXSize / (double) anSrcWin[2] );
            pasGCPs[i].dfGCPLine  *= (nOYSize / (double) anSrcWin[3] );
        }

        poVDS->SetGCPs( nGCPs, pasGCPs, pDset->GetGCPProjection() );
        GDALDeinitGCPs( nGCPs, pasGCPs );
        CPLFree( pasGCPs );
    }

/* -------------------------------------------------------------------- */
/*      Transfer generally applicable metadata.                         */
/* -------------------------------------------------------------------- */
    if(!bSkipSourceMetadata)
        poVDS->SetMetadata( pDset->GetMetadata() );
//.........这里部分代码省略.........