//.........这里部分代码省略.........
                              pdfSrcNoData, nTargetValues, panTargetValues );

        // Write out results.
        eErr = 
            GDALRasterIO( hWorkProximityBand, GF_Write, 0, iLine, nXSize, 1, 
                          pafProximity, nXSize, 1, GDT_Float32, 0, 0 );

        if( eErr != CE_None )
            break;

        if( !pfnProgress( 0.5 * (iLine+1) / (double) nYSize, 
                          "", pProgressArg ) )
        {
            CPLError( CE_Failure, CPLE_UserInterrupt, "User terminated" );
            eErr = CE_Failure;
        }
    }

/* -------------------------------------------------------------------- */
/*      Loop from bottom to top of the image.                           */
/* -------------------------------------------------------------------- */
    for( i = 0; i < nXSize; i++ )
        panNearX[i] = panNearY[i] = -1;

    for( iLine = nYSize-1; eErr == CE_None && iLine >= 0; iLine-- )
    {
        // Read first pass proximity
        eErr = 
            GDALRasterIO( hWorkProximityBand, GF_Read, 0, iLine, nXSize, 1, 
                          pafProximity, nXSize, 1, GDT_Float32, 0, 0 );

        if( eErr != CE_None )
            break;

        // Read pixel values.

        eErr = GDALRasterIO( hSrcBand, GF_Read, 0, iLine, nXSize, 1, 
                             panSrcScanline, nXSize, 1, GDT_Int32, 0, 0 );
        if( eErr != CE_None )
            break;

        // Right to left
        ProcessProximityLine( panSrcScanline, panNearX, panNearY, 
                              FALSE, iLine, nXSize, dfMaxDist, pafProximity,
                              pdfSrcNoData, nTargetValues, panTargetValues );

        // Left to right
        ProcessProximityLine( panSrcScanline, panNearX, panNearY, 
                              TRUE, iLine, nXSize, dfMaxDist, pafProximity,
                              pdfSrcNoData, nTargetValues, panTargetValues );

        // Final post processing of distances. 
        for( i = 0; i < nXSize; i++ )
        {
            if( pafProximity[i] < 0.0 )
                pafProximity[i] = fNoDataValue;
            else if( pafProximity[i] > 0.0 )
            {
                if( bFixedBufVal )
                    pafProximity[i] = (float) dfFixedBufVal;
                else 
                    pafProximity[i] = (float)(pafProximity[i] * dfDistMult);
            }
        }
  
        // Write out results.
        eErr = 
            GDALRasterIO( hProximityBand, GF_Write, 0, iLine, nXSize, 1, 
                          pafProximity, nXSize, 1, GDT_Float32, 0, 0 );

        if( eErr != CE_None )
            break;

        if( !pfnProgress( 0.5 + 0.5 * (nYSize-iLine) / (double) nYSize, 
                          "", pProgressArg ) )
        {
            CPLError( CE_Failure, CPLE_UserInterrupt, "User terminated" );
            eErr = CE_Failure;
        }
    }

/* -------------------------------------------------------------------- */
/*      Cleanup                                                         */
/* -------------------------------------------------------------------- */
end:
    CPLFree( panNearX );
    CPLFree( panNearY );
    CPLFree( panSrcScanline );
    CPLFree( pafProximity );
    CPLFree(panTargetValues);

    if( hWorkProximityDS != NULL )
    {
        CPLString osProxFile = GDALGetDescription( hWorkProximityDS );
        GDALClose( hWorkProximityDS );
        GDALDeleteDataset( GDALGetDriverByName( "GTiff" ), osProxFile );
    }

    return eErr;
}

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

        CSLDestroy( argv );
        CSLDestroy( papszExtraMDDomains );
    
        GDALDumpOpenDatasets( stderr );

        GDALDestroyDriverManager();

        CPLDumpSharedList( NULL );

        exit( 1 );
    }
    
/* -------------------------------------------------------------------- */
/*      Read specified subdataset if requested.                         */
/* -------------------------------------------------------------------- */
    if ( nSubdataset > 0 )
    {
        char **papszSubdatasets = GDALGetMetadata( hDataset, "SUBDATASETS" );
        int nSubdatasets = CSLCount( papszSubdatasets );

        if ( nSubdatasets > 0 && nSubdataset <= nSubdatasets )
        {
            char szKeyName[1024];
            char *pszSubdatasetName;

            snprintf( szKeyName, sizeof(szKeyName),
                      "SUBDATASET_%d_NAME", nSubdataset );
            szKeyName[sizeof(szKeyName) - 1] = '\0';
            pszSubdatasetName =
                CPLStrdup( CSLFetchNameValue( papszSubdatasets, szKeyName ) );
            GDALClose( hDataset );
            hDataset = GDALOpen( pszSubdatasetName, GA_ReadOnly );
            CPLFree( pszSubdatasetName );
        }
        else
        {
            fprintf( stderr,
                     "gdalinfo warning: subdataset %d of %d requested. "
                     "Reading the main dataset.\n",
                     nSubdataset, nSubdatasets );

        }
    }

