void object::test<7>()
    {
        OGRErr err = OGRERR_NONE;

        // Read feature without geometry
        std::string tmp(data_tmp_);
        tmp += SEP;
        tmp += "tpoly.shp";
        OGRDataSourceH ds = OGR_Dr_Open(drv_, tmp.c_str(), false);
        ensure("Can't open layer", NULL != ds);

        OGRLayerH lyr = OGR_DS_GetLayer(ds, 0);
        ensure("Can't get layer", NULL != lyr);

        err = OGR_L_SetAttributeFilter(lyr, "PRFEDEA = 'nulled'");
        ensure_equals("Can't set attribute filter", OGRERR_NONE, err);

        // Fetch feature without geometry
        OGRFeatureH featNonSpatial = OGR_L_GetNextFeature(lyr);
        ensure("Didnt get feature with null geometry back", NULL != featNonSpatial);

        // Null geometry is expected
        OGRGeometryH nonGeom = OGR_F_GetGeometryRef(featNonSpatial);
        ensure("Didnt get null geometry as expected", NULL == nonGeom);

        OGR_F_Destroy(featNonSpatial);
        OGR_DS_Destroy(ds);
    }

void GDALWPrintRecords(GDALWConnection conn) {
	char * wkt;
	int i;
	OGRFeatureH feature;
	OGRGeometryH geometry;
	OGRFeatureDefnH featureDefn;
	featureDefn = OGR_L_GetLayerDefn(conn.layer);
	OGR_L_ResetReading(conn.layer);
	while( (feature = OGR_L_GetNextFeature(conn.layer)) != NULL ) {
		for(i = 0; i < OGR_FD_GetFieldCount(featureDefn); i++ ) {
			OGRFieldDefnH hFieldDefn = OGR_FD_GetFieldDefn( featureDefn, i );
		    if( OGR_Fld_GetType(hFieldDefn) == OFTInteger )
		    	printf( "%d,", OGR_F_GetFieldAsInteger( feature, i ) );
		    else if( OGR_Fld_GetType(hFieldDefn) == OFTReal )
		        printf( "%.3f,", OGR_F_GetFieldAsDouble( feature, i) );
		    else
		    	printf( "%s,", OGR_F_GetFieldAsString( feature, i) );

		}
		geometry = OGR_F_GetGeometryRef(feature);
		OGR_G_ExportToWkt(geometry, &wkt);
		printf("%s", wkt);
		printf("\n");
		CPLFree(wkt);
		OGR_F_Destroy(feature);
	}
}

    void object::test<5>()
    {
        // Original shapefile
        std::string orig(data_);
        orig += SEP;
        orig += "poly.shp";
        OGRDataSourceH dsOrig = OGR_Dr_Open(drv_, orig.c_str(), false);
        ensure("Can't open layer", NULL != dsOrig);

        OGRLayerH lyrOrig = OGR_DS_GetLayer(dsOrig, 0);
        ensure("Can't get layer", NULL != lyrOrig);

        // Copied shapefile
        std::string tmp(data_tmp_);
        tmp += SEP;
        tmp += "tpoly.shp";
        OGRDataSourceH dsTmp = OGR_Dr_Open(drv_, tmp.c_str(), false);
        ensure("Can't open layer", NULL != dsTmp);

        OGRLayerH lyrTmp = OGR_DS_GetLayer(dsTmp, 0);
        ensure("Can't get layer", NULL != lyrTmp);

        // Iterate through features and compare geometries
        OGRFeatureH featOrig = OGR_L_GetNextFeature(lyrOrig);
        OGRFeatureH featTmp = OGR_L_GetNextFeature(lyrTmp);

        while (NULL != featOrig && NULL != featTmp)
        {
            OGRGeometryH lhs = OGR_F_GetGeometryRef(featOrig);
            OGRGeometryH rhs = OGR_F_GetGeometryRef(featTmp);

            ensure_equal_geometries(lhs, rhs, 0.000000001);

            // TODO: add ensure_equal_attributes()

            OGR_F_Destroy(featOrig);
            OGR_F_Destroy(featTmp);

            // Move to next feature
            featOrig = OGR_L_GetNextFeature(lyrOrig);
            featTmp = OGR_L_GetNextFeature(lyrTmp);
        }

        OGR_DS_Destroy(dsOrig);
        OGR_DS_Destroy(dsTmp);
    }

/**
 * \brief Convenience function to check if a geometry is contained in a OGR
 *        datasource for a given layer.
 *
 * The passed geometry is a wkt representation of a geometry of type GeomType.
 * pszFile is opened, and the passed geometry is queried against all
 * geometries in pszLayer.  If the passed geometry is contained in *any* of the
 * geomtries in the layer, TRUE is returned.  FALSE is returned otherwise,
 * including errors.  The SRS of all geometries is assumed to be the same.
 *
 * \param pszWkt Well-known text representation of a geometry.
 * \param pszFile File to open
 * \param pszLayer Layer to extract geometry from, if NULL, use layer 0.
 * \return TRUE if pszWkt is contained in any geometry in pszLayer, FALSE
 *         otherwise, include errors
 */
