C++ GDALOpen函数代码示例

// -------------------------------------------------------------------
int main(int argc, char **argv)
{
	
   char *ERROR_PRM = "Ungueltiger numerischer Wert fuer %s: %s\n!";
	
   // Alle GDAL Treiber registrieren   
   GDALAllRegister();
	
   // Dateiname der Geotiff-Datei     
   const char *format = "GTiff";
   GDALDriverH h_drv = GDALGetDriverByName( format );
   if( h_drv == NULL ) {
	   error_exit(10,"Treiber %s nicht vorhanden!" ,format);
   }

   // Test ob Geotiffdateien erzeugt werden koennen
   char **test_meta;
   test_meta = GDALGetMetadata( h_drv, NULL );
   if( ! CSLFetchBoolean( test_meta, GDAL_DCAP_CREATE, FALSE ) ) {
	   error_exit(10,"Das Format %s kann nicht erzeugt werden" ,format);
	}
      
   // 3 Kommandozeilenparameter
   if (argc<6) {
	   error_exit(10,
	   "Fehlende Parameter\nUsage %s IN OUT EXT SZ ID X Y ID X Y ID X Y....!\n",
	   argv[0]);
   }
   
   // Eingabemuster einlesen
   char *ifile = argv[1];

   // Ausgabemuster einlesen
   char *ofile = argv[2];

   // zusammengesetzte Ausgabedatei 
   char cfile[512];

   // Dateierweiterung setzen
   char *ext   = argv[3];

   // Fenstergroesse
   int size = 64; 
   if (! sscanf(argv[4],"%d",&size) ) {
        error_exit(1000+3,ERROR_PRM,"SZ",argv[4]);
   }
   double trfm[] ={0,0,0,0,0,0};

   // Vektoren fuer die Positionen und ID
   int_vector_t id;
   int_vector_init(&id, 10);
   dbl_vector_t pos_x;
   dbl_vector_init(&pos_x, 10);
   dbl_vector_t pos_y;
   dbl_vector_init(&pos_y, 10);
   
   // Positionen einlesen 
   int a = 5; double dbl; int pk;
   while( a < argc-2 ) {

	   // X Koordinate parsen
       if (! sscanf(argv[a],"%d",&pk) ) {
          error_exit(1000+a,ERROR_PRM,"ID", argv[a]);
       }
	   int_vector_add(&id,pk);
       
	   // X Koordinate parsen
       if (! sscanf(argv[a+1],"%lf",&dbl) ) {
          error_exit(1000+a+1,ERROR_PRM,"X", argv[a+1]);
       }
	   dbl_vector_add(&pos_x,dbl);

	   // Y Koordinate parsen
       if (! sscanf(argv[a+2],"%lf",&dbl) ) {
          error_exit(1000+a+2,ERROR_PRM,"Y", argv[a+2]);
       }
	   dbl_vector_add(&pos_y,dbl);

	   a+=3;
   }
      
   // Alle GDAL Treiber registrieren   
   GDALAllRegister();
   
   // Geotiff oeffnen
   printf("# IN FILE:  %s\n", ifile);
   GDALDatasetH h_dset = GDALOpen( ifile, GA_ReadOnly);

   printf("# OUT FILE: %s\n", ofile);

   printf("# EXTENTION: %s\n", ext);
   
   // Transformation holen 
   if( GDALGetGeoTransform( h_dset, trfm ) == CE_None ) {
        printf("# TRANSFORM: \n");
        printf("#  X = %.6f + %.6f * COL + %.6f * ROW\n",
                   trfm[0], trfm[1], trfm[2] );
        printf("#  Y = %.6f + %.6f * COL + %.6f * ROW\n# EOF:\n",
                   trfm[3], trfm[4], trfm[5] );
//.........这里部分代码省略.........

//.........这里部分代码省略.........
                                     char ** papszOptions )

{
    if( eType != GDT_Float32 )
    {
        CPLError(CE_Failure, CPLE_AppDefined,
                 "Attempt to create CTable2 file with unsupported data type '%s'.",
                 GDALGetDataTypeName( eType ) );
        return NULL;
    }

/* -------------------------------------------------------------------- */
/*      Try to open or create file.                                     */
/* -------------------------------------------------------------------- */
    VSILFILE	*fp;

    fp = VSIFOpenL( pszFilename, "wb" );
    
    if( fp == NULL )
    {
        CPLError( CE_Failure, CPLE_OpenFailed,
                  "Attempt to create file `%s' failed.\n",
                  pszFilename );
        return NULL;
    }

/* -------------------------------------------------------------------- */
/*      Create a file header, with a defaulted georeferencing.          */
/* -------------------------------------------------------------------- */
    char achHeader[160];
    int nValue32;
    double dfValue;

    memset( achHeader, 0, sizeof(achHeader));

    memcpy( achHeader+0, "CTABLE V2.0     ", 16 );
    
    if( CSLFetchNameValue( papszOptions, "DESCRIPTION" ) != NULL )
        strncpy( achHeader + 16, 
                 CSLFetchNameValue( papszOptions, "DESCRIPTION" ), 
                 80 );
    
    // lower left origin (longitude, center of pixel, radians)
    dfValue = 0;
    CPL_LSBPTR64( &dfValue );
    memcpy( achHeader + 96, &dfValue, 8 );

    // lower left origin (latitude, center of pixel, radians)
    dfValue = 0;
    CPL_LSBPTR64( &dfValue );
    memcpy( achHeader + 104, &dfValue, 8 );

    // pixel width (radians)
    dfValue = 0.01 * M_PI / 180.0;
    CPL_LSBPTR64( &dfValue );
    memcpy( achHeader + 112, &dfValue, 8 );

    // pixel height (radians)
    dfValue = 0.01 * M_PI / 180.0;
    CPL_LSBPTR64( &dfValue );
    memcpy( achHeader + 120, &dfValue, 8 );

    // raster width in pixels
    nValue32 = nXSize;
    CPL_LSBPTR32( &nValue32 );
    memcpy( achHeader + 128, &nValue32, 4 );

    // raster width in pixels
    nValue32 = nYSize;
    CPL_LSBPTR32( &nValue32 );
    memcpy( achHeader + 132, &nValue32, 4 );

    VSIFWriteL( achHeader, 1, sizeof(achHeader), fp );

/* -------------------------------------------------------------------- */
/*      Write zeroed grid data.                                         */
/* -------------------------------------------------------------------- */
    float *pafLine = (float *) CPLCalloc(sizeof(float)*2,nXSize);
    int i;

    for( i = 0; i < nYSize; i++ )
    {
        if( (int)VSIFWriteL( pafLine, sizeof(float)*2, nXSize, fp ) != nXSize ) 
        {
            CPLError( CE_Failure, CPLE_FileIO, 
                      "Write failed at line %d, perhaps the disk is full?",
                      i );
            return NULL;
        }
    }
    
/* -------------------------------------------------------------------- */
/*      Cleanup and return.                                             */
/* -------------------------------------------------------------------- */
    CPLFree( pafLine );

    VSIFCloseL( fp );

    return (GDALDataset *) GDALOpen( pszFilename, GA_Update );
}

int main( int argc, char ** argv )

{
    /* Check that we are running against at least GDAL 1.4 (probably older in fact !) */
    /* Note to developers : if we use newer API, please change the requirement */
    if (atoi(GDALVersionInfo("VERSION_NUM")) < 1400)
    {
        fprintf(stderr, "At least, GDAL >= 1.4.0 is required for this version of %s, "
                "which was compiled against GDAL %s\n", argv[0], GDAL_RELEASE_NAME);
        exit(1);
    }

/* -------------------------------------------------------------------- */
/*      Generic arg processing.                                         */
/* -------------------------------------------------------------------- */
    GDALAllRegister();
    GDALSetCacheMax( 100000000 );
    argc = GDALGeneralCmdLineProcessor( argc, &argv, 0 );
    if( argc < 1 )
        exit( -argc );
    
/* -------------------------------------------------------------------- */
/*      Parse arguments.                                                */
/* -------------------------------------------------------------------- */
    int i;
    const char *pszOutFile = NULL;
    const char *pszInFile = NULL;
    int nMaxNonBlack = 2;
    int nNearDist = 15;
    int bNearWhite = FALSE;
    
    for( i = 1; i < argc; i++ )
    {
        if( EQUAL(argv[i], "--utility_version") )
        {
            printf("%s was compiled against GDAL %s and is running against GDAL %s\n",
                   argv[0], GDAL_RELEASE_NAME, GDALVersionInfo("RELEASE_NAME"));
            return 0;
        }
        else if( EQUAL(argv[i], "-o") && i < argc-1 )
            pszOutFile = argv[++i];
        else if( EQUAL(argv[i], "-white") )
            bNearWhite = TRUE;
        else if( EQUAL(argv[i], "-nb") && i < argc-1 )
            nMaxNonBlack = atoi(argv[++i]);
        else if( EQUAL(argv[i], "-near") && i < argc-1 )
            nNearDist = atoi(argv[++i]);
        else if( argv[i][0] == '-' )
            Usage();
        else if( pszInFile == NULL )
            pszInFile = argv[i];
        else
            Usage();
    }

    if( pszInFile == NULL )
        Usage();

    if( pszOutFile == NULL )
        pszOutFile = pszInFile;

/* -------------------------------------------------------------------- */
/*      Open input file.                                                */
/* -------------------------------------------------------------------- */
    GDALDatasetH hInDS, hOutDS = NULL;
    int nXSize, nYSize, nBands;

    if( pszOutFile == pszInFile )
        hInDS = hOutDS = GDALOpen( pszInFile, GA_Update );
    else
        hInDS = GDALOpen( pszInFile, GA_ReadOnly );

    if( hInDS == NULL )
        exit( 1 );

    nXSize = GDALGetRasterXSize( hInDS );
    nYSize = GDALGetRasterYSize( hInDS );
    nBands = GDALGetRasterCount( hInDS );

/* -------------------------------------------------------------------- */
/*      Do we need to create output file?                               */
/* -------------------------------------------------------------------- */
    if( hOutDS == NULL )
    {
        GDALDriverH hDriver = GDALGetDriverByName( "HFA" );
        
        hOutDS = GDALCreate( hDriver, pszOutFile, 
                             nXSize, nYSize, nBands, GDT_Byte, 
                             NULL );
        if( hOutDS == NULL )
            exit( 1 );

        double adfGeoTransform[6];

        if( GDALGetGeoTransform( hInDS, adfGeoTransform ) == CE_None )
        {
            GDALSetGeoTransform( hOutDS, adfGeoTransform );
            GDALSetProjection( hOutDS, GDALGetProjectionRef( hInDS ) );
        }
    }
//.........这里部分代码省略.........

//.........这里部分代码省略.........
    double dfMinY = adfGeoTransform[5] * (nYSize - 0.5) + adfGeoTransform[3];
    double dfMaxY = adfGeoTransform[3] + adfGeoTransform[5] / 2;
    CPLErr eErr = WriteHeader( fp, nXSize, nYSize,
			       dfMinX, dfMaxX, dfMinY, dfMaxY, 0.0, 0.0 );

    if( eErr != CE_None )
    {
	VSIFCloseL( fp );
        return NULL;
    }

/* -------------------------------------------------------------------- */
/*      Copy band data.							*/
/* -------------------------------------------------------------------- */
    float *pfData = (float *)VSIMalloc2( nXSize, sizeof( float ) );
    if( pfData == NULL )
    {
	VSIFCloseL( fp );
	CPLError( CE_Failure, CPLE_OutOfMemory,
		  "Unable to create copy, unable to allocate line buffer.\n" );
	return NULL;
    }

    int     bSrcHasNDValue;
    float   fSrcNoDataValue = (float) poSrcBand->GetNoDataValue( &bSrcHasNDValue );
    double  dfMinZ = DBL_MAX;
    double  dfMaxZ = -DBL_MAX;
    for( GInt16 iRow = nYSize - 1; iRow >= 0; iRow-- )
    {
	eErr = poSrcBand->RasterIO( GF_Read, 0, iRow,
				    nXSize, 1, pfData,
				    nXSize, 1, GDT_Float32, 0, 0 );

	if( eErr != CE_None )
	{
	    VSIFCloseL( fp );
	    VSIFree( pfData );
	    return NULL;
	}

	for( int iCol=0; iCol<nXSize; iCol++ )
	{
	    if( bSrcHasNDValue && pfData[iCol] == fSrcNoDataValue )
	    {
		pfData[iCol] = fNODATA_VALUE;
	    }
	    else
	    {
		if( pfData[iCol] > dfMaxZ )
		    dfMaxZ = pfData[iCol];

		if( pfData[iCol] < dfMinZ )
		    dfMinZ = pfData[iCol];
	    }

	    CPL_LSBPTR32( pfData+iCol );
	}

	if( VSIFWriteL( (void *)pfData, 4, nXSize,
			fp ) != static_cast<unsigned>(nXSize) )
	{
	    VSIFCloseL( fp );
	    VSIFree( pfData );
	    CPLError( CE_Failure, CPLE_FileIO,
		      "Unable to write grid row. Disk full?\n" );
	    return NULL;
	}

	if( !pfnProgress( static_cast<double>(nYSize - iRow)/nYSize,
			  NULL, pProgressData ) )
	{
	    VSIFCloseL( fp );
	    VSIFree( pfData );
	    CPLError( CE_Failure, CPLE_UserInterrupt, "User terminated" );
	    return NULL;
	}
    }

    VSIFree( pfData );

    /* write out the min and max values */
    eErr = WriteHeader( fp, nXSize, nYSize,
			dfMinX, dfMaxX, dfMinY, dfMaxY, dfMinZ, dfMaxZ );

    if( eErr != CE_None )
    {
	VSIFCloseL( fp );
        return NULL;
    }

    VSIFCloseL( fp );

    GDALPamDataset *poDS = (GDALPamDataset *)GDALOpen( pszFilename,
                                                GA_Update );
    if (poDS)
    {
        poDS->CloneInfo( poSrcDS, GCIF_PAM_DEFAULT );
    }
    return poDS;
}

//.........这里部分代码省略.........
            {
                char** papszMEMDSOptions = NULL;
                char szTmp[64];
                memset(szTmp, 0, sizeof(szTmp));
                CPLPrintPointer(szTmp,
                                pabyMEMDSBuffer + iBand * nDataTypeSize *
                                nReqXSize * nReqYSize, sizeof(szTmp));
                papszMEMDSOptions = CSLSetNameValue(papszMEMDSOptions, "DATAPOINTER", szTmp);
                GDALAddBand(hMemDS, eDataType, papszMEMDSOptions);
                CSLDestroy(papszMEMDSOptions);
            }
            
            GDALDatasetH hOutDS = GDALCreateCopy(hTileDriver,
                                        osTempFileName.c_str(), hMemDS, FALSE,
                                        papszTileDriverOptions, NULL, NULL);

            GDALClose(hMemDS);
            if (hOutDS)
                GDALClose(hOutDS);
            else
            {
                eErr = CE_Failure;
                break;
            }

/* -------------------------------------------------------------------- */
/*      Insert new entry into raster table                              */
/* -------------------------------------------------------------------- */

            vsi_l_offset nDataLength;
            GByte *pabyData = VSIGetMemFileBuffer( osTempFileName.c_str(),
                                                   &nDataLength, FALSE);

            OGRFeatureH hFeat = OGR_F_Create( OGR_L_GetLayerDefn(hRasterLayer) );
            OGR_F_SetFieldBinary(hFeat, 0, (int)nDataLength, pabyData);
            
            OGR_L_CreateFeature(hRasterLayer, hFeat);
            /* Query raster ID to set it as the ID of the associated metadata */
            int nRasterID = (int)OGR_F_GetFID(hFeat);
            
            OGR_F_Destroy(hFeat);
            
            VSIUnlink(osTempFileName.c_str());
            
/* -------------------------------------------------------------------- */
/*      Insert new entry into metadata table                            */
/* -------------------------------------------------------------------- */
            
            hFeat = OGR_F_Create( OGR_L_GetLayerDefn(hMetadataLayer) );
            OGR_F_SetFID(hFeat, nRasterID);
            OGR_F_SetFieldString(hFeat, 0, GDALGetDescription(poSrcDS));
            OGR_F_SetFieldInteger(hFeat, 1, nTileId ++);
            OGR_F_SetFieldInteger(hFeat, 2, nReqXSize);
            OGR_F_SetFieldInteger(hFeat, 3, nReqYSize);
            OGR_F_SetFieldDouble(hFeat, 4, adfGeoTransform[1]);
            OGR_F_SetFieldDouble(hFeat, 5, -adfGeoTransform[5]);
            
            minx = adfGeoTransform[0] +
                (nBlockXSize * nBlockXOff) * adfGeoTransform[1];
            maxx = adfGeoTransform[0] +
                (nBlockXSize * nBlockXOff + nReqXSize) * adfGeoTransform[1];
            maxy = adfGeoTransform[3] +
                (nBlockYSize * nBlockYOff) * adfGeoTransform[5];
            miny = adfGeoTransform[3] +
                (nBlockYSize * nBlockYOff + nReqYSize) * adfGeoTransform[5];
            
            OGRGeometryH hRectangle = OGR_G_CreateGeometry(wkbPolygon);
            OGRGeometryH hLinearRing = OGR_G_CreateGeometry(wkbLinearRing);
            OGR_G_AddPoint_2D(hLinearRing, minx, miny);
            OGR_G_AddPoint_2D(hLinearRing, minx, maxy);
            OGR_G_AddPoint_2D(hLinearRing, maxx, maxy);
            OGR_G_AddPoint_2D(hLinearRing, maxx, miny);
            OGR_G_AddPoint_2D(hLinearRing, minx, miny);
            OGR_G_AddGeometryDirectly(hRectangle, hLinearRing);
            
            OGR_F_SetGeometryDirectly(hFeat, hRectangle);
            
            OGR_L_CreateFeature(hMetadataLayer, hFeat);
            OGR_F_Destroy(hFeat);
            
            nBlocks++;
            if (pfnProgress && !pfnProgress(1.0 * nBlocks / nTotalBlocks,
                                            NULL, pProgressData))
                eErr = CE_Failure;
        }
    }
    
    if (eErr == CE_None)
        OGR_DS_ExecuteSQL(hDS, "COMMIT", NULL, NULL);
    else
        OGR_DS_ExecuteSQL(hDS, "ROLLBACK", NULL, NULL);
    
    CSLDestroy(papszTileDriverOptions);
    
    VSIFree(pabyMEMDSBuffer);
    
    OGRReleaseDataSource(hDS);
        
    return (GDALDataset*) GDALOpen(pszFilename, GA_Update);
}