/* -------------------------------------------------------------------- */
/*      Report general info.                                            */
/* -------------------------------------------------------------------- */
    hDriver = GDALGetDatasetDriver( hDataset );
    printf( "Driver: %s/%s\n",
            GDALGetDriverShortName( hDriver ),
            GDALGetDriverLongName( hDriver ) );

    papszFileList = GDALGetFileList( hDataset );
    if( CSLCount(papszFileList) == 0 )
    {
        printf( "Files: none associated\n" );
    }
    else
    {
        printf( "Files: %s\n", papszFileList[0] );
        if( bShowFileList )
        {
            for( i = 1; papszFileList[i] != NULL; i++ )

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);
//.........这里部分代码省略.........

//.........这里部分代码省略.........
                           iLayer);
                    exit(1);
                }

                if( !bAllLayers )
                {
                    printf("%d: %s", iLayer + 1, poLayer->GetName());

                    const int nGeomFieldCount =
                        poLayer->GetLayerDefn()->GetGeomFieldCount();
                    if( nGeomFieldCount > 1 )
                    {
                        printf(" (");
                        for( int iGeom = 0; iGeom < nGeomFieldCount; iGeom++ )
                        {
                            if( iGeom > 0 )
                                printf(", ");
                            OGRGeomFieldDefn* poGFldDefn =
                                poLayer->GetLayerDefn()->
                                    GetGeomFieldDefn(iGeom);
                            printf(
                                "%s",
                                OGRGeometryTypeToName(
                                    poGFldDefn->GetType()));
                        }
                        printf(")");
                    }
                    else if( poLayer->GetGeomType() != wkbUnknown )
                        printf(" (%s)",
                               OGRGeometryTypeToName(
                                   poLayer->GetGeomType()));

                    printf("\n");
                }
                else
                {
                    if( iRepeat != 0 )
                        poLayer->ResetReading();

                    ReportOnLayer(poLayer, pszWHERE,
                                  pszGeomField, poSpatialFilter,
                                  bListMDD, bShowMetadata, papszExtraMDDomains,
                                  bFeatureCount, bExtent, pszWKTFormat);
                }
            }
        }
        else
        {
/* -------------------------------------------------------------------- */
/*      Process specified data source layers.                           */
/* -------------------------------------------------------------------- */

            for( char** papszIter = papszLayers;
                 *papszIter != nullptr;
                 ++papszIter )
            {
                OGRLayer *poLayer = poDS->GetLayerByName(*papszIter);

                if( poLayer == nullptr )
                {
                    printf("FAILURE: Couldn't fetch requested layer %s!\n",
                           *papszIter);
                    exit(1);
                }

                if( iRepeat != 0 )
                    poLayer->ResetReading();

                ReportOnLayer(poLayer, pszWHERE, pszGeomField, poSpatialFilter,
                              bListMDD, bShowMetadata, papszExtraMDDomains,
                              bFeatureCount, bExtent, pszWKTFormat);
            }
        }
    }

/* -------------------------------------------------------------------- */
/*      Close down.                                                     */
/* -------------------------------------------------------------------- */
    GDALClose(poDS);

#ifdef __AFL_HAVE_MANUAL_CONTROL
    }
#else
end:
#endif

    CSLDestroy(papszArgv);
    CSLDestroy(papszLayers);
    CSLDestroy(papszOptions);
    CSLDestroy(papszOpenOptions);
    CSLDestroy(papszExtraMDDomains);
    if( poSpatialFilter )
        OGRGeometryFactory::destroyGeometry(poSpatialFilter);
    CPLFree(pszSQLStatement);
    CPLFree(pszWHERE);

    OGRCleanupAll();

    return nRet;
}

//.........这里部分代码省略.........
    // if crs transform needed
    if ( it->raster->crs() != mOutputCrs )
    {
      QgsRasterProjector proj;
      proj.setCrs( it->raster->crs(), mOutputCrs );
      proj.setInput( it->raster->dataProvider() );
      proj.setPrecision( QgsRasterProjector::Exact );

      block = proj.block( it->bandNumber, mOutputRectangle, mNumOutputColumns, mNumOutputRows );
    }
    else
    {
      block = it->raster->dataProvider()->block( it->bandNumber, mOutputRectangle, mNumOutputColumns, mNumOutputRows );
    }
    if ( block->isEmpty() )
    {
      delete block;
      delete calcNode;
      qDeleteAll( inputBlocks );
      return static_cast<int>( MemoryError );
    }
    inputBlocks.insert( it->ref, block );
  }

  //open output dataset for writing
  GDALDriverH outputDriver = openOutputDriver();
  if ( !outputDriver )
  {
    return static_cast< int >( CreateOutputError );
  }

  GDALDatasetH outputDataset = openOutputFile( outputDriver );
  GDALSetProjection( outputDataset, mOutputCrs.toWkt().toLocal8Bit().data() );
  GDALRasterBandH outputRasterBand = GDALGetRasterBand( outputDataset, 1 );

  float outputNodataValue = -FLT_MAX;
  GDALSetRasterNoDataValue( outputRasterBand, outputNodataValue );

  if ( p )
  {
    p->setMaximum( mNumOutputRows );
  }

  QgsRasterMatrix resultMatrix;
  resultMatrix.setNodataValue( outputNodataValue );

  //read / write line by line
  for ( int i = 0; i < mNumOutputRows; ++i )
  {
    if ( p )
    {
      p->setValue( i );
    }

    if ( p && p->wasCanceled() )
    {
      break;
    }

    if ( calcNode->calculate( inputBlocks, resultMatrix, i ) )
    {
      bool resultIsNumber = resultMatrix.isNumber();
      float* calcData = new float[mNumOutputColumns];

      for ( int j = 0; j < mNumOutputColumns; ++j )
      {
        calcData[j] = ( float )( resultIsNumber ? resultMatrix.number() : resultMatrix.data()[j] );
      }

      //write scanline to the dataset
      if ( GDALRasterIO( outputRasterBand, GF_Write, 0, i, mNumOutputColumns, 1, calcData, mNumOutputColumns, 1, GDT_Float32, 0, 0 ) != CE_None )
      {
        QgsDebugMsg( "RasterIO error!" );
      }

      delete[] calcData;
    }

  }

  if ( p )
  {
    p->setValue( mNumOutputRows );
  }

  //close datasets and release memory
  delete calcNode;
  qDeleteAll( inputBlocks );
  inputBlocks.clear();

  if ( p && p->wasCanceled() )
  {
    //delete the dataset without closing (because it is faster)
    GDALDeleteDataset( outputDriver, TO8F( mOutputFile ) );
    return static_cast< int >( Cancelled );
  }
  GDALClose( outputDataset );

  return static_cast< int >( Success );
}

