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<6>()
    {
        // Create feature without geometry
        std::string tmp(data_tmp_);
        tmp += SEP;
        tmp += "tpoly.shp";
        OGRDataSourceH ds = OGR_Dr_Open(drv_, tmp.c_str(), true);
        ensure("Can't open layer", NULL != ds);

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

        OGRFeatureDefnH featDefn = OGR_L_GetLayerDefn(lyr);
        ensure("Layer schema is NULL", NULL != featDefn);

        OGRFeatureH featNonSpatial = OGR_F_Create(featDefn);
        ensure("Can't create non-spatial feature", NULL != featNonSpatial);

        int fldIndex = OGR_FD_GetFieldIndex(featDefn, "PRFEDEA");
        ensure("Can't find field 'PRFEDEA'", fldIndex >= 0);

        OGR_F_SetFieldString(featNonSpatial, fldIndex, "nulled");
       
        OGRErr err = OGR_L_CreateFeature(lyr, featNonSpatial);
        ensure_equals("Can't write non-spatial feature to layer", OGRERR_NONE, err);

        OGR_F_Destroy(featNonSpatial);
        OGR_DS_Destroy(ds);
    }

int            S57_ogrLoadLayer(const char *layername, void *ogrlayer, S52_loadObject_cb loadObject_cb)
{
    if (NULL==layername || NULL==ogrlayer) {
        PRINTF("ERROR: layername || ogrlayer || S52_loadLayer_cb is NULL\n");
        g_assert(0);
    }

    if (NULL == loadObject_cb) {
        static int  silent  = FALSE;
        if (FALSE == silent) {
            PRINTF("NOTE: using default S52_loadObject() callback\n");
            PRINTF("       (this msg will not repeat)\n");
            silent = TRUE;
        }
        loadObject_cb = S52_loadObject;
    }

    OGRFeatureH feature = NULL;
    while ( NULL != (feature = OGR_L_GetNextFeature((OGRLayerH)ogrlayer))) {
        // debug
        //PRINTF("layer:feature %X:%X\n",  ogrlayer, feature);

#ifdef _MINGW
        // on Windows 32 the callback is broken
        S52_loadObject(layername, feature);
#else
        loadObject_cb(layername, (void*)feature);
#endif

        OGR_F_Destroy(feature);
    }

    return TRUE;
}

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;
}

    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);
    }

bool QgsOgrFeatureIterator::readFeature( OGRFeatureH fet, QgsFeature& feature )
{
  feature.setFeatureId( OGR_F_GetFID( fet ) );
  feature.initAttributes( mSource->mFields.count() );
  feature.setFields( &mSource->mFields ); // allow name-based attribute lookups

  bool useIntersect = mRequest.flags() & QgsFeatureRequest::ExactIntersect;
  bool geometryTypeFilter = mSource->mOgrGeometryTypeFilter != wkbUnknown;
  if ( mFetchGeometry || useIntersect || geometryTypeFilter )
  {
    OGRGeometryH geom = OGR_F_GetGeometryRef( fet );

    if ( geom )
    {
      if ( mGeometrySimplifier )
        mGeometrySimplifier->simplifyGeometry( geom );

      // get the wkb representation
      int memorySize = OGR_G_WkbSize( geom );
      unsigned char *wkb = new unsigned char[memorySize];
      OGR_G_ExportToWkb( geom, ( OGRwkbByteOrder ) QgsApplication::endian(), wkb );

      QgsGeometry* geometry = feature.geometry();
      if ( !geometry ) feature.setGeometryAndOwnership( wkb, memorySize ); else geometry->fromWkb( wkb, memorySize );
    }
    if (( useIntersect && ( !feature.geometry() || !feature.geometry()->intersects( mRequest.filterRect() ) ) )
        || ( geometryTypeFilter && ( !feature.geometry() || QgsOgrProvider::ogrWkbSingleFlatten(( OGRwkbGeometryType )feature.geometry()->wkbType() ) != mSource->mOgrGeometryTypeFilter ) ) )
    {
      OGR_F_Destroy( fet );
      return false;
    }
  }

  if ( !mFetchGeometry )
  {
    feature.setGeometry( 0 );
  }

  // fetch attributes
  if ( mRequest.flags() & QgsFeatureRequest::SubsetOfAttributes )
  {
    const QgsAttributeList& attrs = mRequest.subsetOfAttributes();
    for ( QgsAttributeList::const_iterator it = attrs.begin(); it != attrs.end(); ++it )
    {
      getFeatureAttribute( fet, feature, *it );
    }
  }
  else
  {
    // all attributes
    for ( int idx = 0; idx < mSource->mFields.count(); ++idx )
    {
      getFeatureAttribute( fet, feature, idx );
    }
  }

  return true;
}