//.........这里部分代码省略.........
		    if( dfValue > dfMax )
			dfMax = dfValue;

		    if( dfValue < dfMin )
			dfMin = dfValue;
		}

		std::ostringstream ssOut;
		ssOut.precision(nFIELD_PRECISION);
		ssOut.setf( std::ios::uppercase );
		ssOut << dfValue << " ";
		CPLString sOut = ssOut.str();

		if( VSIFWriteL( sOut.c_str(), 1, sOut.length(), fp )
		    != sOut.length() )
		{
		    VSIFCloseL( fp );
		    VSIFree( pdfData );
		    CPLError( CE_Failure, CPLE_FileIO,
			      "Unable to write grid cell.  Disk full?\n" );
		    return NULL;
		}
	    }

	    if( VSIFWriteL( (void *)"\x0D\x0A", 1, 2, fp ) != 2 )
	    {
		VSIFCloseL( fp );
		VSIFree( pdfData );
		CPLError( CE_Failure, CPLE_FileIO,
			  "Unable to finish write of grid line. Disk full?\n" );
		return NULL;
	    }
	}

	if( VSIFWriteL( (void *)"\x0D\x0A", 1, 2, fp ) != 2 )
	{
	    VSIFCloseL( fp );
	    VSIFree( pdfData );
	    CPLError( CE_Failure, CPLE_FileIO,
		      "Unable to finish write of grid row. Disk full?\n" );
	    return NULL;
	}

	if( !pfnProgress( static_cast<double>(iRow + 1)/nYSize,
			  NULL, pProgressData ) )
	{
	    VSIFCloseL( fp );
	    VSIFree( pdfData );
	    CPLError( CE_Failure, CPLE_UserInterrupt, "User terminated" );
	    return NULL;
	}
    }

    VSIFree( pdfData );

    /* write out the min and max values */
    std::ostringstream ssRange;
    ssRange.precision( nFIELD_PRECISION );
    ssRange.setf( std::ios::uppercase );
    ssRange << dfMin << " " << dfMax << "\x0D\x0A";
    if( ssRange.str().length() != nDummyRangeLen )
    {
	int nShiftSize = ssRange.str().length() - nDummyRangeLen;
	if( ShiftFileContents( fp, nRangeStart + nDummyRangeLen,
			       nShiftSize, "\x0D\x0A" ) != CE_None )
	{
	    VSIFCloseL( fp );
	    CPLError( CE_Failure, CPLE_FileIO,
		      "Unable to shift file contents.\n" );
	    return NULL;
	}
    }

    if( VSIFSeekL( fp, nRangeStart, SEEK_SET ) != 0 )
    {
	VSIFCloseL( fp );
	CPLError( CE_Failure, CPLE_FileIO,
		  "Unable to seek to start of grid file copy.\n" );
	return NULL;
    }

    if( VSIFWriteL( (void *)ssRange.str().c_str(), 1, ssRange.str().length(),
		    fp ) != ssRange.str().length() )
    {
	VSIFCloseL( fp );
        CPLError( CE_Failure, CPLE_FileIO,
                  "Unable to write range information.\n" );
        return NULL;
    }

    VSIFCloseL( fp );

    GDALPamDataset *poDS = (GDALPamDataset *)GDALOpen( pszFilename,
                                                GA_Update );
    if (poDS)
    {
        poDS->CloneInfo( poSrcDS, GCIF_PAM_DEFAULT );
    }
    return poDS;
}