// Slot called when the menu item is triggered
// If you created more menu items / toolbar buttons in initiGui, you should
// create a separate handler for each action - this single run() method will
// not be enough
void Heatmap::run()
{
  HeatmapGui d( mQGisIface->mainWindow(), QgisGui::ModalDialogFlags, &mSessionSettings );

  if ( d.exec() == QDialog::Accepted )
  {
    // everything runs here

    // Get the required data from the dialog
    QgsRectangle myBBox = d.bbox();
    int columns = d.columns();
    int rows = d.rows();
    double cellsize = d.cellSizeX(); // or d.cellSizeY();  both have the same value
    mDecay = d.decayRatio();
    int kernelShape = d.kernelShape();

    // Start working on the input vector
    QgsVectorLayer* inputLayer = d.inputVectorLayer();

    // Getting the rasterdataset in place
    GDALAllRegister();

    GDALDataset *emptyDataset;
    GDALDriver *myDriver;

    myDriver = GetGDALDriverManager()->GetDriverByName( d.outputFormat().toUtf8() );
    if ( myDriver == NULL )
    {
      QMessageBox::information( 0, tr( "GDAL driver error" ), tr( "Cannot open the driver for the specified format" ) );
      return;
    }

    double geoTransform[6] = { myBBox.xMinimum(), cellsize, 0, myBBox.yMinimum(), 0, cellsize };
    emptyDataset = myDriver->Create( d.outputFilename().toUtf8(), columns, rows, 1, GDT_Float32, NULL );
    emptyDataset->SetGeoTransform( geoTransform );
    // Set the projection on the raster destination to match the input layer
    emptyDataset->SetProjection( inputLayer->crs().toWkt().toLocal8Bit().data() );

    GDALRasterBand *poBand;
    poBand = emptyDataset->GetRasterBand( 1 );
    poBand->SetNoDataValue( NO_DATA );

    float* line = ( float * ) CPLMalloc( sizeof( float ) * columns );
    for ( int i = 0; i < columns ; i++ )
    {
      line[i] = NO_DATA;
    }
    // Write the empty raster
    for ( int i = 0; i < rows ; i++ )
    {
      poBand->RasterIO( GF_Write, 0, i, columns, 1, line, columns, 1, GDT_Float32, 0, 0 );
    }

    CPLFree( line );
    //close the dataset
    GDALClose(( GDALDatasetH ) emptyDataset );

    // open the raster in GA_Update mode
    GDALDataset *heatmapDS;
    heatmapDS = ( GDALDataset * ) GDALOpen( TO8F( d.outputFilename() ), GA_Update );
    if ( !heatmapDS )
    {
      QMessageBox::information( 0, tr( "Raster update error" ), tr( "Could not open the created raster for updating. The heatmap was not generated." ) );
      return;
    }
    poBand = heatmapDS->GetRasterBand( 1 );

    QgsAttributeList myAttrList;
    int rField = 0;
    int wField = 0;

    // Handle different radius options
    double radius;
    double radiusToMapUnits = 1;
    int myBuffer = 0;
    if ( d.variableRadius() )
    {
      rField = d.radiusField();
      myAttrList.append( rField );
      QgsDebugMsg( QString( "Radius Field index received: %1" ).arg( rField ) );

      // If not using map units, then calculate a conversion factor to convert the radii to map units
      if ( d.radiusUnit() == HeatmapGui::Meters )
      {
        radiusToMapUnits = mapUnitsOf( 1, inputLayer->crs() );
      }
    }
    else
    {
      radius = d.radius(); // radius returned by d.radius() is already in map units
      myBuffer = bufferSize( radius, cellsize );
    }

    if ( d.weighted() )
    {
      wField = d.weightField();
//.........这里部分代码省略.........

//.........这里部分代码省略.........
   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] );
   } else {
     error_exit(10, "Keine Transformation im TIFF vorhanden!\n");
   }  
   
   // Geotiff Fehler abfangen 
   if( h_dset == NULL ) {
     error_exit(10, "Datensatz %s kann nicht geoeffnet werden!\n",
                ifile);
   }

  // Bilddimensionen ermitteln 
  int img_width  = GDALGetRasterXSize( h_dset ); 
  int img_height = GDALGetRasterYSize( h_dset );
  int num_bands  = GDALGetRasterCount (h_dset );   

  // Bilddimensionen ermitteln 
  GDALRasterBandH h_band[num_bands];
  GDALDataType    h_type[num_bands];
  int             h_tsize[num_bands];

  for(int b=0 ; b<num_bands; b++) {
    h_band[b]  = GDALGetRasterBand( h_dset, b+1 );   
    h_type[b]  = GDALGetRasterDataType(h_band[b]);   
    h_tsize[b] = GDALGetDataTypeSize(h_type[b]);
  }  

  // Erzeuge Bildschnitte
  for (int c=0; c< pos_x.length; c++ ) {
	
	 // Welt zu Pixeltransformation
	 long icol = -1; long irow = -1;
     trfm_geo_pix(trfm, pos_x.data[c], pos_y.data[c],
                &icol , &irow);

	 // Dateinamen erzeugen
	 sprintf(cfile,"%s.%d%s",ofile, id.data[c],ext);

	 // Test ob das Schnittfenstrer passt
	 if (icol-size/2<=0 || 
	      irow-size/2<=0 || 
	      icol+size/2>=img_width || 
	      irow+size/2>=img_height){
		printf ("IGN %d %s\n",id.data[c],ofile);
		continue;
	 }		  	

	 // Schneiden
  	 printf ("ADD %d %s\n",id.data[c],ofile);
	      
	 int ioffs_col = icol-size/2;  
     int ioffs_row = irow-size/2; 
	
	 // Tiff-Datei anlegen 
     // char **options = NULL;
     GDALDatasetH h_dset_out = GDALCreate( h_drv, 
         cfile, size, size, num_bands, h_type[0], NULL);
	
	 // Iteration ueber alle Baender 
     for (int b=0; b< num_bands; b++) {

   	    // IO Buffer allozieren @todo static
	    void *io_buffer = CPLMalloc(h_tsize[0] * size * size);

		// Pixel lesen 
        GDALRasterIO( h_band[b], GF_Read, 
                      ioffs_col, ioffs_row, size, size, 
                      io_buffer, size, size, h_type[b], 0, 0 );
		// Pixel schreiben  
		GDALRasterBandH h_band_out = GDALGetRasterBand(h_dset_out, b+1);
		GDALRasterIO( h_band_out, GF_Write, 0, 0, size, size, 
                      io_buffer, size, size, h_type[b], 0, 0 );  
	    // IO Buffer freigeben
        CPLFree(io_buffer);
     }
     // resultierendes Tiff flushen und schliessen
     GDALClose( h_dset_out );         

  } // EOF Positions

  // Eingangsbild
  GDALClose( h_dset);         
  return 0;
}