int NinjaOGRContain(const char *pszWkt, const char *pszFile,
                    const char *pszLayer)
{
    int bContains = FALSE;
    if( pszWkt == NULL || pszFile == NULL )
    {
        return FALSE;
    }
    CPLDebug( "WINDNINJA", "Checking for containment of %s in %s:%s",
              pszWkt, pszFile, pszLayer ? pszLayer : "" );
    OGRGeometryH hTestGeometry = NULL;
    int err = OGR_G_CreateFromWkt( (char**)&pszWkt, NULL, &hTestGeometry );
    if( hTestGeometry == NULL || err != CE_None )
    {
        return FALSE;
    }
    OGRDataSourceH hDS = OGROpen( pszFile, 0, NULL );
    if( hDS == NULL )
    {
        CPLError( CE_Failure, CPLE_AppDefined,
                  "Failed to open datasource: %s", pszFile );
        OGR_G_DestroyGeometry( hTestGeometry );
        bContains = FALSE;
        return bContains;
    }
    OGRLayerH hLayer;
    if( pszLayer == NULL )
    {
        hLayer = OGR_DS_GetLayer( hDS, 0 );
    }
    else
    {
        hLayer = OGR_DS_GetLayerByName( hDS, pszLayer );
    }
    OGRFeatureH hFeature;
    if( hLayer != NULL )
    {
        OGRGeometryH hGeometry;
        OGR_L_ResetReading( hLayer );
        while( ( hFeature = OGR_L_GetNextFeature( hLayer ) ) != NULL )
        {
            hGeometry = OGR_F_GetGeometryRef( hFeature );
            if( OGR_G_Contains( hGeometry, hTestGeometry ) )
            {
                bContains = TRUE;
                OGR_F_Destroy( hFeature );
                break;
            }
            OGR_F_Destroy( hFeature );
        }
    }
    OGR_G_DestroyGeometry( hTestGeometry );
    OGR_DS_Destroy( hDS );
    return bContains;
}

bool QgsOgrFeatureIterator::fetchFeature( QgsFeature& feature )
{
  feature.setValid( false );

  if ( mClosed || !ogrLayer )
    return false;

  if ( mRequest.filterType() == QgsFeatureRequest::FilterFid )
  {
    bool result = fetchFeatureWithId( mRequest.filterFid(), feature );
    close(); // the feature has been read or was not found: we have finished here
    return result;
  }
  else if ( mRequest.filterType() == QgsFeatureRequest::FilterFids )
  {
    while ( mFilterFidsIt != mFilterFids.constEnd() )
    {
      QgsFeatureId nextId = *mFilterFidsIt;
      mFilterFidsIt++;

      if ( fetchFeatureWithId( nextId, feature ) )
        return true;
    }
    close();
    return false;
  }

  OGRFeatureH fet;

  while (( fet = OGR_L_GetNextFeature( ogrLayer ) ) )
  {
    if ( !readFeature( fet, feature ) )
      continue;
    else
      OGR_F_Destroy( fet );

    if ( !mRequest.filterRect().isNull() && !feature.hasGeometry() )
      continue;

    // we have a feature, end this cycle
    feature.setValid( true );
    return true;

  } // while

  close();
  return false;
}

bool QgsOgrFeatureIterator::fetchFeature( QgsFeature& feature )
{
  feature.setValid( false );

  if ( mClosed )
    return false;

  if ( !P->mRelevantFieldsForNextFeature )
    ensureRelevantFields();

  if ( mRequest.filterType() == QgsFeatureRequest::FilterFid )
  {
    OGRFeatureH fet = OGR_L_GetFeature( ogrLayer, FID_TO_NUMBER( mRequest.filterFid() ) );
    if ( !fet )
    {
      close();
      return false;
    }

    if ( readFeature( fet, feature ) )
      OGR_F_Destroy( fet );

    feature.setValid( true );
    close(); // the feature has been read: we have finished here
    return true;
  }

  OGRFeatureH fet;

  while (( fet = OGR_L_GetNextFeature( ogrLayer ) ) )
  {
    if ( !readFeature( fet, feature ) )
      continue;

    // we have a feature, end this cycle
    feature.setValid( true );
    OGR_F_Destroy( fet );
    return true;

  } // while

  QgsDebugMsg( "Feature is null" );

  close();
  return false;
}

bool QgsOgrFeatureIterator::fetchFeature( QgsFeature& feature )
{
  feature.setValid( false );

  if ( mClosed )
    return false;

  if ( mRequest.filterType() == QgsFeatureRequest::FilterFid )
  {
    OGRFeatureH fet = OGR_L_GetFeature( ogrLayer, FID_TO_NUMBER( mRequest.filterFid() ) );
    if ( !fet )
    {
      close();
      return false;
    }

    if ( readFeature( fet, feature ) )
      OGR_F_Destroy( fet );

    feature.setValid( true );
    close(); // the feature has been read: we have finished here
    return true;
  }

  OGRFeatureH fet;

  while (( fet = OGR_L_GetNextFeature( ogrLayer ) ) )
  {
    if ( !readFeature( fet, feature ) )
      continue;

    if ( !mRequest.filterRect().isNull() && !feature.constGeometry() )
      continue;

    // we have a feature, end this cycle
    feature.setValid( true );
    OGR_F_Destroy( fet );
    return true;

  } // while

  close();
  return false;
}