int main( int nArgc, char ** papszArgv )

{
    GDALDatasetH     hDataset;
    const char      *pszResampling = "nearest";
    const char      *pszFilename = NULL;
    int              anLevels[1024];
    int              nLevelCount = 0;
    int              nResultStatus = 0;
    int              bReadOnly = FALSE;
    int              bClean = FALSE;
    GDALProgressFunc pfnProgress = GDALTermProgress; 

    /* Check that we are running against at least GDAL 1.7 */
    /* Note to developers : if we use newer API, please change the requirement */
    if (atoi(GDALVersionInfo("VERSION_NUM")) < 1700)
    {
        fprintf(stderr, "At least, GDAL >= 1.7.0 is required for this version of %s, "
                        "which was compiled against GDAL %s\n", papszArgv[0], GDAL_RELEASE_NAME);
        exit(1);
    }

    GDALAllRegister();

    nArgc = GDALGeneralCmdLineProcessor( nArgc, &papszArgv, 0 );
    if( nArgc < 1 )
        exit( -nArgc );

/* -------------------------------------------------------------------- */
/*      Parse commandline.                                              */
/* -------------------------------------------------------------------- */
    for( int iArg = 1; iArg < nArgc; iArg++ )
    {
        if( EQUAL(papszArgv[iArg], "--utility_version") )
        {
            printf("%s was compiled against GDAL %s and is running against GDAL %s\n",
                   papszArgv[0], GDAL_RELEASE_NAME, GDALVersionInfo("RELEASE_NAME"));
            return 0;
        }
        else if( EQUAL(papszArgv[iArg],"-r") && iArg < nArgc-1 )
            pszResampling = papszArgv[++iArg];
        else if( EQUAL(papszArgv[iArg],"-ro"))
            bReadOnly = TRUE;
        else if( EQUAL(papszArgv[iArg],"-clean"))
            bClean = TRUE;
        else if( EQUAL(papszArgv[iArg],"-q") || EQUAL(papszArgv[iArg],"-quiet") ) 
            pfnProgress = GDALDummyProgress; 
        else if( pszFilename == NULL )
            pszFilename = papszArgv[iArg];
        else if( atoi(papszArgv[iArg]) > 0 )
            anLevels[nLevelCount++] = atoi(papszArgv[iArg]);
        else
            Usage();
    }

    if( pszFilename == NULL || (nLevelCount == 0 && !bClean) )
        Usage();

/* -------------------------------------------------------------------- */
/*      Open data file.                                                 */
/* -------------------------------------------------------------------- */
    if (bReadOnly)
        hDataset = NULL;
    else
    {
        CPLPushErrorHandler( CPLQuietErrorHandler );
        hDataset = GDALOpen( pszFilename, GA_Update );
        CPLPopErrorHandler();
    }

    if( hDataset == NULL )
        hDataset = GDALOpen( pszFilename, GA_ReadOnly );

    if( hDataset == NULL )
        exit( 2 );

/* -------------------------------------------------------------------- */
/*      Clean overviews.                                                */
/* -------------------------------------------------------------------- */
    if ( bClean &&
        GDALBuildOverviews( hDataset,pszResampling, 0, 0, 
                             0, NULL, pfnProgress, NULL ) != CE_None )
    {
        printf( "Cleaning overviews failed.\n" );
        nResultStatus = 200;
    }

/* -------------------------------------------------------------------- */
/*      Generate overviews.                                             */
/* -------------------------------------------------------------------- */
    if (nLevelCount > 0 && nResultStatus == 0 &&
        GDALBuildOverviews( hDataset,pszResampling, nLevelCount, anLevels,
                             0, NULL, pfnProgress, NULL ) != CE_None )
    {
        printf( "Overview building failed.\n" );
        nResultStatus = 100;
    }

/* -------------------------------------------------------------------- */
/*      Cleanup                                                         */
//.........这里部分代码省略.........

CPLErr GDALWMSRasterBand::ReadBlockFromFile(int x, int y, const char *file_name, int to_buffer_band, void *buffer, int advise_read) {
    CPLErr ret = CE_None;
    GDALDataset *ds = 0;
    GByte *color_table = NULL;
    int i;

    //CPLDebug("WMS", "ReadBlockFromFile: to_buffer_band=%d, (x,y)=(%d, %d)", to_buffer_band, x, y);

    /* expected size */
    const int esx = MIN(MAX(0, (x + 1) * nBlockXSize), nRasterXSize) - MIN(MAX(0, x * nBlockXSize), nRasterXSize);
    const int esy = MIN(MAX(0, (y + 1) * nBlockYSize), nRasterYSize) - MIN(MAX(0, y * nBlockYSize), nRasterYSize);
    ds = reinterpret_cast<GDALDataset*>(GDALOpen(file_name, GA_ReadOnly));
    if (ds != NULL) {
        int sx = ds->GetRasterXSize();
        int sy = ds->GetRasterYSize();
        bool accepted_as_no_alpha = false;  // if the request is for 4 bands but the wms returns 3  
        /* Allow bigger than expected so pre-tiled constant size images work on corners */
        if ((sx > nBlockXSize) || (sy > nBlockYSize) || (sx < esx) || (sy < esy)) {
            CPLError(CE_Failure, CPLE_AppDefined, "GDALWMS: Incorrect size %d x %d of downloaded block, expected %d x %d, max %d x %d.",
                sx, sy, esx, esy, nBlockXSize, nBlockYSize);
            ret = CE_Failure;
        }
        if (ret == CE_None) {
            int nDSRasterCount = ds->GetRasterCount();
            if (nDSRasterCount != m_parent_dataset->nBands) {
                /* Maybe its an image with color table */
                bool accepted_as_ct = false;
                if ((eDataType == GDT_Byte) && (ds->GetRasterCount() == 1)) {
                    GDALRasterBand *rb = ds->GetRasterBand(1);
                    if (rb->GetRasterDataType() == GDT_Byte) {
                        GDALColorTable *ct = rb->GetColorTable();
                        if (ct != NULL) {
                            accepted_as_ct = true;
                            if (!advise_read) {
                                color_table = new GByte[256 * 4];
                                const int count = MIN(256, ct->GetColorEntryCount());
                                for (i = 0; i < count; ++i) {
                                    GDALColorEntry ce;
                                    ct->GetColorEntryAsRGB(i, &ce);
                                    color_table[i] = static_cast<GByte>(ce.c1);
                                    color_table[i + 256] = static_cast<GByte>(ce.c2);
                                    color_table[i + 512] = static_cast<GByte>(ce.c3);
                                    color_table[i + 768] = static_cast<GByte>(ce.c4);
                                }
                                for (i = count; i < 256; ++i) {
                                    color_table[i] = 0;
                                    color_table[i + 256] = 0;
                                    color_table[i + 512] = 0;
                                    color_table[i + 768] = 0;
                                }
                            }
                        }
                    }
                }

                if (nDSRasterCount == 4 && m_parent_dataset->nBands == 3)
                {
                    /* metacarta TMS service sometimes return a 4 band PNG instead of the expected 3 band... */
                }
                else if (!accepted_as_ct) {
                   if (ds->GetRasterCount()==3 && m_parent_dataset->nBands == 4 && (eDataType == GDT_Byte))
                   { // WMS returned a file with no alpha so we will fill the alpha band with "opaque" 
                      accepted_as_no_alpha = true;
                   }
                   else
                   {
                      CPLError(CE_Failure, CPLE_AppDefined, "GDALWMS: Incorrect bands count %d in downloaded block, expected %d.",
                         nDSRasterCount, m_parent_dataset->nBands);
                      ret = CE_Failure;
                   }
                }
            }
        }
        if (!advise_read) {
            for (int ib = 1; ib <= m_parent_dataset->nBands; ++ib) {
                if (ret == CE_None) {
                    void *p = NULL;
                    GDALRasterBlock *b = NULL;
                    if ((buffer != NULL) && (ib == to_buffer_band)) {
                        p = buffer;
                    } else {
                        GDALWMSRasterBand *band = static_cast<GDALWMSRasterBand *>(m_parent_dataset->GetRasterBand(ib));
                        if (m_overview >= 0) band = static_cast<GDALWMSRasterBand *>(band->GetOverview(m_overview));
                        if (!band->IsBlockInCache(x, y)) {
                            b = band->GetLockedBlockRef(x, y, true);
                            if (b != NULL) {
                                p = b->GetDataRef();
                                if (p == NULL) {
                                  CPLError(CE_Failure, CPLE_AppDefined, "GDALWMS: GetDataRef returned NULL.");
                                  ret = CE_Failure;
                                }
                            }
                        }
                        else
                        {
                            //CPLDebug("WMS", "Band %d, block (x,y)=(%d, %d) already in cache", band->GetBand(), x, y);
                        }
                    }
                    if (p != NULL) {
                        int pixel_space = GDALGetDataTypeSize(eDataType) / 8;
//.........这里部分代码省略.........

//		GDALRegister_ADRG();
		GDALRegister_COASP();
		GDALRegister_BLX();
		GDALRegister_LCP();
//		GDALRegister_PGCHIP();
//		GDALRegister_TMS();
		GDALRegister_EIR();
//		GDALRegister_GEOR();

        MemoryIgnoreLeaksEnd();
    }
 
	Close();
	
	MemoryIgnoreLeaksBegin();
	dataset = (GDALDataset *) GDALOpen(fn, GA_ReadOnly);
	MemoryIgnoreLeaksEnd();
	if(!dataset)
		return false;
	int nbands = dataset->GetRasterCount();
	for(int i = 1; i <= nbands; i++)
		bands.Add(new Band(*this, dataset->GetRasterBand(i)));
	double affine[6];
	dataset->GetGeoTransform(affine);
	transform.a.x = affine[0];
	transform.x.x = affine[1];
	transform.y.x = affine[2];
	transform.a.y = affine[3];
	transform.x.y = affine[4];
	transform.y.y = affine[5];
	if(fabs(Determinant(transform)) <= 1e-10)

static void msTransformToGeospatialPDF(imageObj *img, mapObj *map, cairo_renderer *r)
{
  /* We need a GDAL 1.10 PDF driver at runtime, but as far as the C API is concerned, GDAL 1.9 is */
  /* largely sufficient. */
#if defined(USE_GDAL) && defined(GDAL_VERSION_NUM) && GDAL_VERSION_NUM >= 1900
  GDALDatasetH hDS = NULL;
  const char* pszGEO_ENCODING = NULL;
  GDALDriverH hPDFDriver = NULL;
  const char* pszVirtualIO = NULL;
  int bVirtualIO = FALSE;
  char* pszTmpFilename = NULL;
  VSILFILE* fp = NULL;

  if (map == NULL)
    return;

  pszGEO_ENCODING = msGetOutputFormatOption(img->format, "GEO_ENCODING", NULL);
  if (pszGEO_ENCODING == NULL)
    return;

  msGDALInitialize();

  hPDFDriver = GDALGetDriverByName("PDF");
  if (hPDFDriver == NULL)
    return;

  /* When compiled against libpoppler, the PDF driver is VirtualIO capable */
  /* but not, when it is compiled against libpodofo. */
  pszVirtualIO = GDALGetMetadataItem( hPDFDriver, GDAL_DCAP_VIRTUALIO, NULL );
  if (pszVirtualIO)
    bVirtualIO = CSLTestBoolean(pszVirtualIO);

  if (bVirtualIO)
    pszTmpFilename = msTmpFile(map, NULL, "/vsimem/mscairopdf/", "pdf");
  else
    pszTmpFilename = msTmpFile(map, map->mappath, NULL, "pdf");

  /* Copy content of outputStream buffer into file */
  fp = VSIFOpenL(pszTmpFilename, "wb");
  if (fp == NULL) {
    msFree(pszTmpFilename);
    return;
  }
  VSIFWriteL(r->outputStream->data, 1, r->outputStream->size, fp);
  VSIFCloseL(fp);
  fp = NULL;

  hDS = GDALOpen(pszTmpFilename, GA_Update);
  if ( hDS != NULL ) {
    char* pszWKT = msProjectionObj2OGCWKT( &(map->projection) );
    if( pszWKT != NULL ) {
      double adfGeoTransform[6];
      int i;

      /* Add user-specified options */
      for( i = 0; i < img->format->numformatoptions; i++ ) {
        const char* pszOption = img->format->formatoptions[i];
        if( strncasecmp(pszOption,"METADATA_ITEM:",14) == 0 ) {
          char* pszKey = NULL;
          const char* pszValue = CPLParseNameValue(pszOption + 14,
                                 &pszKey);
          if( pszKey != NULL ) {
            GDALSetMetadataItem(hDS, pszKey, pszValue, NULL);
            CPLFree(pszKey);
          }
        }
      }

      /* We need to rescale the geotransform because GDAL will not necessary */
      /* open the PDF with the DPI that was used to generate it */
      memcpy(adfGeoTransform, map->gt.geotransform, 6 * sizeof(double));
      adfGeoTransform[1] = adfGeoTransform[1] * map->width / GDALGetRasterXSize(hDS);
      adfGeoTransform[5] = adfGeoTransform[5] * map->height / GDALGetRasterYSize(hDS);
      GDALSetGeoTransform(hDS, adfGeoTransform);
      GDALSetProjection(hDS, pszWKT);

      msFree( pszWKT );
      pszWKT = NULL;

      CPLSetThreadLocalConfigOption("GDAL_PDF_GEO_ENCODING", pszGEO_ENCODING);

      GDALClose(hDS);
      hDS = NULL;

      CPLSetThreadLocalConfigOption("GDAL_PDF_GEO_ENCODING", NULL);

      /* We need to replace the buffer with the content of the GDAL file */
      fp = VSIFOpenL(pszTmpFilename, "rb");
      if( fp != NULL ) {
        int nFileSize;

        VSIFSeekL(fp, 0, SEEK_END);
        nFileSize = (int)VSIFTellL(fp);

        msBufferResize(r->outputStream, nFileSize);

        VSIFSeekL(fp, 0, SEEK_SET);
        r->outputStream->size = VSIFReadL(r->outputStream->data, 1, nFileSize, fp);

        VSIFCloseL(fp);
//.........这里部分代码省略.........

GDALDataset *PNMDataset::Create( const char * pszFilename,
                                 int nXSize, int nYSize, int nBands,
                                 GDALDataType eType,
                                 char ** papszOptions )

{
/* -------------------------------------------------------------------- */
/*      Verify input options.                                           */
/* -------------------------------------------------------------------- */
    if( eType != GDT_Byte && eType != GDT_UInt16 )
    {
        CPLError( CE_Failure, CPLE_AppDefined,
              "Attempt to create PNM dataset with an illegal\n"
              "data type (%s), only Byte and UInt16 supported.\n",
              GDALGetDataTypeName(eType) );

        return NULL;
    }

    if( nBands != 1 && nBands != 3 )
    {
        CPLError( CE_Failure, CPLE_AppDefined,
                  "Attempt to create PNM dataset with an illegal number\n"
                  "of bands (%d).  Must be 1 (greyscale) or 3 (RGB).\n",
                  nBands );

        return NULL;
    }

/* -------------------------------------------------------------------- */
/*      Try to create the file.                                         */
/* -------------------------------------------------------------------- */
    VSILFILE        *fp;

    fp = VSIFOpenL( pszFilename, "wb" );

    if( fp == NULL )
    {
        CPLError( CE_Failure, CPLE_OpenFailed,
                  "Attempt to create file `%s' failed.\n",
                  pszFilename );
        return NULL;
    }

/* -------------------------------------------------------------------- */
/*      Write out the header.                                           */
/* -------------------------------------------------------------------- */
    char        szHeader[500];
    const char  *pszMaxValue = NULL;
    int         nMaxValue = 0;

    pszMaxValue = CSLFetchNameValue( papszOptions, "MAXVAL" );
    if ( pszMaxValue )
    {
        nMaxValue = atoi( pszMaxValue );
        if ( eType == GDT_Byte && (nMaxValue > 255 || nMaxValue < 0) )
            nMaxValue = 255;
        else if ( nMaxValue > 65535 || nMaxValue < 0 )
            nMaxValue = 65535;
    }
    else
    {
        if ( eType == GDT_Byte )
            nMaxValue = 255;
        else
            nMaxValue = 65535;
    }


    memset( szHeader, 0, sizeof(szHeader) );

    if( nBands == 3 )
        sprintf( szHeader, "P6\n%d %d\n%d\n", nXSize, nYSize, nMaxValue );
    else
        sprintf( szHeader, "P5\n%d %d\n%d\n", nXSize, nYSize, nMaxValue );

    VSIFWriteL( (void *) szHeader, strlen(szHeader) + 2, 1, fp );
    VSIFCloseL( fp );

    return (GDALDataset *) GDALOpen( pszFilename, GA_Update );
}

int LLVMFuzzerTestOneInput(const uint8_t *buf, size_t len)
{
    VSILFILE* fp = VSIFileFromMemBuffer( "/vsimem/test.tar",
            reinterpret_cast<GByte*>(const_cast<uint8_t*>(buf)), len, FALSE );
    VSIFCloseL(fp);

    CPLPushErrorHandler(CPLQuietErrorHandler);

    char** papszArgv = nullptr;

    // Prevent generating too big output raster. Make sure they are set at
    // the beginning to avoid being accidentally eaten by invalid arguments
    // afterwards.
    papszArgv = CSLAddString(papszArgv, "-limit_outsize");
    papszArgv = CSLAddString(papszArgv, "1000000");

    fp = VSIFOpenL("/vsitar//vsimem/test.tar/cmd.txt", "rb");
    if( fp != nullptr )
    {
        const char* pszLine = nullptr;
        while( (pszLine = CPLReadLineL(fp)) != nullptr )
        {
            if( !EQUAL(pszLine, "-limit_outsize") )
                papszArgv = CSLAddString(papszArgv, pszLine);
        }
        VSIFCloseL(fp);
    }

    int nXDim = -1;
    int nYDim = -1;
    bool bXDimPct = false;
    bool bYDimPct = false;
    bool bNonNearestResampling = false;
    int nBlockXSize = 0;
    int nBlockYSize = 0;
    bool bStatsEnabled = false;
    bool bHFA = false;
    if( papszArgv != nullptr )
    {
        int nCount = CSLCount(papszArgv);
        for( int i = 0; i < nCount; i++ )
        {
            if( EQUAL(papszArgv[i], "-outsize") && i + 2 < nCount )
            {
                nXDim = atoi(papszArgv[i+1]);
                bXDimPct = (papszArgv[i+1][0] != '\0' &&
                            papszArgv[i+1][strlen(papszArgv[i+1])-1] == '%');
                nYDim = atoi(papszArgv[i+2]);
                bYDimPct = (papszArgv[i+2][0] != '\0' &&
                            papszArgv[i+2][strlen(papszArgv[i+2])-1] == '%');
            }
            else if( EQUAL(papszArgv[i], "-r") && i + 1 < nCount )
            {
                bNonNearestResampling = !STARTS_WITH_CI(papszArgv[i+1], "NEAR");
            }
            else if( EQUAL(papszArgv[i], "-co") && i + 1 < nCount )
            {
                if( STARTS_WITH_CI(papszArgv[i+1], "BLOCKSIZE=") )
                {
                    nBlockXSize = std::max(nBlockXSize,
                                atoi(papszArgv[i+1]+strlen("BLOCKSIZE=")));
                    nBlockYSize = std::max(nBlockYSize,
                                atoi(papszArgv[i+1]+strlen("BLOCKSIZE=")));
                }
                else if( STARTS_WITH_CI(papszArgv[i+1], "BLOCKXSIZE=") )
                {
                    nBlockXSize = std::max(nBlockXSize,
                                atoi(papszArgv[i+1]+strlen("BLOCKXSIZE=")));
                }
                else if( STARTS_WITH_CI(papszArgv[i+1], "BLOCKYSIZE=") )
                {
                    nBlockYSize = std::max(nBlockYSize,
                                atoi(papszArgv[i+1]+strlen("BLOCKYSIZE=")));
                }
            }
            else if( EQUAL(papszArgv[i], "-stats") )
            {
                bStatsEnabled = true;
            }
            else if( EQUAL(papszArgv[i], "-of") && i + 1 < nCount )
            {
                bHFA = EQUAL( papszArgv[i+1], "HFA" );
            }
        }
        if( bHFA )
        {
            // Disable statistics computation for HFA, as it can be time
            // consuming.
            // See https://bugs.chromium.org/p/oss-fuzz/issues/detail?id=10067
            papszArgv = CSLInsertString(papszArgv, 0, "-co");
            papszArgv = CSLInsertString(papszArgv, 1, "STATISTICS=NO");
        }
    }

    if( papszArgv != nullptr )
    {
        GDALTranslateOptions* psOptions = GDALTranslateOptionsNew(papszArgv, nullptr);
        if( psOptions )
        {
            GDALDatasetH hSrcDS = GDALOpen( "/vsitar//vsimem/test.tar/in", GA_ReadOnly );
//.........这里部分代码省略.........

//.........这里部分代码省略.........
        else if( EQUAL(argv[i], "-nomd") )
            bShowMetadata = FALSE;
        else if( EQUAL(argv[i], "-norat") )
            bShowRAT = FALSE;
        else if( EQUAL(argv[i], "-noct") )
            bShowColorTable = FALSE;
        else if( EQUAL(argv[i], "-mdd") )
        {
            CHECK_HAS_ENOUGH_ADDITIONAL_ARGS(1);
            papszExtraMDDomains = CSLAddString( papszExtraMDDomains,
                                                argv[++i] );
        }
        else if( EQUAL(argv[i], "-nofl") )
            bShowFileList = FALSE;
        else if( EQUAL(argv[i], "-sd") )
        {
            CHECK_HAS_ENOUGH_ADDITIONAL_ARGS(1);
            nSubdataset = atoi(argv[++i]);
        }
        else if( argv[i][0] == '-' )
            Usage(CPLSPrintf("Unkown option name '%s'", argv[i]));
        else if( pszFilename == NULL )
            pszFilename = argv[i];
        else
            Usage("Too many command options.");
    }

    if( pszFilename == NULL )
        Usage("No datasource specified.");

/* -------------------------------------------------------------------- */
/*      Open dataset.                                                   */
/* -------------------------------------------------------------------- */
    hDataset = GDALOpen( pszFilename, GA_ReadOnly );

    if( hDataset == NULL )
    {
        fprintf( stderr,
                 "gdalinfo failed - unable to open '%s'.\n",
                 pszFilename );

/* -------------------------------------------------------------------- */
/*      If argument is a VSIFILE, then print its contents               */
/* -------------------------------------------------------------------- */
        if ( strncmp( pszFilename, "/vsizip/", 8 ) == 0 || 
             strncmp( pszFilename, "/vsitar/", 8 ) == 0 ) 
        {
            papszFileList = VSIReadDirRecursive( pszFilename );
            if ( papszFileList )
            {
                int nCount = CSLCount( papszFileList );
                fprintf( stdout, 
                         "Unable to open source `%s' directly.\n"
                         "The archive contains %d files:\n", 
                         pszFilename, nCount );
                for ( i = 0; i < nCount; i++ )
                {
                    fprintf( stdout, "       %s/%s\n", pszFilename, papszFileList[i] );
                }
                CSLDestroy( papszFileList );
                papszFileList = NULL;
            }
        }

        CSLDestroy( argv );
        CSLDestroy( papszExtraMDDomains );

int equalize_density(char *infile, char *outfile, int fast, int accurate) {
	
	int xsize, ysize;				// Size of the density grid.
	double *gridx, *gridy;			// Array for grid	
	double **rho;					// Initial population density
	GDALDatasetH hDataset;			// The input density raster file.
	GDALRasterBandH hBand;			// The raster band we are going to use.
	FILE *outfp;					// The morphing file (a text file).
	double adfGeoTransform[6];		// For the georeference of the raster.
	
	
	// Register all GDAL drivers.
    GDALAllRegister();
	
	
#if defined (_OPENMP)
	omp_set_num_threads(omp_get_num_procs());
#endif
	
	
	hDataset = GDALOpen(infile, GA_ReadOnly);
    if (hDataset == NULL) {
		fprintf(stderr,"Error. Unable to open file `%s'\n", infile);
		exit(1);
	}
	
	outfp = fopen(outfile, "w");
	if (outfp == NULL) {
		fprintf(stderr,"Error. Unable to open file `%s'\n", outfile);
		exit(1);
	}
	
	
	// Get the raster band for the dataset; we are using the first band.
	hBand = GDALGetRasterBand(hDataset, 1);
	if (hBand == NULL) {
		fprintf(stderr, "Error. Unable to read band 1 in file `%s'\n", infile);
		exit(1);
	}
	
	// Determine the raster size
	xsize = GDALGetRasterBandXSize(hBand);
	ysize = GDALGetRasterBandYSize(hBand);
	
	
	
	// Allocate space for the cartogram code to use
	cart_makews(xsize, ysize);
	
	
	// Read in the population data, transform it, then destroy it again
	rho = cart_dmalloc(xsize, ysize);
	if (readpop(hBand, rho, xsize, ysize)) {
		fprintf(stderr,"Error. Density file contains too few or incorrect data\n");
		exit(1);
	}
	cart_transform(rho, xsize, ysize);
	cart_dfree(rho);
	
	
	// Create the grid of points
	gridx = malloc((xsize+1)*(ysize+1)*sizeof(double));
	gridy = malloc((xsize+1)*(ysize+1)*sizeof(double));
	creategrid(gridx, gridy, xsize, ysize);
	
	
	// Compute the cartogram
	cart_makecart(gridx, gridy, (xsize+1)*(ysize+1), xsize, ysize, 0.0);
	
	
	// Write out the final positions of the grid points
	GDALGetGeoTransform(hDataset, adfGeoTransform);
	writepoints(outfp, gridx, gridy, xsize, ysize, adfGeoTransform);
	//writepoints(outfp, gridx, gridy, (xsize+1)*(ysize+1));
	
	
	// Free up the allocated memory
	cart_freews(xsize, ysize);
	free(gridx);
	free(gridy);
	
	
	// Close the input and output files
	GDALClose(hDataset);
	fclose(outfp);
	
	return 0;
}

//.........这里部分代码省略.........
/* -------------------------------------------------------------------- */
    VSILFILE *fp = VSIFOpenL( pszFilename, "w" );
    if( fp == NULL )
    {
        CPLError( CE_Failure, CPLE_OpenFailed,
                  "Failed to create file '%s': %s",
                  pszFilename, VSIStrerror( errno ) );
        return NULL;
    }

/* -------------------------------------------------------------------- */
/*      Prepare and write 512 byte header.                              */
/* -------------------------------------------------------------------- */
    GByte abyHeader[512];

    memset( abyHeader, 0, 512 );

    abyHeader[0] = 1;
    abyHeader[1] = 218;
    abyHeader[2] = 1; // RLE
    abyHeader[3] = 1;  // 8bit

    GInt16 nShortValue;
    if( nBands == 1 )
        nShortValue = CPL_MSBWORD16(2);
    else
        nShortValue = CPL_MSBWORD16(3);
    memcpy( abyHeader + 4, &nShortValue, 2 );

    nShortValue = CPL_MSBWORD16(nXSize);
    memcpy( abyHeader + 6, &nShortValue, 2 );

    nShortValue = CPL_MSBWORD16(nYSize);
    memcpy( abyHeader + 8, &nShortValue, 2 );

    nShortValue = CPL_MSBWORD16(nBands);
    memcpy( abyHeader + 10, &nShortValue, 2 );

    GInt32 nIntValue = CPL_MSBWORD32(0);
    memcpy( abyHeader + 12, &nIntValue, 4 );

    GUInt32 nUIntValue = CPL_MSBWORD32(255);
    memcpy( abyHeader + 16, &nUIntValue, 4 );

    VSIFWriteL( abyHeader, 1, 512, fp );

/* -------------------------------------------------------------------- */
/*      Create our RLE compressed zero-ed dummy line.                   */
/* -------------------------------------------------------------------- */
    GByte *pabyRLELine = reinterpret_cast<GByte *>(
        CPLMalloc( ( nXSize / 127 ) * 2 + 4 ) );

    int nPixelsRemaining = nXSize;
    GInt32 nRLEBytes = 0;
    while( nPixelsRemaining > 0 )
    {
        pabyRLELine[nRLEBytes] = static_cast<GByte>(
            std::min( 127, nPixelsRemaining  ) );
        pabyRLELine[nRLEBytes+1] = 0;
        nPixelsRemaining -= pabyRLELine[nRLEBytes];

        nRLEBytes += 2;
    }

/* -------------------------------------------------------------------- */
/*      Prepare and write RLE offset/size tables with everything        */
/*      zeroed indicating dummy lines.                                  */
/* -------------------------------------------------------------------- */
    const int nTableLen = nYSize * nBands;
    GInt32 nDummyRLEOffset = 512 + 4 * nTableLen * 2;

    CPL_MSBPTR32( &nRLEBytes );
    CPL_MSBPTR32( &nDummyRLEOffset );

    for( int i = 0; i < nTableLen; i++ )
        VSIFWriteL( &nDummyRLEOffset, 1, 4, fp );

    for( int i = 0; i < nTableLen; i++ )
        VSIFWriteL( &nRLEBytes, 1, 4, fp );

/* -------------------------------------------------------------------- */
/*      write the dummy RLE blank line.                                 */
/* -------------------------------------------------------------------- */
    CPL_MSBPTR32( &nRLEBytes );
    if( static_cast<GInt32>( VSIFWriteL( pabyRLELine, 1, nRLEBytes, fp ) )
        != nRLEBytes )
    {
        CPLError( CE_Failure, CPLE_FileIO,
                  "Failure writing SGI file '%s'.\n%s",
                  pszFilename,
                  VSIStrerror( errno ) );
        return NULL;
    }

    VSIFCloseL( fp );
    CPLFree( pabyRLELine );

    return reinterpret_cast<GDALDataset *>(
        GDALOpen( pszFilename, GA_Update ) );
}