//.........这里部分代码省略.........
        if (pszWkt && pszWkt[0])
        {
            hSRS = OSRNewSpatialReference(pszWkt);
        }

        if (pszLayerName == NULL)
            pszLayerName = CPLStrdup(CPLGetBasename(pszDstFilename));

        hOGRLayer = OGR_DS_CreateLayer( hOGRDS, pszLayerName,
                                        hSRS, eType, NULL);

        if (hSRS)
            OSRDestroySpatialReference(hSRS);

        if (hOGRLayer == NULL)
        {
            fprintf(stderr, "Can't create layer %s\n", pszLayerName);
            exit(1);
        }

        if (eType != wkbPoint25D)
        {
            OGRFieldDefnH hFieldDefn =  OGR_Fld_Create( "z", OFTReal );
            OGR_L_CreateField(hOGRLayer, hFieldDefn, 0);
            OGR_Fld_Destroy( hFieldDefn );
        }
    }


    padfBuffer = (double*)CPLMalloc(nXSize * sizeof(double));

#define GET_X(j, i) adfGeotransform[0] + (j) * adfGeotransform[1] + (i) * adfGeotransform[2]
#define GET_Y(j, i) adfGeotransform[3] + (j) * adfGeotransform[4] + (i) * adfGeotransform[5]
#define GET_XY(j, i) GET_X(j, i), GET_Y(j, i)