static void
feature_destroy(ErlNifEnv *env, void *obj)
{
    // Feature is owned by the caller
    EnvFeature_t **feature = (EnvFeature_t**)obj;
    OGR_F_Destroy((**feature).obj);
    enif_free_env((**feature).env);
    enif_free(*feature);
}

bool QgsOgrFeatureIterator::readFeature( OGRFeatureH fet, QgsFeature& feature ) const
{
  feature.setFeatureId( OGR_F_GetFID( fet ) );
  feature.initAttributes( mSource->mFields.count() );
  feature.setFields( mSource->mFields ); // allow name-based attribute lookups

  bool useIntersect = mRequest.flags() & QgsFeatureRequest::ExactIntersect;
  bool geometryTypeFilter = mSource->mOgrGeometryTypeFilter != wkbUnknown;
  if ( mFetchGeometry || useIntersect || geometryTypeFilter )
  {
    OGRGeometryH geom = OGR_F_GetGeometryRef( fet );

    if ( geom )
    {
      feature.setGeometry( QgsOgrUtils::ogrGeometryToQgsGeometry( geom ) );
    }
    else
      feature.clearGeometry();

    if ( mSource->mOgrGeometryTypeFilter == wkbGeometryCollection &&
         geom && wkbFlatten( OGR_G_GetGeometryType( geom ) ) == wkbGeometryCollection )
    {
      // OK
    }
    else if (( useIntersect && ( !feature.hasGeometry() || !feature.geometry().intersects( mRequest.filterRect() ) ) )
             || ( geometryTypeFilter && ( !feature.hasGeometry() || QgsOgrProvider::ogrWkbSingleFlatten(( OGRwkbGeometryType )feature.geometry().wkbType() ) != mSource->mOgrGeometryTypeFilter ) ) )
    {
      OGR_F_Destroy( fet );
      return false;
    }
  }

  if ( !mFetchGeometry )
  {
    feature.clearGeometry();
  }

  // fetch attributes
  if ( mRequest.flags() & QgsFeatureRequest::SubsetOfAttributes )
  {
    QgsAttributeList attrs = mRequest.subsetOfAttributes();
    for ( QgsAttributeList::const_iterator it = attrs.begin(); it != attrs.end(); ++it )
    {
      getFeatureAttribute( fet, feature, *it );
    }
  }
  else
  {
    // all attributes
    for ( int idx = 0; idx < mSource->mFields.count(); ++idx )
    {
      getFeatureAttribute( fet, feature, idx );
    }
  }

  return true;
}

/*!
  \brief Read feature from OGR layer at given offset (level 1)
  
  This function implements random access on level 1.

  \param Map pointer to Map_info structure 
  \param[out] line_p container used to store line points within
  \param[out] line_c container used to store line categories within
  \param offset given offset 
  
  \return line type
  \return 0 dead line
  \return -2 no more features
  \return -1 out of memory
*/
int V1_read_line_ogr(struct Map_info *Map,
		     struct line_pnts *line_p, struct line_cats *line_c, off_t offset)
{
    long FID;
    int type;
    OGRGeometryH hGeom;

    G_debug(4, "V1_read_line_ogr() offset = %lu offset_num = %lu",
	    (long) offset, (long) Map->fInfo.ogr.offset_num);

    if (offset >= Map->fInfo.ogr.offset_num)
	return -2;
    
    if (line_p != NULL)
	Vect_reset_line(line_p);
    if (line_c != NULL)
	Vect_reset_cats(line_c);

    FID = Map->fInfo.ogr.offset[offset];
    G_debug(4, "  FID = %ld", FID);
    
    /* coordinates */
    if (line_p != NULL) {
	/* Read feature to cache if necessary */
	if (Map->fInfo.ogr.feature_cache_id != FID) {
	    G_debug(4, "Read feature (FID = %ld) to cache", FID);
	    if (Map->fInfo.ogr.feature_cache) {
		OGR_F_Destroy(Map->fInfo.ogr.feature_cache);
	    }
	    Map->fInfo.ogr.feature_cache =
		OGR_L_GetFeature(Map->fInfo.ogr.layer, FID);
	    if (Map->fInfo.ogr.feature_cache == NULL) {
		G_fatal_error(_("Unable to get feature geometry, FID %ld"),
			      FID);
	    }
	    Map->fInfo.ogr.feature_cache_id = FID;
	}
	
	hGeom = OGR_F_GetGeometryRef(Map->fInfo.ogr.feature_cache);
	if (hGeom == NULL) {
	    G_fatal_error(_("Unable to get feature geometry, FID %ld"),
			  FID);
	}
	
	type = read_line(Map, hGeom, offset + 1, line_p);
    }
    else {
	type = get_line_type(Map, FID);
    }