QgsFeatureList QgsOgrUtils::stringToFeatureList( const QString& string, const QgsFields& fields, QTextCodec* encoding )
{
  QgsFeatureList features;
  if ( string.isEmpty() )
    return features;

  QString randomFileName = QString( "/vsimem/%1" ).arg( QUuid::createUuid().toString() );

  // create memory file system object from string buffer
  QByteArray ba = string.toUtf8();
  VSIFCloseL( VSIFileFromMemBuffer( TO8( randomFileName ), reinterpret_cast< GByte* >( ba.data() ),
                                    static_cast< vsi_l_offset >( ba.size() ), FALSE ) );

  OGRDataSourceH hDS = OGROpen( TO8( randomFileName ), false, nullptr );
  if ( !hDS )
  {
    VSIUnlink( TO8( randomFileName ) );
    return features;
  }

  OGRLayerH ogrLayer = OGR_DS_GetLayer( hDS, 0 );
  if ( !ogrLayer )
  {
    OGR_DS_Destroy( hDS );
    VSIUnlink( TO8( randomFileName ) );
    return features;
  }

  OGRFeatureH oFeat;
  while (( oFeat = OGR_L_GetNextFeature( ogrLayer ) ) )
  {
    QgsFeature feat = readOgrFeature( oFeat, fields, encoding );
    if ( feat.isValid() )
      features << feat;

    OGR_F_Destroy( oFeat );
  }

  OGR_DS_Destroy( hDS );
  VSIUnlink( TO8( randomFileName ) );

  return features;
}

void OgrFileImport::importLayer(MapPart* map_part, OGRLayerH layer)
{
	Q_ASSERT(map_part);
	
	auto feature_definition = OGR_L_GetLayerDefn(layer);
	
	OGR_L_ResetReading(layer);
	while (auto feature = ogr::unique_feature(OGR_L_GetNextFeature(layer)))
	{
		auto geometry = OGR_F_GetGeometryRef(feature.get());
		if (!geometry || OGR_G_IsEmpty(geometry))
		{
			++empty_geometries;
			continue;
		}
		
		OGR_G_FlattenTo2D(geometry);
		importFeature(map_part, feature_definition, feature.get(), geometry);
	}
}

QgsFeatureList QgsOgrUtils::stringToFeatureList( const QString &string, const QgsFields &fields, QTextCodec *encoding )
{
  QgsFeatureList features;
  if ( string.isEmpty() )
    return features;

  QString randomFileName = QStringLiteral( "/vsimem/%1" ).arg( QUuid::createUuid().toString() );

  // create memory file system object from string buffer
  QByteArray ba = string.toUtf8();
  VSIFCloseL( VSIFileFromMemBuffer( randomFileName.toUtf8().constData(), reinterpret_cast< GByte * >( ba.data() ),
                                    static_cast< vsi_l_offset >( ba.size() ), FALSE ) );

  gdal::ogr_datasource_unique_ptr hDS( OGROpen( randomFileName.toUtf8().constData(), false, nullptr ) );
  if ( !hDS )
  {
    VSIUnlink( randomFileName.toUtf8().constData() );
    return features;
  }

  OGRLayerH ogrLayer = OGR_DS_GetLayer( hDS.get(), 0 );
  if ( !ogrLayer )
  {
    hDS.reset();
    VSIUnlink( randomFileName.toUtf8().constData() );
    return features;
  }

  gdal::ogr_feature_unique_ptr oFeat;
  while ( oFeat.reset( OGR_L_GetNextFeature( ogrLayer ) ), oFeat )
  {
    QgsFeature feat = readOgrFeature( oFeat.get(), fields, encoding );
    if ( feat.isValid() )
      features << feat;
  }

  hDS.reset();
  VSIUnlink( randomFileName.toUtf8().constData() );

  return features;
}