/* -------------------------------------------------------------------- */
/*     "Translate" the source dataset                                   */
/* -------------------------------------------------------------------- */
    for(i=0;i<nYSize;i+=yStep)
    {
        GDALRasterIO( hBand, GF_Read, 0, i, nXSize, 1,
                      padfBuffer, nXSize, 1, GDT_Float64, 0, 0);
        for(j=0;j<nXSize;j+=xStep)
        {
            if (fOut)
            {
                fprintf(fOut, "%f,%f,%f\n",
                        GET_XY(j + .5, i + .5), padfBuffer[j]);
            }
            else
            {
                OGRFeatureH hFeature = OGR_F_Create(OGR_L_GetLayerDefn(hOGRLayer));
                OGRGeometryH hGeometry = OGR_G_CreateGeometry(eType);
                if (eType == wkbPoint25D)
                {
                    OGR_G_SetPoint(hGeometry, 0, GET_XY(j + .5, i + .5),
                                   padfBuffer[j]);
                }
                else if (eType == wkbPoint)
                {
                    OGR_G_SetPoint_2D(hGeometry, 0, GET_XY(j + .5, i + .5));
                    OGR_F_SetFieldDouble(hFeature, 0, padfBuffer[j]);
                }
                else
                {
                    OGRGeometryH hLinearRing = OGR_G_CreateGeometry(wkbLinearRing);
                    OGR_G_SetPoint_2D(hLinearRing, 0, GET_XY(j + 0, i + 0));
                    OGR_G_SetPoint_2D(hLinearRing, 1, GET_XY(j + 1, i + 0));
                    OGR_G_SetPoint_2D(hLinearRing, 2, GET_XY(j + 1, i + 1));
                    OGR_G_SetPoint_2D(hLinearRing, 3, GET_XY(j + 0, i + 1));
                    OGR_G_SetPoint_2D(hLinearRing, 4, GET_XY(j + 0, i + 0));
                    OGR_G_AddGeometryDirectly(hGeometry, hLinearRing);
                    OGR_F_SetFieldDouble(hFeature, 0, padfBuffer[j]);
                }
                OGR_F_SetGeometryDirectly(hFeature, hGeometry);
                OGR_L_CreateFeature(hOGRLayer, hFeature);
                OGR_F_Destroy(hFeature);
            }
        }
    }

/* -------------------------------------------------------------------- */
/*     Cleanup                                                          */
/* -------------------------------------------------------------------- */
    if (fOut)
        fclose(fOut);
    else
        OGR_DS_Destroy(hOGRDS);
    GDALClose(hDS);

    CPLFree(padfBuffer);
    CPLFree(pszLayerName);

    GDALDumpOpenDatasets( stderr );
    GDALDestroyDriverManager();
    OGRCleanupAll();
    CSLDestroy( argv );

    return 0;
}

//.........这里部分代码省略.........
    int          bHasNoData;

/* -------------------------------------------------------------------- */
/*      Identify arguments.                                             */
/* -------------------------------------------------------------------- */

    for( iArg = 1; iArg < argc; iArg++ )
    {
        if( EQUAL(argv[iArg],"-trim") )
            bEnableTrim = TRUE;

        else if( EQUAL(argv[iArg],"-fill") )
            nFillDist = atoi(argv[++iArg]);

        else if( EQUAL(argv[iArg],"-level") )
            nOutLevel = atoi(argv[++iArg]);
        else
        {
            if( pszFilename != NULL )
                Usage();
            pszFilename = argv[iArg];
        }
    }		

    if( pszFilename == NULL )
        Usage();

/* -------------------------------------------------------------------- */
/*      Open input file.                                                */
/* -------------------------------------------------------------------- */
    GDALAllRegister();
    hSrcDS = GDALOpen( pszFilename, GA_ReadOnly );
    if( hSrcDS == NULL )
        exit(1);

    hSrcBand = GDALGetRasterBand( hSrcDS, 1 );

    noDataValue = (GInt16)GDALGetRasterNoDataValue(hSrcBand, &bHasNoData);

    nXSize = GDALGetRasterXSize( hSrcDS );
    nYSize = GDALGetRasterYSize( hSrcDS );

    GDALGetGeoTransform( hSrcDS, adfGeoTransform );

/* -------------------------------------------------------------------- */
/*      Create output stream.                                           */
/* -------------------------------------------------------------------- */
    pStream = DTEDCreatePtStream( ".", nOutLevel );

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

/* -------------------------------------------------------------------- */
/*      Process all the profiles.                                       */
/* -------------------------------------------------------------------- */
    panData = (GInt16 *) malloc(sizeof(GInt16) * nXSize);

    for( iY = 0; iY < nYSize; iY++ )
    {
        GDALRasterIO( hSrcBand, GF_Read, 0, iY, nXSize, 1, 
                      panData, nXSize, 1, GDT_Int16, 0, 0 );

        if (bHasNoData)
        {
            for( iX = 0; iX < nXSize; iX++ )
            {
                if (panData[iX] == noDataValue)
                    panData[iX] = DTED_NODATA_VALUE;
            }
        }

        for( iX = 0; iX < nXSize; iX++ )
        {
            DTEDWritePt( pStream, 
                         adfGeoTransform[0] 
                         + adfGeoTransform[1] * (iX + 0.5)
                         + adfGeoTransform[2] * (iY + 0.5),
                         adfGeoTransform[3] 
                         + adfGeoTransform[4] * (iX + 0.5)
                         + adfGeoTransform[5] * (iY + 0.5),
                         panData[iX] );
        }
    }

    free( panData );