    /* category */
    if (line_c != NULL) {
	Vect_cat_set(line_c, 1, (int) FID);
    }

    return type;
}

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;
}

	bool getNextFeature()
	{
		if (m_current_feature)
			OGR_F_Destroy(m_current_feature);
		m_current_feature = OGR_L_GetFeature(m_layer, (long)m_index);
		if (!m_current_feature)
			return false;
		
		m_current_geometry = OGR_F_GetGeometryRef(m_current_feature);
		return true;
	}

bool QgsOgrFeatureIterator::readFeature( OGRFeatureH fet, QgsFeature& feature )
{
  feature.setFeatureId( OGR_F_GetFID( fet ) );
  feature.initAttributes( P->fields().count() );
  feature.setFields( &P->mAttributeFields ); // allow name-based attribute lookups

  bool fetchGeom    = !( mRequest.flags() & QgsFeatureRequest::NoGeometry );
  bool useIntersect = mRequest.flags() & QgsFeatureRequest::ExactIntersect;
  bool geometryTypeFilter = P->mOgrGeometryTypeFilter != wkbUnknown;
  if ( fetchGeom || useIntersect || geometryTypeFilter )
  {
    OGRGeometryH geom = OGR_F_GetGeometryRef( fet );

    if ( geom )
    {
      // get the wkb representation
      unsigned char *wkb = new unsigned char[OGR_G_WkbSize( geom )];
      OGR_G_ExportToWkb( geom, ( OGRwkbByteOrder ) QgsApplication::endian(), wkb );

      feature.setGeometryAndOwnership( wkb, OGR_G_WkbSize( geom ) );
    }
    if (( useIntersect && ( !feature.geometry() || !feature.geometry()->intersects( mRequest.filterRect() ) ) )
        || ( geometryTypeFilter && ( !feature.geometry() || wkbFlatten(( OGRwkbGeometryType )feature.geometry()->wkbType() ) != wkbFlatten( P->mOgrGeometryTypeFilter ) ) ) )
    {
      OGR_F_Destroy( fet );
      return false;
    }
  }

  if ( !fetchGeom )
  {
    feature.setGeometry( 0 );
  }

  // fetch attributes
  if ( mRequest.flags() & QgsFeatureRequest::SubsetOfAttributes )
  {
    const QgsAttributeList& attrs = mRequest.subsetOfAttributes();
    for ( QgsAttributeList::const_iterator it = attrs.begin(); it != attrs.end(); ++it )
    {
      getFeatureAttribute( fet, feature, *it );
    }
  }
  else
  {
    // all attributes
    for ( int idx = 0; idx < P->mAttributeFields.count(); ++idx )
    {
      getFeatureAttribute( fet, feature, idx );
    }
  }

  return true;
}

bool QgsOgrFeatureIterator::fetchFeatureWithId( QgsFeatureId id, QgsFeature& feature ) const
{
  feature.setValid( false );
  OGRFeatureH fet = OGR_L_GetFeature( ogrLayer, FID_TO_NUMBER( id ) );
  if ( !fet )
  {
    return false;
  }

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

  feature.setValid( true );
  return true;
}

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;
}

/*!
  \brief Recursively descend to feature and read the part
  
  \param Map pointer to Map_info structure
  \param hGeom OGR geometry
  \param offset given offset
  \param[out] Points container used to store line pointes within
  
  \return feature type
  \return -1 on error
*/
static int get_line_type(const struct Map_info *Map, long FID)
{
    int eType;
    OGRFeatureH hFeat;
    OGRGeometryH hGeom;

    G_debug(4, "get_line_type() fid = %ld", FID);

    hFeat = OGR_L_GetFeature(Map->fInfo.ogr.layer, FID);
    if (hFeat == NULL)
	return -1;

    hGeom = OGR_F_GetGeometryRef(hFeat);
    if (hGeom == NULL)
	return -1;
    
    eType = wkbFlatten(OGR_G_GetGeometryType(hGeom));

    OGR_F_Destroy(hFeat);

    G_debug(4, "OGR Geometry of type: %d", eType);

    switch (eType) {
    case wkbPoint:
    case wkbMultiPoint:
	return GV_POINT;
	break;
	
    case wkbLineString:
    case wkbMultiLineString:
	return GV_LINE;
	break;

    case wkbPolygon:
    case wkbMultiPolygon:
    case wkbGeometryCollection:
	return GV_BOUNDARY;
	break;

    default:
	G_warning(_("OGR feature type %d not supported"), eType);
	break;
    }

    return -1;
}

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 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);
    }