bool CUtils::insideInPolygons(OGRDataSourceH poDS, double x, double y)
{
	bool res = false;
	OGRGeometryH pt = OGR_G_CreateGeometry(wkbPoint);
	OGR_G_AddPoint_2D(pt, x, y);

	for(int iLayer = 0; iLayer < OGR_DS_GetLayerCount(poDS); iLayer++)	
	{
		OGRLayerH poLayer = OGR_DS_GetLayer(poDS, iLayer);
		if(poLayer!=NULL)
		{
			OGREnvelope layerBounds;
			OGR_L_GetExtent(poLayer, &layerBounds, 1);
			
			if(	(layerBounds.MinX <= x) && (layerBounds.MinY <= y) && 
				(layerBounds.MaxX >= x) && (layerBounds.MaxY >= y) )
			{
				OGR_L_ResetReading(poLayer);
				if(OGR_FD_GetGeomType( OGR_L_GetLayerDefn(poLayer) ) == wkbPolygon)
				{
					OGRFeatureH poFeat;
					while((poFeat = OGR_L_GetNextFeature(poLayer))!= NULL)
					{
						OGRGeometryH hGeom = OGR_F_GetGeometryRef(poFeat);
						if(OGR_G_Within(pt, hGeom))
						{
							res = true;
							break;
						}
					}
					if(res) { OGR_L_ResetReading(poLayer); break; }
				}
			}
		}
	}	
	OGR_G_DestroyGeometry(pt);
	return res;
}

    void object::test<9>()
    {
        // Open directory as a datasource
        OGRDataSourceH ds = OGR_Dr_Open(drv_, data_tmp_ .c_str(), false);
        ensure("Can't open datasource", NULL != ds);

        std::string sql("select * from tpoly where prfedea = '35043413'");
        OGRLayerH lyr = OGR_DS_ExecuteSQL(ds, sql.c_str(), NULL, NULL);
        ensure("Can't create layer from query", NULL != lyr);

        // Prepare tester collection
        std::vector<std::string> list;
        list.push_back("35043413");
       
        // Test attributes
        ensure_equal_attributes(lyr, "prfedea", list);

        // Test geometry
        const char* wkt = "POLYGON ((479750.688 4764702.000,479658.594 4764670.000,"
            "479640.094 4764721.000,479735.906 4764752.000,"
            "479750.688 4764702.000))";

        OGRGeometryH testGeom = NULL;
        OGRErr err = OGR_G_CreateFromWkt((char**) &wkt, NULL, &testGeom);
        ensure_equals("Can't create geometry from WKT", OGRERR_NONE, err);

        OGR_L_ResetReading(lyr);
        OGRFeatureH feat = OGR_L_GetNextFeature(lyr);
        ensure("Can't featch feature", NULL != feat);

        ensure_equal_geometries(OGR_F_GetGeometryRef(feat), testGeom, 0.001);

        OGR_F_Destroy(feat);
        OGR_G_DestroyGeometry(testGeom);
        OGR_DS_ReleaseResultSet(ds, lyr);
        OGR_DS_Destroy(ds);
    }

/* OGRFeatureH OGR_L_GetNextFeature(OGRLayerH hLayer)

{ok, DataSource} = lgeo_ogr:open("test/polygon.shp"),
{ok, Layer} = lgeo_ogr:ds_get_layer(DataSource, 0),
{ok, Feature} = lgeo_ogr:l_get_next_feature(Layer).

*/
static ERL_NIF_TERM
l_get_next_feature(ErlNifEnv* env, int argc, const ERL_NIF_TERM argv[])
{
    EnvLayer_t **layer;
    ERL_NIF_TERM eterm;

    if(argc != 1) {
        return enif_make_badarg(env);
    }

    if(!enif_get_resource(env, argv[0], OGR_L_RESOURCE, (void**)&layer)) {
        return enif_make_badarg(env);
    }

    OGRFeatureH feat = OGR_L_GetNextFeature((**layer).obj);
    if(feat == NULL) {
        eterm = enif_make_string(env, "No more features", ERL_NIF_LATIN1);
        return enif_make_tuple2(env, enif_make_atom(env, "error"), eterm);
    }

    EnvFeature_t **feature = \
        enif_alloc_resource(OGR_F_RESOURCE, sizeof(EnvFeature_t*));

    ErlNifEnv *feature_env = enif_alloc_env();

    *feature = (EnvFeature_t*) enif_alloc(sizeof(EnvFeature_t));
    (**feature).env = feature_env;
    (**feature).obj = feat;

    // Save copy of layer so is not garbage collected
    enif_make_copy(feature_env, argv[0]);

    eterm = enif_make_resource(env, feature);
    enif_release_resource(feature);
    return enif_make_tuple2(env, enif_make_atom(env, "ok"), eterm);
}