/* -------------------------------------------------------------------- */
/*      Cleanup.                                                        */
/* -------------------------------------------------------------------- */
    if( bEnableTrim )
        DTEDPtStreamTrimEdgeOnlyTiles( pStream );

    if( nFillDist > 0 )
        DTEDFillPtStream( pStream, nFillDist );

    DTEDClosePtStream( pStream );
    GDALClose( hSrcDS );

    exit( 0 );
}

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

            GDALTranslateOptionsFree(psOptions);
            GDALTranslateOptionsForBinaryFree(psOptionsForBinary);

            GDALDestroyDriverManager();
            exit(1);
        }
    }

/* -------------------------------------------------------------------- */
/*      Attempt to open source file.                                    */
/* -------------------------------------------------------------------- */

    GDALDatasetH hDataset =
        GDALOpenEx(psOptionsForBinary->pszSource,
                   GDAL_OF_RASTER | GDAL_OF_VERBOSE_ERROR, nullptr,
                   psOptionsForBinary->papszOpenOptions, nullptr);

    if( hDataset == nullptr )
    {
        GDALDestroyDriverManager();
        exit( 1 );
    }

/* -------------------------------------------------------------------- */
/*      Handle subdatasets.                                             */
/* -------------------------------------------------------------------- */
    if( !psOptionsForBinary->bCopySubDatasets
        && GDALGetRasterCount(hDataset) == 0
        && CSLCount(GDALGetMetadata( hDataset, "SUBDATASETS" )) > 0 )
    {
        fprintf( stderr,
                 "Input file contains subdatasets. Please, select one of them for reading.\n" );
        GDALClose( hDataset );
        GDALDestroyDriverManager();
        exit( 1 );
    }

    int bUsageError = FALSE;
    GDALDatasetH hOutDS = nullptr;
    if( psOptionsForBinary->bCopySubDatasets &&
        CSLCount(GDALGetMetadata( hDataset, "SUBDATASETS" )) > 0 )
    {
        char **papszSubdatasets = GDALGetMetadata(hDataset,"SUBDATASETS");
        char *pszSubDest = static_cast<char *>(
            CPLMalloc(strlen(psOptionsForBinary->pszDest) + 32));

        CPLString osPath = CPLGetPath(psOptionsForBinary->pszDest);
        CPLString osBasename = CPLGetBasename(psOptionsForBinary->pszDest);
        CPLString osExtension = CPLGetExtension(psOptionsForBinary->pszDest);
        CPLString osTemp;

        const char* pszFormat = nullptr;
        if ( CSLCount(papszSubdatasets)/2 < 10 )
        {
            pszFormat = "%s_%d";
        }
        else if ( CSLCount(papszSubdatasets)/2 < 100 )
        {
            pszFormat = "%s_%002d";
        }
        else
        {
            pszFormat = "%s_%003d";
        }

int FindSRS( const char *pszInput, OGRSpatialReference &oSRS )

{
    int            bGotSRS = FALSE;
    VSILFILE      *fp = NULL;
    GDALDataset	  *poGDALDS = NULL; 
    OGRLayer      *poLayer = NULL;
    const char    *pszProjection = NULL;
    CPLErrorHandler oErrorHandler = NULL;
    int bIsFile = FALSE;
    OGRErr eErr = OGRERR_NONE;
    int bDebug  = FALSE;

    /* temporarily suppress error messages we may get from xOpen() */
    bDebug = CSLTestBoolean(CPLGetConfigOption("CPL_DEBUG", "OFF"));
    if ( ! bDebug )
        oErrorHandler = CPLSetErrorHandler ( CPLQuietErrorHandler );

    /* Test if argument is a file */
    fp = VSIFOpenL( pszInput, "r" );
    if ( fp )  {
        bIsFile = TRUE;
        VSIFCloseL( fp );
        CPLDebug( "gdalsrsinfo", "argument is a file" );
    } 
       
    /* try to open with GDAL */
    if( strncmp(pszInput, "http://spatialreference.org/",
                strlen("http://spatialreference.org/")) != 0 )
    {
        CPLDebug( "gdalsrsinfo", "trying to open with GDAL" );
        poGDALDS = (GDALDataset *) GDALOpenEx( pszInput, 0, NULL, NULL, NULL );
    }
    if ( poGDALDS != NULL ) {
        pszProjection = poGDALDS->GetProjectionRef( );
        if( pszProjection != NULL && pszProjection[0] != '\0' )
        {
            char* pszProjectionTmp = (char*) pszProjection;
            if( oSRS.importFromWkt( &pszProjectionTmp ) == OGRERR_NONE ) {
                CPLDebug( "gdalsrsinfo", "got SRS from GDAL" );
                bGotSRS = TRUE;
            }
        }
        else if( poGDALDS->GetLayerCount() > 0 )
        {
            poLayer = poGDALDS->GetLayer( 0 );
            if ( poLayer != NULL ) {
                OGRSpatialReference *poSRS = poLayer->GetSpatialRef( );
                if ( poSRS != NULL ) {
                    CPLDebug( "gdalsrsinfo", "got SRS from OGR" );
                    bGotSRS = TRUE;
                    OGRSpatialReference* poSRSClone = poSRS->Clone();
                    oSRS = *poSRSClone;
                    OGRSpatialReference::DestroySpatialReference( poSRSClone );
                }
            }
        }
        GDALClose( (GDALDatasetH) poGDALDS );
        if ( ! bGotSRS ) 
            CPLDebug( "gdalsrsinfo", "did not open with GDAL" );
    }    
    
    /* Try ESRI file */
    if ( ! bGotSRS && bIsFile && (strstr(pszInput,".prj") != NULL) ) {
        CPLDebug( "gdalsrsinfo", 
                  "trying to get SRS from ESRI .prj file [%s]", pszInput );

        char **pszTemp;
        if ( strstr(pszInput,"ESRI::") != NULL )
            pszTemp = CSLLoad( pszInput+6 );
        else 
            pszTemp = CSLLoad( pszInput );

        if( pszTemp ) {
            eErr = oSRS.importFromESRI( pszTemp );
            CSLDestroy( pszTemp );
        }
        else 
            eErr = OGRERR_UNSUPPORTED_SRS;

        if( eErr != OGRERR_NONE ) {
            CPLDebug( "gdalsrsinfo", "did not get SRS from ESRI .prj file" );
        }
        else {
            CPLDebug( "gdalsrsinfo", "got SRS from ESRI .prj file" );
            bGotSRS = TRUE;
        }
    }

    /* Last resort, try OSRSetFromUserInput() */
    if ( ! bGotSRS ) {
        CPLDebug( "gdalsrsinfo", 
                  "trying to get SRS from user input [%s]", pszInput );

        eErr = oSRS.SetFromUserInput( pszInput );
 
       if(  eErr != OGRERR_NONE ) {
            CPLDebug( "gdalsrsinfo", "did not get SRS from user input" );
        }
        else {
//.........这里部分代码省略.........

//.........这里部分代码省略.........
        exit(1);
    }

    if (strcmp(pszDest, "/vsistdout/") == 0)
    {
        bQuiet      = TRUE;
        pfnProgress = GDALDummyProgress;
    }

/* -------------------------------------------------------------------- */
/*      Attempt to open source file.                                    */
/* -------------------------------------------------------------------- */

    hDataset = GDALOpenShared(pszSource, GA_ReadOnly);

    if (hDataset == NULL)
    {
        fprintf(stderr,
                "GDALOpen failed - %d\n%s\n",
                CPLGetLastErrorNo(), CPLGetLastErrorMsg());
        GDALDestroyDriverManager();
        exit(1);
    }

/* -------------------------------------------------------------------- */
/*      Handle subdatasets.                                             */
/* -------------------------------------------------------------------- */
    if (!bCopySubDatasets
        && CSLCount(GDALGetMetadata(hDataset, "SUBDATASETS")) > 0
        && GDALGetRasterCount(hDataset) == 0)
    {
        fprintf(stderr,
                "Input file contains subdatasets. Please, select one of them for reading.\n");
        GDALClose(hDataset);
        GDALDestroyDriverManager();
        exit(1);
    }

    if (CSLCount(GDALGetMetadata(hDataset, "SUBDATASETS")) > 0
        && bCopySubDatasets)
    {
        char **papszSubdatasets = GDALGetMetadata(hDataset, "SUBDATASETS");
        char *pszSubDest        = (char*) CPLMalloc(strlen(pszDest) + 32);
        int  i;
        int  bOldSubCall    = bSubCall;
        char **papszDupArgv = CSLDuplicate(argv);
        int  nRet           = 0;

        CPLFree(papszDupArgv[iDstFileArg]);
        papszDupArgv[iDstFileArg] = pszSubDest;
        bSubCall                  = TRUE;

        for (i = 0; papszSubdatasets[i] != NULL; i += 2)
        {
            CPLFree(papszDupArgv[iSrcFileArg]);
            papszDupArgv[iSrcFileArg] = CPLStrdup(strstr(papszSubdatasets[i], "=") + 1);
            sprintf(pszSubDest, "%s%d", pszDest, i / 2 + 1);
            nRet = ProxyMain(argc, papszDupArgv);
            if (nRet != 0)
                break;
        }

        CSLDestroy(papszDupArgv);

        bSubCall = bOldSubCall;
        CSLDestroy(argv);

//.........这里部分代码省略.........
                ahLayers.push_back(hLayer);

                hDstDS = CreateOutputDataset(ahLayers, hSRS,
                                 bGotBounds, sEnvelop,
                                 hDriver, pszDstFilename,
                                 nXSize, nYSize, dfXRes, dfYRes,
                                 bTargetAlignedPixels,
                                 anBandList.size(), eOutputType,
                                 papszCreateOptions, adfInitVals,
                                 bNoDataSet, dfNoData);
            }

            ProcessLayer( hLayer, hSRS != NULL, hDstDS, anBandList, 
                          adfBurnValues, b3D, bInverse, pszBurnAttribute,
                          papszRasterizeOptions, pfnProgress, NULL );

            OGR_DS_ReleaseResultSet( hSrcDS, hLayer );
        }
    }

/* -------------------------------------------------------------------- */
/*      Create output file if necessary.                                */
/* -------------------------------------------------------------------- */
    int nLayerCount = CSLCount(papszLayers);

    if (bCreateOutput && hDstDS == NULL)
    {
        std::vector<OGRLayerH> ahLayers;

        for( i = 0; i < nLayerCount; i++ )
        {
            OGRLayerH hLayer = OGR_DS_GetLayerByName( hSrcDS, papszLayers[i] );
            if( hLayer == NULL )
            {
                continue;
            }
            ahLayers.push_back(hLayer);
        }

        hDstDS = CreateOutputDataset(ahLayers, hSRS,
                                bGotBounds, sEnvelop,
                                hDriver, pszDstFilename,
                                nXSize, nYSize, dfXRes, dfYRes,
                                bTargetAlignedPixels,
                                anBandList.size(), eOutputType,
                                papszCreateOptions, adfInitVals,
                                bNoDataSet, dfNoData);
    }

/* -------------------------------------------------------------------- */
/*      Process each layer.                                             */
/* -------------------------------------------------------------------- */

    for( i = 0; i < nLayerCount; i++ )
    {
        OGRLayerH hLayer = OGR_DS_GetLayerByName( hSrcDS, papszLayers[i] );
        if( hLayer == NULL )
        {
            fprintf( stderr, "Unable to find layer %s, skipping.\n", 
                      papszLayers[i] );
            continue;
        }

        if( pszWHERE )
        {
            if( OGR_L_SetAttributeFilter( hLayer, pszWHERE ) != OGRERR_NONE )
                break;
        }

        void *pScaledProgress;
        pScaledProgress =
            GDALCreateScaledProgress( 0.0, 1.0 * (i + 1) / nLayerCount,
                                      pfnProgress, NULL );

        ProcessLayer( hLayer, hSRS != NULL, hDstDS, anBandList, 
                      adfBurnValues, b3D, bInverse, pszBurnAttribute,
                      papszRasterizeOptions, GDALScaledProgress, pScaledProgress );

        GDALDestroyScaledProgress( pScaledProgress );
    }

/* -------------------------------------------------------------------- */
/*      Cleanup                                                         */
/* -------------------------------------------------------------------- */

    OGR_DS_Destroy( hSrcDS );
    GDALClose( hDstDS );

    OSRDestroySpatialReference(hSRS);

    CSLDestroy( argv );
    CSLDestroy( papszRasterizeOptions );
    CSLDestroy( papszLayers );
    CSLDestroy( papszCreateOptions );
    
    GDALDestroyDriverManager();
    OGRCleanupAll();

    return 0;
}