static CPLErr
EmitPolygonToLayer( OGRLayerH hOutLayer, int iPixValField,
                    RPolygonF *poRPoly, double *padfGeoTransform )

{
    OGRFeatureH hFeat;
    OGRGeometryH hPolygon;

/* -------------------------------------------------------------------- */
/*      Turn bits of lines into coherent rings.                         */
/* -------------------------------------------------------------------- */
    poRPoly->Coalesce();

/* -------------------------------------------------------------------- */
/*      Create the polygon geometry.                                    */
/* -------------------------------------------------------------------- */
    size_t iString;

    hPolygon = OGR_G_CreateGeometry( wkbPolygon );
    
    for( iString = 0; iString < poRPoly->aanXY.size(); iString++ )
    {
        std::vector<int> &anString = poRPoly->aanXY[iString];
        OGRGeometryH hRing = OGR_G_CreateGeometry( wkbLinearRing );

        int iVert;

        // we go last to first to ensure the linestring is allocated to 
        // the proper size on the first try.
        for( iVert = anString.size()/2 - 1; iVert >= 0; iVert-- )
        {
            double dfX, dfY;
            int    nPixelX, nPixelY;
            
            nPixelX = anString[iVert*2];
            nPixelY = anString[iVert*2+1];

            dfX = padfGeoTransform[0] 
                + nPixelX * padfGeoTransform[1]
                + nPixelY * padfGeoTransform[2];
            dfY = padfGeoTransform[3] 
                + nPixelX * padfGeoTransform[4]
                + nPixelY * padfGeoTransform[5];

            OGR_G_SetPoint_2D( hRing, iVert, dfX, dfY );
        }

        OGR_G_AddGeometryDirectly( hPolygon, hRing );
    }
    
/* -------------------------------------------------------------------- */
/*      Create the feature object.                                      */
/* -------------------------------------------------------------------- */
    hFeat = OGR_F_Create( OGR_L_GetLayerDefn( hOutLayer ) );

    OGR_F_SetGeometryDirectly( hFeat, hPolygon );

    if( iPixValField >= 0 )
        OGR_F_SetFieldDouble( hFeat, iPixValField, (double)poRPoly->fPolyValue );

/* -------------------------------------------------------------------- */
/*      Write the to the layer.                                         */
/* -------------------------------------------------------------------- */
    CPLErr eErr = CE_None;

    if( OGR_L_CreateFeature( hOutLayer, hFeat ) != OGRERR_NONE )
        eErr = CE_Failure;

    OGR_F_Destroy( hFeat );

    return eErr;
}

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

  if ( mClosed )
    return false;

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

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

    // skip features without geometry
    if ( !OGR_F_GetGeometryRef( fet ) && !P->mFetchFeaturesWithoutGeom )
    {
      OGR_F_Destroy( fet );
      close();
      return false;
    }

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

  OGRFeatureH fet;
  QgsRectangle selectionRect;

  while (( fet = OGR_L_GetNextFeature( P->ogrLayer ) ) )
  {
    // skip features without geometry
    if ( !P->mFetchFeaturesWithoutGeom && !OGR_F_GetGeometryRef( fet ) )
    {
      OGR_F_Destroy( fet );
      continue;
    }

    readFeature( fet, feature );

    if ( mRequest.flags() & QgsFeatureRequest::ExactIntersect )
    {
      //precise test for intersection with search rectangle
      //first make QgsRectangle from OGRPolygon
      OGREnvelope env;
      memset( &env, 0, sizeof( env ) );
      if ( mSelectionRectangle )
        OGR_G_GetEnvelope( mSelectionRectangle, &env );
      if ( env.MinX != 0 || env.MinY != 0 || env.MaxX != 0 || env.MaxY != 0 ) //if envelope is invalid, skip the precise intersection test
      {
        selectionRect.set( env.MinX, env.MinY, env.MaxX, env.MaxY );
        if ( !feature.geometry() || !feature.geometry()->intersects( selectionRect ) )
        {
          OGR_F_Destroy( fet );
          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;
}

//.........这里部分代码省略.........
	else { hLayer1 = OGR_DS_GetLayer(hDS1,outletslyr);} // get layerinfo from layername

	if(hLayer1 == NULL)getlayerfail(hDS1,outletsds,outletslyr);
	OGRwkbGeometryType gtype;
	gtype=OGR_L_GetGeomType(hLayer1);

	// Test that the type is a point
	if(gtype != wkbPoint)getlayerfail(hDS1,outletsds,outletslyr);

	const char* RasterProjectionName;
	const char* sprs;
	const char* sprso;
	const char* OutletProjectionName;
	int pj_raster,pj_outlet;

	// Spatial reference of outlet
	hRSOutlet = OGR_L_GetSpatialRef(hLayer1);
	if(hSRSRaster!=NULL){
	  pj_raster=OSRIsProjected(hSRSRaster); // find if projected or not
	  if(pj_raster==0) {sprs="GEOGCS";} else { sprs="PROJCS"; }
	  RasterProjectionName = OSRGetAttrValue(hSRSRaster,sprs,0); // get projection name
	}
	if(hRSOutlet!=NULL){
	  pj_outlet=OSRIsProjected(hRSOutlet);
	  if(pj_outlet==0) {sprso="GEOGCS";} else { sprso="PROJCS"; }
	  OutletProjectionName = OSRGetAttrValue(hRSOutlet,sprso,0);
	}

	//Write warnings where projections may not match
	if(hRSOutlet!=NULL && hSRSRaster!=NULL){
	
		if (pj_raster==pj_outlet){
			  
			 int rc= strcmp(RasterProjectionName,OutletProjectionName); // compare string
			 if(rc!=0){
				printf( "Warning: Projection of Outlet feature and Raster data may be different.\n" );
				printf("Projection of Raster datasource %s.\n",RasterProjectionName);
                printf("Projection of Outlet feature %s.\n",OutletProjectionName);
			}
		}
    
		else {
			  printf( "Warning: Spatial References of Outlet feature and Raster data are different.\n" );
			  printf("Projection of Raster datasource %s.\n",RasterProjectionName);
              printf("Projection of Outlet feature %s.\n",OutletProjectionName);
		}
	}
	
	else if(hSRSRaster==NULL && hRSOutlet!=NULL) {
		      printf( "Warning: Spatial Reference of Raster is missing.\n" );
              printf("Projection of Outlet feature %s.\n",OutletProjectionName);

	}
	else if(hSRSRaster!=NULL && hRSOutlet==NULL) {
	          printf( "Warning: Spatial Reference of Outlet feature is missing.\n" );
			  printf("Projection of Raster datasource %s.\n",RasterProjectionName);
	}
	else {
	          printf( "Warning: Spatial References of Outlet feature and Raster data are missing.\n" );
	}



	long countPts=0;
	// count number of feature
	countPts=OGR_L_GetFeatureCount(hLayer1,0); 
	// get schema i.e geometry, properties (e.g. ID)
	hFDefn1 = OGR_L_GetLayerDefn(hLayer1); 
	x = new double[countPts];
	y = new double[countPts];
	int iField;
	int nxy=0;
	id = new int[countPts];
	// loop through each feature and get lat,lon and id information

    OGR_L_ResetReading(hLayer1);
    while( (hFeature1 = OGR_L_GetNextFeature(hLayer1)) != NULL ) {

		 //hFeature1=OGR_L_GetFeature(hLayer1,j); // get feature info
		 geometry = OGR_F_GetGeometryRef(hFeature1); // get geometry
         x[nxy] = OGR_G_GetX(geometry, 0);
		 y[nxy] =  OGR_G_GetY(geometry, 0);
		 int idfld =OGR_F_GetFieldIndex(hFeature1,"id");
		 if (idfld >= 0)
		   {
			 
			hFieldDefn1 = OGR_FD_GetFieldDefn( hFDefn1,idfld); // get field definiton based on index
			if( OGR_Fld_GetType(hFieldDefn1) == OFTInteger ) {
					id[nxy] =OGR_F_GetFieldAsInteger( hFeature1, idfld );} // get id value 
		    }
		 else {
		      id[nxy]=1;// if there is no id field         
		 } 
		 nxy++; // count number of outlets point
		 OGR_F_Destroy( hFeature1 ); // destroy feature
	}
	*noutlets=nxy;
	OGR_DS_Destroy( hDS1); // destroy data source
	return 0;
}