//.........这里部分代码省略.........
    const char *pszTmpFilename =CPLFormFilename( pszDataPath, "irs/FWA", ".dat" );
    std::vector<CFWA>fwas = LoadFwas( pszTmpFilename );
    
    int nFeatures = OGR_L_GetFeatureCount( hInputLayer, TRUE );
    FILE *fout = fopen( pszOutputfile, "w" );

    //CFWA *fwa;

    Random random;

    char pszDb[8192];
    sprintf( pszDb, "%s/omffr.sqlite", pszDataPath );

    IRSDataAccess *poDA = IRSDataAccess::Create( 0, pszDb );
    int rc;
    sqlite3 *db;
    rc = sqlite3_open_v2( pszDb, &db, SQLITE_OPEN_READONLY, NULL );
    rc = sqlite3_enable_load_extension( db, 1 );
    rc = sqlite3_load_extension( db, "/usr/local/lib/libspatialite.so", 0, NULL );
    sqlite3_stmt *stmt;
    rc = sqlite3_prepare_v2( db, "SELECT * from fwa join fwa_bndry USING(fwa_gis_id) " \
                                 "WHERE ST_Contains(fwa_bndry.geometry, MakePoint(?, ?, 4269))",
                                 -1, &stmt, NULL );

    if(rc)
    {
        CPLError( CE_Failure, CPLE_AppDefined, "Could not open DB");
    }

    GDALTermProgress( 0.0, NULL, NULL );
    OGR_L_ResetReading( hInputLayer );
    const char *pszFwaName;
    int nDone = 0;
    while( ( hInputFeature = OGR_L_GetNextFeature( hInputLayer ) ) != NULL )
    {
        /*
        fwaid = OGR_F_GetFieldAsString( hInputFeature,
                                        OGR_FD_GetFieldIndex( hInputFeatureDefn,
                                                              "fwa_name" ) );
        abyFwa = CPLStrdup( fwaid );
        LaunderFwaName( abyFwa );
        fwa = FindFwa( fwas, abyFwa );
        if( fwa == NULL )
        {
            CPLError( CE_Warning, CPLE_FileIO,
                    "Could not load fwa (%s)from file, missing\n", abyFwa );
            continue;
        }
        */
        /* Get fwa by point */
        hGeometry = OGR_F_GetGeometryRef( hInputFeature );
        /* Try to handle non-geometry types (csv) */
        if( hGeometry != NULL )
        {
            lat = OGR_G_GetY( hGeometry, 0 );
            lon = OGR_G_GetX( hGeometry, 0 );
        }
        else
        {
            lat = OGR_F_GetFieldAsDouble( hInputFeature, OGR_FD_GetFieldIndex( hInputFeatureDefn, 
                                                            "Y") );
            lon = OGR_F_GetFieldAsDouble( hInputFeature, OGR_FD_GetFieldIndex( hInputFeatureDefn, 
                                                            "X") );
        }

        std::string oFwaName = poDA->PointQuery( "fwa_bndry", "fwa_lndr_name",

static int RasterliteInsertSRID(OGRDataSourceH hDS, const char* pszWKT)
{
    CPLString osSQL;
    
    int nAuthorityCode = 0;
    CPLString osAuthorityName, osProjCS, osProj4;
    if (pszWKT != NULL && strlen(pszWKT) != 0)
    {
        OGRSpatialReferenceH hSRS = OSRNewSpatialReference(pszWKT);
        if (hSRS)
        {
            const char* pszAuthorityName = OSRGetAuthorityName(hSRS, NULL);
            if (pszAuthorityName) osAuthorityName = pszAuthorityName;
            
            const char* pszProjCS = OSRGetAttrValue(hSRS, "PROJCS", 0);
            if (pszProjCS) osProjCS = pszProjCS;
            
            const char* pszAuthorityCode = OSRGetAuthorityCode(hSRS, NULL);
            if (pszAuthorityCode) nAuthorityCode = atoi(pszAuthorityCode);
            
            char    *pszProj4 = NULL;
            if( OSRExportToProj4( hSRS, &pszProj4 ) != OGRERR_NONE )
                pszProj4 = CPLStrdup("");
            osProj4 = pszProj4;
            CPLFree(pszProj4);
        }
        OSRDestroySpatialReference(hSRS);
    }
        
    int nSRSId = -1;
    if (nAuthorityCode != 0 && osAuthorityName.size() != 0)
    {
        osSQL.Printf   ("SELECT srid FROM spatial_ref_sys WHERE auth_srid = %d", nAuthorityCode);
        OGRLayerH hLyr = OGR_DS_ExecuteSQL(hDS, osSQL.c_str(), NULL, NULL);
        if (hLyr == NULL)
        {
            nSRSId = nAuthorityCode;
            
            if ( osProjCS.size() != 0 )
                osSQL.Printf(
                    "INSERT INTO spatial_ref_sys "
                    "(srid, auth_name, auth_srid, ref_sys_name, proj4text) "
                    "VALUES (%d, '%s', '%d', '%s', '%s')",
                    nSRSId, osAuthorityName.c_str(),
                    nAuthorityCode, osProjCS.c_str(), osProj4.c_str() );
            else
                osSQL.Printf(
                    "INSERT INTO spatial_ref_sys "
                    "(srid, auth_name, auth_srid, proj4text) "
                    "VALUES (%d, '%s', '%d', '%s')",
                    nSRSId, osAuthorityName.c_str(),
                    nAuthorityCode, osProj4.c_str() );

            
            OGR_DS_ExecuteSQL(hDS, osSQL.c_str(), NULL, NULL);
        }
        else
        {
            OGRFeatureH hFeat = OGR_L_GetNextFeature(hLyr);
            if (hFeat)
            {
                nSRSId = OGR_F_GetFieldAsInteger(hFeat, 0);
                OGR_F_Destroy(hFeat);
            }
            OGR_DS_ReleaseResultSet(hDS, hLyr);
        }
    }
    
    return nSRSId;
}

OGRDataSourceH RasterliteCreateTables(OGRDataSourceH hDS, const char* pszTableName,
                                      int nSRSId, int bWipeExistingData)
{
    CPLString osSQL;
    
    CPLString osOldVal = CPLGetConfigOption("SQLITE_LIST_ALL_TABLES", "FALSE");
    CPLString osDBName = OGR_DS_GetName(hDS);
    
    CPLString osRasterLayer;
    osRasterLayer.Printf("%s_rasters", pszTableName);
    
    CPLString osMetatadataLayer;
    osMetatadataLayer.Printf("%s_metadata", pszTableName);

    OGRLayerH hLyr;

    if (OGR_DS_GetLayerByName(hDS, osRasterLayer.c_str()) == NULL)
    {
/* -------------------------------------------------------------------- */
/*      The table don't exist. Create them                              */
/* -------------------------------------------------------------------- */

        /* Create _rasters table */
        osSQL.Printf   ("CREATE TABLE \"%s\" ("
                        "id INTEGER NOT NULL PRIMARY KEY AUTOINCREMENT,"
                        "raster BLOB NOT NULL)", osRasterLayer.c_str());
        OGR_DS_ExecuteSQL(hDS, osSQL.c_str(), NULL, NULL);
            
        /* Create _metadata table */
        osSQL.Printf   ("CREATE TABLE \"%s\" ("
                        "id INTEGER NOT NULL PRIMARY KEY,"
                        "source_name TEXT NOT NULL,"
                        "tile_id INTEGER NOT NULL,"
                        "width INTEGER NOT NULL,"
                        "height INTEGER NOT NULL,"
                        "pixel_x_size DOUBLE NOT NULL,"
                        "pixel_y_size DOUBLE NOT NULL)",
                        osMetatadataLayer.c_str());
        OGR_DS_ExecuteSQL(hDS, osSQL.c_str(), NULL, NULL);

        /* Add geometry column to _metadata table */
        osSQL.Printf("SELECT AddGeometryColumn('%s', 'geometry', %d, 'POLYGON', 2)",
                      osMetatadataLayer.c_str(), nSRSId);
        if ((hLyr = OGR_DS_ExecuteSQL(hDS, osSQL.c_str(), NULL, NULL)) == NULL)
        {
            CPLError(CE_Failure, CPLE_AppDefined,
                     "Check that the OGR SQLite driver has Spatialite support");
            OGRReleaseDataSource(hDS);
            return NULL;
        }
        OGR_DS_ReleaseResultSet(hDS, hLyr);
                
        /* Create spatial index on _metadata table */
        osSQL.Printf("SELECT CreateSpatialIndex('%s', 'geometry')",
                      osMetatadataLayer.c_str());
        if ((hLyr = OGR_DS_ExecuteSQL(hDS, osSQL.c_str(), NULL, NULL)) == NULL)
        {
            OGRReleaseDataSource(hDS);
            return NULL;
        }
        OGR_DS_ReleaseResultSet(hDS, hLyr);
        
        /* Re-open the DB to take into account the new tables*/
        OGRReleaseDataSource(hDS);
        
        CPLSetConfigOption("SQLITE_LIST_ALL_TABLES", "TRUE");
        hDS = OGROpen(osDBName.c_str(), TRUE, NULL);
        CPLSetConfigOption("SQLITE_LIST_ALL_TABLES", osOldVal.c_str());
    }
    else
    {
        /* Check that the existing SRS is consistent with the one of the new */
        /* data to be inserted */
        osSQL.Printf("SELECT srid FROM geometry_columns WHERE f_table_name = '%s'",
                     osMetatadataLayer.c_str());
        hLyr = OGR_DS_ExecuteSQL(hDS, osSQL.c_str(), NULL, NULL);
        if (hLyr)
        {
            int nExistingSRID = -1;
            OGRFeatureH hFeat = OGR_L_GetNextFeature(hLyr);
            if (hFeat)
            {
                nExistingSRID = OGR_F_GetFieldAsInteger(hFeat, 0);
                OGR_F_Destroy(hFeat);
            }
            OGR_DS_ReleaseResultSet(hDS, hLyr);

            if (nExistingSRID != nSRSId)
            {
                if (bWipeExistingData)
                {
                    osSQL.Printf("UPDATE geometry_columns SET srid = %d "
                                 "WHERE f_table_name = \"%s\"",
                                 nSRSId, osMetatadataLayer.c_str());
                    OGR_DS_ExecuteSQL(hDS, osSQL.c_str(), NULL, NULL);
        
                    /* Re-open the DB to take into account the change of SRS */
                    OGRReleaseDataSource(hDS);
                    
                    CPLSetConfigOption("SQLITE_LIST_ALL_TABLES", "TRUE");
//.........这里部分代码省略.........

//.........这里部分代码省略.........
    if (hRasterLayer == NULL || hMetadataLayer == NULL)
    {
        CPLError(CE_Failure, CPLE_AppDefined,
                 "Cannot find metadata and/or raster tables");
        OGRReleaseDataSource(hDS);
        return NULL;
    }

/* -------------------------------------------------------------------- */
/*      Check if there is overlapping data and warn the user            */
/* -------------------------------------------------------------------- */
    double minx = adfGeoTransform[0];
    double maxx = adfGeoTransform[0] + nXSize * adfGeoTransform[1];
    double maxy = adfGeoTransform[3];
    double miny = adfGeoTransform[3] + nYSize * adfGeoTransform[5];
    
    osSQL.Printf("SELECT COUNT(geometry) FROM \"%s\" "
                 "WHERE rowid IN "
                 "(SELECT pkid FROM \"idx_%s_metadata_geometry\" "
                  "WHERE xmin < %.15f AND xmax > %.15f "
                  "AND ymin < %.15f  AND ymax > %.15f) "
                 "AND pixel_x_size >= %.15f AND pixel_x_size <= %.15f AND "
                 "pixel_y_size >= %.15f AND pixel_y_size <= %.15f",
                  osMetatadataLayer.c_str(),
                  osTableName.c_str(),
                  maxx, minx, maxy, miny,
                  adfGeoTransform[1] - 1e-15, adfGeoTransform[1] + 1e-15,
                  - adfGeoTransform[5] - 1e-15, - adfGeoTransform[5] + 1e-15);
    
    int nOverlappingGeoms = 0;
    OGRLayerH hCountLyr = OGR_DS_ExecuteSQL(hDS, osSQL.c_str(), NULL, NULL);
    if (hCountLyr)
    {
        OGRFeatureH hFeat = OGR_L_GetNextFeature(hCountLyr);
        if (hFeat)
        {
            nOverlappingGeoms = OGR_F_GetFieldAsInteger(hFeat, 0);
            OGR_F_Destroy(hFeat);
        }
        OGR_DS_ReleaseResultSet(hDS, hCountLyr);
    }
    
    if (nOverlappingGeoms != 0)
    {
        CPLError(CE_Warning, CPLE_AppDefined,
                 "Raster tiles already exist in the %s table within "
                 "the extent of the data to be inserted in",
                 osTableName.c_str());
    }
   
/* -------------------------------------------------------------------- */
/*      Iterate over blocks to add data into raster and metadata tables */
/* -------------------------------------------------------------------- */
    int nXBlocks = (nXSize + nBlockXSize - 1) / nBlockXSize;
    int nYBlocks = (nYSize + nBlockYSize - 1) / nBlockYSize;

    GDALDataType eDataType = poSrcDS->GetRasterBand(1)->GetRasterDataType();
    int nDataTypeSize = GDALGetDataTypeSize(eDataType) / 8;
    GByte* pabyMEMDSBuffer =
        (GByte*)VSIMalloc3(nBlockXSize, nBlockYSize, nBands * nDataTypeSize);
    if (pabyMEMDSBuffer == NULL)
    {
        OGRReleaseDataSource(hDS);
        return NULL;
    }
    

//.........这里部分代码省略.........
                    nOffset += 11 + nThisTRESize;

                    NITFDESFreeTREData(pabyTREData);
                }
            }

            /* Report info from location table, if found. */
            if( !bHasFoundLocationTable && nRPFDESOffset >= 0 )
            {
                int i;
                int nLocCount = 0;
                NITFLocation* pasLocations;

                VSIFSeekL(psFile->fp, psSegInfo->nSegmentStart + nRPFDESOffset, SEEK_SET);
                pasLocations = NITFReadRPFLocationTable(psFile->fp, &nLocCount);
                if (pasLocations)
                {
                    printf( "  Location Table\n" );
                    for( i = 0; i < nLocCount; i++ )
                    {
                        printf( "    LocName=%s, LocId=%d, Offset=%d, Size=%d\n",
                                GetLocationNameFromId(pasLocations[i].nLocId),
                                pasLocations[i].nLocId,
                                pasLocations[i].nLocOffset,
                                pasLocations[i].nLocSize );
                    }

                    CPLFree(pasLocations);
                    printf( "\n" );
                }
            }

            DumpMetadata( "  DES Metadata:", "    ", psDES->papszMetadata );

            if ( bExtractSHP && CSLFetchNameValue(psDES->papszMetadata, "NITF_SHAPE_USE") != NULL )
            {
                char szFilename[32];
                char szRadix[32];
                if (bExtractSHPInMem)
                    snprintf(szRadix, sizeof(szRadix), "/vsimem/nitf_segment_%d", iSegment + 1);
                else
                    snprintf(szRadix, sizeof(szRadix), "nitf_segment_%d", iSegment + 1);

                if (NITFDESExtractShapefile(psDES, szRadix))
                {
                    OGRDataSourceH hDS;
                    OGRRegisterAll();
                    snprintf(szFilename, sizeof(szFilename), "%s.SHP", szRadix);
                    hDS = OGROpen(szFilename, FALSE, NULL);
                    if (hDS)
                    {
                        int nGeom = 0;
                        OGRLayerH hLayer = OGR_DS_GetLayer(hDS, 0);
                        if (hLayer)
                        {
                            OGRFeatureH hFeat;
                            printf("\n");
                            while ( (hFeat = OGR_L_GetNextFeature(hLayer)) != NULL )
                            {
                                OGRGeometryH hGeom = OGR_F_GetGeometryRef(hFeat);
                                if (hGeom)
                                {
                                    char* pszWKT = NULL;
                                    OGR_G_ExportToWkt(hGeom, &pszWKT);
                                    if (pszWKT)
                                        printf("    Geometry %d : %s\n", nGeom ++, pszWKT);
                                    CPLFree(pszWKT);
                                }
                                OGR_F_Destroy(hFeat);
                            }
                        }
                        OGR_DS_Destroy(hDS);
                    }
                }

                if (bExtractSHPInMem)
                {
                    snprintf(szFilename, sizeof(szFilename), "%s.SHP", szRadix);
                    VSIUnlink(szFilename);
                    snprintf(szFilename, sizeof(szFilename), "%s.SHX", szRadix);
                    VSIUnlink(szFilename);
                    snprintf(szFilename, sizeof(szFilename), "%s.DBF", szRadix);
                    VSIUnlink(szFilename);
                }
            }
        }

/* -------------------------------------------------------------------- */
/*      Close.                                                          */
/* -------------------------------------------------------------------- */
        NITFClose( psFile );
    }

    CPLFinderClean();
    CPLCleanupTLS();
    VSICleanupFileManager();
    OGRCleanupAll();

    exit( 0 );
}

static CPLErr ProcessLayer(
    OGRLayerH hSrcLayer, int bSRSIsSet,
    GDALDatasetH hDstDS, std::vector<int> anBandList,
    const std::vector<double> &adfBurnValues, int b3D, int bInverse,
    const char *pszBurnAttribute, char **papszRasterizeOptions,
    GDALProgressFunc pfnProgress, void* pProgressData )

{
/* -------------------------------------------------------------------- */
/*      Checkout that SRS are the same.                                 */
/*      If -a_srs is specified, skip the test                           */
/* -------------------------------------------------------------------- */
    OGRCoordinateTransformationH hCT = NULL;
    if (!bSRSIsSet)
    {
        OGRSpatialReferenceH  hDstSRS = NULL;
        if( GDALGetProjectionRef( hDstDS ) != NULL )
        {
            char *pszProjection;

            pszProjection = (char *) GDALGetProjectionRef( hDstDS );

            hDstSRS = OSRNewSpatialReference(NULL);
            if( OSRImportFromWkt( hDstSRS, &pszProjection ) != OGRERR_NONE )
            {
                OSRDestroySpatialReference(hDstSRS);
                hDstSRS = NULL;
            }
        }

        OGRSpatialReferenceH hSrcSRS = OGR_L_GetSpatialRef(hSrcLayer);
        if( hDstSRS != NULL && hSrcSRS != NULL )
        {
            if( OSRIsSame(hSrcSRS, hDstSRS) == FALSE )
            {
                hCT = OCTNewCoordinateTransformation(hSrcSRS, hDstSRS);
                if( hCT == NULL )
                {
                    CPLError(CE_Warning, CPLE_AppDefined,
                        "The output raster dataset and the input vector layer do not have the same SRS.\n"
                        "And reprojection of input data did not work. Results might be incorrect.");
                }
            }
        }
        else if( hDstSRS != NULL && hSrcSRS == NULL )
        {
            CPLError(CE_Warning, CPLE_AppDefined,
                    "The output raster dataset has a SRS, but the input vector layer SRS is unknown.\n"
                    "Ensure input vector has the same SRS, otherwise results might be incorrect.");
        }
        else if( hDstSRS == NULL && hSrcSRS != NULL )
        {
            CPLError(CE_Warning, CPLE_AppDefined,
                    "The input vector layer has a SRS, but the output raster dataset SRS is unknown.\n"
                    "Ensure output raster dataset has the same SRS, otherwise results might be incorrect.");
        }

        if( hDstSRS != NULL )
        {
            OSRDestroySpatialReference(hDstSRS);
        }
    }

/* -------------------------------------------------------------------- */
/*      Get field index, and check.                                     */
/* -------------------------------------------------------------------- */
    int iBurnField = -1;

    if( pszBurnAttribute )
    {
        iBurnField = OGR_FD_GetFieldIndex( OGR_L_GetLayerDefn( hSrcLayer ),
                                           pszBurnAttribute );
        if( iBurnField == -1 )
        {
            CPLError(CE_Failure, CPLE_AppDefined, "Failed to find field %s on layer %s, skipping.",
                    pszBurnAttribute,
                    OGR_FD_GetName( OGR_L_GetLayerDefn( hSrcLayer ) ) );
            if( hCT != NULL )
                OCTDestroyCoordinateTransformation(hCT);
            return CE_Failure;
        }
    }

/* -------------------------------------------------------------------- */
/*      Collect the geometries from this layer, and build list of       */
/*      burn values.                                                    */
/* -------------------------------------------------------------------- */
    OGRFeatureH hFeat;
    std::vector<OGRGeometryH> ahGeometries;
    std::vector<double> adfFullBurnValues;

    OGR_L_ResetReading( hSrcLayer );

    while( (hFeat = OGR_L_GetNextFeature( hSrcLayer )) != NULL )
    {
        OGRGeometryH hGeom;

        if( OGR_F_GetGeometryRef( hFeat ) == NULL )
        {
            OGR_F_Destroy( hFeat );
//.........这里部分代码省略.........