//.........这里部分代码省略.........
        if (hTileIndexDS != NULL)
	  {
            hLayer = OGR_DS_GetLayer(hTileIndexDS, 0);
	  }
      }
    else
      {
        OGRSFDriverH hDriver;
        const char* pszDriverName = "ESRI Shapefile";

        fprintf( stderr,"Creating new index file...\n" );
        hDriver = OGRGetDriverByName( pszDriverName );
        if( hDriver == NULL )
	  {
	    char msg[1024];
            sprintf( msg, "%s driver not available.", pszDriverName );
	    setMapInMaps(conf,"lenv","message",msg);
	    return SERVICE_FAILED;
	  }

        hTileIndexDS = OGR_Dr_CreateDataSource( hDriver, index_filename, NULL );
        if (hTileIndexDS)
	  {
            char* pszLayerName = CPLStrdup(CPLGetBasename(index_filename));
            OGRSpatialReferenceH hSpatialRef = NULL;
            GDALDatasetH hDS = GDALOpen( index_filename, GA_ReadOnly );
            if (hDS)
	      {
                const char* pszWKT = GDALGetProjectionRef(hDS);
                if (pszWKT != NULL && pszWKT[0] != '\0')
		  {
                    hSpatialRef = OSRNewSpatialReference(pszWKT);
		  }
                GDALClose(hDS);
	      }

	    DIR *dirp = opendir(pszDataDir);
	    if(dirp==NULL){
	      char tmp1[1024];
	      sprintf(tmp1,_ss("The specified path %s doesn't exist."),pszDataDir);
	      setMapInMaps(conf,"lenv","message",tmp1);
	      return SERVICE_FAILED;
	    }
	    char *res=NULL;
	    struct dirent *dp;
	    tmpMap=getMapFromMaps(inputs,"ext","value");
	    char *ext;
	    if(tmpMap!=NULL)
	      ext=tmpMap->value;
	    while ((dp = readdir(dirp)) != NULL){
	      if(strncmp(dp->d_name,".",1)!=0&&strncmp(dp->d_name,"..",2)!=0){
		char* argv=(char*) malloc((strlen(dp->d_name)+strlen(pszDataDir)+2)*sizeof(char));
		sprintf(argv,"%s/%s",pszDataDir,dp->d_name);
		GDALDatasetH hDS0 = GDALOpen( argv, GA_ReadOnly );
		if (hDS0)
		  {
		    const char* pszWKT = GDALGetProjectionRef(hDS0);
		    fprintf(stderr,"SRS %s \n",pszWKT);
		    if (pszWKT != NULL && pszWKT[0] != '\0')
		      {
			if (hSpatialRef)
			  OSRRelease(hSpatialRef);
			hSpatialRef = OSRNewSpatialReference(pszWKT);
		      }
		    GDALClose(hDS);
		  }